$title Traveling Salesman Problem - One (TSP1,SEQ=177)
$onText
This is the first problem in a series of traveling salesman
problems. In this problem we first solve an assignment
problem as a relaxation of the TSP. Subtours of this solution
are detected and printed. The subtours are then eliminated via
cuts (constraints that eliminate solution with subtours).
Note: we deal here with an unsymmetric TSP. If symmetric
one can add 2 cuts in each cycle: forward and
backward.
Additional information can be found at:
http://www.gams.com/modlib/adddocs/tsp1doc.htm
Kalvelagen, E, Model Building with GAMS. forthcoming
de Wetering, A V, private communication.
Keywords: mixed integer linear programming, traveling salesman problem, iterative
subtour elimination, cut generation
$offText
$eolCom //
$include br17.inc
* For this simple algorithm the problem is too difficult
* so we consider only the first 6 cities.
Set i(ii) / i1*i6 /;
* options. Make sure MIP solver finds global optima.
option optCr = 0;
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";
Parameter tour(i,j,t) 'subtours';
Singleton Set
from(i) 'contains always one element: the from city'
next(j) 'contains always one element: the to city'
tt(t) 'contains always one element: the current subtour';
Set visited(i) 'flag whether a city is already visited';
* initialize
from(i)$(ord(i) = 1) = yes; // turn first element on
tt(t)$( ord(t) = 1) = yes; // turn first element on
loop(i,
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(k$(not visited(k)), // find starting point of new subtour
from(k) = yes;
);
);
);
display tour;
* subtour elimination by adding cuts
* the logic to detect if there are subtours is similar
* to the code above
Set cc / c1*c100 /; // we allow up to 100 cuts
Alias (cc,ccc);
Set
curcut(cc) 'current cut always one element'
allcuts(cc) 'total cuts';
Parameter cutcoeff(cc, i, j);
Equation cut(cc) 'dynamic cuts';
cut(allcuts).. sum((i,j), cutcoeff(allcuts,i,j)*x(i,j)) =l= card(i) - 1;
Model tspcut / objective, rowsum, colsum, cut /;
curcut(cc)$(ord(cc) = 1) = yes;
Scalar ok;
loop(ccc,
from(i) = ord(i) = 1; // initialize from to first city
visited(i) = no;
ok = 1;
loop(i$((ord(i) < card(i)) and ok), // last city can be ignored
next(j) = x.l(from,j) > 0.5; // find next city
visited(from) = yes;
from(j) = next(j);
ok$sum(visited(next),1) = 0; // we have detected a subtour
);
break$(ok = 1); // done: no subtours
// introduce cut
cutcoeff(curcut, i, j)$(x.l(i,j) > 0.5) = 1;
// next one is needed in the general case but not for TSP
// cutcoeff(curcut, i, j)$(x.l(i,j) < 0.5) = -1;
allcuts(curcut) = yes; // include this cut in set
curcut(cc) = curcut(cc-1); // get next element
solve tspcut using mip minimizing z;
tspcut.solPrint = %solPrint.Quiet%;
);
* print solution in proper order
Set xtour 'ordered tour';
from(i) = ord(i) = 1; // initialize from to first city
visited(i) = no;
ok = 1;
loop(t$ok,
next(j) = x.l(from,j) > 0.5; // find next city
xtour(t,from,next) = yes;
visited(from) = yes;
from(j) = next(j);
ok$sum(visited(next),1) = 0; // we have detected a subtour
);
option xtour:0:0:1;
display xtour,x.l;
abort$(card(allcuts) = card(cc)) "Too many cuts needed";