AddDisjunction function does not work like it supposed

[teemusanteri]

Hi,

Adddisjunction does not work like it is documented.

See the forum topic:

https://groups.google.com/forum/m/?utm_medium=email&utm_source=footer#!msg/or-tools-discuss/703V05No26s/oQbpsPUCAgAJ

[jmarca]

To repeat the forum posts here a little bit.

The disjunctions by themselves aren't the problem. The problem comes when there is some sort of binary flip flop between one state and another. Here we want to choose one of two vehicles. The idea is a truck by itself, or a truck with a trailer---obviously you can't drive both.

This is accomplished by the lines:

    for veh_pair in range(0, num_veh//2):
        combo = veh_pair*2
        single = veh_pair*2 + 1
        index_combo = routing.End(combo)
        index_single = routing.End(single)

        end_time_combo = time_dimension.CumulVar(index_combo)
        end_time_single = time_dimension.CumulVar(index_single)

        combo_on = end_time_combo > 0
        single_on = end_time_single > 0
        solver.Add(combo_on * single_on == 0)

This works, but only without disjunctions. So if the problem has more destinations than can be served by the vehicles, you need disjunctions to allow dropping of nodes. However, once you turn on disjunctions, the solver fails to find the correct solution.

Attached code---disjunction_fail.zip---demonstrates the problem. See also my repo at https://github.com/jmarca/disjunction_tests.

Running without disjunctions (the default) and with just 5 destinations, the solver does the right thing---uses two combo trucks, one with three visits, the other with two. Logging turned on to show the solver's progress.

python disjunction_fail.py    -l 1
WARNING: Logging before InitGoogleLogging() is written to STDERR
I0611 19:59:20.584077     8 search.cc:252] Start search (memory used = 42.73 MB)
I0611 19:59:20.584228     8 search.cc:252] Root node processed (time = 0 ms, constraints = 110, memory used = 42.73 MB)
I0611 19:59:20.584597     8 search.cc:252] Solution #0 (75, time = 0 ms, branches = 52, failures = 8, depth = 33, memory used = 42.73 MB)
I0611 19:59:20.611366     8 search.cc:252] Finished search tree (time = 27 ms, branches = 2599, failures = 1518, neighbors = 338, filtered neighbors = 133, accepted neigbors = 0, memory used = 42.73 MB)
I0611 19:59:20.611438     8 search.cc:252] End search (time = 27 ms, branches = 2599, failures = 1518, memory used = 42.73 MB, speed = 96259 branches/s)
The Objective Value is 75
Truck 0 travel time
     combo: 40 
     single: 0
Truck 1 travel time
     combo: 35 
     single: 0
(etc)

Run with disjunctions turned on, again with just 5 demand nodes, the solver fails. Instead of assigning all five as before, it decides to drop one node, no matter how large the penalty is.

 python disjunction_fail.py  -d 1 --singlepenalty 500  -l 1
single node disjunction penalty is 500
added 5 disjunctions, one per node
WARNING: Logging before InitGoogleLogging() is written to STDERR
I0611 20:01:10.471397     9 routing.cc:2833] All Unperformed Solution (2500, time = 0 ms, memory used = 42.73 MB)
I0611 20:01:10.471446     9 search.cc:252] Start search (memory used = 42.73 MB)
I0611 20:01:10.471575     9 search.cc:252] Root node processed (time = 0 ms, constraints = 125, memory used = 42.73 MB)
I0611 20:01:10.471832     9 search.cc:252] Solution #0 (2030, time = 0 ms, branches = 38, failures = 0, depth = 33, memory used = 42.77 MB)
I0611 20:01:10.472105     9 search.cc:252] Solution #1 (2006, objective maximum = 2030, time = 0 ms, branches = 41, failures = 2, depth = 33, Relocate<1>, neighbors = 1, filtered neighbors = 1, accepted neighbors = 1, memory used = 42.77 MB)
I0611 20:01:10.472709     9 search.cc:252] Solution #2 (1536, objective maximum = 2030, time = 1 ms, branches = 45, failures = 5, depth = 33, MakeActiveOperator, neighbors = 18, filtered neighbors = 3, accepted neighbors = 2, memory used = 42.77 MB)
I0611 20:01:10.472937     9 search.cc:252] Solution #3 (1041, objective maximum = 2030, time = 1 ms, branches = 48, failures = 7, depth = 33, MakeActiveOperator, neighbors = 19, filtered neighbors = 4, accepted neighbors = 3, memory used = 42.77 MB)
I0611 20:01:10.473183     9 search.cc:252] Solution #4 (546, objective maximum = 2030, time = 1 ms, branches = 53, failures = 9, depth = 33, MakeActiveOperator, neighbors = 20, filtered neighbors = 5, accepted neighbors = 4, memory used = 42.77 MB)
I0611 20:01:10.487522     9 search.cc:252] Finished search tree (time = 16 ms, branches = 827, failures = 622, neighbors = 300, filtered neighbors = 147, accepted neigbors = 4, memory used = 42.77 MB)
I0611 20:01:10.487587     9 search.cc:252] End search (time = 16 ms, branches = 827, failures = 622, memory used = 42.77 MB, speed = 51687 branches/s)
The Objective Value is 546
Truck 0 travel time
     combo: 0 
     single: 6
Truck 1 travel time
     combo: 40 
     single: 0
Route for vehicle 0:
 0 Load(0) Cost(0) ->  0 Load(0) Cost(0)
Distance of the route: 0m
Load of the route: 0
(etc)

[jmarca]

Any thoughts or comments on this? Is there something more I can do to clarify the issue?

[jmarca]

I can't speak for the original poster's bug, but my confirmation of it seems to be fixed.

My possible solution is that I correctly added the "guided_local_search" flag, and it seems to be working now.

Updated the python in my repo: https://github.com/jmarca/disjunction_tests/blob/master/disjunction_fail.py

I have not tested all of the failing cases I found, but so far things seem to work.

The fix:

Previously I was using

    search_parameters.first_solution_strategy = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)

Which is nonsense, but passes for some reason (I think because GUIDED_LOCAL_SEARCH is just 2, which translates to a solution strategy that makes sense but certainly not the local search metaheuristic of guided local search...

So, fixing that,

    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.GLOBAL_CHEAPEST_ARC) 
    search_parameters.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
    search_parameters.time_limit.seconds =  60  # timelimit

And then the proper solution pops out

[kevindalmeijer]

Would it be possible to reopen this issue to update the documentation? I agree with @teemusanteri that the behavior of AddDisjunction() does not match the documentation when len(indices) > 1. The Python documentation and the C++ documentation both state:

This is equivalent to adding the constraint: p + Sum(i)active[i] == max_cardinality where p is an integer variable, and the following cost to the cost function: p * penalty.

The code below sets up a small CVRP with four customers and one vehicle of capacity 2. When I add a disjunction on [1, 2, 3, 4] with penalty=1 and max_capacity=4. I get the route 0 Load(0) -> 4 Load(1) -> 3 Load(2) -> 0 Load(2) with penalty 1. However, based on the documentation I expected to see a penalty of 2, because max_cardinality - Sum(i)active[i] = 2.

To be clear, I do not suggest to change this behavior, as the functionality is great. However, it would be helpful to update the documentation to avoid confusion.

Example code tested on ortools 9.10.4067 (simplified version of the CVRP example):

"""Capacited Vehicles Routing Problem (CVRP)."""

from ortools.constraint_solver import pywrapcp

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data["distance_matrix"] = [
        # fmt: off
      [0, 548, 776, 696, 582],
      [548, 0, 684, 308, 194],
      [776, 684, 0, 992, 878],
      [696, 308, 992, 0, 114],
      [582, 194, 878, 114, 0],
        # fmt: on
    ]
    data["demands"] = [0, 1, 1, 1, 1]
    data["vehicle_capacities"] = [2]
    data["num_vehicles"] = 1
    data["depot"] = 0
    return data

def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    print(f"Objective: {solution.ObjectiveValue()}")
    for vehicle_id in range(data["num_vehicles"]):
        index = routing.Start(vehicle_id)
        plan_output = f"Route for vehicle {vehicle_id}:\n"
        route_distance = 0
        route_load = 0
        while not routing.IsEnd(index):
            node_index = manager.IndexToNode(index)
            route_load += data["demands"][node_index]
            plan_output += f" {node_index} Load({route_load}) -> "
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            route_distance += routing.GetArcCostForVehicle(
                previous_index, index, vehicle_id
            )
        plan_output += f" {manager.IndexToNode(index)} Load({route_load})\n"
        print(plan_output)

def main():
    """Solve the CVRP problem."""
    # Instantiate the data problem.
    data = create_data_model()

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

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

    # Add 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
        "Capacity",
    )

    ### RELEVANT CODE ###
    # add disjunction for all nodes with penalty 1
    penalty = 1
    max_cardinality = 4
    routing.AddDisjunction([manager.NodeToIndex(node) for node in range(1, 5)], penalty, max_cardinality)
    ### RELEVANT CODE ###

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.time_limit.FromSeconds(1)

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

    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)
        print("Total penalty: {}".format(solution.ObjectiveValue()))

if __name__ == "__main__":
    main()

[jmarca]

I agree the docs seem wrong.

From the C++ I see:

  IntVar* no_active_var = solver_->MakeBoolVar();
  IntVar* number_active_vars = solver_->MakeIntVar(0, max_cardinality);
  solver_->AddConstraint(
      solver_->MakeSumEquality(disjunction_vars, number_active_vars));
  solver_->AddConstraint(solver_->MakeIsDifferentCstCt(
      number_active_vars, max_cardinality, no_active_var));
  const int64_t penalty = disjunctions_[disjunction].value.penalty;
  if (penalty < 0) {
    no_active_var->SetMax(0);
    return nullptr;
  } else {
    return solver_->MakeProd(no_active_var, penalty)->Var();
  }

To me that looks like no_active_var is a 0/1 boolean that flags if the max cardinality is matched or not. So the penalty is multiplied by 1 or 0.

Regards,

James

On Tue, May 21, 2024 at 07:39:49PM -0700, Kevin Dalmeijer wrote:

Would it be possible to reopen this issue to update the documentation? I agree with @teemusanteri that the behavior of AddDisjunction() does not match the documentation when len(indices) > 1. The Python documentation and the C++ documentation both state:

This is equivalent to adding the constraint: p + Sum(i)active[i] == max_cardinality where p is an integer variable, and the following cost to the cost function: p * penalty.

The code below sets up a small CVRP with four customers and one vehicle of capacity 2. When I add a disjunction on [1, 2, 3, 4] with penalty=1 and max_capacity=4. I get the route 0 Load(0) -> 4 Load(1) -> 3 Load(2) -> 0 Load(2) with penalty 1. However, based on the documentation I expected to see a penalty of 2, because max_cardinality - Sum(i)active[i] = 2.

To be clear, I do not suggest to change this behavior, as the functionality is great. However, it would be helpful to update the documentation to avoid confusion.

Example code tested on ortools 9.10.4067 (simplified version of the CVRP example):

"""Capacited Vehicles Routing Problem (CVRP)."""

from ortools.constraint_solver import pywrapcp

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data["distance_matrix"] = [
        # fmt: off
      [0, 548, 776, 696, 582],
      [548, 0, 684, 308, 194],
      [776, 684, 0, 992, 878],
      [696, 308, 992, 0, 114],
      [582, 194, 878, 114, 0],
        # fmt: on
    ]
    data["demands"] = [0, 1, 1, 1, 1]
    data["vehicle_capacities"] = [2]
    data["num_vehicles"] = 1
    data["depot"] = 0
    return data

def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    print(f"Objective: {solution.ObjectiveValue()}")
    for vehicle_id in range(data["num_vehicles"]):
        index = routing.Start(vehicle_id)
        plan_output = f"Route for vehicle {vehicle_id}:\n"
        route_distance = 0
        route_load = 0
        while not routing.IsEnd(index):
            node_index = manager.IndexToNode(index)
            route_load += data["demands"][node_index]
            plan_output += f" {node_index} Load({route_load}) -> "
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            route_distance += routing.GetArcCostForVehicle(
                previous_index, index, vehicle_id
            )
        plan_output += f" {manager.IndexToNode(index)} Load({route_load})\n"
        print(plan_output)

def main():
    """Solve the CVRP problem."""
    # Instantiate the data problem.
    data = create_data_model()

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

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

    # Add 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
        "Capacity",
    )

    ### RELEVANT CODE ###
    # add disjunction for all nodes with penalty 1
    penalty = 1
    max_cardinality = 4
    routing.AddDisjunction([manager.NodeToIndex(node) for node in range(1, 5)], penalty, max_cardinality)
    ### RELEVANT CODE ###

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.time_limit.FromSeconds(1)

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

    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)
        print("Total penalty: {}".format(solution.ObjectiveValue()))

if __name__ == "__main__":
    main()

-- Reply to this email directly or view it on GitHub: https://github.com/google/or-tools/issues/1348#issuecomment-2123763486 You are receiving this because you commented.

Message ID: @.***>

--

James E. Marca Activimetrics LLC

[Mizux]

source of the ref doc: https://github.com/google/or-tools/blob/c72e8bbb05d299adf19c199af9718141e7a357d9/ortools/constraint_solver/routing.h#L861-L876

and just to get syntax coloring of jmarca post (not working on email reply)...

https://github.com/google/or-tools/blob/c72e8bbb05d299adf19c199af9718141e7a357d9/ortools/constraint_solver/routing.cc#L1757-L1771

[jmarca]

(I never noticed email replies do not support markdown. I'll keep that in mind in the future)