0
$\begingroup$

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.

$\endgroup$

2 Answers 2

6
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ 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 !! $\endgroup$
    – satk0
    Oct 26, 2023 at 21:57
  • 1
    $\begingroup$ Glad to help. Notice that the WY* constraints and subtour57 together yield an irreducible infeasible subset (IIS). $\endgroup$
    – RobPratt
    Oct 26, 2023 at 22:19
  • $\begingroup$ Oh, That is right! $6$ could never be smaller than $5$ haha. Brilliant! $\endgroup$
    – satk0
    Oct 26, 2023 at 22:58
1
$\begingroup$

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])}")
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.