Description

The model writes an AWK script on the fly to process the an 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

The subtour elimination constraints are supplied on the fly while
GAMS/CPLEX is running. The lazy constraint BCH call checks if
the integer solution contains subtours, supplies the corresponding
cuts, and rejects the solution.

Keywords: mixed integer linear programming, traveling salesman problem, branch
          and cut and heuristic facility, iterative subtour elimination

Small Model of Type : MIP

Category : GAMS Model library

Main file : bchtsp.gms   includes :  p43.inc

$title Traveling Salesman Problem Instance with BCH (BCHTSP,SEQ=348)

$onText
The model writes an AWK script on the fly to process the an 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

The subtour elimination constraints are supplied on the fly while
GAMS/CPLEX is running. The lazy constraint BCH call checks if
the integer solution contains subtours, supplies the corresponding
cuts, and rejects the solution.

Keywords: mixed integer linear programming, traveling salesman problem, branch
          and cut and heuristic facility, iterative subtour elimination
$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;

option optCr = 0;

Model assign / all /;

execute_unload 'tspdata', i;
option mip = cplex;
assign.optFile = 1;

solve assign using mip minimizing z;
abort$(assign.modelStat <> 1 and abs(78 - z.l) < 1e-6) 'wrong solution';

$onEcho > cplex.opt
preInd 0
userLazyConCall subtourcheck
$offEcho

$onEchoV > subtourcheck.gms
$title TSP Subtour Emlinitaion - Check

Set i cities;
$gdxIn tspdata
$load i
$eval cardi card(i)

Alias (i,j,k);

Binary Variable x(i,j) 'decision variables - leg of trip';

$gdxIn bchout_i.gdx
$load x

Set
   t           'cuts tours'  / 1*%cardi% /
   tour(t,i,j) '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';

$eolCom //

from(i)$(ord(i) = 1) = yes;    // turn first element on
tt(t)$(ord(t)   = 1) = yes;    // turn first element on
tour(t,i,j) = no;
visited(i)  = no;

while(card(visited) < card(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(tt,from,next) = 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);           // advance tt index
      loop(k$(not visited(k)),   // find starting point of new subtour
         from(j) =  no;          // make sure we only have one element turned on
         from(k) = yes;
      );
   );
);
* display tour;

* In order to communicate with the MIP solver we need certain conventions
* 1. Cut matrix interface requirement with fixed names
Parameter
   rhs_c(t)   'cut rhs'
   sense_c(t) 'the sense of the cuts 1 =l=, 2 =e=, 3 =g='
   numcuts    'number of cuts passed back'

* 2. Nonzeros in cut matrix using the original variable name plus _c
   x_c(t,i,j) 'coefficient of the x variables in the cut';

option tt < tour;  // Get active tours
abort$(card(tt) = 1) 'Accept solution';

numcuts = card(tt);
loop(tt,
   from(i)           = sum(tour(tt,i,j), 1); // all cities of the tour;
   next(i)           = from(i);
   x_c(tt,from,next) = 1;
   rhs_c(tt)         = sum(tour(tt,i,j), 1) - 1;
   sense_c(tt)       = 1;
);
display numcuts, x_c, rhs_c, sense_c;
$offEcho