Repeated vehicle-breaks: "SetBreakDistanceDurationOfVehicle" sounds fitting but doesn't behave as expected

[sschnug]

(This is borderline issue/discussion: even if not a bug, there might be an issue here in regards to documentation / example-availability/ test-availability | feel free to convert)

I'm trying to model repeated vehicle-breaks (e.g. 540 min break after 1080 min of driving).

"Repeated" meaning, that the tour-duration is long enough to imply a need for more than one break (to keep feasibility).

The feasible solution-space should be reduced as much as possibe implying, that we cannot a-priori fix those breaks but should schedule them in clever fashion:

• Taking a break EARLY (e.g. after 800min of driving) might have positive effects

Now taking an break early (e.g. the first of many breaks) implies a shift of break-starts towards the root (time=0) for all following breaks as breaks are defined in "relative" manner.

There are two "customized" approaches (not outlined in detail here) which come to mind:

• dummy-nodes (like recharging) -> powerful concept but probably less good in terms of scalability
• cp-model based constraints which link break-start- and and break-ends -> probably rather disruptive for the solver

Skimming through the sources, i found:

/// With breaks supposed to be consecutive, this forces the distance between
/// breaks of size at least minimum_break_duration to be at most distance.
/// This supposes that the time until route start and after route end are
/// infinite breaks.
void SetBreakDistanceDurationOfVehicle(int64_t distance, int64_t duration,
                                       int vehicle);

which reads as this is exactly what i'm looking for in regards to out of the box support for repeated vehicle-breaks.

With some understanding of CP-technology, i would also assume, that this feature of dynamically linking breaks together (with max-duration in-between) is pretty naturally covered by whatever is done during propagation and value-section. Meaning: I would assume out of the box support and this is the closest function i found.

Howewer, some issues appeared for me:

• A: There is no documentation (besides the header-docstring) afaik
• B: There is no example nor test afaik (github code-search finds no caller)
• C: Using it following my interpretation / intuition is not successfull

Let's talk about issue C, my naive approach of using it.

We define a break as:

• at least 540 min of break
• at most 1080 min of driving

A perfect span-minimizing solution with 5 repeated breaks will look like:

Interval occupation	Interval duration	Accumulated total	Accumulated driving
Driving	1080	1080	1080
Break	540	1620	1080
Driving	1080	2700	2160
Break	540	3240	2160
Driving	1080	4320	3240
Break	540	4860	3240
Driving	1080	5940	4320
Break	540	6480	4320
Driving	1080	7560	5400
Break	540	8100	5400
Driving	1080	9180	6480

We can try to model this by mutating the example or-tools/ortools/constraint_solver/samples /vrp_breaks.py

• We simply introduce a single non-depot node and transit to+from costs the same as accumulated driving colum span (3240 * 2)
• We also remove the original service-times
• We then introduce multiple cases:
  • A: Simply 5 non-optional breaks -> feasible + all breaks after depot-return ("modelling-defect")
  • B: Like A, but with SetBreakDistanceDurationOfVehicle -> infeasible -> WHY?
  • C: Like B, but 6>5 breaks available
  • Initial idea: fix first break to be 80 minutes early -> one more break needed
  • But issues of B are blocking this experiment

# [START import]
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
# [END import]

# [START data_model]
def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data["num_vehicles"] = 1
    data["depot"] = 0
    data["time_matrix"] = [
        [0, 3240],
        [3240, 0],
    ]
    data["service_time"] = [0] * len(data["time_matrix"])
    data["service_time"][data["depot"]] = 0
    assert len(data["time_matrix"]) == len(data["service_time"])
    return data
    # [END data_model]

# [START solution_printer]
def print_solution(manager, routing, solution):
    """Prints solution on console."""
    print(f"Objective: {solution.ObjectiveValue()}")

    print("Breaks:")
    intervals = solution.IntervalVarContainer()
    for i in range(intervals.Size()):
        brk = intervals.Element(i)
        if brk.PerformedValue():
            print(
                f"{brk.Var().Name()}: "
                + f"Start({brk.StartValue()}) Duration({brk.DurationValue()})"
            )
        else:
            print(f"{brk.Var().Name()}: Unperformed")

    time_dimension = routing.GetDimensionOrDie("Time")
    total_time = 0
    for vehicle_id in range(manager.GetNumberOfVehicles()):
        index = routing.Start(vehicle_id)
        plan_output = f"Route for vehicle {vehicle_id}:\n"
        while not routing.IsEnd(index):
            time_var = time_dimension.CumulVar(index)
            plan_output += f"{manager.IndexToNode(index)} "
            plan_output += f"Time({solution.Value(time_var)}) -> "
            index = solution.Value(routing.NextVar(index))
        time_var = time_dimension.CumulVar(index)
        plan_output += f"{manager.IndexToNode(index)} "
        plan_output += f"Time({solution.Value(time_var)})\n"
        plan_output += f"Time of the route: {solution.Value(time_var)}min\n"
        print(plan_output)
        total_time += solution.Value(time_var)
    print(f"Total time of all routes: {total_time}min")
    # [END solution_printer]

def main():
    """Solve the VRP with time windows."""
    # 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_matrix"]), data["num_vehicles"], data["depot"]
    )
    # [END index_manager]

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

    # Create and register a transit callback.
    # [START transit_callback]
    def time_callback(from_index, to_index):
        """Returns the travel time + service 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] + data["service_time"][from_node]

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

    # Define cost of each arc.
    # [START arc_cost]
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
    # [END arc_cost]

    # Add Time Windows constraint.
    time = "Time"
    routing.AddDimension(
        transit_callback_index,
        10000,  # needed optional waiting time to place break
        10000,  # maximum time per vehicle
        True,  # Force start cumul to zero.
        time,
    )
    time_dimension = routing.GetDimensionOrDie(time)
    time_dimension.SetGlobalSpanCostCoefficient(10)

    # Breaks
    # [START break_constraint]
    # warning: Need a pre-travel array using the solver's index order.
    node_visit_transit = [0] * routing.Size()
    for index in range(routing.Size()):
        node = manager.IndexToNode(index)
        node_visit_transit[index] = data["service_time"][node]

    eval_case = 'B'                                      # CASE SELECTION

    break_intervals = {}
    for v in range(manager.GetNumberOfVehicles()):

        n_breaks = 5 if eval_case in ('A', 'B') else 6   # CASE DEPENDENCY

        break_intervals[v] = [
            routing.solver().FixedDurationIntervalVar(
                0,  # start min
                10000,  # start max
                540,  
                False,  # optional: no
                f"Break for vehicle {v}",
            ) for i in range(n_breaks)                   # CASE DEPENDENCY
        ]

        time_dimension.SetBreakIntervalsOfVehicle(
            break_intervals[v], v, node_visit_transit)  # breaks  # vehicle index

        # CASE DEPENDENCY
        if eval_case in ['B', 'C']:
            time_dimension.SetBreakDistanceDurationOfVehicle(
              1080, 540, v
            )

    # [END break_constraint]

    # Setting first solution heuristic.
    # [START parameters]
    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(2)
    # [END parameters]

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

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

if __name__ == "__main__":
    main()
# [END program]

Output with case B and or-tools = 9.6.2534

x@y:~/Dev/repeated_breaks$ python3 own_example.py 
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1692622053.891369   93861 search.cc:265] Start search (memory used = 13.86 MB)
I0000 00:00:1692622053.891457   93861 search.cc:265] Root node processed (time = 0 ms, constraints = 31, memory used = 13.89 MB)
I0000 00:00:1692622055.890628   93861 search.cc:265] Finished search tree (time = 1999 ms, branches = 487703, failures = 243878, memory used = 14.01 MB)
I0000 00:00:1692622055.890706   93861 search.cc:265] End search (time = 1999 ms, branches = 487703, failures = 243878, memory used = 14.01 MB, speed = 243973 branches/s)
No solution found !

Questions related to the issue:

• Why is B infeasible?
  • Is my interpretation of breaks wrong?
  • Is my interpretation of SetBreakDistanceDurationOfVehicle wrong?
  • Is this a result of incomplete / non-exhaustive search? (howewer it should be simple enough)
  • Is this a bug?

Additional questions:

• Can we get an example making use of this functionality (which hopefully is the one i interpreted)?

Future:

Ultimately, i want case C to result in the solver finding a span-minimized solution AFTER fixing the start of the first break-start = 1000 (not shown in sources above) < 1080 (break too early compared to perfect sol) which leads to missing out on the perfect solution as shown in the table: we would need 1 more break to cover the transits -> increasing the span.

[sschnug]

I guess my single non-depot node introduces problems due to break-definition. TODO: Adapt the example to have non-depot node cardinality >= breaks.

EDIT

I now use 5 non-depot nodes which should be theory-compatible. I also check this with case A by fixing the breaks leading to the perfect solution table. Still, case B is infeasible which i expected to be feasible!

# [START program]
"""Vehicle Routing Problem (VRP) with breaks.

   This is a sample using the routing library python wrapper to solve a VRP
   problem.
   A description of the problem can be found here:
   http://en.wikipedia.org/wiki/Vehicle_routing_problem.

   Durations are in minutes.
"""

# [START import]
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
# [END import]

# [START data_model]
def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data["num_vehicles"] = 1
    data["depot"] = 0
    data["time_matrix"] = [
        [0, 1080, 1080, 1080, 1080, 1080],
        [1080, 0, 1080, 1080, 1080, 1080],
        [1080, 1080, 0, 1080, 1080, 1080],
        [1080, 1080, 1080, 0, 1080, 1080],
        [1080, 1080, 1080, 1080, 0, 1080],
        [1080, 1080, 1080, 1080, 1080, 0],
    ]
    data["service_time"] = [0] * len(data["time_matrix"])
    data["service_time"][data["depot"]] = 0
    assert len(data["time_matrix"]) == len(data["service_time"])
    return data
    # [END data_model]

# [START solution_printer]
def print_solution(manager, routing, solution):
    """Prints solution on console."""
    print(f"Objective: {solution.ObjectiveValue()}")

    print("Breaks:")
    intervals = solution.IntervalVarContainer()
    for i in range(intervals.Size()):
        brk = intervals.Element(i)
        if brk.PerformedValue():
            print(
                f"{brk.Var().Name()}: "
                + f"Start({brk.StartValue()}) Duration({brk.DurationValue()})"
            )
        else:
            print(f"{brk.Var().Name()}: Unperformed")

    time_dimension = routing.GetDimensionOrDie("Time")
    total_time = 0
    for vehicle_id in range(manager.GetNumberOfVehicles()):
        index = routing.Start(vehicle_id)
        plan_output = f"Route for vehicle {vehicle_id}:\n"
        while not routing.IsEnd(index):
            time_var = time_dimension.CumulVar(index)
            plan_output += f"{manager.IndexToNode(index)} "
            plan_output += f"Time({solution.Value(time_var)}) -> "
            index = solution.Value(routing.NextVar(index))
        time_var = time_dimension.CumulVar(index)
        plan_output += f"{manager.IndexToNode(index)} "
        plan_output += f"Time({solution.Value(time_var)})\n"
        plan_output += f"Time of the route: {solution.Value(time_var)}min\n"
        print(plan_output)
        total_time += solution.Value(time_var)
    print(f"Total time of all routes: {total_time}min")
    # [END solution_printer]

def main():
    """Solve the VRP with time windows."""
    # 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_matrix"]), data["num_vehicles"], data["depot"]
    )
    # [END index_manager]

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

    # Create and register a transit callback.
    # [START transit_callback]
    def time_callback(from_index, to_index):
        """Returns the travel time + service 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] + data["service_time"][from_node]

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

    # Define cost of each arc.
    # [START arc_cost]
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
    # [END arc_cost]

    # Add Time Windows constraint.
    time = "Time"
    routing.AddDimension(
        transit_callback_index,
        10000,  # needed optional waiting time to place break
        10000,  # maximum time per vehicle
        True,  # Force start cumul to zero.
        time,
    )
    time_dimension = routing.GetDimensionOrDie(time)
    time_dimension.SetGlobalSpanCostCoefficient(10)

    # Breaks
    # [START break_constraint]
    # warning: Need a pre-travel array using the solver's index order.
    node_visit_transit = [0] * routing.Size()
    for index in range(routing.Size()):
        node = manager.IndexToNode(index)
        node_visit_transit[index] = data["service_time"][node]

    eval_case = 'A'                                      # CASE SELECTION

    break_intervals = {}
    for v in range(manager.GetNumberOfVehicles()):

        n_breaks = 5 if eval_case in ('A', 'B') else 6   # CASE DEPENDENCY

        break_intervals[v] = [
            routing.solver().FixedDurationIntervalVar(
                0,  # start min
                10000,  # start max
                540,  
                False,  # optional: no
                f"Break for vehicle {v}",
            ) for i in range(n_breaks)                   # CASE DEPENDENCY
        ]

        time_dimension.SetBreakIntervalsOfVehicle(
            break_intervals[v], v, node_visit_transit)  # breaks  # vehicle index

        # CASE DEPENDENCY
        if eval_case in ['B', 'C']:
            time_dimension.SetBreakDistanceDurationOfVehicle(
              1080, 540, v
            )

        # CASE DEPENDENCY
        if eval_case in ['A']:
            for ind, i in enumerate(break_intervals[v]):
              i.SetStartRange(1080 + ind * (540 + 1080), 1080 + ind * (540 + 1080))

    # [END break_constraint]

    # Setting first solution heuristic.
    # [START parameters]
    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(2)
    # [END parameters]

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

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

if __name__ == "__main__":
    main()
# [END program]

Output case A:

Objective: 98280
Breaks:
Break for vehicle 0: Start(1080) Duration(540)
Break for vehicle 0: Start(2700) Duration(540)
Break for vehicle 0: Start(4320) Duration(540)
Break for vehicle 0: Start(5940) Duration(540)
Break for vehicle 0: Start(7560) Duration(540)
Route for vehicle 0:
0 Time(0) -> 5 Time(1080) -> 4 Time(2700) -> 3 Time(4320) -> 2 Time(5940) -> 1 Time(7560) -> 0 Time(9180)
Time of the route: 9180min

Total time of all routes: 9180min

Output case B:

No solution found !

[sschnug]

Next observations:

• A: Taking above and do the (perfect-solution) fixing on case B (with SetBreakDistanceDurationOfVehicle ) leads to a solution

  • Implication: SetBreakDistanceDurationOfVehicle has side-effects / interferes with solution-construction solving

• B: BUT mutating the arguments such that it should not be feasible anymore => not respected -> the feasible solution (from the fixings) is not feasible imho

        if eval_case in ['B', 'C']:
            bigM = 1000
            time_dimension.SetBreakDistanceDurationOfVehicle(
              1080 + bigM, 540, v
            )

leads to (imho invalid) solution:

Objective: 98280
Breaks:
Break for vehicle 0: Start(1080) Duration(540)
Break for vehicle 0: Start(2700) Duration(540)
Break for vehicle 0: Start(4320) Duration(540)
Break for vehicle 0: Start(5940) Duration(540)
Break for vehicle 0: Start(7560) Duration(540)
Route for vehicle 0:
0 Time(0) -> 5 Time(1080) -> 4 Time(2700) -> 3 Time(4320) -> 2 Time(5940) -> 1 Time(7560) -> 0 Time(9180)
Time of the route: 9180min

Total time of all routes: 9180min