# "GLPK: No Primal Feasible Solution" for Travelling Salesman Problem (Dantzig–Fulkerson–Johnson formulation)

I have a problem with implementing the Dantzig–Fulkerson–Johnson formulation to solve the following Travelling Salesman Problem:

$$\begin{bmatrix} M & 21 & 17 & 19 & 21 & 22\\ 26 & M & 18 & 27 & 27 & 26\\ 20 & 25 & M & 27 & 26 & 24\\ 18 & 20 & 15 & M & 17 & 19\\ 24 & 25 & 15 & 26 & M & 25\\ 18 & 18 & 27 & 21 & 25 & M \end{bmatrix}$$

I managed to create the following code (sorry for not optimizing it, I'm a total newbie):

var x12 >= 0 binary;
var x13 >= 0 binary;
var x14 >= 0 binary;
var x15 >= 0 binary;
var x16 >= 0 binary;

var x21 >= 0 binary;
var x23 >= 0 binary;
var x24 >= 0 binary;
var x25 >= 0 binary;
var x26 >= 0 binary;

var x31 >= 0 binary;
var x32 >= 0 binary;
var x34 >= 0 binary;
var x35 >= 0 binary;
var x36 >= 0 binary;

var x41 >= 0 binary;
var x42 >= 0 binary;
var x43 >= 0 binary;
var x45 >= 0 binary;
var x46 >= 0 binary;

var x51 >= 0 binary;
var x52 >= 0 binary;
var x53 >= 0 binary;
var x54 >= 0 binary;
var x56 >= 0 binary;

var x61 >= 0 binary;
var x62 >= 0 binary;
var x63 >= 0 binary;
var x64 >= 0 binary;
var x65 >= 0 binary;

minimize z: 21*x12 + 17*x13 + 19*x14 + 21*x15 + 22*x16 + 26*x21 +18*x23 + 27*x24+ 27*x25
+ 26*x26+ 20*x31 + 25*x32 + 27*x34 + 26*x35 + 24*x36 + 18*x41 + 20*x42 + 15*x43 + 17*x45 + 19*x46
+ 24*x51 + 25*x52 + 15*x53 + 26*x54 + 25*x56 + 18*x61 + 18*x62 + 27*x63 + 21*x64 + 25*x65;

subject to WY1:   x12+x13+x14+x15+x16=1;
subject to WY2:   x21+x23+x24+x25+x26=1;
subject to WY3:   x31+x32+x34+x35+x36=1;
subject to WY4:   x41+x42+x43+x45+x46=1;
subject to WY5:   x51+x52+x53+x54+x56=1;
subject to WY6:   x61+x62+x63+x64+x65=1;

subject to WE1:   x21+x31+x41+x51+x61=1;
subject to WE2:   x12+x32+x42+x52+x62=1;
subject to WE3:   x13+x23+x43+x53+x63=1;
subject to WE4:   x14+x24+x34+x54+x64=1;
subject to WE5:   x15+x25+x35+x45+x65=1;
subject to WE6:   x16+x26+x36+x46+x56=1;

# subtour for each subset of vertices of length >= 2:
subject to subtour1: x12 + x21 <= 1;
subject to subtour2: x13 + x31 <= 1;
subject to subtour3: x14 + x41 <= 1;
subject to subtour4: x15 + x51 <= 1;
subject to subtour5: x16 + x61 <= 1;
subject to subtour6: x23 + x32 <= 1;
subject to subtour7: x24 + x42 <= 1;
subject to subtour8: x25 + x52 <= 1;
subject to subtour9: x26 + x62 <= 1;
subject to subtour10: x34 + x43 <= 1;
subject to subtour11: x35 + x53 <= 1;
subject to subtour12: x36 + x63 <= 1;
subject to subtour13: x45 + x54 <= 1;
subject to subtour14: x46 + x64 <= 1;
subject to subtour15: x56 + x65 <= 1;
subject to subtour16: x12 + x13 + x21 + x23 + x31 + x32 <= 2;
subject to subtour17: x12 + x14 + x21 + x24 + x41 + x42 <= 2;
subject to subtour18: x12 + x15 + x21 + x25 + x51 + x52 <= 2;
subject to subtour19: x12 + x16 + x21 + x26 + x61 + x62 <= 2;
subject to subtour20: x13 + x14 + x31 + x34 + x41 + x43 <= 2;
subject to subtour21: x13 + x15 + x31 + x35 + x51 + x53 <= 2;
subject to subtour22: x13 + x16 + x31 + x36 + x61 + x63 <= 2;
subject to subtour23: x14 + x15 + x41 + x45 + x51 + x54 <= 2;
subject to subtour24: x14 + x16 + x41 + x46 + x61 + x64 <= 2;
subject to subtour25: x15 + x16 + x51 + x56 + x61 + x65 <= 2;
subject to subtour26: x23 + x24 + x32 + x34 + x42 + x43 <= 2;
subject to subtour27: x23 + x25 + x32 + x35 + x52 + x53 <= 2;
subject to subtour28: x23 + x26 + x32 + x36 + x62 + x63 <= 2;
subject to subtour29: x24 + x25 + x42 + x45 + x52 + x54 <= 2;
subject to subtour30: x24 + x26 + x42 + x46 + x62 + x64 <= 2;
subject to subtour31: x25 + x26 + x52 + x56 + x62 + x65 <= 2;
subject to subtour32: x34 + x35 + x43 + x45 + x53 + x54 <= 2;
subject to subtour33: x34 + x36 + x43 + x46 + x63 + x64 <= 2;
subject to subtour34: x35 + x36 + x53 + x56 + x63 + x65 <= 2;
subject to subtour35: x45 + x46 + x54 + x56 + x64 + x65 <= 2;
subject to subtour36: x12 + x13 + x14 + x21 + x23 + x24 + x31 + x32 + x34 + x41 + x42 + x43 <= 3;
subject to subtour37: x12 + x13 + x15 + x21 + x23 + x25 + x31 + x32 + x35 + x51 + x52 + x53 <= 3;
subject to subtour38: x12 + x13 + x16 + x21 + x23 + x26 + x31 + x32 + x36 + x61 + x62 + x63 <= 3;
subject to subtour39: x12 + x14 + x15 + x21 + x24 + x25 + x41 + x42 + x45 + x51 + x52 + x54 <= 3;
subject to subtour40: x12 + x14 + x16 + x21 + x24 + x26 + x41 + x42 + x46 + x61 + x62 + x64 <= 3;
subject to subtour41: x12 + x15 + x16 + x21 + x25 + x26 + x51 + x52 + x56 + x61 + x62 + x65 <= 3;
subject to subtour42: x13 + x14 + x15 + x31 + x34 + x35 + x41 + x43 + x45 + x51 + x53 + x54 <= 3;
subject to subtour43: x13 + x14 + x16 + x31 + x34 + x36 + x41 + x43 + x46 + x61 + x63 + x64 <= 3;
subject to subtour44: x13 + x15 + x16 + x31 + x35 + x36 + x51 + x53 + x56 + x61 + x63 + x65 <= 3;
subject to subtour45: x14 + x15 + x16 + x41 + x45 + x46 + x51 + x54 + x56 + x61 + x64 + x65 <= 3;
subject to subtour46: x23 + x24 + x25 + x32 + x34 + x35 + x42 + x43 + x45 + x52 + x53 + x54 <= 3;
subject to subtour47: x23 + x24 + x26 + x32 + x34 + x36 + x42 + x43 + x46 + x62 + x63 + x64 <= 3;
subject to subtour48: x23 + x25 + x26 + x32 + x35 + x36 + x52 + x53 + x56 + x62 + x63 + x65 <= 3;
subject to subtour49: x24 + x25 + x26 + x42 + x45 + x46 + x52 + x54 + x56 + x62 + x64 + x65 <= 3;
subject to subtour50: x34 + x35 + x36 + x43 + x45 + x46 + x53 + x54 + x56 + x63 + x64 + x65 <= 3;
subject to subtour51: x12 + x13 + x14 + x15 + x21 + x23 + x24 + x25 + x31 + x32 + x34 + x35 + x41 + x42 + x43 + x45 + x51 + x52 + x53 + x54 <= 4;
subject to subtour52: x12 + x13 + x14 + x16 + x21 + x23 + x24 + x26 + x31 + x32 + x34 + x36 + x41 + x42 + x43 + x46 + x61 + x62 + x63 + x64 <= 4;
subject to subtour53: x12 + x13 + x15 + x16 + x21 + x23 + x25 + x26 + x31 + x32 + x35 + x36 + x51 + x52 + x53 + x56 + x61 + x62 + x63 + x65 <= 4;
subject to subtour54: x12 + x14 + x15 + x16 + x21 + x24 + x25 + x26 + x41 + x42 + x45 + x46 + x51 + x52 + x54 + x56 + x61 + x62 + x64 + x65 <= 4;
subject to subtour55: x13 + x14 + x15 + x16 + x31 + x34 + x35 + x36 + x41 + x43 + x45 + x46 + x51 + x53 + x54 + x56 + x61 + x63 + x64 + x65 <= 4;
subject to subtour56: x23 + x24 + x25 + x26 + x32 + x34 + x35 + x36 + x42 + x43 + x45 + x46 + x52 + x53 + x54 + x56 + x62 + x63 + x64 + x65 <= 4;
subject to subtour57: x12 + x13 + x14 + x15 + x16 + x21 + x23 + x24 + x25 + x26 + x31 + x32 + x34 + x35 + x36 + x41 + x42 + x43 + x45 + x46 + x51 + x52 + x53 + x54 + x56 + x61 + x62 + x63 + x64 + x65 <= 5;

end;

I managed to generate all the constraints in Python as I don't know how to implement it in the model language.

Unfortunately, the code is not working:

glpsol --model code.mod
Model has been successfully generated
GLPK Integer Optimizer 5.0
70 rows, 30 columns, 570 non-zeros
30 integer variables, all of which are binary
Preprocessing...
69 rows, 30 columns, 540 non-zeros
30 integer variables, all of which are binary
Scaling...
A: min|aij| =  1.000e+00  max|aij| =  1.000e+00  ratio =  1.000e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part is 68
Solving LP relaxation...
GLPK Simplex Optimizer 5.0
69 rows, 30 columns, 540 non-zeros
0: obj =   1.450000000e+02 inf =   3.400e+01 (18)
7: obj =   1.280000000e+02 inf =   1.000e+00 (1)
LP HAS NO PRIMAL FEASIBLE SOLUTION
Time used:   0.0 secs
Memory used: 0.4 Mb (433124 bytes)

I know it is possible to limit the constraints thanks to the lazy constraint but a teacher wants us to write all the constraints.

Would be great if anyone could help me solve this issue or even simplify my code.

### EDIT -> SOLUTION:

First of all, big thanks to @RobPratt and @ytsao for the answers!

Turns out that, @RobPratt's answer was the solution.

Expanding on his solution, I forgot about the fact that $$Q$$ must be a proper subset of vertices set:

$$Q\subset \{1, 2, 3, 4, 5, 6\} \implies Q \ne \{1, 2, 3, 4, 5, 6\}$$

, that means it cannot be equal to vertices set.

Therefore, the error is the subtour57 constraint that comes from the $$\{1, 2, 3, 4, 5, 6\}$$ subset, deleting it solved the problem.

The problem is indeed infeasible but becomes feasible if you omit the subtour57 constraint. For a tour on $$6$$ nodes, you need to use $$6$$ arcs, but that constraint allows only $$5$$.

• Exactly, works perfectly, I forgot about the fact that the $Q$ must be a proper subset, as I explained in an edit to the question. Thanks a lot !! Commented Oct 26, 2023 at 21:57
• Glad to help. Notice that the WY* constraints and subtour57 together yield an irreducible infeasible subset (IIS). Commented Oct 26, 2023 at 22:19
• Oh, That is right! $6$ could never be smaller than $5$ haha. Brilliant! Commented Oct 26, 2023 at 22:58

For python implementation using pyscipopt, you can refer as following:

def Combination(n:int):
Q = []

for i in range(2, n):
Q.append(list(itertools.combinations(cities, i)))

result = []
for i in Q:
for j in i:
result.append(list(j))

return result

n = 6
cities = list(range(n))

distance_matrix = [
[9999, 21, 17, 19, 21, 22],
[26, 9999, 18, 27, 27, 26],
[20, 25, 9999, 27, 26, 24],
[18, 20, 15, 9999, 17, 19],
[24, 25, 15, 26, 9999, 25],
[18, 18, 27, 21, 25, 9999]
]

model = pyscipopt.Model("TSP")

# Create variables
x = {}
for i in range(n):
for j in range(n):
if i == j:
x[i,j] = model.addVar(lb=0, ub=0, vtype="B", name="x(%s,%s)"%(i,j))
else:
x[i,j] = model.addVar(lb=0, ub=1, vtype="B", name="x(%s,%s)"%(i,j))

# Create objective function
model.setObjective(pyscipopt.quicksum(distance_matrix[i][j]*x[i,j] for i in range(n) for j in range(n)), sense="minimize")

# Create constraints
# Constraints 1
for i in range(n):
model.addCons(pyscipopt.quicksum(x[i,j] for j in range(n)) == 1)

# Constraints 2
for j in range(n):
model.addCons(pyscipopt.quicksum(x[i,j] for i in range(n)) == 1)

# Subtour-Elimination constraints (DFJ formulation)
Q = Combination(n)

for each_q in Q:
model.addCons(pyscipopt.quicksum(x[i,j]for i in each_q for j in each_q if i != j) <= len(each_q)-1)

model.writeProblem("TSP.lp")
model.optimize()
status = model.getStatus()
print(f"Solution status: {status}")
if status == "optimal":
print(f"Objective value = {model.getObjVal()}")
for i in range(n):
for j in range(n):
print(f"{x[i,j]} = {model.getVal(x[i,j])}")