tsp5.gms : TSP solution with Miller et al subtour elimination

Description

```This is the fifth problem in a series of traveling salesman
problems. Here we use a compact formulation that ensures no subtours
due to Miller, Tucker and Zemlin. We either run MTZ for a limited time
or apply a greedy heuristic with 2-opt improvement swaps and use this
as backup solution in case the subtour elimination runs out of time.

Moreover, we use a dynamic subtour elimination using even stronger
cuts than in the previous models. Also the GAMS programming is more
compact than in the previous examples.
```

Large Model of Type : MIP

Category : GAMS Model library

Main file : tsp5.gms   includes :  br17.inc

``````\$title Traveling Salesman Problem - Five (TSP5,SEQ=345)

\$onText
This is the fifth problem in a series of traveling salesman
problems. Here we use a compact formulation that ensures no subtours
due to Miller, Tucker and Zemlin. We either run MTZ for a limited time
or apply a greedy heuristic with 2-opt improvement swaps and use this
as backup solution in case the subtour elimination runs out of time.

Moreover, we use a dynamic subtour elimination using even stronger
cuts than in the previous models. Also the GAMS programming is more
compact than in the previous examples.

Miller, C E, Tucker, A W, and Zemlin, R A, Integer Programming
Formulation of Traveling Salesman Problems. J. ACM 7, 4 (1960),
326-329.

Keywords: mixed integer linear programming, traveling salesman problem, Miller-
Tucker-Zemlin subtour elimination, 2-opt heuristic, iterative subtour elimination
\$offText

\$include br17.inc

Set
quicktour(ii,jj) 'a tour found by MTZ model or 2-opt tour'
i(ii) / #ii /;

\$if not set BACKUP \$set BACKUP 2-opt
\$ifthen "%BACKUP%"=="2-opt"
EmbeddedCode Python:
import gams.transfer as gt
import numpy as np
def swap(p,i,j):
tmp = p[i]
p[i] = p[j]
p[j] = tmp

m = gt.Container(gams.db, system_directory=r"%gams.sysdir% ".rstrip())
n = len(gams.db['ii'])
c = m.data['c'].toDense()
# Quick greedy tour
p = np.zeros(n, dtype=int) # permutation vector
next = 0
for i in range(n-1):
c[:,next] = np.nan
next = np.nanargmin(c[next])
p[i+1] = next
# 2-opt swaps
c = m.data['c'].toDense()
tlen = sum(c[p[i],p[(i+1) % n]] for i in range(n))
gams.printLog(f'Greedy tour length: {tlen:.2f}')
nswaps = 1
while nswaps>0:
nswaps = 0
for i in range(n):
for j in range(i+1,n):
swap(p,i,j)
newlen = sum(c[p[i],p[(i+1) % n]] for i in range(n))
if newlen < tlen:
gams.printLog(f'Swap of {p[j]} and {p[i]} improved tour by {tlen-newlen:.2f}. New length: {newlen:.2f}')
tlen = newlen
nswaps = nswaps + 1
else:
swap(p,i,j) # swap back

gams.set('quicktour', [ (int(p[i]+1), int(p[(i+1) % n]+1)) for i in range(n)])
endEmbeddedCode quicktour
\$else

* Build compact TSP model
Positive Variable p(ii) 'position in tour';

Equations defMTZ(ii,jj) 'Miller, Tucker and Zemlin subtour elimination';

defMTZ(i,j).. p(i) - p(j) =l= card(j) - card(j)*x(i,j) - 1 + card(j)\$(ord(j) = 1);

Model MTZ / all /;

p.fx(j)\$(ord(j) = 1) = 0;
p.up(j) = card(j) - 1;

MTZ.resLim = 30;
solve MTZ min z using mip;
if (MTZ.modelStat= %modelStat.optimal% or MTZ.modelStat = %ModelStat.Integer Solution%,
quicktour(i,j) = round(x.l(i,j));
);
\$endIf

* Dynamic subtour elimination
Set
ste         'possible subtour elimination cuts' / c1*c1000 /
a(ste)      'active cuts'
tour(ii,jj) 'possible subtour'
n(jj)       'nodes visited by subtour';
Singleton Set
cur(ste)    'current cut'                       /c1/;

Parameter
cc(ste,ii,jj) 'cut coefficients'
rhs(ste)      'right hand side of cut';

Equation defste(ste) 'subtour elimination cut';

defste(a).. sum((i,j), cc(a,i,j)*x(i,j)) =l= rhs(a);

Model DSE / rowsum, colsum, objective, defste /;

a(ste)    = no;
cc(a,i,j) =  0;
rhs(a)    =  0;

option limRow = 0, limCol = 0, solPrint = silent, solveLink = 5;

repeat
solve DSE min z using mip;
abort\$(DSE.modelStat <> %modelStat.optimal% and DSE.modelStat <> %ModelStat.Integer Solution% ) 'problems with MIP solver';
x.l(i,j) = round(x.l(i,j));

*  Check for subtours
tour(i,j) = no;
n(i) = ord(i)=1;

while(card(n),
*  Construct a single subtour and remove it by setting x.l=0 for its edges
while(sum((n,j), x.l(n,j)),
loop((i,j)\$(n(i) and x.l(i,j)),
tour(i,j) = yes;
x.l(i,j)  =   0;
n(j)      = yes;
);
);
* Subtour
if(card(n) < card(j),
a(cur)   = yes;
rhs(cur) = card(n) - 1;
cc(cur,i,j)\$(n(i) and n(j)) = 1;
cur(ste+1) = cur(ste);
abort\$(card(cur)=0) 'Out of subtour cuts, enlarge set ste';
else
* Optimal
break 2;
);
n(j) = no;
loop((i,j)\$(card(n) = 0 and x.l(i,j)), n(i) = yes);
);
put_utility 'log' / 'Added ' (card(a)-card(defste)):0:0 ' subtour cuts'
until timeElapsed>30;

if(card(n) = card(j),
option tour:0:0:1;
display 'Optimal tour found', tour;
put_utility 'log' / 'Optimal Tour';
else
put_utility 'log' / 'Out of time';
put_utility\$card(quicktour) 'log' / 'Report Heuristic Tour (2-opt or time-limited MTZ)';
tour(i,j) = quicktour(i,j);
);

put_utility\$card(tour) 'log' / 'Tour length = ' (sum(tour(i,j), c(i,j))):0 ':';
n(i) = ord(i)=1;
while(1,
loop(tour(i,j)\$n(i),
put_utility 'log' / i.te(i) ' --> ' j.te(j) ' (' c(i,j):6:2 ')';
n(i) = no;
n(j) = yes;
break\$(ord(j)=1) 2;
);
);
``````
GAMS Development Corp.
GAMS Software GmbH

General Information and Sales
U.S. (+1) 202 342-0180
Europe: (+49) 221 949-9170