How do I use Google OR Tools to add disjunctions, set penalties, and prevent certain locations from being dropped in VRPTW?

Question

I have a working Vehicle Routing Problem (with Time Windows) solution implemented using Google's OR Tools python library. I have a time matrix of 15 locations and time windows for each location. I also have the duration of a visit to each location baked into the cost of each visit. All values are in units of seconds. I am intentionally solving this with just one vehicle (essentially solving a Traveling Salesperson Problem).

I am not trying to add the ability to drop locations from a solution if they are preventing a valid solution from being created. Right now, if I have the duration for each visit set to 3600 seconds, it won't be possible to visit all 15 locations. However, if I instead set the duration for each visit to 900 seconds, then a solution will be possible for all 15 locations. I want to add a disjunction to allow for a solution to be created with these large durations in place, and instead of failing, simply drop a location from the solution.

I have some locations that I do not want to be dropped from the solution, so I have given them absurdly large penalties to ensure that they won't be dropped, while others I assign a penalty of zero. But now, all the other locations are being dropped because their penalty is zero - I assume that this is because the penalty is less than the cost of transit, but I am not completely sure if this is indeed the reason. How should I allow locations to be dropped from the solution, but prevent other locations from being droppable?

Right now the only code I have added is:

# Allow to drop nodes.
for node in range(1, len(Penalties)):
  routing.AddDisjunction([manager.NodeToIndex(node)], Penalties[node])

Source

from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import routing_enums_pb2

Matrix = [
  [0, 557, 763, 813, 618, 822, 700, 1527, 112, 1011, 734, 551, 604, 1156, 732],             # Depot
  [523, 0, 598, 934, 607, 658, 535, 1529, 589, 857, 424, 475, 725, 1107, 899],              # 1 - Location
  [631, 480, 0, 960, 570, 451, 135, 1698, 582, 1075, 642, 743, 751, 968, 925],              # 2 - Location
  [746, 1000, 1135, 0, 1168, 1186, 1071, 1012, 776, 1358, 1162, 947, 594, 1283, 426],       # 3 - Location
  [685, 627, 810, 990, 0, 712, 709, 1649, 550, 1221, 788, 726, 734, 1227, 653],             # 4 - Location
  [869, 718, 558, 1105, 650, 0, 384, 1328, 821, 1313, 880, 981, 989, 732, 993],             # 5 - Location
  [679, 528, 202, 1008, 618, 412, 0, 1740, 630, 1123, 690, 791, 799, 878, 973],             # 6 - Location
  [1378, 1553, 1766, 1031, 1595, 1355, 1731, 0, 1408, 1990, 1715, 1579, 1226, 1452, 1098],  # 7 - Location
  [149, 626, 762, 696, 556, 821, 698, 1410, 0, 999, 803, 536, 487, 1156, 614],              # 8 - Location
  [918, 943, 1184, 1296, 1193, 1244, 1121, 2010, 1030, 0, 509, 870, 1063, 1693, 1250],      # 9 - Location
  [763, 541, 782, 1118, 791, 842, 719, 1713, 719, 558, 0, 567, 909, 1291, 1083],            # 10 - Location
  [449, 543, 843, 887, 769, 902, 780, 1601, 561, 876, 573, 0, 678, 1352, 806],              # 11 - Location
  [628, 873, 1008, 491, 933, 1068, 945, 1205, 657, 1014, 1035, 825, 0, 1276, 444],          # 12 - Location
  [1343, 1247, 1367, 1270, 1289, 809, 1193, 1335, 1253, 1842, 1409, 1509, 1287, 0, 1158],   # 13 - Location
  [520, 774, 909, 385, 875, 1052, 845, 1078, 550, 1132, 936, 721, 368, 1149, 0]             # 14 - Location
]

Windows = [
  [28800, 28800],   # Depot
  [43200, 43200],   # 1 - Location
  [50400, 50400],   # 2 - Location
  [21600, 79200],   # 3 - Location
  [21600, 79200],   # 4 - Location
  [21600, 79200],   # 5 - Location
  [21600, 79200],   # 6 - Location
  [21600, 79200],   # 7 - Location
  [21600, 79200],   # 8 - Location
  [21600, 79200],   # 9 - Location
  [21600, 79200],   # 10 - Location
  [21600, 79200],   # 11 - Location
  [21600, 79200],   # 12 - Location
  [21600, 79200],   # 13 - Location
  [21600, 79200]    # 14 - Location
]

Durations = [0, 1800, 3600, 3600, 7200, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800, 1800]

Penalties = [99999999, 99999999, 99999999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

# The inputs to RoutingIndexManager are:
#   1. The number of rows of the time matrix, which is the number of locations (including the depot).
#   2. The number of vehicles in the problem.
#   3. The node corresponding to the depot.

# Create the routing index manager.
manager = pywrapcp.RoutingIndexManager(len(Matrix), 1, 0)

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

# Create and register a transit callback.
def transit_callback(from_index, to_index):
  # Returns the travel time between the two nodes.
  # Convert from routing variable Index to time matrix NodeIndex.
  from_node = manager.IndexToNode(from_index)
  to_node = manager.IndexToNode(to_index)
  return Matrix[from_node][to_node] + Durations[from_node]

transit_callback_index = routing.RegisterTransitCallback(transit_callback)

# Define cost of each arc.
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

# Add Time Windows constraint.
routing.AddDimension(
    transit_callback_index,
    64800,  # An upper bound for slack (the wait times at the locations).
    64800,  # An upper bound for the total time over each vehicle's route.
    False,  # Determine whether the cumulative variable is set to zero at the start of the vehicle's route.
    'Time')
time_dimension = routing.GetDimensionOrDie('Time')

# Allow to drop nodes.
for node in range(1, len(Penalties)):
  routing.AddDisjunction([manager.NodeToIndex(node)], Penalties[node])

# Add time window constraints for each location except depot.
for location_idx, time_window in enumerate(Windows):
  if location_idx == 0:
    continue
  index = manager.NodeToIndex(location_idx)
  time_dimension.CumulVar(index).SetRange(time_window[0], time_window[1])

# Add time window constraints for each vehicle start node.
index = routing.Start(0)
time_dimension.CumulVar(index).SetRange(Windows[0][0],Windows[0][1])

# Instantiate route start and end times to produce feasible times.
routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.Start(0)))
routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.End(0)))

# Setting first solution heuristic. 
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

# Setting local search metaheuristics:
search_parameters.local_search_metaheuristic = (routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
search_parameters.time_limit.seconds = 5
search_parameters.log_search = False

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

# Display dropped nodes.
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)

# Return the solution.
time = 0
index = routing.Start(0)
print('Solution:')
while not routing.IsEnd(index):
  time = time_dimension.CumulVar(index)
  print('{0} ({1},{2})'.format(manager.IndexToNode(index), solution.Min(time), solution.Max(time)))
  index = solution.Value(routing.NextVar(index))
time = time_dimension.CumulVar(index)
print('{0} ({1},{2})'.format(manager.IndexToNode(index), solution.Min(time), solution.Max(time)))

Output

Dropped nodes: 3 4 5 6 7 8 9 10 11 12 13 14
Solution:
0 (28800,28800)
1 (43200,43200)
2 (50400,50400)
0 (54631,54631)

Again, if I remove those 3 lines of code that I added, and set all of the values in the Duration array to be equal to 900 seconds, a solution is create just fine. But when I increase the durations, a solution is not able to be created. And when I add the Disjunction and set all penalties equal to zero, the solution omits every location except for the depot. Are there any blatant errors in my usage of OR Tools? Any help will be greatly appreciate!

Answer 1

1. To make location "mandatory" you must use the max int64 value ( 0x7FFFFFFFFFFFFFF) since solver can't overflow it will forbid it to drop this location.

2. Since you are using your time matrix to feed the arc cost evaluator you should set your penalty around 10k otherwise solver has incentive to drop node and paying the cost than paying the transit cost.

3. All your time window ranges should be in [0, max dimension value], here you set it to 64800 BUT you have some Time Windows with 79200 as upper bound, this can make the problem infeasible from the solver point of view so you should set the time dimension upper bound to at least 79200 IMHO.