VRP example: [SCIP] Exception: SCIP: method cannot be called at this time in solution process!

[leoneifler]

Describe the bug When running the "VRP.py" example, the problem occurs

To Reproduce Run the "VRP.py" example

Expected behavior No error should be triggered

Screenshots

System

• no information, bug reported by other user

Additional context Add any other context about the problem here.

[Joao-Dionisio]

Thank you for the report, @leoneifler! vrp.py is in the unfinished examples folder, and from a quick glance, the problem comes from running line 44 (which tries to add a constraint) after a model is already solved. Something like model.freeTransform() right before line 75 is missing. But doing so leads to an infinite loop, it seems.

Looking at this properly would take some time, it will not be a priority for now, on my end.

[leoneifler]

Yes, I know. I just wanted to open an issue so it does not get lost.

[abb-omidi]

Dear support team,

I am trying to play around with the code and I can see the following points. (Also at least as I could try)

def make_data(n):

  customers = [*range(1, n + 1)]
  V = [depot] + customers
  x = dict([(i,random.randint(10,20)) for i in V])
  y = dict([(i,random.randint(10,20)) for i in V])
  c= {}
  Q = 150

  q = {i: random.randint(50, 60) for i in V}
  q[0] = 0  # depot has no demand

  for i in V:
    for j in V:
      if i != j:
        c[i,j] = round(distance(x[i],y[i],x[j],y[j]))
  return customers, V, x, y, c ,q , Q

depot = 0
n = 7
seed = 0
random.seed(seed)
customers, V, x, y, c, q, Q = make_data(n)

from pyscipopt import *
model = Model("vrp")

x = {}
for i in V:
  for j in V:
    if j > i and i == V[0]: # depot
      x[i,j] = model.addVar(ub=2, vtype="I", name="x(%s,%s)"%(i,j))
    elif j > i:
      x[i,j] = model.addVar(ub=1, vtype="I", name="x(%s,%s)"%(i,j))

model.setObjective(quicksum(c[i,j]*x[i,j] for i in V for j in V if j>i), "minimize")

model.addCons(quicksum(x[V[0],j] for j in V[1:]) == 2*m, "DegreeDepot")

for i in V[1:]:
    model.addCons(quicksum(x[j,i] for j in V if j < i) + quicksum(x[i,j] for j in V if j > i) == 2, "Degree(%s)"%i)

model.redirectOutput()
model.printVersion()
print("--------------------------------------")
model.optimize()
model.printSol()

The second part is to add the capacity cuts of the form ($\sum{i,j \in S} x{ij} \le |S| - \left\lceil \frac{\sum_{c \in S} d_c}{Q} \right\rceil \quad \forall S \subseteq C$).

edges = []
for (i,j) in x:
  if model.getVal(x[i,j]) > 0.9:
    if i != V[0] and j != V[0]:
      edges.append((i,j))

model.freeTransform()

for S in S:
  S_card = len(S)
  q_sum = sum(q[i] for i in S)
  NS = int(math.ceil(float(q_sum)/Q))
  S_edges = [(i,j) for i in S for j in S if i<j and (i,j) in edges]

  if S_card >= 2 and (len(S_edges) >= S_card or NS > 1):
    model.addCons(quicksum(x[i,j] for i in S for j in S if j > i) <= 2-NS)
    print(f"cut: len_{S} <= {len(S)-NS}")

    model.redirectOutput()
    model.printVersion()

    print()
    print("--------------------------------------") 
    model.optimize()

    print()
    print("--------------------------------------")
    model.printSol()

It seems the code works well without throwing any issues and the solver output is as follows:

SCIP version 8.0.4 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver: Soplex 6.0.4] [GitHash: a8e51afd1e]
Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB)

--------------------------------------
  [linear] <c9>: <x(3,5)>[B] (+0) +<x(3,7)>[B] (+0) +<x(5,7)>[B] (+1) <= 0;
;
violation: right hand side is violated by 1
  [linear] <c9>: <x(3,5)>[B] (+0) +<x(3,7)>[B] (+0) +<x(5,7)>[B] (+1) <= 0;
;
violation: right hand side is violated by 1
  [linear] <c9>: <x(3,5)>[B] (+0) +<x(3,7)>[B] (+1) +<x(5,7)>[B] (+1) <= 0;
;
violation: right hand side is violated by 2
1/4 feasible solution given by solution candidate storage, new primal bound 5.300000e+01

presolving:
(round 1, fast)       3 del vars, 1 del conss, 0 add conss, 3 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs
   (0.0s) probing cycle finished: starting next cycle
   (0.0s) symmetry computation started: requiring (bin +, int -, cont +), (fixed: bin -, int +, cont -)
   (0.0s) no symmetry present
   (0.0s) symmetry computation started: requiring (bin +, int +, cont +), (fixed: bin -, int -, cont -)
   (0.0s) no symmetry present
presolving (2 rounds: 2 fast, 1 medium, 1 exhaustive):
 3 deleted vars, 1 deleted constraints, 0 added constraints, 3 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
 36 implications, 0 cliques
presolved problem has 25 variables (18 bin, 7 int, 0 impl, 0 cont) and 8 constraints
      8 constraints of type <linear>
transformed objective value is always integral (scale: 1)
Presolving Time: 0.00
transformed 1/1 original solutions to the transformed problem space

 time | node  | left  |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr|  dualbound   | primalbound  |  gap   | compl. 
p 0.0s|     1 |     0 |     0 |     - | vbounds|   0 |  25 |   9 |   8 |   0 |  0 |   2 |   0 | 0.000000e+00 | 4.300000e+01 |    Inf | unknown
  0.0s|     1 |     0 |     7 |     - |  1172k |   0 |  25 |  16 |   8 |   0 |  0 |   9 |   0 | 4.300000e+01 | 4.300000e+01 |   0.00%| unknown

SCIP Status        : problem is solved [optimal solution found]
Solving Time (sec) : 0.01
Solving Nodes      : 1
Primal Bound       : +4.30000000000000e+01 (2 solutions)
Dual Bound         : +4.30000000000000e+01
Gap                : 0.00 %

Although the solver added the capacity cuts, it seems these kinds of cuts are not sufficient to achive a feasible solution and it does need to add other cuts. As I am not well familiar with VRP types, maybe someone with more experience would help to add some comments. (If I can find more details, I will add).

Please, see and check that. All the best Abbas