# Subtour elimination implementation (DFJ) using Python (PULP). While loop never exits

The simple test problem I'm trying to implement is

\begin{align} \min &\quad c_{ij}x_{ij} \\ \text{s.t} &\quad \\ &\quad \sum_{j\in N}x_{ij} = 1, \quad i\in N\\ &\quad \sum_{i\in N}x_{ij} = 1, \quad j\in N\\ &\quad x_{ii} = 0, \quad i\in N\\ &\quad \sum_{i\in S}\sum_{j\in S, \ j\neq i} x_{ij} \leq |S|-1, \quad \forall S \subset N, 2\leq |S| \leq n-1 \\ & \quad x_{ij}\in\{0,1\}, \quad i,j\in{N} \end{align}

where $$N={1,...,n}$$ is number of locations and $$S$$ is the set of sub-tours. I have the following locations with their coordinates and the corresponding distance matrix for each pair of locations named locations_df and dist_mat respectively.

I've followed this article (github-link) and I managed to correctly implement the MTZ version, however I'm running into troubles when trying to implement the DFJ method of sub-tour elimination. More specifically, the while loop (NOT any of the while loops in the function get_plan but further below, the last one) below never exits and I can't figure out why the size of sub-tour list never goes to 1 in order to exit the while loop. I've spent quite a lot of time trying to debug this and I'd appreciate any help.

The code below should be completely reproducible, just copy and paste. Note that pip install pulp is required.

import pulp
import pandas as pd
import numpy as np
import copy

location_df = pd.DataFrame({'Location': ['Depot','LL716','LL384','LR59','LL701','LL866','LR830','LL1034','LR80','LR220','LL804'],
'x': [0.00,140.21,76.48,6.37,133.84,172.07,159.33,203.94,12.75,38.24,159.33],
'y': [0.00,30.62,0.00,68.90,74.00,5.10,76.55,25.52,40.83,71.45,10.21]})

N = len(location_df)

dist_mat = np.array([[  0.  , 170.83,  76.48,  75.27, 207.84, 177.17, 235.88, 229.46,
53.58, 109.69, 169.54],
[170.83,   0.  ,  94.35, 172.12,  49.75,  67.58,  65.05,  86.69,
137.67, 142.8 ,  57.39],
[ 76.48,  94.35,   0.  , 139.01, 131.36, 100.69, 159.4 , 152.98,
104.56, 109.69,  93.06],
[ 75.27, 172.12, 139.01,   0.  , 147.72, 229.5 , 170.66, 240.95,
37.01,  54.67, 211.65],
[207.84,  49.75, 131.36, 147.72,   0.  , 107.13,  38.09, 118.58,
156.82, 113.3 ,  89.28],
[177.17,  67.58, 100.69, 229.5 , 107.13,   0.  ,  84.19,  62.49,
195.05, 200.18,  28.05],
[235.88,  65.05, 159.4 , 170.66,  38.09,  84.19,   0.  ,  95.64,
184.86, 136.24,  66.34],
[229.46,  86.69, 152.98, 240.95, 118.58,  62.49,  95.64,   0.  ,
206.5 , 211.63,  80.34],
[ 53.58, 137.67, 104.56,  37.01, 156.82, 195.05, 184.86, 206.5 ,
0.  ,  58.67, 177.2 ],
[109.69, 142.8 , 109.69,  54.67, 113.3 , 200.18, 136.24, 211.63,
58.67,   0.  , 182.33],
[169.54,  57.39,  93.06, 211.65,  89.28,  28.05,  66.34,  80.34,
177.2 , 182.33,   0.  ]])

##################### Solve model using the DFJ subtour elimination

# find all sub-tours in the solution
def get_plan(r0):
r = copy.copy(r0)
route = []
while len(r) != 0:
plan = [r]
del (r)
l = 0
while len(plan) > l:
l = len(plan)
for i, j in enumerate(r):
if plan[-1] == j:
plan.append(j)
del (r[i])
route.append(plan)
return(route)

model = pulp.LpProblem('tspDFJ',pulp.LpMinimize)
#define variable
x = pulp.LpVariable.dicts("x",((i,j) for i in range(N) \
for j in range(N)), \
cat='Binary')

#set objective
model += pulp.lpSum(dist_mat[i][j] * x[i,j] for i in range(N) \
for j in range(N))
# st constraints
for i in range(len(location_df)):
model += x[i,i] == 0
model += pulp.lpSum(x[i,j] for j in range(N)) == 1
model += pulp.lpSum(x[j,i] for j in range(N)) == 1

status = model.solve()

route = [(i,j) for i in range(N) \
for j in range(N) if pulp.value(x[i,j]) == 1]

S = get_plan(route)
subtour = []

#Check if we got subtours, if we do, we

while len(S) != 1:
for i in range(len(S)):
#print(S[i])
model += pulp.lpSum(x[S[i][j], S[i][j]] \
for j in range(len(S[i])) if j!=i) <= len(S[i]) - 1

status = model.solve()
route = [(i,j) for i in range(N) \
for j in range(N) if pulp.value(x[i,j]) == 1]

S = get_plan(route)
subtour.append(len(S))

print("-----------------")
print(status,pulp.LpStatus[status],pulp.value(model.objective))
print(S)
print("no. of times LP model is solved = ", len(subtour))
print("subtour log (no. of subtours in each solution))", subtour)

• What is the meaning of if j!=i if the pulp.lpSum which adds the subtour elimination constraint? Jun 17, 2022 at 11:49
• @fontanf - I've tried it without it as well. I got it from this post: or.stackexchange.com/questions/6153/…. There, under one of the sums the condition is $j\neq i$. Jun 17, 2022 at 11:59
• In your code i is not the index of a vertex but the index of the subtour. i \neq j doesn't mean anything Jun 17, 2022 at 13:17
• @fontanf - Yes that's correct. Desperation makes you try anything. Removing if j!=i still does not change anything, the while loop never exits. I've implemented a lot more complex models than this with no issues, however this one just does not work. Jun 17, 2022 at 13:25
• It's not a good idea to add multiple times the same constraint. If I add if plan not in route: before route.append(plan), I don't know if it's good, but it terminates Jun 17, 2022 at 14:25

1. In constraint
while len(S) != 1:
for i in range(len(S)):
#print(S[i])
model += pulp.lpSum(x[S[i][j], S[i][j]] \
for j in range(len(S[i])) if j!=i) <= len(S[i]) - 1


if j!=i doesn't make sense since i is the index of the subtour, and not an index of a vertex

1. Multiple identical subtour are computed in the get_plan function. Replace
      route.append(plan)


by

      if plan not in route:
route.append(plan)

• Concerning point 1: The real mistake is that route.append(plan) is in the while loop but it should not be. So instead of this if statement, you could just move it out of the while loop by "re-indenting" it. Jun 18, 2022 at 6:23
• @PeterD maybe you're right, I didn't go that much into the code Jun 18, 2022 at 7:58