Google OR-Tools VRP - Nodes getting dropped unnecessarily

Question

I am trying to solve a vehicle routing problem with time window constraints where the solver is allowed to drop nodes if no feasible solution is found. However, I have found that nodes that getting dropped unnecessarily after adding disjunctions, even after imposing a large penalty.

Below is a simple example program so illustrate the issue. The solver drops node 1 and returns the solution of 0 -> 3 -> 2 -> 0. The correct route of 0 -> 1 -> 2 -> 3 -> 0 is returned if the routing.AddDisjunction([manager.NodeToIndex(node)], penalty) code is commented out.

Am I going about this the wrong way? Any help would be much appreciated.

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

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data['time_matrix'] = [
        [0, 3, 2, 1], # depot
        [3, 0, 3, 3],
        [2, 3, 0, 2],
        [1, 3, 2, 0],
    ]
    data['time_windows'] = [
        (0, 0),  # depot
        (0, 3),  
        (0, 6),  
        (0, 9),  
    ]
    data['num_vehicles'] = 1
    data['depot'] = 0

    return data


def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    # Display dropped nodes.
    dropped_nodes = '\nDropped 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 + '\n')

    time_dimension = routing.GetDimensionOrDie('Time')
    for vehicle_id in range(data['num_vehicles']):
        index = routing.Start(vehicle_id)
        plan_output = 'Route for vehicle {}:\n'.format(vehicle_id)
        while not routing.IsEnd(index):
            time_var = time_dimension.CumulVar(index)
            plan_output += '{0} Time ({1}, {2}) -> '.format(
                manager.IndexToNode(index), 
                solution.Min(time_var),
                solution.Max(time_var))
            index = solution.Value(routing.NextVar(index))
        time_var = time_dimension.CumulVar(index)
        plan_output += '{0} Time({1},{2})\n'.format(manager.IndexToNode(index),
                                                    solution.Min(time_var),
                                                    solution.Max(time_var))
        plan_output += '\nTime of the route: {}min\n'.format(
            solution.Min(time_var))
        print(plan_output)



def main():
    """Solve the VRP with time windows."""
    # Instantiate the data problem.
    data = create_data_model()

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(len(data['time_matrix']),
                                           data['num_vehicles'], data['depot'])

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

    # Create and register a transit callback.
    def time_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 data['time_matrix'][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(time_callback)

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

    # Add Time Windows constraint.
    time = 'Time'
    routing.AddDimension(
        transit_callback_index,
        9999999,  # maximum waiting time
        9999999,  # maximum travel time per vehicle
        False,  # Don't force start cumul to zero.
        time)
    time_dimension = routing.GetDimensionOrDie(time)

    # Add time window constraints for each location except depot.
    for location_idx, time_window in enumerate(data['time_windows']):
        if location_idx == data['depot']:
            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.
    depot_idx = data['depot']
    for vehicle_id in range(data['num_vehicles']):
        index = routing.Start(vehicle_id)
        time_dimension.CumulVar(index).SetRange(
            data['time_windows'][depot_idx][0],
            data['time_windows'][depot_idx][1])

    # Instantiate route start and end times to produce feasible times.
    for i in range(data['num_vehicles']):
        routing.AddVariableMinimizedByFinalizer(
            time_dimension.CumulVar(routing.Start(i)))
        routing.AddVariableMinimizedByFinalizer(
            time_dimension.CumulVar(routing.End(i)))

    # Allow to drop nodes.
    penalty = 100000000000
    for node in range(1, len(data['time_matrix'])):
        routing.AddDisjunction([manager.NodeToIndex(node)], penalty)

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

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

    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)
    else:
        print("\nNo solutions found")


if __name__ == '__main__':
    main()

Answer 1

You were on the good direction with adding the disjunction you just miss this lines at the end

    # 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 # to get some logs
-    search_parameters.time_limit.seconds = 10
+    search_parameters.time_limit.seconds = 1 # 1s is large enough ;)

I.e. you forget to enable the Guided Local Search (GLS) so you end up having only the first solution (with the node 1 dropped) and you didn't run the GLS so solver has no chance to bring it back to the solution...