None formulations_sol

DM872 - MILP Formulations for the Traveling Salesman Problem

MTZ

The following formulation (C. E. Miller, A. W. Tucker, and R. A. Zemlin, “Integer programming formulations and traveling salesman problems,” J. ACM, 7 (1960), pp. 326–329) which also eliminates subtours has only a polynomial number of constraints. Let $x_{ij}$ be the variables for directed graphs and let $u_{i}$, $i=0,1,\ldots,n-1$ be the number of cities visited already when we arrive at city $i$. $$\begin{aligned} \min & \sum_{(ij) \in A} c_{ij} x_{ij}\\ &\sum_{i:i\neq j} x_{ij}=1 & \forall j=0,\ldots,n-1 \\ &\sum_{j:i\neq j} x_{ij}=1 &\forall i=0,\ldots,n-1 \\ &u_i - u_j + n x_{ij} \leq n - 1, & \forall i, j = 1,2,\ldots, n-1, i\neq j\\ &x_{ij}\in \mathbb{B} &\forall i,j,i\neq j\\ &u_i\in \mathbb{R} & \forall i = 1,\ldots,n-1 \end{aligned}$$

Note that the third set of constraints can be equivalently written as: $$u_i +1 \leq u_j + n (1-x_{ij}), \qquad i, j = 1, 2,\ldots, n-1, i\neq j $$

Let's fix a vertex, say $0$, to be the home base and for each other vertex $i$ let $u_i$ be an arbitrary real number. The $u_i$ variables play a role similar to node potentials in a network and the inequalities involving them serve to eliminate tours that do not begin and end at city $0$ and tours that visit more than $n-1$ cities. Consider a feasible solution to the TSP: it can be expressed as a permutation of cities, $\pi=(\pi_1=0,\pi_2,\ldots\pi_n,\pi_1=0)$. Hence $x_{\pi_1,\pi_2}=1$. Unless $\pi_2=0$ then there exists $\pi_3$ such that $x_{\pi_2\pi_3}=1$. We proceed in this way until some $\pi_i=0$, which if the solution is feasible happens only at the end of the permutation. The constraint forces $u_{\pi_j}\geq u_{\pi_i}+1$ when $x_{\pi_i,\pi_j}=1$, except when $\pi_{i+1}=0$. Hence, an assignment of $\vec u$ can be found when the solution is feasible. On the contrary, if there were subtours or a node was visited more than once, then an assignment for the variables $u_i$ that satisfies $u_{\pi_j}\geq u_{\pi_i}+1$ when $x_{\pi_i,\pi_j}=1$ could not be found (the constraint would be violated where the subtour closes).

You find an implementation in the git repository: mtz.py. Let's try it on the 20-cities random instance.

In [4]:
%run src/mtz.py
%matplotlib inline
GLPSOL: GLPK LP/MIP Solver, v4.65
Parameter(s) specified in the command line:
 --write /var/folders/w9/sxcp2ljj4wq3fdhy78rf5djh0000gn/T/tmp09x8lgry.glpk.raw
 --wglp /var/folders/w9/sxcp2ljj4wq3fdhy78rf5djh0000gn/T/tmpnmcqmj0o.glpk.glp
 --cpxlp /var/folders/w9/sxcp2ljj4wq3fdhy78rf5djh0000gn/T/tmpipypypsc.pyomo.lp
Reading problem data from '/var/folders/w9/sxcp2ljj4wq3fdhy78rf5djh0000gn/T/tmpipypypsc.pyomo.lp'...
/var/folders/w9/sxcp2ljj4wq3fdhy78rf5djh0000gn/T/tmpipypypsc.pyomo.lp:3724: warning: lower bound of variable 'x1' redefined
/var/folders/w9/sxcp2ljj4wq3fdhy78rf5djh0000gn/T/tmpipypypsc.pyomo.lp:3724: warning: upper bound of variable 'x1' redefined
383 rows, 400 columns, 1787 non-zeros
380 integer variables, all of which are binary
4104 lines were read
Writing problem data to '/var/folders/w9/sxcp2ljj4wq3fdhy78rf5djh0000gn/T/tmpnmcqmj0o.glpk.glp'...
3356 lines were written
GLPK Integer Optimizer, v4.65
383 rows, 400 columns, 1787 non-zeros
380 integer variables, all of which are binary
Preprocessing...
342 constraint coefficient(s) were reduced
382 rows, 399 columns, 1786 non-zeros
380 integer variables, all of which are binary
Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  1.900e+01  ratio =  1.900e+01
GM: min|aij| =  9.069e-01  max|aij| =  1.103e+00  ratio =  1.216e+00
EQ: min|aij| =  8.537e-01  max|aij| =  1.000e+00  ratio =  1.171e+00
2N: min|aij| =  1.000e+00  max|aij| =  1.188e+00  ratio =  1.188e+00
Constructing initial basis...
Size of triangular part is 381
Solving LP relaxation...
GLPK Simplex Optimizer, v4.65
382 rows, 399 columns, 1786 non-zeros
      0: obj =   9.499897478e+03 inf =   4.475e+01 (35)
     58: obj =   7.622956718e+03 inf =   2.887e-14 (0)
*   159: obj =   2.538879799e+03 inf =   2.503e-14 (0) 1
OPTIMAL LP SOLUTION FOUND
Integer optimization begins...
Long-step dual simplex will be used
+   159: mip =     not found yet >=              -inf        (1; 0)
+   474: >>>>>   4.200212637e+03 >=   2.544703204e+03  39.4% (32; 0)
+   800: >>>>>   3.564047418e+03 >=   2.625299150e+03  26.3% (65; 1)
+  1701: >>>>>   3.557525827e+03 >=   2.632103682e+03  26.0% (138; 36)
+  2223: >>>>>   3.415775367e+03 >=   2.671406077e+03  21.8% (188; 40)
+  4616: >>>>>   3.286400141e+03 >=   2.775354912e+03  15.6% (369; 105)
+ 16433: >>>>>   3.276447950e+03 >=   2.959466979e+03   9.7% (1304; 436)
+ 16594: >>>>>   3.211743952e+03 >=   2.959721857e+03   7.8% (1289; 500)
+ 27097: >>>>>   3.182275143e+03 >=   3.007513974e+03   5.5% (1832; 1212)
+ 70003: mip =   3.182275143e+03 >=   3.109710054e+03   2.3% (2355; 4414)
+ 93480: mip =   3.182275143e+03 >=     tree is empty   0.0% (0; 12885)
INTEGER OPTIMAL SOLUTION FOUND
Time used:   11.2 secs
Memory used: 10.1 Mb (10613183 bytes)
Writing MIP solution to '/var/folders/w9/sxcp2ljj4wq3fdhy78rf5djh0000gn/T/tmp09x8lgry.glpk.raw'...
792 lines were written
The optimal objective is  3182.275143072144
Computation time 11.36

Svestka

In JA Svestka, "A continuous variable representation of the TSP," Math Prog, v15, 1978, pp 211-213 a similar flow formulation is presented. The idea is the same as in the MTZ formulation but instead of working with potentials at nodes we have now flow $y$ on the arcs. The gain at every node except node 0 can be any value $\epsilon>0$ (differently from the MTZ formulation where it has been fixed to 1). The formulation is as follows:

$$\begin{aligned} \min & \sum_{ij \in A} c_{ij} x_{ij}\\ &\sum_{j:j\neq 0} y_{0j}=1 \\ &\sum_{j:j\neq i} y_{ji}\geq 1 & i=1..n-1\\ &\sum_{j:j\neq i} y_{ij}-\sum_{j:j\neq i} y_{ji}=\epsilon & i=1..n-1 \\ &y_{ij} \leq (1+n\cdot \epsilon)x_{ij} & i,j = 0..n-1\\ &\sum_{ij \in A} x_{ij}\leq n \\ &x_{ij}\in \mathbb{B} &i,j=0..n-1\\ &y_{ij} \geq 0 & i,j = 0..n-1 \end{aligned} $$

The fourth set of constraints links the non negative variables $y_{ij}$ to the binary variables $x_{ij}$ in such a way that $x_{ij}$ is one when the variable $y_{ij}$ are strictly positive.

With each solution we associate a directed graph $H\subset D$ in which an arc goes from $i$ to $j$ if and only if $x_{ij}$ is 1. To show that $H$ is an Hamiltonian cycle we have to show that:

(i) $H$ is connected

(ii) every node in $H$ is entered by at least one arc and left by at least one arc

(iii) $H$ has at most $n$ arcs.

If each $y_{ij}$ is thought of as a flow from $i$ to $j$ then each node except node $0$ has a gain of $\epsilon$ and so the node $0$ must have a loss of $(n - 1)\epsilon$. Now (i) follows since each connected component of $H$ must have a zero total gain. To verify (ii), let us first consider the node 0. By the first constraint, at least one arc leaves this node; to account for the positive loss, another arc must enter. In case $i\neq 0$, the second set of conditions guarantee that at least one arc enters $i$; another arc must leave it to account for the positive gain. Finally, the property (iii) is equivalent to the fifth constraint.

You find an implementation in the git repository: svestka.py. GLPK is unable to find an upper bound to this formulation within a couple of hours on the 20-cities random instance. Gurobi is instead very fast. We use gurobi therefore.

In [15]:
%run src/svestka.py
%matplotlib inline
Using license file /Library/gurobi900/gurobi.lic
Academic license - for non-commercial use only
Read LP format model from file /var/folders/w9/sxcp2ljj4wq3fdhy78rf5djh0000gn/T/tmpmbg4eiev.pyomo.lp
Reading time = 0.01 seconds
x761: 421 rows, 761 columns, 2243 nonzeros
Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (mac64)
Optimize a model with 421 rows, 761 columns and 2243 nonzeros
Model fingerprint: 0xe0756b41
Variable types: 381 continuous, 380 integer (380 binary)
Coefficient statistics:
  Matrix range     [1e+00, 3e+00]
  Objective range  [4e+01, 9e+02]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e-01, 2e+01]
Presolve removed 1 rows and 1 columns
Presolve time: 0.01s
Presolved: 420 rows, 760 columns, 2242 nonzeros
Variable types: 380 continuous, 380 integer (380 binary)

Root relaxation: objective 1.032800e+03, 545 iterations, 0.01 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 1032.80036    0   36          - 1032.80036      -     -    0s
     0     0 2582.13008    0   35          - 2582.13008      -     -    0s
     0     0 2586.08618    0   32          - 2586.08618      -     -    0s
     0     0 2748.29088    0   51          - 2748.29088      -     -    0s
     0     0 2845.37188    0   52          - 2845.37188      -     -    0s
     0     0 2855.33498    0   49          - 2855.33498      -     -    0s
     0     0 2855.36226    0   51          - 2855.36226      -     -    0s
     0     0 2932.48116    0   52          - 2932.48116      -     -    0s
     0     0 2933.00377    0   52          - 2933.00377      -     -    0s
     0     0 2966.79737    0   45          - 2966.79737      -     -    0s
     0     0 3012.78870    0   47          - 3012.78870      -     -    0s
     0     0 3065.13490    0   47          - 3065.13490      -     -    0s
     0     0 3067.26543    0   56          - 3067.26543      -     -    0s
     0     0 3071.86283    0   56          - 3071.86283      -     -    0s
     0     0 3073.90272    0   53          - 3073.90272      -     -    0s
     0     0 3075.36792    0   52          - 3075.36792      -     -    0s
     0     0 3075.36792    0   44          - 3075.36792      -     -    0s
H    0     0                    10282.571940 3075.36792  70.1%     -    0s
H    0     0                    3440.4862307 3095.78474  10.0%     -    0s
     0     0 3095.78474    0   17 3440.48623 3095.78474  10.0%     -    0s
     0     0 3095.78474    0   33 3440.48623 3095.78474  10.0%     -    0s
     0     0 3095.78474    0   33 3440.48623 3095.78474  10.0%     -    0s
     0     0 3095.78474    0   36 3440.48623 3095.78474  10.0%     -    0s
     0     0 3095.78474    0   45 3440.48623 3095.78474  10.0%     -    1s
     0     0 3095.78474    0   53 3440.48623 3095.78474  10.0%     -    1s
     0     0 3095.78474    0   52 3440.48623 3095.78474  10.0%     -    1s
     0     0 3095.78474    0   59 3440.48623 3095.78474  10.0%     -    1s
     0     0 3095.78474    0   56 3440.48623 3095.78474  10.0%     -    1s
     0     0 3095.78474    0   59 3440.48623 3095.78474  10.0%     -    1s
     0     0 3095.78474    0   55 3440.48623 3095.78474  10.0%     -    1s
     0     0 3095.78474    0   55 3440.48623 3095.78474  10.0%     -    1s
     0     0 3095.78474    0   51 3440.48623 3095.78474  10.0%     -    1s
     0     0 3095.78474    0   51 3440.48623 3095.78474  10.0%     -    1s
H    0     0                    3440.4858495 3105.16247  9.75%     -    1s
     0     2 3105.16247    0   50 3440.48585 3105.16247  9.75%     -    1s
*   19    22               7    3182.2751431 3139.84107  1.33%  49.0    1s

Cutting planes:
  Gomory: 3
  Cover: 14
  MIR: 49
  Flow cover: 84
  Flow path: 1
  Zero half: 2
  Mod-K: 1
  Network: 2

Explored 73 nodes (5422 simplex iterations) in 1.52 seconds
Thread count was 4 (of 4 available processors)

Solution count 3: 3182.28 3440.49 10282.6 

Optimal solution found (tolerance 1.00e-04)
Best objective 3.182275143072e+03, best bound 3.182275143072e+03, gap 0.0000%
The optimal objective is  3182.275143072144
Computation time 3.10

Dantzig

In "Linear Programming and Extensions" 1963, Dantzig G.B. gave a three index formulation for the asymmetric TSP. It considers $n$ different levels indexed by $t$ and looks for tours that visit exactly one node at each level and step up in a coordinated way between levels.

$$ \begin{aligned} \min & \sum_{ijt} c_{ij} x_{ijt}\\ &\sum_{t=0}^{n-1}\sum_{j:j\neq i} x_{i,j,t}=1 & i=0..n-1 \\ &\sum_{i:i\neq j} x_{i,j,t}=\sum_{k:k\neq j} x_{j,k,t+1} & j=0..n-1, t=0..n-2\\ &\sum_{i:i\neq j} x_{i,j,n-1}=\sum_{k:k\neq j} x_{j,k,0} & j=0..n-1\\ &x_{ijt}\in \mathbb{B} &i,j=0..n-1\\ \end{aligned} $$

Here the first set of constraints ensure that every node is visited, the second that the steps are done from the nodes we entered to and the third closes the cycle among the levels.

You find an implementation in the git repository: dantzig.py.

In [14]:
%run src/dantzig.py
%matplotlib inline
Using license file /Library/gurobi900/gurobi.lic
Academic license - for non-commercial use only
Read LP format model from file /var/folders/w9/sxcp2ljj4wq3fdhy78rf5djh0000gn/T/tmpslxvti6l.pyomo.lp
Reading time = 0.03 seconds
x7601: 421 rows, 7601 columns, 22801 nonzeros
Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (mac64)
Optimize a model with 421 rows, 7601 columns and 22801 nonzeros
Model fingerprint: 0x72b3f570
Variable types: 1 continuous, 7600 integer (7600 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [4e+01, 9e+02]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Presolve removed 1 rows and 1 columns
Presolve time: 0.03s
Presolved: 420 rows, 7600 columns, 22800 nonzeros
Variable types: 0 continuous, 7600 integer (7600 binary)

Root relaxation: objective 2.504182e+03, 888 iterations, 0.04 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 2504.18187    0  190          - 2504.18187      -     -    0s
H    0     0                    3624.6250828 2504.18187  30.9%     -    0s
     0     0 2504.18187    0  213 3624.62508 2504.18187  30.9%     -    0s
     0     0 2504.18187    0  190 3624.62508 2504.18187  30.9%     -    0s
     0     2 2536.19642    0  188 3624.62508 2536.19642  30.0%     -    0s
H   28    31                    3352.0600286 2608.62139  22.2%   187    1s
H  144   149                    3295.4767738 2608.62139  20.8%  72.0    2s
   618   565 3158.41699   55  105 3295.47677 2785.31393  15.5%  44.5    5s
H  672   573                    3277.1644177 2813.44583  14.1%  44.7    5s
H  688   576                    3269.4973350 2813.99972  13.9%  45.1    5s
   752   612 3083.00927   27  160 3269.49734 2884.53938  11.8%  47.1   10s
H  752   581                    3239.8213064 2885.83983  10.9%  47.1   10s
   785   603 3074.15719   40  109 3239.82131 2939.14658  9.28%  71.2   15s
*  788   572              40    3232.9779809 2939.14658  9.09%  71.3   15s
H  792   544                    3203.5087657 2939.14658  8.25%  72.2   15s
H  795   516                    3189.4144781 2939.14658  7.85%  72.2   15s
*  812   484              42    3182.2751431 2939.14658  7.64%  76.7   15s
  1026   472 2965.50829   56  274 3182.27514 2965.50829  6.81%  98.2   22s
H 1040   444                    3182.2751283 2965.50829  6.81%  99.5   22s
  1295   394     cutoff   68      3182.27513 2969.62002  6.68%   113   25s
  1665   378 3080.59141   79  104 3182.27513 2978.35594  6.41%   122   30s
  1963   374     cutoff  133      3182.27513 3000.46403  5.71%   130   38s
H 1964   351                    3182.2751248 3000.46403  5.71%   130   38s
  2022   346 3101.52181  143  155 3182.27512 3006.63593  5.52%   133   40s
  2138   344 infeasible   81      3182.27512 3026.60789  4.89%   150   45s
  2291   347 infeasible   80      3182.27512 3044.89760  4.32%   172   50s
  2446   265 infeasible  324      3182.27512 3076.05505  3.34%   198   56s
  2764     0     cutoff   94      3182.27512 3127.48587  1.72%   208   60s

Cutting planes:
  Gomory: 10
  Cover: 3
  MIR: 7
  StrongCG: 1
  Inf proof: 2
  Zero half: 53

Explored 2823 nodes (580194 simplex iterations) in 60.16 seconds
Thread count was 4 (of 4 available processors)

Solution count 10: 3182.28 3189.41 3203.51 ... 3624.63

Optimal solution found (tolerance 1.00e-04)
Best objective 3.182275124785e+03, best bound 3.182275124785e+03, gap 0.0000%
The optimal objective is  3182.2751247851697
Computation time 62.19

Comparison

Let's compare the different formulations encountered on the basis of their linear relaxation:

Instance DFJ MTZ Svestka Dantzig
ran20points
dantzig42.dat 2538.8 1032.8 2504.2
berlin52.dat
bier127.dat

Keep now the integrality constraints of the main variables in the formulations so that an optimal solution can be found. Which is the fastest formulation to solve the problem? How many branch and bound nodes are needed in the formulations to find the optimal solution?

Experiment ideas. For example, try combining the two formulations DFJ and MTZ.