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

Question

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[0]]
    del (r[0])
    l = 0
    while len(plan) > l:
      l = len(plan)
      for i, j in enumerate(r):
        if plan[-1][1] == j[0]:
          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][0], S[i][j][1]] \
                          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)

Answer 1

1. In constraint

while len(S) != 1:
  for i in range(len(S)):
    #print(S[i])
    model += pulp.lpSum(x[S[i][j][0], S[i][j][1]] \
                          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

2. Multiple identical subtour are computed in the get_plan function. Replace

      route.append(plan)

by

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