TSP with refuel stations

tsp_refuel.py

#!/usr/bin/env python3
from __future__ import print_function
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
import math
import numpy as np
import matplotlib.pyplot as plt


def create_data_model():
    data = {}
    fuel_capacity = 250
    _locations = [
        (0, 2),  # start
        (1, 2),
        (2, 2),
        (3, 2),
        (4, 2),
        (5, 2), # end
        (1, 3),  # locations to visit
        (2, 4), (3, 1),
        (4, 0)]

    data["locations"] = [(l[0] * 50, l[1] * 50) for l in _locations]
    data["num_locations"] = len(data["locations"])
    # data["demands"] = [0, fuel_capacity, fuel_capacity, fuel_capacity, fuel_capacity, fuel_capacity, -10, -15, -10, -15] # changing this caused to give a sensible solution
    data["time_windows"] = [(0, 5),
                            (0, 5),
                            (4, 10),
                            (6, 15),
                            (7, 18),
                            (8, 13),
                            (0, 3),
                            (7, 9),
                            (10, 13),
                            (14, 15)]
    data["num_vehicles"] = 1
    data["fuel_capacity"] = fuel_capacity
    data["vehicle_speed"] = 35
    data["starts"] = [0]
    data["ends"] = [5]
    distance_matrix = np.zeros((data["num_locations"], data["num_locations"]), dtype=int)
    for i in range(data["num_locations"]):
        for j in range(data["num_locations"]):
            if i == j:
                distance_matrix[i][j] = 0
            else:
                distance_matrix[i][j] = euclidean_distance(data["locations"][i], data["locations"][j])
    dist_matrix = distance_matrix.tolist()
    #print(dist_matrix)
    data["distance_matrix"] = dist_matrix
    assert len(data['distance_matrix']) == len(data['locations'])
    assert len(data['distance_matrix']) == len(data['time_windows'])
    assert len(data['starts']) == len(data['ends'])
    assert data['num_vehicles'] == len(data['starts'])
    return data


def euclidean_distance(position_1, position_2):
    return int(math.hypot((position_1[0] - position_2[0]), (position_1[1] - position_2[1])))


def print_solution(data, manager, routing, solution):
    print("Objective: {}".format(solution.ObjectiveValue()))
    total_distance = 0
    total_load = 0
    total_time = 0
    fuel_dimension = routing.GetDimensionOrDie("Fuel")
    time_dimension = routing.GetDimensionOrDie("Time")
    dropped_nodes = "Dropped nodes:"
    for node in range(routing.Size()):
        if routing.IsStart(node) or routing.IsEnd(node):
            continue
        if solution.Value(routing.NextVar(node)) == node:
            dropped_nodes += " {}".format(manager.IndexToNode(node))
    print(dropped_nodes)
    dum_list = [(0, 100)]
    for vehicle_id in range(data["num_vehicles"]):
        index = routing.Start(vehicle_id)
        plan_output = "Route for vehicle {}:\n".format(vehicle_id)
        distance = 0
        while not routing.IsEnd(index):
            fuel_var = fuel_dimension.CumulVar(index)
            time_var = time_dimension.CumulVar(index)
            plan_output += "{0} Fuel({1}) Time({2},{3}) -> ".format(
                    manager.IndexToNode(index),
                    solution.Value(fuel_var),
                    solution.Min(time_var),
                    solution.Max(time_var))
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            if index <= 8:
                if index == 1 or index == 2 or index == 3 or index == 4:
                    dum_list.append(data["locations"][index])
                else:
                    dum_list.append(data["locations"][index + 1])
            #print(f"plop {routing.GetArcCostForVehicle(previous_index, index, vehicle_id)}\n")
            distance += routing.GetArcCostForVehicle(previous_index, index, vehicle_id)
        fuel_var = fuel_dimension.CumulVar(index)
        time_var = time_dimension.CumulVar(index)
        plan_output += "{0} Fuel({1}) Time({2},{3})\n".format(
                manager.IndexToNode(index),
                solution.Value(fuel_var),
                solution.Min(time_var),
                solution.Max(time_var))
        plan_output += "Distance of the route: {}units\n".format(distance)
        plan_output += "Remaining Fuel of the route: {}\n".format(solution.Value(fuel_var))
        plan_output += "Total Time of the route: {}\n".format(solution.Value(time_var))
        print(plan_output)
        total_distance += distance
        total_load += solution.Value(fuel_var)
        total_time += solution.Value(time_var)
    print('Total Distance of all routes: {}units'.format(total_distance))
    print('Total Fuel consumed of all routes: {}units'.format(total_load))
    print('Total Time of all routes: {}units'.format(total_time))

    #dum_list.append((250, 100))
    #x1 = []
    #y1 = []
    #for i in dum_list:
    #    x1.append(i[0])
    #    y1.append(i[1])
    #fig = plt.figure()
    #ax = fig.add_subplot(111)
    #for x, y in data["locations"]:
    #    ax.plot(x, y, 'r.', markersize=15)
    ## plot_vehicle_routes(dum_list, ax, data['locations'], data['num_vehicles'])
    #for i in range(len(x1)):
    #    if i <= len(x1) - 2:
    #        ax.arrow(x1[i], y1[i], (x1[i+1] - x1[i]), (y1[i+1] - y1[i]), width=3)
    #    # else:
    #    #     ax.arrow(x1[6], y1[6], (x1[0] - x1[6]), (y1[0] - y1[6]), width=2)
    #plt.show()


def main():
    # Instantiate the data problem.
    data = create_data_model()

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(
            len(data["distance_matrix"]),
            data["num_vehicles"],
            data["starts"],
            data["ends"])

    # Create Routing Model.
    routing = pywrapcp.RoutingModel(manager)

    # Distance
    def distance_callback(from_index, to_index):
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data["distance_matrix"][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Fuel constraints
    def fuel_callback(from_index, to_index):
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return -euclidean_distance(data["locations"][from_node], data["locations"][to_node])

    fuel_callback_index = routing.RegisterTransitCallback(fuel_callback)
    #fuel_callback_index = routing.RegisterUnaryTransitCallback(fuel_callback)
    routing.AddDimension(
            fuel_callback_index,
            data["fuel_capacity"],
            data["fuel_capacity"],
            False,
            'Fuel')

    penalty = 0
    fuel_dimension = routing.GetDimensionOrDie('Fuel')
    fuel_dimension.SlackVar(routing.Start(0)).SetValue(0)
    for node in range(len(data["distance_matrix"])):
        if node == 0 or node == 5:
            continue
        if node > 5:
            index = manager.NodeToIndex(node)
            fuel_dimension.SlackVar(index).SetValue(0)
            routing.AddVariableMinimizedByFinalizer(fuel_dimension.CumulVar(node))
        else:
            index = manager.NodeToIndex(node)
            routing.AddDisjunction([index], penalty)


    # Time
    def time_callback(from_index, to_index):
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return int(data["distance_matrix"][from_node][to_node] / data["vehicle_speed"])

    time_callback_index = routing.RegisterTransitCallback(time_callback)
    routing.AddDimension(
            time_callback_index,
            3,
            300,
            False,
            'Time')

    time_dimension = routing.GetDimensionOrDie('Time')
    for location_idx, time_window in enumerate(data["time_windows"]):
        if location_idx == 0 or location_idx == 5:
            continue
        index = manager.NodeToIndex(location_idx)
        time_dimension.CumulVar(index).SetRange(time_window[0], time_window[1])
        routing.AddToAssignment(time_dimension.SlackVar(index))
    # Add time window constraints for each vehicle start node
    # and "copy" the slack var in the solution object (aka Assignment) to print it
    for vehicle_id in range(data["num_vehicles"]):
        index = routing.Start(vehicle_id)
        time_dimension.CumulVar(index).SetRange(data["time_windows"][0][0], data["time_windows"][0][1])
        routing.AddToAssignment(time_dimension.SlackVar(index))

    for i in range(data['num_vehicles']):
        routing.AddVariableMinimizedByFinalizer(
            time_dimension.CumulVar(routing.Start(i)))
        routing.AddVariableMinimizedByFinalizer(
            time_dimension.CumulVar(routing.End(i)))

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        #routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
        routing_enums_pb2.FirstSolutionStrategy.GLOBAL_CHEAPEST_ARC)
    search_parameters.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
    search_parameters.time_limit.FromSeconds(5)

    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)

    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)

    print("Solver status:", routing.status())


if __name__ == '__main__':
    main()

zz.md

Possible output:

[*v*]─mizux@nuc10i7 %./last.py
Objective: 491
Dropped nodes: 1 3 4
Route for vehicle 0:
0 Fuel(241) Time(0,0) -> 6 Fuel(171) Time(2,3) -> 7 Fuel(101) Time(7,8) -> 2 Fuel(1) Time(9,10) -> 8 Fuel(181) Time(11,12) -> 9 Fuel(111) Time(14,14) -> 5 Fuel(0) Time(17,17)
Distance of the route: 491units
Remaining Fuel of the route: 0
Total Time of the route: 17

Total Distance of all routes: 491units
Total Fuel consumed of all routes: 0units
Total Time of all routes: 17units
Solver status: 1