Constraint on Station TW for Subramanian

problem.py

#!/usr/bin/env python3
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
import math
import numpy as np


def create_data_model():
    data = {}
    fuel_capacity = 40000
    _locations = [
        (2.5, 2.5),  # start
        (2.43, 2.094),
        (2.365, 1.688),
        (2.3, 1.28),
        (2.23, 0.875),
        (2.173, 1.28),
        (2.115, 1.688),
        (2.06, 2.09),
        (2, 2.5),
        (1.695, 2.81),
        (1.39, 3.11),
        (1.085, 3.415),
        (0.78, 3.72),
        (1.57, 3.765),
        (2.36, 3.81),
        (3.15, 3.855),
        (3.94, 3.9),  # end
        (1.3, 1.1),  # locations to visit
        (0.6, 4.2), (4.6, 4.0),
        (1.3, 4.2), (4.3, 0.5), (4.8, 1.3),
        (0.6, 0.2), (5, 2.1), (3.1, 1.3), (2.0, 2.5),
        (3.2, 3.9), (4.1, 4.7), (3.3, 0.6), (3.3, 4.5), (0.7, 2.8),
        (4.1, 2.6), (1.2, 1.0), (0.7, 3.2), (0.3, 0.3), (4.3, 4.7), (2, 0),
        (0.2, 1.7), (0.5, 4.2), (0.7, 0.4), (4, 2.9)]
    data["coordinates"] = _locations

    data["locations"] = [(l[0] * 5280, l[1] * 5280) for l in _locations]
    data["num_locations"] = len(data["locations"])

    data["time_windows"] = [
    (0, 60),  # veh_start_node
    (860, 960),  # refuelstation_1
    (1660, 1760),
    (2460, 2560),
    (3260, 3360),
    (4060, 4160),
    (4860, 4960),
    (5660, 5766),
    (6460, 6560),
    (7260, 7360),
    (8060, 8169),
    (8860, 8960),
    (9660, 9760),
    (10460, 10560),
    (11260, 11360),
    (12060, 12160),  # refuel_station_15
    (12860, 100000),  # veh_end_node
    (0, 100000),
    (7, 100000),
    (10, 100000),
    (14, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
    (0, 100000),
     ]

    data["time_windows_stations"] = [
    (0, 60), # start
    (860, 960), # station 1
    (1660, 1760),
    (2460, 2560),
    (3260, 3360),
    (4060, 4160),
    (4860, 4960),
    (5660, 5766),
    (6460, 6560),
    (7260, 7360),
    (8060, 8169),
    (8860, 8960),
    (9660, 9760),
    (10460, 10560),
    (11260, 11360),
    (12060, 12160), # station 15
    ]

    data["num_vehicles"] = 2
    data["fuel_capacity"] = fuel_capacity
    data["vehicle_speed"] = 33
    data["starts"] = [0, 0]
    data["ends"] = [16, 16]
    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()
    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'])
    data["stations"] = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
    return data

def euclidean_distance(position_1, position_2):
    return round(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_fuel = 0
    total_time = 0
    #distance_dimension = routing.GetDimensionOrDie("Distance")
    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)
    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):
            #dist_var = distance_dimension.CumulVar(index)
            fuel_var = fuel_dimension.CumulVar(index)
            time_var = time_dimension.CumulVar(index)
            slack_var = time_dimension.SlackVar(index)
            plan_output += "{0} Fuel({1}) Time({2},{3}) Slack({4},{5}) -> ".format(
                manager.IndexToNode(index),
                solution.Value(fuel_var),
                solution.Min(time_var), solution.Max(time_var),
                solution.Min(slack_var), solution.Max(slack_var))
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            #print(f"plop {routing.GetArcCostForVehicle(previous_index, index, vehicle_id)}\n")
            # distance += routing.GetArcCostForVehicle(previous_index, index, vehicle_id)
        #dist_var = distance_dimension.CumulVar(index)
        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: {} ft\n".format(solution.Value(dist_var))
        plan_output += "Remaining Fuel of the route: {}\n".format(solution.Value(fuel_var))
        plan_output += "Total Time of the route: {} seconds\n".format(solution.Value(time_var))
        print(plan_output)
        #total_distance += solution.Value(dist_var)
        total_fuel += solution.Value(fuel_var)
        total_time += solution.Value(time_var)
    #print('Total Distance of all routes: {} ft'.format(total_distance))
    print('Total Fuel remaining of all routes: {}'.format(total_fuel))
    print('Total Time of all routes: {} seconds'.format(total_time))

def main():
    """Entry point of the program."""
    # Instantiate the data problem.
    # [START data]
    data = create_data_model()
    # [END data]

    # Create the routing index manager.
    # [START index_manager]
    manager = pywrapcp.RoutingIndexManager(
            len(data['time_windows']),
            data['num_vehicles'],
            data['starts'],
            data['ends'])
    # [END index_manager]

    # Create Routing Model.
    # [START routing_model]
    routing = pywrapcp.RoutingModel(manager)
    # [END routing_model]

    # Distance
    stations = data["stations"]
    #[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]  # depot + refill stations

    # Create and register a transit callback.
    # [START transit_callback]
    def time_callback(from_index, to_index):
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        if from_node in stations and from_node != 0:
            return 1000 + int(data["distance_matrix"][from_node][to_node] / data["vehicle_speed"])
        else:
            return 300 + int(data["distance_matrix"][from_node][to_node] / data["vehicle_speed"])


    time_callback_index = routing.RegisterTransitCallback(time_callback)
    routing.SetArcCostEvaluatorOfAllVehicles(time_callback_index)
    # [END transit_callback]

    routing.AddDimension(
            time_callback_index,
            300, # waiting time
            50_000, # vehicle max time
            False, # don't force cumul to zero
            'Time')

    time_dimension = routing.GetDimensionOrDie('Time')
    time_dimension.SetGlobalSpanCostCoefficient(100)
    for location_idx, time_window in enumerate(data["time_windows"]):
        if location_idx in data["starts"] or location_idx in data["ends"]:
            continue
        elif location_idx in data["stations"] and location_idx != 0:
            index = manager.NodeToIndex(location_idx)
            time_dimension.CumulVar(index).SetRange(time_window[0],
                                                    time_window[0] + 1500)
            routing.AddToAssignment(time_dimension.SlackVar(index))
        else:
            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 from_node in range(len(data["time_windows"])):
        for to_node, time_window2 in enumerate(data["time_windows_stations"]):
            if from_node == to_node or to_node in [0,16] or from_node in [0,16]:
                continue
            from_index = manager.NodeToIndex(from_node)
            to_index = manager.NodeToIndex(to_node)
            expr = (routing.NextVar(from_index) == to_index)
            print(f"if NextVar({from_node}) == {to_node}")
            next_node = to_node + 1
            if next_node in stations:
                next_index = manager.NodeToIndex(next_node)
                routing.solver().Add(
                        (expr * time_dimension.CumulVar(next_index)) >= (expr * (time_window2[0] + 1800)))
                routing.solver().Add(
                        (expr * time_dimension.CumulVar(next_index)) <= (expr * (time_window2[1] + 1800)))

            next_node_2 = to_node + 2
            if next_node_2 in stations:
                next_index_2 = manager.NodeToIndex(next_node_2)
                routing.solver().Add(
                        (expr * time_dimension.CumulVar(next_index_2)) >= (expr * (time_window2[0] + 2600)))
                routing.solver().Add(
                        (expr * time_dimension.CumulVar(next_index_2)) <= (expr * (time_window2[1] + 2600)))

            next_node_3 = to_node + 3
            if next_node_3 in stations:
                next_index_3 = manager.NodeToIndex(next_node_3)
                routing.solver().Add(
                        (expr * time_dimension.CumulVar(next_index_3)) >= (expr * (time_window2[0] + 3400)))
                routing.solver().Add(
                        (expr * time_dimension.CumulVar(next_index_3)) <= (expr * (time_window2[1] + 3400)))

    # Fuel constraints
    def fuel_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]

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

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

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

    # Solve the problem.
    # [START solve]
    solution = routing.SolveWithParameters(search_parameters)
    # [END solve]

    # Print solution on console.
    # [START print_solution]
    if solution:
        print_solution(data, manager, routing, solution)
    else:
        print("No solution found !")
    # [END print_solution]
    print("Solver status:", routing.status())


if __name__ == '__main__':
    main()

zz.md

Possible output:

./problem.py
...
Objective: 839925
Dropped nodes: 1 2 5 6 7 9 10 11 12 13 14 15
Route for vehicle 0:
0 Fuel(39923) Time(0,0) Slack(0,0) -> 27 Fuel(31658) Time(550,550) Slack(0,0) -> 30 Fuel(28446) Time(947,947) Slack(0,0) -> 28 Fuel(24092) Time(1378,1378) Slack(0,0) -> 36 Fuel(23036) Time(1710,1710) Slack(0,0) -> 19 Fuel(19015) Time(2131,2131) Slack(0,0) -> 41 Fuel(12399) Time(2631,2631) Slack(0,0) -> 3 Fuel(0) Time(3306,3306) Slack(0,0) -> 37 Fuel(27844) Time(4516,4516) Slack(0,0) -> 40 Fuel(20662) Time(5033,5033) Slack(0,0) -> 23 Fuel(19481) Time(5368,5368) Slack(0,0) -> 35 Fuel(17811) Time(5718,5718) Slack(0,0) -> 38 Fuel(10400) Time(6242,6242) Slack(0,0) -> 8 Fuel(0) Time(6857,6857) Slack(0,0) -> 16 Fuel(0) Time(8239,8239)
Remaining Fuel of the route: 0
Total Time of the route: 8239 seconds

Route for vehicle 1:
0 Fuel(39731) Time(60,60) Slack(0,0) -> 26 Fuel(37091) Time(440,440) Slack(0,0) -> 20 Fuel(27384) Time(1034,1034) Slack(0,0) -> 18 Fuel(23688) Time(1446,1446) Slack(0,0) -> 39 Fuel(23160) Time(1762,1762) Slack(0,0) -> 34 Fuel(17775) Time(2225,2225) Slack(0,0) -> 31 Fuel(15663) Time(2589,2589) Slack(0,0) -> 33 Fuel(5799) Time(3187,3187) Slack(0,0) -> 17 Fuel(5052) Time(3509,3509) Slack(298,298) -> 4 Fuel(0) Time(4260,4260) Slack(0,0) -> 25 Fuel(30837) Time(5414,5414) Slack(0,0) -> 29 Fuel(26993) Time(5830,5830) Slack(0,0) -> 21 Fuel(21687) Time(6290,6290) Slack(0,0) -> 22 Fuel(16706) Time(6740,6740) Slack(0,0) -> 24 Fuel(12352) Time(7171,7171) Slack(0,0) -> 32 Fuel(6916) Time(7635,7635) Slack(0,0) -> 16 Fuel(0) Time(8144,8144)
Remaining Fuel of the route: 0
Total Time of the route: 8144 seconds

Total Fuel remaining of all routes: 0
Total Time of all routes: 16383 seconds
Solver status: 1