Prefered driver penalty

[Isaksjoberg]

I'm working on a VRPTW for the home care services.
Each visit has a set of prefered caretakers (drivers).
I want to implement a penalty if the solver choses a driver that is not prefered for that visit. I have constructed a code below and OR-tools gives no error;

    penalty_pref = 36000
    penalty_vars = []

    for node, preferred_vehicles in data["pref_driver"].items():
        if not preferred_vehicles:
            continue  

        index = manager.NodeToIndex(node)
        vehicle_var = routing.VehicleVar(index)

        is_preferred_bools = [
            routing.solver().IsEqualCstVar(vehicle_var, v)
            for v in preferred_vehicles
        ]

        is_preferred = routing.solver().Max(is_preferred_bools)
        is_not_preferred = routing.solver().IntVar(0, 1, f'is_not_preferred_{node}')
        routing.solver().Add(is_not_preferred + is_preferred == 1)

        penalty_var = routing.solver().IntVar(0, penalty_pref, f'penalty_{node}')
        routing.solver().Add(penalty_var == is_not_preferred * (penalty_pref-1))
        penalty_vars.append(penalty_var)

    penalty_sum = routing.solver().Sum(penalty_vars)
    penalty_total = routing.solver().IntVar(0, 3600000000, "penalty_total")
    routing.solver().Add(penalty_total == penalty_sum)
    routing.AddVariableMinimizedByFinalizer(penalty_total)

This was my only finalizer before adding the prefered drivers, which produces good results. It looks like this;

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

The model is not affected no matter what penalty i give in the prefered driver section. It seems like the model doesn't take this variable into account. What am I doing wrong?

Thanks in advance

[Mizux]

Few questions:
why using two variables (1 IntExpr and 1 IntVar) here

 penalty_sum = routing.solver().Sum(penalty_vars)
 penalty_total = routing.solver().IntVar(0, 3600000000, "penalty_total")
 routing.solver().Add(penalty_total == penalty_sum)
 routing.AddVariableMinimizedByFinalizer(penalty_total)|

You should be able to use:

penalty_sum = routing.solver().Sum(penalty_vars).Var()
routing.AddVariableMinimizedByFinalizer(penalty_sum)

ref:

or-tools/ortools/constraint_solver/constraint_solveri.h

Line 115 in cae6cd6

IntVar* Var() override;

After Finalizer is kind of post processing mostyl use when there is a SlackVar e.g. waiting time to force the solver to minimize/maximize the value of the cumulvar instead of returning a range (kind of meta solution of all possibilities entangled)
see: #1051 (comment)

In your case you want to have the prefered vehicle cost embeded in the objective function not in post processing
Still ned to think how to model it, from top of my head (i.e. not tested) I would try:

1. Create a capacity dimension, forced to zero with
2. An unary transit function per vehicle returning:
   0: if vehicle prefered for visiting the node
   1: if not a prefered vehicle
   note: in python you can use partial to have a function with vehicle index and wrap it multiple time...
   see: https://github.com/google/or-tools/blob/main/ortools/routing/samples/cvrptw_break.py
3. Add a GlobalSpanCost of this dimension to be 36000 following your example...

Now the objective function should have an extra cost equal to number_of_not_prefered_vehicle_visit x 36000

[Isaksjoberg]

Thanks for the fast reply!

I tried creating another dimension but ended up with a new problem. My code looks like this;

def capacity_callback(data):
        vehicles = data["pref_driver"]
        def capacity_evaluator(manager, node):
            indexes = manager.NodeToIndex(node)
            vehicle_vares = routing.VehicleVar(indexes)
            if not vehicles[indexes]:
                pass
            else:
                if vehicle_vares in vehicles[indexes]:
                    return 0
                else:
                    return 1
        return capacity_evaluator
capacity_callback_index = routing.RegisterUnaryTransitCallback(partial(capacity_callback(data), manager))

capacity = "Capacity"
routing.AddDimension(
    capacity_callback_index,
    0,
    100000000,
    True,
    capacity,
    )
capacity_dimension = routing.GetDimensionOrDie(capacity)
capacity_dimension.SetGlobalSpanCostCoefficient(360)

The problem I'm running into is that all visits are assigned to vehichle 0 even though the visits overlap. No matter what SpanCostCoefficient I set it ends up the same. The rest of the model works as intended when I remove the capacity dimension. Clearly I did something wrong here.

This is a sample of how the data["pref_driver"] looks: {0: [], 1: [], 2: [], 3: [0], 4: [19], 5: [18], 6: [14], 7: [18], 8: [21, 6], 9: [14], 10: [8], 11: [15], 12: [7, 6], 13: [21, 12], 14: [0], 15: [21, 12]}

[Mizux]

few problems here imho:

• callback use "Index" that you need to convert to "Node" (you do the contrary)
• instead of pass you should return 0 otherwise it is UB
• didn't know you can probe the VehicleVar from the callback itself -> I need to check
• you should also pass routing
• SetGlobalSpanCostCoefficient() will only take into consideration the route with the most unprefered vehicle visit...
  while you want the sum of all CumulVar to be taken into account -> prefer to use SetSpanCostCoefficientForAllVehicles(int64_t coefficient)

ref:

or-tools/ortools/routing/routing.h

Lines 3323 to 3330 in 1c5daab

	/// Sets a cost proportional to the dimension span on a given vehicle,
	/// or on all vehicles at once. "coefficient" must be nonnegative.
	/// This is handy to model costs proportional to idle time when the dimension
	/// represents time.
	/// The cost for a vehicle is
	/// span_cost = coefficient * (dimension end value - dimension start value).
	void SetSpanCostCoefficientForVehicle(int64_t coefficient, int vehicle);
	void SetSpanCostCoefficientForAllVehicles(int64_t coefficient);

blind fix:

def capacity_callback(data):
        vehicles = data["pref_driver"]
        def capacity_evaluator(manager, routing, index):
            node = manager.IndexToNode(index)
            vehicle_var = routing.VehicleVar(index)
            if not vehicles[node]:
                return 0 # no pref -> no penalty
            else:
                if vehicle_var in vehicles[node]:
                    return 0 # prefered vehicle used -> no penalty
                else:
                    return 1 # non prefered vehicle used -> penalty
        return capacity_evaluator
capacity_callback_index = routing.RegisterUnaryTransitCallback(partial(capacity_callback(data), manager, routing))

[Mizux]

As expected vehicle_var = routing.VehicleVar(index) only return an IntVar object so if vehicle_var in vehicles[node]: won't work and always return true....

On my way to fix, testing on a pet project....

[Mizux]

(based on vrp_capacity.py....)

vrp_capacity_with_prefered_vehicle.py:

#!/usr/bin/env python3
# [START program]
"""Capacited Vehicles Routing Problem (CVRP) With Prefered Vehicles."""

# [START import]
from functools import partial
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["distance_matrix"] = [
        # fmt: off
      [0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354, 468, 776, 662],
      [548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674, 1016, 868, 1210],
      [776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164, 1130, 788, 1552, 754],
      [696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822, 1164, 560, 1358],
      [582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708, 1050, 674, 1244],
      [274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628, 514, 1050, 708],
      [502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856, 514, 1278, 480],
      [194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320, 662, 742, 856],
      [308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662, 320, 1084, 514],
      [194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388, 274, 810, 468],
      [536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764, 730, 388, 1152, 354],
      [502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114, 308, 650, 274, 844],
      [388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194, 536, 388, 730],
      [354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0, 342, 422, 536],
      [468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536, 342, 0, 764, 194],
      [776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274, 388, 422, 764, 0, 798],
      [662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730, 536, 194, 798, 0],
        # fmt: on
    ]
    # [START demands_capacities]
    data["demands"] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    data["vehicle_capacities"] = [15, 15, 15, 15]
    #########################
    #########################
    data["prefered_vehicle"] = {
        0: [],
        1: [2, 3],
        2: [],
        3: [0],
        4: [0],
        5: [0],
        6: [0],
        7: [0],
        8: [0, 1],
        9: [1],
        10: [2],
        11: [0],
        12: [0, 4],
        13: [0, 2],
        14: [0],
        15: [0, 2],
        16: [],
    }
    #########################
    #########################
    # [END demands_capacities]
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data
    # [END data_model]


# [START solution_printer]
def print_solution(manager, routing, solution):
    """Prints solution on console."""
    print(f"Objective: {solution.ObjectiveValue()}")
    total_distance = 0
    total_load = 0
    load_dimension = routing.GetDimensionOrDie("Load")
    pref_dimension = routing.GetDimensionOrDie("Prefered")
    for vehicle_id in range(manager.GetNumberOfVehicles()):
        if not routing.IsVehicleUsed(solution, vehicle_id):
            continue
        index = routing.Start(vehicle_id)
        plan_output = f"Route for vehicle {vehicle_id}:\n"
        route_distance = 0
        while not routing.IsEnd(index):
            load_var = load_dimension.CumulVar(index)
            load_value = solution.Value(load_var)
            pref_var = pref_dimension.CumulVar(index)
            pref_value = solution.Value(pref_var)
            node_index = manager.IndexToNode(index)
            plan_output += f" {node_index} Load({load_value}) Pref({pref_value}) -> "
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            route_distance += routing.GetArcCostForVehicle(
                previous_index, index, vehicle_id
            )
        load_var = load_dimension.CumulVar(index)
        load_value = solution.Value(load_var)
        pref_var = pref_dimension.CumulVar(index)
        pref_value = solution.Value(pref_var)
        plan_output += (
            f" {manager.IndexToNode(index)} Load({load_value}) Pref({pref_value})\n"
        )
        route_distance -= pref_value * 100_000
        plan_output += f"Distance of the route: {route_distance}m\n"
        plan_output += f"Load of the route: {load_value}\n"
        print(plan_output)
        total_distance += route_distance
        total_load += load_value
    print(f"Total distance of all routes: {total_distance}m")
    print(f"Total load of all routes: {total_load}")
    # [END solution_printer]


def main():
    """Solve the CVRP problem."""
    # 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["distance_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 distance_callback(from_index, to_index):
        """Returns the distance between the two nodes."""
        # Convert from routing variable Index to distance matrix NodeIndex.
        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)
    # [END transit_callback]

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

    # Add Capacity constraint.
    # [START capacity_constraint]
    def demand_callback(from_index):
        """Returns the demand of the node."""
        # Convert from routing variable Index to demands NodeIndex.
        from_node = manager.IndexToNode(from_index)
        return data["demands"][from_node]

    demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)
    routing.AddDimensionWithVehicleCapacity(
        demand_callback_index,
        0,  # null capacity slack
        data["vehicle_capacities"],  # vehicle maximum capacities
        True,  # start cumul to zero
        "Load",
    )
    # [END capacity_constraint]

    #########################
    #########################
    # Prefered Vehicles
    def prefered_callback(data):
        vehicles = data["prefered_vehicle"]

        def prefered_evaluator(manager, vehicle_id, index):
            node = manager.IndexToNode(index)
            if not vehicles[node]:
                print(f"v:{vehicle_id}, node:{node} -> pass")
                return 0  # no pref -> no penalty
            if vehicle_id in vehicles[node]:
                print(f"v:{vehicle_id}, node:{node} -> prefered")
                return 0  # prefered vehicle used -> no penalty
            else:
                print(f"v:{vehicle_id}, node:{node} -> not prefered")
                return 1  # non prefered vehicle used -> penalty

        return prefered_evaluator

    prefered_callback_indices = []
    for vehicle_id in range(manager.GetNumberOfVehicles()):
        prefered_callback_indices.append(
            routing.RegisterUnaryTransitCallback(
                partial(prefered_callback(data), manager, vehicle_id)
            )
        )

    name = "Prefered"
    routing.AddDimensionWithVehicleTransits(
        prefered_callback_indices,
        0,  # no slack
        100_000_000,
        True,  # start forced to zero
        name,
    )
    prefered_dimension = routing.GetDimensionOrDie(name)
    prefered_dimension.SetSpanCostCoefficientForAllVehicles(100_000)
    #########################
    #########################

    # 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.time_limit.FromSeconds(1)
    # [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)
    # [END print_solution]


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

[Mizux]

potential output:

$ ./vrp_capacity_with_prefered_vehicle.py
...
Objective: 408628
v:0, node:0 -> pass
v:0, node:11 -> prefered
v:0, node:4 -> prefered
v:0, node:3 -> prefered
v:0, node:7 -> prefered
Route for vehicle 0:
 0 Load(0) Pref(0) ->  11 Load(0) Pref(0) ->  4 Load(1) Pref(0) ->  3 Load(5) Pref(0) ->  7 Load(7) Pref(0) ->  0 Load(15) Pref(0)
Distance of the route: 1712m
Load of the route: 15

v:1, node:0 -> pass
v:1, node:8 -> prefered
v:1, node:6 -> not prefered
v:1, node:5 -> not prefered
v:1, node:9 -> prefered
Route for vehicle 1:
 0 Load(0) Pref(0) ->  8 Load(0) Pref(0) ->  6 Load(8) Pref(0) ->  5 Load(12) Pref(1) ->  9 Load(14) Pref(2) ->  0 Load(15) Pref(2)
Distance of the route: 1164m
Load of the route: 15

v:2, node:0 -> pass
v:2, node:10 -> prefered
v:2, node:13 -> prefered
v:2, node:15 -> prefered
v:2, node:1 -> prefered
Route for vehicle 2:
 0 Load(0) Pref(0) ->  10 Load(0) Pref(0) ->  13 Load(2) Pref(0) ->  15 Load(6) Pref(0) ->  1 Load(14) Pref(0) ->  0 Load(15) Pref(0)
Distance of the route: 3104m
Load of the route: 15

v:3, node:0 -> pass
v:3, node:2 -> pass
v:3, node:16 -> pass
v:3, node:14 -> not prefered
v:3, node:12 -> not prefered
Route for vehicle 3:
 0 Load(0) Pref(0) ->  2 Load(0) Pref(0) ->  16 Load(1) Pref(0) ->  14 Load(9) Pref(0) ->  12 Load(13) Pref(1) ->  0 Load(15) Pref(2)
Distance of the route: 2648m
Load of the route: 15

Total distance of all routes: 8628m
Total load of all routes: 60

[Isaksjoberg]

Thank you so much for the quick and extensive reply! I've now implemented your code and got the prefered vehicle printing correctly but unfortunately I'm running into the same issue of all visits being assigned to one vehicle. I've tried to do some debugging and I think I've narrowed it down to how I define my time dimension which I use as the ArcCostEvaluator. I define my time dimension like this:

def time_callback(manager, 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)
        travel_time = data["time_matrix"][from_node][to_node]
        service_time = data["service_times"].get(from_node, 0)
        change_time = 300


        return travel_time + service_time + change_time

transit_callback_index = routing.RegisterTransitCallback(partial(time_callback, manager))

routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
time = "Time"
    routing.AddDimension(
        transit_callback_index,
        300000,  # allow waiting time
        86400,  # maximum time per vehicle
        False,  # Don't force start cumul to zero.
        time,
    )
time_dimension = routing.GetDimensionOrDie(time)

Because workers have different schedules it's important to define the time dimension using partial. But the issue seems to lie in that it is defined as a partial and some clash between the "Time" dimension and the "prefered" dimension occurs. If I change the definition of the distance dimension to partial in your code, all distances are set to 0. This mirrors the behaviour of my time dimension not working when i include the prefered dimension. Time dimension and prefered dimension are the only dimensions in my program. Do you have any idea why this occurs and how I could solve it?

Thanks again for all your time and help.

[Mizux]

as usual:

• verify your data["time_matrix"] and data["service_times"] always return int values and no floating point values.
• having a 300000, # allow waiting time aka Max SlackVar greater than the 86400 (CumulVar max value) bug me.

Because workers have different schedules it's important to define the time dimension using partial.

I don't see in your snippet code a "worker" (vehicle?) specific -> don't understand your partial justification

if I change the definition of the distance dimension to partial in your code, all distances are set to 0.

Do you have a diff to share, so I can reproduce it ? not sure to fully understand

EDIT:
Reproduce the "0" distance using this modification

--- vrp_capacity_with_prefered_vehicle.py.bckp	2025-04-24 15:52:32.748021068 +0200
+++ vrp_capacity_with_prefered_vehicle.py	2025-04-24 15:50:48.388530417 +0200
@@ -140,7 +140,8 @@
         to_node = manager.IndexToNode(to_index)
         return data["distance_matrix"][from_node][to_node]
 
-    transit_callback_index = routing.RegisterTransitCallback(distance_callback)
+    transit_callback_index = routing.RegisterTransitCallback(
+            partial(distance_callback, manager))
     # [END transit_callback]
 
     # Define cost of each arc.

Fixed by using:

--- vrp_capacity_with_prefered_vehicle.py.bckp	2025-04-24 15:52:32.748021068 +0200
+++ vrp_capacity_with_prefered_vehicle.py	2025-04-24 15:57:32.534557703 +0200
@@ -133,14 +133,15 @@
 
     # Create and register a transit callback.
     # [START transit_callback]
-    def distance_callback(from_index, to_index):
+    def distance_callback(manager, from_index, to_index):
         """Returns the distance between the two nodes."""
         # Convert from routing variable Index to distance matrix NodeIndex.
         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)
+    transit_callback_index = routing.RegisterTransitCallback(
+            partial(distance_callback, manager))
     # [END transit_callback]
 
     # Define cost of each arc.

Please notice the change in the distance_callback prototype having manager as first arg...
note: using partial(distance_callback, manager=manager) do not work too

[Isaksjoberg]

I've verified that the matrix and data returns integer values. Changed the Max SlackVar to 30000. My code works when not including the prefered dimension. This is the code snippet for why I need to use partial definition for the time dimension:

for vehicle_id in range(data["num_vehicles"]):
        start_index = routing.Start(vehicle_id)
        end_index = routing.End(vehicle_id)
        start_time, end_time = data["vehicle_time_windows"][vehicle_id]
        time_dimension.CumulVar(start_index).SetMin(start_time)
        time_dimension.CumulVar(end_index).SetMin(start_time)
        time_dimension.SetCumulVarSoftUpperBound(start_index, end_time, 10)
        time_dimension.SetCumulVarSoftUpperBound(end_index, end_time, 10)

Each worker starts and ends their shift at different times. Your EDIT mirrors what I did with your program to make it not work in the same way that my program does. Would you like me to include full script and some visual results before and after adding the prefered dimension?