Subtour elimination in SDVRP

Question

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[0]
            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[1])
K = int(sys.argv[2])
'''
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 = [0]
demand += [random.randint(1,10) for i in range(n)]

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

volume = [0]
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 = m.addVars(
    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):
            m.addConstr(
                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[0])
                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
                    )
                    print("added lazy cut")


# 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[0]
            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 = [0]
demand += [random.randint(1,10) for i in range(n)]

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

volume = [0]
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] == 0:
        i = "Source"
    else:
        i = k[0]
    if k[1] == 0:
        j = "Sink"
    else:
        j = k[1]
    G.add_edge(i, j, dist=v)

m = gp.Model()

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

x = m.addVars(
    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 = m.addVars(
    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"]:
        m.addConstrs(
            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
m.addConstrs(
    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
m.addConstrs(
    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):
            m.addConstr(
                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}")

Answer 1

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] == 0:
        i = "Source"
    else:
        i = k[0]
    if k[1] == 0:
        j = "Sink"
    else:
        j = k[1]
    G.add_edge(i, j, dist=v)

m = gp.Model()

# Create variables:
x_keys = {(e[0], e[1], k): e[2]["dist"] for e in G.edges(data=True) for k in range(K)}
x = m.addVars(
    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"]:
        m.addConstrs(
            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
m.addConstrs(
    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
m.addConstrs(
    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
z = m.addVars(
    [(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
m.addConstrs(
    (
        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")
m.addConstrs(
    (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}")

Answer 2

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