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).