$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;
   );
);