# Subtour elimination in SDVRP

I have the model below, based on this paper. It's a vehicle routing problem with split deliveries, i.e the locations/customers can be visited by multiple vehicles that share the demand at that vehicle. However, in this case, the "customers" are in fact pick locations and the vehicles are forklifts picking items in a warehouse. The cost is simply the distance traveled, which is to be minimized.

1. $$x_{ijk}$$ is the indicator variable of fork-lift $$k$$ traversing from $$i$$ to $$j$$, if so $$x_{ijk}=1$$, else $$0$$.
2. $$y_{ik}$$ is the positive integer variable denoting the quantity fork-lift $$k$$ picks at location $$i$$.
3. $$c_{ij}$$ is the cost of traversing location $$i$$ to $$j$$.
4. $$d_i$$ is the daily demand to pe picked at location $$i$$.
5. $$w_i$$ is the weight of item $$i$$, placed at location $$i$$ in $$kg$$.
6. $$v_i$$ is the volume of item $$i$$, placed at location $$i$$ in $$m^3$$.
7. $$W$$ is the maximum weight capacity of a fork-lift in $$kg$$.
8. $$V$$ is the maximum volume capacity of a fork-lift in $$m^3$$.
9. $$\mathcal{K}$$ is the set of all fork-lifts.
10. $$\mathcal{N}_0$$ is the set of all pick-locations, including depot.
11. $$\mathcal{N}$$ is the set of all pick-locations, excluding depot.
12. $$\mathcal{S}$$ is the set of all cycles on the set $$\mathcal{N}$$ which include the depot. The point $$0$$ denots the depot.

# SDVRP model

\begin{align} % objective \min_{\boldsymbol{x}, \ (i,j)\in \mathcal{E}, \ k\in\mathcal{K}} \quad & z=\displaystyle\sum_{k\in\mathcal{K}}\displaystyle\sum_{i\in\mathcal{N}_0}\displaystyle\sum_{j\in\mathcal{N}_0}c_{ij}x_{ijk} \\ \nonumber\\ \textrm{subject to} \quad & \displaystyle\sum_{i\in\mathcal{N}_0}\sum_{k\in\mathcal{K}} x_{ijk} \ge 1, & j\in \mathcal{N}_0 \tag{1} \\ &\displaystyle\sum_{i\in\mathcal{N}_0}x_{ipk} - \displaystyle\sum_{j\in\mathcal{N}_0}x_{pjk} = 0, & p\in \mathcal{N}_0,\quad k\in\mathcal{K} \tag{2}\\ &\sum_{i\in\mathcal{S}}\sum_{j\in\mathcal{S}} x_{ijk} \le |\mathcal{S}|-1, & k \in \mathcal{K}, \quad \mathcal{S}\subseteq \mathcal{N}_0\tag{3}\\ &d_i\sum_{j\in\mathcal{N}_0}x_{ijk}\ge y_{ik}, & i \in \mathcal{N}, \quad k\in\mathcal{K}\tag{4}\\ & \displaystyle\sum_{k\in\mathcal{K}} y_{ik} = d_i, & i\in\mathcal{N}\tag{5}\\ & \displaystyle\sum_{i\in\mathcal{N}} y_{ik}w_i \le W, & k\in\mathcal{K}\tag{6}\\ & \displaystyle\sum_{i\in\mathcal{N}} y_{ik}v_i \le V, & k\in\mathcal{K}\tag{7}\\ & x_{ijk} \in \{0,1\}, & i\in \mathcal{N}_0,\quad j\in \mathcal{N}_0, \quad k\in\mathcal{K}\tag{8}\\ & y_{ik} \in \mathbb{Z}^{+}, & i\in\mathcal{N}_0, \quad k\in\mathcal{K}\tag{9} \end{align}

I've implemented this using lazy constraints (callback) where i modified this code in Gurobis website. I tested my program using 3 fork lifts (=K) and 10 pick-locations (=n) but the solution is not what I expect. See end of this post for explanation.

Gurobi and I wonder how one would debug a situation like this when there is no error at all.

Below is a completely reproducible example of my implementation and the output of the variables is shown at the end. I appreciate any help!

Below is a completely reproducible example of my implementation and the output.

# Gurobipy implementation

import sys
import math
import random
from itertools import permutations
import gurobipy as gp
from gurobipy import GRB

# Callback - use lazy constraints to eliminate sub-tours
def subtourelim(model, where):
if where == GRB.Callback.MIPSOL:
# make a list of edges selected in the solution
vals = model.cbGetSolution(model._x)
selected = gp.tuplelist((i,j,k) for i, j, k in model._x.keys() if vals[i, j, k] > 0.5)
# find the shortest cycle in the selected edge list
tour = subtour(selected)
if len(tour) < n:
for k in range(1,K+1): # subtouring for each fork lift
model.cbLazy(gp.quicksum(model._x[i, j, k] for i,j,k in permutations(tour, 3)) <= len(tour)-1)

# Given a tuplelist of edges, find the shortest subtour not containing depot
def subtour(edges):
unvisited = list(range(1, n))
cycle = range(n+1)  # initial length has 1 more location
while unvisited:
thiscycle = []
neighbors = unvisited
while neighbors:
current = neighbors
thiscycle.append(current)
if current != 0:
unvisited.remove(current)
neighbors = [j for i,j,k in edges.select(current,'*','*') if j == 0 or j in unvisited]
if 0 not in thiscycle and len(cycle) > len(thiscycle):
cycle = thiscycle
return cycle

# Parse argument
'''
if len(sys.argv) < 3:
print('Usage: vrp.py npoints nforklifts')
print('{:7}npoints includes depot'.format(""))
sys.exit(1)
n = int(sys.argv)
K = int(sys.argv)
'''
n = 10 # number of items to pick, equivalen to number of locations to visit
K = 3 # number of fork-lifts to use

w_capacity = 1000 # kg
v_capacity = 10000 # dm^3

# Create n random points
points = [(0, 0)]
points += [(random.randint(0, 100), random.randint(0, 100))
for i in range(n-1)]

# Dictionary of Manhattan distance between each pair of points
dist = {(i,j,k):
math.sqrt(sum((points[i][p]-points[j][p])**2 for p in range(2)))
for i in range(n) for j in range(n) for k in range(1,K+1) if i != j}
#random.seed(1)
demand = 
demand += [random.randint(1,10) for i in range(n)]

weight = 
weight += [random.randint(0,100) for i in range(n)]

volume = 
volume += [random.randint(0,1000) for i in range(n)]

print('Number of items ', n)
print('Number of pickers ', K)
print('Total demand: ', sum(demand))
print('Single picker weight capacity: ', w_capacity)
print('Single picker volume capacity: ', v_capacity)
print('Total pick weight capacity: ', K*w_capacity)
print('Total pick volume capacity: ', K*v_capacity)
print('Weights: ', weight)
print('Volumes: ', volume)
print('Demands: ', demand)
print('Total weight demand: ', sum(weight[i]*demand[i] for i in range(len(weight))))
print('Total volume demand: ', sum(volume[i]*demand[i] for i in range(len(volume))))

m = gp.Model('SDVRP')

# Create variables:

# x[i,j,k] = 1, if forklift k visits and goes directly from location i to location j
x = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='x')

# y[i,k] = number of items of type i that forklift k picks
y = m.addVars(gp.tuplelist([(i,k) for i in range(1,n) for k in range(1,K+1)]), vtype=gp.GRB.INTEGER, lb=0, name="y")

# Constraints (1)
m.addConstr(gp.quicksum(x[i,j,k] for i,j,k in dist.keys()) >= 1)

# Constraints (2)
m.addConstrs(gp.quicksum(x[i,p,k] for i in range(n) if i!=p) - gp.quicksum(x[p,j,k] for j in range(n) if j!=p) == 0 for p in range(n) for k in range(1,K+1))

# Constraints (4)
m.addConstrs(demand[i]*gp.quicksum(x[i,j,k] for j in range(n) if i!=j) >= y[i,k] for i in range(1,n) for k in range(1,K+1))

# Constraints (5)
m.addConstrs(gp.quicksum(y[i,k] for k in range(1,K+1)) == demand[i] for i in range(1,n))

# Constraints (6)
m.addConstrs(gp.quicksum(y[i,k]*weight[i] for i in range(1,n)) <= w_capacity for k in range(1,K+1))

# Constraints (7)
m.addConstrs(gp.quicksum(y[i,k]*volume[i] for i in range(1,n)) <= v_capacity for k in range(1,K+1))

# Optimize model

m._x = x
m.Params.LazyConstraints = 1
m.setParam('MIPGap', 0.05) # finish running once 5% gap is reached
m.setParam('Timelimit', 60) # finish running once 1 minute has passed
m.optimize(subtourelim)


# Output

Set parameter LazyConstraints to value 1
Set parameter MIPGap to value 0.05
Set parameter TimeLimit to value 60
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 73 rows, 297 columns and 1161 nonzeros
Model fingerprint: 0x49013b49
Variable types: 0 continuous, 297 integer (270 binary)
Coefficient statistics:
Matrix range     [1e+00, 1e+03]
Objective range  [1e+01, 1e+02]
Bounds range     [1e+00, 1e+00]
RHS range        [1e+00, 1e+04]
Presolve time: 0.00s
Presolved: 73 rows, 297 columns, 1161 nonzeros
Variable types: 0 continuous, 297 integer (276 binary)
Found heuristic solution: objective 3588.1169986
Found heuristic solution: objective 3426.2294702
Found heuristic solution: objective 1246.6673887
Found heuristic solution: objective 647.2200397

Root relaxation: objective 2.256603e+02, 100 iterations, 0.00 seconds (0.00 work units)

Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
...
Optimal solution found (tolerance 5.00e-02)
Best objective 2.896915417605e+02, best bound 2.801074130851e+02, gap 3.3084%

User-callback calls 878, time in user-callback 0.01 sec


# Variable result

    Variable            X
-------------------------
x[0,6,3]            1
x[1,9,1]            1
x[1,9,2]            1
x[2,7,2]            1
x[3,8,3]            1
x[4,5,2]            1
x[5,4,2]            1
x[6,0,3]            1
x[7,2,2]            1
x[8,3,3]            1
x[9,1,1]            1
x[9,1,2]            1
y[1,1]            5
y[2,2]            3
y[3,3]            1
y[4,2]            9
y[5,2]            1
y[6,3]           10
y[7,2]            4
y[8,3]           10
y[9,1]            7
y[9,2]            1


The y seems fine since they all correctly sum up to the demand for each location. But the x variables building the tour are way off. Firstly, only forklift 3 starts at zero and none of the forklifts return to zero, despite the flow constraints. Secondly forklift 3 goes 0 -> 6 ->0 and then also 8 -> 3, this tells me that there must be something wrong with the subtour elimination.

# EDIT - Following David Torres answer.

To declare the positive integer variable y I did the following:

y_keys = {(i, k) for i in list(G.predecessors("Sink")) for k in range(K)}

y_keys,
vtype=GRB.INTEGER,
name='y'
)


To implement constraints (4),(5),(6) and (7) in the mathematical model I did the following:

# (4)

for i in list(G.predecessors("Sink")):
for k in range(K):
y[i,k] <= demand[i]*gp.quicksum(x[i,j,k] for j in range(1,n) if j!=i)
)

# (5)

m.addConstrs(gp.quicksum(y[i,k] for k in range(K)) == demand[i] for i in list(G.predecessors("Sink")))

# (6) & (7) Note that instead of using weights and volumes, I just use demand constraint Q = 20, i.e. each forklift can only carry 20 items.

m.addConstrs(gp.quicksum(y[i,k] for i in list(G.predecessors("Sink"))) <= Q for k in range(K))


When I run this, the subtours are not eliminate. This is what I get:

k=0
Source->2
Source->9
1->2
1->6
2->Sink
2->5
5->1
6->5
6->9
9->Sink
k=1
Source->9
1->2
2->7
3->1
3->6
4->3
6->8
7->4
8->7
9->Sink
k=2
Source->9
2->5
5->1
6->Sink
8->6
9->8


All of the trucks go from source to sink immediately, and then there are still subtours.

# Full new implementation below:

import math
import random

import gurobipy as gp
from gurobipy import GRB
import networkx as nx
import matplotlib.pyplot as plt

# Callback - use lazy constraints to eliminate sub-tours
def subtourelim(model, where):
if where == GRB.Callback.MIPSOL:
# make a list of edges selected in the solution
vals = model.cbGetSolution(model._x)
for k in range(K):
selected = gp.tuplelist(
(i, j, k) for i, j in G.edges() if vals[i, j, k] > 0.0
)
tours = allsubtours(selected)
for tour in tours:
tourlen = len(tour)
tour.append(tour)
if tourlen < n:
model.cbLazy(
gp.quicksum(
model._x[i, j, k]
for i, j in zip(tour, tour[1:])
if (i, j) in G.edges()
)
<= tourlen - 1
)

# Given a tuplelist of edges, find all subtours
def allsubtours(edges):
unvisited = list(range(n))
cycles = []
while unvisited:
thiscycle = []
neighbors = unvisited
while neighbors:
current = neighbors
thiscycle.append(current)
unvisited.remove(current)
neighbors = [
j for i, j, _ in edges.select(current, "*", "*") if j in unvisited
]
if len(thiscycle) > 1:
cycles.append(thiscycle)
print(f"{cycles=}")
return cycles

n = 10  # number of items to pick, equivalent to number of locations to visit
K = 3  # number of fork-lifts to use

# Create n random points
points = [(0, 0)]
points += [(random.randint(0, 100), random.randint(0, 100)) for i in range(n - 1)]

# Dictionary of Manhattan distance between each pair of points
dist = {
(i, j): math.sqrt(sum((points[i][p] - points[j][p]) ** 2 for p in range(2)))
for i in range(n)
for j in range(n)
if i != j
}

Q = 20
random.seed(1)
demand = 
demand += [random.randint(1,10) for i in range(n)]

'''
weight = 
weight += [random.randint(0,100) for i in range(n)]

volume = 
volume += [random.randint(0,1000) for i in range(n)]

print('Number of items ', n)
print('Number of pickers ', K)
print('Total demand: ', sum(demand))
print('Single picker weight capacity: ', w_capacity)
print('Single picker volume capacity: ', v_capacity)
print('Total pick weight capacity: ', K*w_capacity)
print('Total pick volume capacity: ', K*v_capacity)
print('Weights: ', weight)
print('Volumes: ', volume)
print('Demands: ', demand)
print('Total weight demand: ', sum(weight[i]*demand[i] for i in range(len(weight))))
print('Total volume demand: ', sum(volume[i]*demand[i] for i in range(len(volume))))'''

# Create graph
G = nx.DiGraph()
for k, v in dist.items():
if k == 0:
i = "Source"
else:
i = k
if k == 0:
j = "Sink"
else:
j = k

m = gp.Model()

# Create variables:
x_keys = {(e, e, k): e["dist"] for e in G.edges(data=True) for k in range(K)}

x_keys,
obj=x_keys,
vtype=GRB.BINARY,
name="x",
)

y_keys = {(i, k) for i in list(G.predecessors("Sink")) for k in range(K)}

y_keys,
vtype=GRB.INTEGER,
name='y'
)

# Visit all nodes
for j in G.nodes():
if j not in ["Sink"]:
pred = list(G.predecessors(j))
if len(pred) > 0:
m.addConstr(gp.quicksum(x[i, j, k] for i in pred for k in range(K)) >= 1)

# Flow-balance
for v in G.nodes():
if v not in ["Source", "Sink"]:
gp.quicksum(x[i, v, k] for i in G.predecessors(v))
- gp.quicksum(x[v, j, k] for j in G.successors(v))
== 0
for k in range(K)
)

# All k's must start at Source
gp.quicksum(x["Source", j, k] for j in G.successors("Source")) == 1
for k in range(K)
)

# All k's must end at Sink
gp.quicksum(x[i, "Sink", k] for i in G.predecessors("Sink")) == 1 for k in range(K)
)

# (4)
for i in list(G.predecessors("Sink")):
for k in range(K):
y[i,k] <= demand[i]*gp.quicksum(x[i,j,k] for j in range(1,n) if j!=i)
)

# (5)

m.addConstrs(gp.quicksum(y[i,k] for k in range(K)) == demand[i] for i in list(G.predecessors("Sink")))

# (6) & (7) Note that instead of using weights and volumes, I just use demand constraint Q = 20, i.e. each forklift can only carry 20 items.

m.addConstrs(gp.quicksum(y[i,k] for i in list(G.predecessors("Sink"))) <= Q for k in range(K))

m._x = x
m.Params.LazyConstraints = 1
m.optimize(subtourelim)

for k in range(K):
print(f"{k=}")
for i, j in G.edges():
if x[i, j, k].X > 0.0:
print(f"\t{i}->{j}")


The subtour elimination constraints are not correct, they are mostly empty. Try to debug the callback calls, as you can see also from the log, none of lazy constraints are applied.

You can:

• use callbacks to enumerate and forbid all cycles for a given forklift. This is what the last comment in the thread where you got the subtour elimination code suggested. This is very expensive for large graphs.
• Use MTZ constraints.

For more in-depth study of the SDVRP, please see Dror et al. (1994).

Here it is using MTZ constraints (disregarding the y constraints): code modified from Kuifje02/vrpy/:vrpy/subproblem_lp.py#L227:

import math
import random

import gurobipy as gp
from gurobipy import GRB
import networkx as nx
import matplotlib.pyplot as plt

n = 100  # number of items to pick, equivalent to number of locations to visit
K = 3  # number of fork-lifts to use

# Create n random points
points = [(0, 0)]
points += [(random.randint(0, 100), random.randint(0, 100)) for i in range(n - 1)]

# Dictionary of Manhattan distance between each pair of points
dist = {
(i, j): math.sqrt(sum((points[i][p] - points[j][p]) ** 2 for p in range(2)))
for i in range(n)
for j in range(n)
if i != j
}

# Create graph
G = nx.DiGraph()
for k, v in dist.items():
if k == 0:
i = "Source"
else:
i = k
if k == 0:
j = "Sink"
else:
j = k

m = gp.Model()

# Create variables:
x_keys = {(e, e, k): e["dist"] for e in G.edges(data=True) for k in range(K)}
x_keys,
obj=x_keys,
vtype=GRB.BINARY,
name="x",
)

# Visit all nodes
for j in G.nodes():
if j not in ["Source", "Sink"]:
pred = list(G.predecessors(j))
m.addConstr(gp.quicksum(x[i, j, k] for i in pred for k in range(K)) >= 1)

# Flow-balance
for v in G.nodes():
if v not in ["Source", "Sink"]:
gp.quicksum(x[i, v, k] for i in G.predecessors(v))
- gp.quicksum(x[v, j, k] for j in G.successors(v))
== 0
for k in range(K)
)

# All k's must start at Source
gp.quicksum(x["Source", j, k] for j in G.successors("Source")) == 1
for k in range(K)
)

# All k's must end at Sink
gp.quicksum(x[i, "Sink", k] for i in G.predecessors("Sink")) == 1 for k in range(K)
)

# MTZ
M = len(G.nodes())

# Add rank varibles
[(v, k) for v in G.nodes() for k in range(K)],
name="z",
lb=0,
ub=len(G.nodes()),
vtype=GRB.INTEGER,
)
# Add big-M constraints
(
z[i, k] + 1 <= z[j, k] + M * (1 - x[i, j, k])
for k in range(K)
for (i, j) in G.edges()
),
name="elementary",
)

# Source is first, Sink is last (optional)
m.addConstrs((z["Source", k] == 0 for k in range(K)), name="Source_is_first")
(z[v, k] <= z["Sink", k] for k in range(K) for v in G.nodes() if v != "Sink"),
name="Sink_after_%s" % v,
)

m.params.timelimit = 60
m.optimize()

for k in range(K):
print(f"{k=}")
for i, j in G.edges():
if x[i, j, k].X > 0.0:
print(f"\t{i}->{j}")

• Hi David - No worries, I appreciate your answer. Just been somewhat frustrated lately because I've been stuck on this issue for quite some time now and nothing I do seems to work. I have a few questions regarding your implementation. (1) When you print lazy cut, what is the 'expr' variable? Or how do I actually obtain it? (2) in constraint 1, each node can be visited many times but at least once, so it should be >= 1. (3) I'm having trouble declaring variable y and implementing constraint 4 in my mathematical model. I'll update the question shortly and show what I've done. See main question. Aug 4 at 13:46
• (1) Sorry forgot to delete that from my draft, updated now. It was to print out the LHS of the constraint to see what variables were actually in the constraint. (2) That is true, that breaks in that case. Aug 4 at 18:03
• Try the updated version where we actually enumerate the cycles, with maximum of max_no_cycles for each k. Aug 4 at 19:01
• Hmm, it works ok for n = 10...13 and K = 3, but when I go n= 14+ it breaks. I get for example for forklift k=1 Source->3 ; 3->Sink ; 5->6 ; 5->13 ; 6->5 ; 13->5 Aug 4 at 19:26
• Have you tried increasing the value of max_no_cycles? Otherwise it is just a heuristic. I have updated my answer with MTZ constraints from VRPy. Please read up on subtour elemination constraints. Aug 5 at 10:12

Your implementations of constraints (1), (6) and (7) do not match the algebraic formulation. (You sum over different variables than what the formulation indicates.) Does fixing that solve the problem?

Addendum: After modification, the new solution uses only variables $$x_{i,j,k}$$ where $$i=j$$ (and consequently $$c_{i,j}=0.$$ Without knowing the full details of the model, I suspect that you want either to constrain $$x_{i,i,k}=0$$ for all $$i$$ or change the index of summation in (1) (and possibly elsewhere) to skip the case $$i=j.$$

• Hello! Yes I fixed that so that it coincides with the mathematical formulation, which changed the output slightly but the objective value is still zero. Also, the variables of x that the model chose are only (7,7,1), (8,8,2) and (9,9,2). I've updated the question to correct the variables. Aug 2 at 17:36
• You're right - I added if j!=i on all the constraints that required it. I also edited the post to include that. However now, it just runs forever, despite very few locations and picker. Is there a way to see a log during execution or limit the execution time? Aug 2 at 19:48
• I'm not a Gurobi user, but I am confident that it (like every other solver) has a parameter that you can use to set a time limit. As far as seeing the log, you printed a log in the output section of your question, so I don't understand what you mean.
– prubin
Aug 2 at 21:07
• its like me hahah, i did vrp with time windows with just 10 nodes and it runs forever. you can use param timelimit to limit the running time in gurobi Aug 3 at 14:02
• @prubin - I've now done some modifications to the code and updated it accordingly. The program now finishes very fast with 0 gap, however the result I get is not what I want. Aug 3 at 17:57