No feasible solution is found on a 5x5 matrix when I prohibit certain edges.

[lampretl]

Details of my setup:
Version of ORtools: 9.11.4210
Language: Python
Solver: Routing
Operating system: Linux Mint 20

Background: I am solving a VRP where there are several locations, for each location I have several assignment (pickup and delivery) nodes and 2 park nodes (one for pickups, one for deliveries). For each vehicle, before visiting any assignment node, I want to force it to first visit the corresponding park node (and hence spend 10min parking, which is done only once at a location). To enforce this, I tried 2 ways: hard constraint of prohibiting certain edges, or soft constraint of adding a large penalty to the objective function for those edges.

All assignment nodes have the option to not be visited, by adding a penalty to the objective function (AddDisjunction). If all assignment nodes at a location are unvisited, it makes no sense to travell all the way to that location just to visit the park node, so I also allow the park nodes to be omitted with 0 penalty.

In the simplest case below, we have a 5x5 matrix, where nodes are 0=depot, 1=pickup, 2=delivery, 3=park_for_pickup, 4=park_for_delivery. Hence a feasible solution is 0 -> 3 -> 1 -> 4 -> 2 -> 0.

Issue: Of the 13 initial solution strategies and 5 metaheuristics, none finds this simple solution. There are only 4! = 24 possible paths starting and ending in 0. On top of that, in special cases, I get an invalid solution, since constraint solver.Add(model.NextVar(_i) != _j) is violated. Self-contained code that demonstrates this:

import numpy as np
from ortools.constraint_solver.pywrapcp import RoutingIndexManager, RoutingModel, DefaultRoutingSearchParameters
from ortools.constraint_solver.routing_enums_pb2 import FirstSolutionStrategy, LocalSearchMetaheuristic

def get_data():
   cost = np.array([
      [  0,   0, 191,   0, 191],
      [  1,   0, 192,   1, 192],
      [192, 192,   0, 192,   0],
      [  2,   2, 193,   0, 193],
      [193, 193,   2, 193,   0]])
   cheat_edge = np.array([
      [0, 1, 1, 0, 0],
      [0, 0, 1, 1, 0],
      [0, 1, 0, 0, 1],
      [1, 0, 1, 0, 1],
      [1, 1, 0, 1, 0]])
   return cost, cheat_edge

def get_paths(solution, model, _node):
   if not solution:
      return None
   paths = []
   for v in range(model.vehicles()):  # iterate over all available vehicles
      _i = model.Start(v)  # _i = index of node where vehicle starts its route
      paths.append([int(_node(_i))])
      while not model.IsEnd(_i):
         _i = solution.Value(model.NextVar(_i))
         paths[v].append(int(_node(_i)))
   return paths

def solve_vrp(strategy:int, heuristic:int, forbidden_edge_penalty=np.inf, omit_node_penalty=[10**4, 0]):
   cost, cheat_edge = get_data()
   n_nodes, n_vehicles, depot_node = cost.shape[0], 1, 0
   manager = RoutingIndexManager(n_nodes, n_vehicles, depot_node)
   model = RoutingModel(manager)
   solver, _node, _idx = model.solver(), manager.IndexToNode, manager.NodeToIndex

   path = [0,3,1,4,2,0]  # this solution is feasible, but not found by the algorithm
   assert all(not cheat_edge[i,j] for i,j in zip(path[:-1],path[1:]))  # verify feasibility

   if np.isinf(forbidden_edge_penalty):  # hard constraint
      for i in range(n_nodes):
         for j in range(n_nodes):
            if cheat_edge[i,j]: 
               _i, _j = _idx(i), _idx(j)
               solver.Add(model.NextVar(_i) != _j)  # forbid certain edges to be used
   elif isinstance(forbidden_edge_penalty, int):  # soft constraint
      cost += cheat_edge * forbidden_edge_penalty
   
   for v in range(n_vehicles):
      cb_cost = model.RegisterTransitCallback(lambda _i, _j: cost[_node(_i), _node(_j)])
      model.SetArcCostEvaluatorOfVehicle(cb_cost, v)  # this creates the objective function

   _i, _j = _idx(1), _idx(2)
   model.AddPickupAndDelivery(_i, _j)  # vehicle must move i->j on same path, not edge

   if not np.isinf(omit_node_penalty[0]):  # allow certain nodes to not be visited, but with a penalty added to the objective function
      # AddDisjunction([i1,...,in], penalty, k) => at most k of the n nodes are used, for l<k we pay (k-l)*penalty
      model.AddDisjunction([_idx(1),_idx(2)], omit_node_penalty[0], 2)  # due to pickup-delivery constraint, 0 or 2 are used
   for i in range(3,5): 
      model.AddDisjunction([_idx(i)], omit_node_penalty[1])  # park nodes 

   search_params = DefaultRoutingSearchParameters()
   search_params.time_limit.FromSeconds(5)
   if strategy == 3: 
      model.SetFirstSolutionEvaluator(lambda _i,_j: 1)
   search_params.first_solution_strategy = [
      FirstSolutionStrategy.AUTOMATIC, 
      FirstSolutionStrategy.PATH_CHEAPEST_ARC,
      FirstSolutionStrategy.PATH_MOST_CONSTRAINED_ARC,
      FirstSolutionStrategy.EVALUATOR_STRATEGY,
      FirstSolutionStrategy.SAVINGS,
      FirstSolutionStrategy.SWEEP,
      FirstSolutionStrategy.CHRISTOFIDES,
      FirstSolutionStrategy.ALL_UNPERFORMED,
      FirstSolutionStrategy.BEST_INSERTION,
      FirstSolutionStrategy.PARALLEL_CHEAPEST_INSERTION,
      FirstSolutionStrategy.LOCAL_CHEAPEST_INSERTION,
      FirstSolutionStrategy.GLOBAL_CHEAPEST_ARC,
      FirstSolutionStrategy.LOCAL_CHEAPEST_ARC,
      FirstSolutionStrategy.FIRST_UNBOUND_MIN_VALUE][strategy]
   search_params.local_search_metaheuristic = [
      LocalSearchMetaheuristic.AUTOMATIC,
      LocalSearchMetaheuristic.GREEDY_DESCENT,
      LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH,
      LocalSearchMetaheuristic.SIMULATED_ANNEALING,
      LocalSearchMetaheuristic.TABU_SEARCH,
      LocalSearchMetaheuristic.GENERIC_TABU_SEARCH][heuristic]
   solution = model.SolveWithParameters(search_params)

   paths = get_paths(solution, model, _node)
   if paths is None: 
      print("No solution found.")
   elif all(len(path)<=2 for path in paths):
      print("No vehicle is used.")
   else: 
      print(f"Solution/paths for s={strategy} and h={heuristic}:",*paths, sep="\n")

for s in range(13):
   for h in range(6):
      solve_vrp(strategy=s, heuristic=h, forbidden_edge_penalty=np.inf, omit_node_penalty=[10_000, 0])

All outputs are "No solution found" or "No vehicle is used". However, if I run

solve_vrp(strategy=11, heuristic=2, forbidden_edge_penalty=np.inf, omit_node_penalty=[10_000, 1000])

i.e. also penalize omitting a park node, then I get a solution [0, 3, 0] which is not even feasible (edge 3->0 is not allowed). On the other hand, if I use a soft constraint

solve_vrp(strategy=11, heuristic=2, forbidden_edge_penalty=1000, omit_node_penalty=[10_000, 0])

I get the desired solution [0, 3, 1, 4, 2, 0]. But in my actual case, where I have many nodes, vehicles and constraints, using a soft constraint for forbidden edges results in many violations (vehicles avoid park nodes and go directly to assignment nodes).

[lperron]

you cannot use _idx on a depot.

[lampretl]

Thank you for a reply! Hmm, was not aware of that. So instead of _idx(i) I should use _idx(i) if i>0 else 0 everywhere where i could be 0? Why is that so? Experimentally, _idx(0) always returned 0.

In any case, I tried the above code with _idx(i) if i>0 else 0 and still no solution. Am I using ORtools wrong or could this be a bug?

[lperron]

The depot is duplicated when num vehicles > 1. So all start and end of each vehicles are duplicated. Le dim. 26 janv. 2025, 21:31, Leon Lampret ***@***.***> a écrit :

…

Hmm, was not aware of that. So instead of _idx(i) I should use _idx(i) if i>0 else 0 everywhere where i could be 0? Why is that so? Experimentally, _idx(0) always returned 0. In any case, I tried the above code with _idx(i) if i>0 else 0 and still no solution. Am I using ORtools wrong or could this be a bug? — Reply to this email directly, view it on GitHub <#4511 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACUPL3JU6KCVEMJ4INY374T2MVAZFAVCNFSM6AAAAABV4VWTFCVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCOJWGE2TCMY> . You are receiving this because you commented.Message ID: ***@***.***>

[Mizux]

i.e. as Laurent Say you must use routing.Start(vehicle_idx), routing.End(vehicle_idx) to get vehicle's start and end index respectively...

[lampretl]

I have edited my MWE that demonstrates the issue, to include your suggestions and also the option of specifying the initial solution.

import numpy as np
from ortools.constraint_solver.pywrapcp import RoutingIndexManager, RoutingModel, DefaultRoutingSearchParameters
from ortools.constraint_solver.routing_enums_pb2 import FirstSolutionStrategy, LocalSearchMetaheuristic

def get_data():
   cost = np.array([
      [  0,   0, 191,   0, 191],
      [  1,   0, 192,   1, 192],
      [192, 192,   0, 192,   0],
      [  2,   2, 193,   0, 193],
      [193, 193,   2, 193,   0]])
   cheat_edge = np.array([
      [0, 1, 1, 0, 0],
      [0, 0, 1, 1, 0],
      [0, 1, 0, 0, 1],
      [1, 0, 1, 0, 1],
      [1, 1, 0, 1, 0]])
   return cost, cheat_edge

def get_paths(solution, model, _node):
   if not solution:
      return None
   paths = []
   for v in range(model.vehicles()):  # iterate over all available vehicles
      _i = model.Start(v)  # _i = index of node where vehicle starts its route
      paths.append([int(_node(_i))])
      while not model.IsEnd(_i):
         _i = solution.Value(model.NextVar(_i))
         paths[v].append(int(_node(_i)))
   return paths

def solve_vrp(strategy:int, heuristic:int, forbidden_edge_penalty=np.inf, omit_node_penalty=[10**4, 0], use_initial_sol=None):
   cost, cheat_edge = get_data()
   n_nodes, n_vehicles, depot_node = cost.shape[0], 1, 0
   manager = RoutingIndexManager(n_nodes, n_vehicles, depot_node)
   model = RoutingModel(manager)
   solver, _node, _idx = model.solver(), manager.IndexToNode, manager.NodeToIndex

   path = [0,3,1,4,2,0]  # this solution is feasible, but not found by the algorithm
   assert all(not cheat_edge[i,j] for i,j in zip(path[:-1],path[1:]))  # verify feasibility

   if np.isinf(forbidden_edge_penalty):  # hard constraint
      for i in range(1,n_nodes):
         for j in range(1,n_nodes):
            if cheat_edge[i,j]:
               solver.Add(model.NextVar(_idx(i)) != _idx(j))  # forbid certain edges to be used
      for v in range(n_vehicles):
         start_idx, end_idx = model.Start(v), model.End(v)
         for j in range(1,n_nodes):
            if cheat_edge[0,j]:
               solver.Add(model.NextVar(start_idx) != _idx(j))
         for i in range(1,n_nodes):
            if cheat_edge[i,0]:
               solver.Add(model.NextVar(_idx(i)) != end_idx)
   elif isinstance(forbidden_edge_penalty, int):  # soft constraint
      cost += cheat_edge * forbidden_edge_penalty
   
   for v in range(n_vehicles):
      cb_cost = model.RegisterTransitCallback(lambda _i, _j: cost[_node(_i), _node(_j)])
      model.SetArcCostEvaluatorOfVehicle(cb_cost, v)  # this creates the objective function

   model.AddDimensionWithVehicleTransitAndCapacity(  # counting dimension used in pickup-delivery constraint
      [model.RegisterTransitCallback(lambda _i, _j: 1) for _ in range(n_vehicles)],
      0, # max slack
      n_vehicles*[10**9], # upper bound
      True, # start at 0
      "VehicleDistance")

   dim = model.GetDimensionOrDie("VehicleDistance")
   _i, _j = _idx(1), _idx(2)
   model.AddPickupAndDelivery(_i, _j)  # vehicle must move i->j on same path, not edge
   solver.Add(model.VehicleVar(_i) == model.VehicleVar(_j))  # same vehicle is at i,j
   solver.Add(dim.CumulVar(_i) <= dim.CumulVar(_j))  # vehicle moves in direction i->j

   if not np.isinf(omit_node_penalty[0]):  # allow certain nodes to not be visited, but with a penalty added to the objective function
      # AddDisjunction([i1,...,in], penalty, k) => at most k of the n nodes are used, for l<k we pay (k-l)*penalty
      model.AddDisjunction([_idx(1),_idx(2)], omit_node_penalty[0], 2)  # due to pickup-delivery constraint, 0 or 2 are used
   for i in range(3,5): 
      model.AddDisjunction([_idx(i)], omit_node_penalty[1])  # park nodes 

   search_params = DefaultRoutingSearchParameters()
   search_params.time_limit.FromSeconds(5)
   if strategy == 3: 
      model.SetFirstSolutionEvaluator(lambda _i,_j: 1)
   search_params.first_solution_strategy = [
      FirstSolutionStrategy.AUTOMATIC, 
      FirstSolutionStrategy.PATH_CHEAPEST_ARC,
      FirstSolutionStrategy.PATH_MOST_CONSTRAINED_ARC,
      FirstSolutionStrategy.EVALUATOR_STRATEGY,
      FirstSolutionStrategy.SAVINGS,
      FirstSolutionStrategy.SWEEP,
      FirstSolutionStrategy.CHRISTOFIDES,
      FirstSolutionStrategy.ALL_UNPERFORMED,
      FirstSolutionStrategy.BEST_INSERTION,
      FirstSolutionStrategy.PARALLEL_CHEAPEST_INSERTION,
      FirstSolutionStrategy.LOCAL_CHEAPEST_INSERTION,
      FirstSolutionStrategy.GLOBAL_CHEAPEST_ARC,
      FirstSolutionStrategy.LOCAL_CHEAPEST_ARC,
      FirstSolutionStrategy.FIRST_UNBOUND_MIN_VALUE][strategy]
   search_params.local_search_metaheuristic = [
      LocalSearchMetaheuristic.AUTOMATIC,
      LocalSearchMetaheuristic.GREEDY_DESCENT,
      LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH,
      LocalSearchMetaheuristic.SIMULATED_ANNEALING,
      LocalSearchMetaheuristic.TABU_SEARCH,
      LocalSearchMetaheuristic.GENERIC_TABU_SEARCH][heuristic]
   if use_initial_sol:
      model.CloseModelWithParameters(search_params)
      sol = model.ReadAssignmentFromRoutes([path[1:-1]], True)
      solution = model.SolveFromAssignmentWithParameters(sol, search_params)
   else:
      solution = model.SolveWithParameters(search_params)

   paths = get_paths(solution, model, _node)
   if paths is None: 
      print("No solution found.", end="   ")
   elif all(len(path)<=2 for path in paths):
      print("No vehicle is used.", end="   ")
   else: 
      print(f"Solution/paths for s={strategy} and h={heuristic}:",*paths)

for s in range(14): 
   for h in range(6):
      solve_vrp(strategy=s, heuristic=h, forbidden_edge_penalty=np.inf, omit_node_penalty=[10_000, 0], use_initial_sol=False)
   print()

Now, a solution is found only for 3 out of 14 strategies: GLOBAL_CHEAPEST_ARC, LOCAL_CHEAPEST_ARC, FIRST_UNBOUND_MIN_VALUE. If, however, I specify use_initial_sol=True, then all 14*6 variations return a solution.

Seems it's still not working as it should. Any info or suggestion would be appreciated.

[lperron]

So working as intended. All initial solution heuristics are incomplete. Le lun. 27 janv. 2025, 20:46, Leon Lampret ***@***.***> a écrit :

…

I have edited my MWE that demonstrates the issue, to include your suggestions and also the option of specifying the initial solution. import numpy as np from ortools.constraint_solver.pywrapcp import RoutingIndexManager, RoutingModel, DefaultRoutingSearchParameters from ortools.constraint_solver.routing_enums_pb2 import FirstSolutionStrategy, LocalSearchMetaheuristic def get_data(): cost = np.array([ [ 0, 0, 191, 0, 191], [ 1, 0, 192, 1, 192], [192, 192, 0, 192, 0], [ 2, 2, 193, 0, 193], [193, 193, 2, 193, 0]]) cheat_edge = np.array([ [0, 1, 1, 0, 0], [0, 0, 1, 1, 0], [0, 1, 0, 0, 1], [1, 0, 1, 0, 1], [1, 1, 0, 1, 0]]) return cost, cheat_edge def get_paths(solution, model, _node): if not solution: return None paths = [] for v in range(model.vehicles()): # iterate over all available vehicles _i = model.Start(v) # _i = index of node where vehicle starts its route paths.append([int(_node(_i))]) while not model.IsEnd(_i): _i = solution.Value(model.NextVar(_i)) paths[v].append(int(_node(_i))) return paths def solve_vrp(strategy:int, heuristic:int, forbidden_edge_penalty=np.inf, omit_node_penalty=[10**4, 0], use_initial_sol=None): cost, cheat_edge = get_data() n_nodes, n_vehicles, depot_node = cost.shape[0], 1, 0 manager = RoutingIndexManager(n_nodes, n_vehicles, depot_node) model = RoutingModel(manager) solver, _node, _idx = model.solver(), manager.IndexToNode, manager.NodeToIndex path = [0,3,1,4,2,0] # this solution is feasible, but not found by the algorithm assert all(not cheat_edge[i,j] for i,j in zip(path[:-1],path[1:])) # verify feasibility if np.isinf(forbidden_edge_penalty): # hard constraint for i in range(1,n_nodes): for j in range(1,n_nodes): if cheat_edge[i,j]: solver.Add(model.NextVar(_idx(i)) != _idx(j)) # forbid certain edges to be used for v in range(n_vehicles): start_idx, end_idx = model.Start(v), model.End(v) for j in range(1,n_nodes): if cheat_edge[0,j]: solver.Add(model.NextVar(start_idx) != _idx(j)) for i in range(1,n_nodes): if cheat_edge[i,0]: solver.Add(model.NextVar(_idx(i)) != end_idx) elif isinstance(forbidden_edge_penalty, int): # soft constraint cost += cheat_edge * forbidden_edge_penalty for v in range(n_vehicles): cb_cost = model.RegisterTransitCallback(lambda _i, _j: cost[_node(_i), _node(_j)]) model.SetArcCostEvaluatorOfVehicle(cb_cost, v) # this creates the objective function model.AddDimensionWithVehicleTransitAndCapacity( # counting dimension used in pickup-delivery constraint [model.RegisterTransitCallback(lambda _i, _j: 1) for _ in range(n_vehicles)], 0, # max slack n_vehicles*[10**9], # upper bound True, # start at 0 "VehicleDistance") dim = model.GetDimensionOrDie("VehicleDistance") _i, _j = _idx(1), _idx(2) model.AddPickupAndDelivery(_i, _j) # vehicle must move i->j on same path, not edge solver.Add(model.VehicleVar(_i) == model.VehicleVar(_j)) # same vehicle is at i,j solver.Add(dim.CumulVar(_i) <= dim.CumulVar(_j)) # vehicle moves in direction i->j if not np.isinf(omit_node_penalty[0]): # allow certain nodes to not be visited, but with a penalty added to the objective function # AddDisjunction([i1,...,in], penalty, k) => at most k of the n nodes are used, for l<k we pay (k-l)*penalty model.AddDisjunction([_idx(1),_idx(2)], omit_node_penalty[0], 2) # due to pickup-delivery constraint, 0 or 2 are used for i in range(3,5): model.AddDisjunction([_idx(i)], omit_node_penalty[1]) # park nodes search_params = DefaultRoutingSearchParameters() search_params.time_limit.FromSeconds(5) if strategy == 3: model.SetFirstSolutionEvaluator(lambda _i,_j: 1) search_params.first_solution_strategy = [ FirstSolutionStrategy.AUTOMATIC, FirstSolutionStrategy.PATH_CHEAPEST_ARC, FirstSolutionStrategy.PATH_MOST_CONSTRAINED_ARC, FirstSolutionStrategy.EVALUATOR_STRATEGY, FirstSolutionStrategy.SAVINGS, FirstSolutionStrategy.SWEEP, FirstSolutionStrategy.CHRISTOFIDES, FirstSolutionStrategy.ALL_UNPERFORMED, FirstSolutionStrategy.BEST_INSERTION, FirstSolutionStrategy.PARALLEL_CHEAPEST_INSERTION, FirstSolutionStrategy.LOCAL_CHEAPEST_INSERTION, FirstSolutionStrategy.GLOBAL_CHEAPEST_ARC, FirstSolutionStrategy.LOCAL_CHEAPEST_ARC, FirstSolutionStrategy.FIRST_UNBOUND_MIN_VALUE][strategy] search_params.local_search_metaheuristic = [ LocalSearchMetaheuristic.AUTOMATIC, LocalSearchMetaheuristic.GREEDY_DESCENT, LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH, LocalSearchMetaheuristic.SIMULATED_ANNEALING, LocalSearchMetaheuristic.TABU_SEARCH, LocalSearchMetaheuristic.GENERIC_TABU_SEARCH][heuristic] if use_initial_sol: model.CloseModelWithParameters(search_params) sol = model.ReadAssignmentFromRoutes([path[1:-1]], True) solution = model.SolveFromAssignmentWithParameters(sol, search_params) else: solution = model.SolveWithParameters(search_params) paths = get_paths(solution, model, _node) if paths is None: print("No solution found.", end=" ") elif all(len(path)<=2 for path in paths): print("No vehicle is used.", end=" ") else: print(f"Solution/paths for s={strategy} and h={heuristic}:",*paths) for s in range(14): for h in range(6): solve_vrp(strategy=s, heuristic=h, forbidden_edge_penalty=np.inf, omit_node_penalty=[10_000, 0], use_initial_sol=False) print() Now, a solution is found only for 3 out of 14 strategies: GLOBAL_CHEAPEST_ARC, LOCAL_CHEAPEST_ARC, FIRST_UNBOUND_MIN_VALUE. If, however, I specify use_initial_sol=True, then all 14*6 variations return a solution. Seems it's still not working as it should. Any info or suggestion would be appreciated. — Reply to this email directly, view it on GitHub <#4511 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACUPL3PUHMEJCGJVNH3UGOD2M2EILAVCNFSM6AAAAABV4VWTFCVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCOJXGQ4DAOI> . You are receiving this because you commented.Message ID: ***@***.***>