$title Traveling Salesman Problem Instance prepared with AWK (AWKTSP,SEQ=298)

$onText
The model writes an AWK script on the fly to process the input file
format defined by the maintainers of the TSPLib. More input instances
are available from

https://github.com/rhgrant10/tsplib95/tree/master/archives/problems

Keywords: mixed integer linear programming, traveling salesman problem,
          iterative subtour elimination, TSPLib, awk script
$offText

$set fn p43.inc
$if not exist %fn% $abort %fn% ist not present

$echoN "$setGlobal n " > %gams.scrdir%n.%gams.scrext%
$call grep DIMENSION %fn% | cut -d: -f2 >> "%gams.scrdir%n.%gams.scrext%"
$include %gams.scrdir%n.%gams.scrext%

$onEcho > %gams.scrdir%x.%gams.scrext%
BEGIN                  { i=1 }
/^EDGE_WEIGHT_TYPE/    { if ( $2 != "EXPLICIT" )    error("Wrong type")   }
/^EDGE_WEIGHT_FORMAT/  { if ( $2 != "FULL_MATRIX" ) error("Wrong format") }
/^EDGE_WEIGHT_SECTION/ { readM = 1; next }
/^EOF/                 { exit }
readM == 1             { for (k=1; k <= NF; k++) print "i" i,"i" j+k,$k;
                         if (j+NF < %n%) j += NF; else { i++; j=0 }       }
function error(s) { print("--- " s ". Exiting.");  exit; }
$offEcho

$call awk -f "%gams.scrdir%x.%gams.scrext%" %fn% > "%gams.scrdir%data.%gams.scrext%"

Set
   ii    'cities'                                       / i1*i%n% /
   i(ii) 'subset of cities to fit GAMS demo model size' / i1*i7   /;

Alias (ii,jj), (i,j,k);

Parameter c(ii,jj) /
$onDelim
$include %gams.scrdir%data.%gams.scrext%
$offDelim
/;

Variable
   x(ii,jj) 'decision variables - leg of trip'
   z        'objective variable';

Binary Variable x;

Equation
   objective  'total cost'
   rowsum(ii) 'leave each city only once'
   colsum(jj) 'arrive at each city only once';

* The assignment problem is a relaxation of the TSP
objective.. z =e= sum((i,j), c(i,j)*x(i,j));

rowsum(i).. sum(j, x(i,j)) =e= 1;

colsum(j).. sum(i, x(i,j)) =e= 1;

* exclude diagonal
x.fx(ii,ii) = 0;

Set ij(ii,jj) 'exclude first row and column';
ij(ii,jj) = ord(ii) > 1 and ord(jj) > 1;

option optCr = 0, optCa = 0.99;

Model assign / objective, rowsum, colsum /;

solve assign using mip minimizing z;

* find and display tours
Set t 'tours' / t1*t17 /;
abort$(card(t) < card(i)) 'Set t is possibly too small';

Set
   tour(i,j,t) 'subtours'
   from(i)     'contains always one element: the from city'
   next(j)     'contains always one element: the to city'
   visited(i)  'flag whether a city is already visited'
   tt(t)       'contains always one element: the current subtour';

Alias (i,ix);

$eolCom //

* initialize
from(i)$(ord(i) = 1) = yes;    // turn first element on
tt(t)$(ord(t)   = 1) = yes;    // turn first element on

* Subtour elimination by adding cuts
Set cc / c1*c1000 /;

Alias (cc,ccc);

Set
   curcut(cc)  'current cut always one element'
   allcuts(cc) 'total cuts';

Parameter
   cutcoeff(cc,i,j)
   rhs(cc)
   nosubtours 'number of subtours';

Equation cut(cc) 'dynamic cuts';

cut(allcuts).. sum((i,j), cutcoeff(allcuts,i,j)*x(i,j)) =l= rhs(allcuts);

Model tspcut / objective, rowsum, colsum, cut /;

curcut(cc)$(ord(cc) = 1) = yes;
Scalar ok;

loop(ccc,
*  initialize
   from(i) = no;
   tt(t)   = no;
   from(i)$(ord(i) = 1) = yes;  // turn first element on
   tt(t)$(ord(t)   = 1) = yes;  // turn first element on
   tour(i,j,t) = no;
   visited(i)  = no;
   loop(i,
      next(j) = no;                  // find city visited after city 'from'
      loop((from,j),next(j)$(x.l(from,j) > 0.5) = yes);
                                     // check x.l(from,j) = 1 would be dangerous
      tour(from,next,tt) = yes;      // store in table
      visited(from)      = yes;      // mark city 'from' as visited
      from(j) = next(j);
      if(sum(visited(next),1) > 0,   // if already visited...
         tt(t) = tt(t-1);
         loop(ix$(not visited(ix)),  // find starting point of new subtour
            from(j)  = no;           // make sure we only have one element turned on
            from(ix) = yes;
          );
      );
   );
   display tour;
   nosubtours = sum(t, max(0,smax(tour(i,j,t),1)));
   display nosubtours;

   break$(nosubtours = 1); // done: no subtours

   // else: introduce cut
   loop(t$(ord(t) <= nosubtours),
      rhs(curcut) = -1;
      loop(tour(i,j,t),
         cutcoeff(curcut,i,j)$(x.l(i,j) > 0.5) = 1;
*  not needed due to nature of assignment constraints
*  cutcoeff(curcut,i,j)$(x.l(i,j) < 0.5) = -1;
         rhs(curcut) = rhs(curcut) + 1;
      );
      allcuts(curcut) = yes;   // include this cut in set
      curcut(cc) = curcut(cc-1);
   );
   solve tspcut using mip minimizing z;
   tspcut.solPrint = %solPrint.quiet%;
   tspcut.limRow = 0;
   tspcut.limCol = 0;
);

display x.l;
abort$(nosubtours <> 1) "Too many cuts needed";