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.integerSolution%,
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.integerSolution% ) '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;
);
);