# largetext

largetext:
$ontext

  Minimum Spanning Tree using an  LP formulation

  Original model: Erwin Kalvelagen, april, 2002
  Model solved from Access: Erwin Kalvelagen, November 2007

$offtext

$if not set db $error db is not set

$onecho >  mdb2gms.txt
i=%db%
x=db.gdx
q1=select distinct code from cities
s1=i
q2=select city1,city2,distance from distances
p2=d
$offecho

$call =mdb2gms.exe @mdb2gms.txt

set i 'cities';
alias (i,j,k);
parameter d(i,j) 'distances';

$gdxin db.gdx
$load i,d

display i,d;

alias(i,t);

set s(i) 'source node (can be any node)' /c40/;

set term(i) 'terminal nodes';
term(i)$(not s(i)) = yes;

parameter uarc(i,j) 'undirected arcs';
uarc(i,j)=d(i,j);


parameter c(i,j) 'cost directed arcs';
c(i,j) = uarc(i,j) + uarc(j,i);

set links(i,j);
links(i,j)$c(i,j) = yes;

parameter b(i,t) 'right-hand side';
b(s,term) = -1;
b(term,term) = 1;

variables
   z         'objective'
   x(i,j)    '0-1 variable indicating tree'
   y(i,j,t)  'flow variables (last index is terminal node)'
;

binary variable x,y;

equations
   objective      'minimize total cost'
   nodebal(i,t)   'node balance'
   compare(i,j,t) 'force x(i,j)>=y(i,j,t)'
;


objective..    Z =e= sum(links, c(links)*x(links));
nodebal(i,term(t))..  Sum(links(i,j), y(i,j,t)) - sum(links(j,i), y(j,i,t)) =e= b(i,t);
compare(links(i,j),term(t)).. X(i,j) =g= y(i,j,t);


option iterlim=100000;
model m /all/;
m.optfile=1;
solve m minimizing z using rmip;

display y.l, x.l;

file log /''/;
if (m.modelstat<>1,
    put log,  "**** ERROR: Model not solved to optimality"/;
    abort$1  "Model not solved to optimality";
);

$if exist solution.csv $call del solution.csv
file sol /solution.csv/;
sol.pc = 5;
put sol;
loop((i,j)$(x.l(i,j)>0.5),
   put i.tl,j.tl/;
);


$onecho > cplex.opt
lpmethod 4
$offecho


largetext:
$ontext

This model solves a Symmetric TSP using a simple algorithm
that adds cuts to exclude subtours found in the previous
solution.

Reference: Dantzig, Fulkerson, Johnson, "Solution of a Large-Scale
           Traveling-Salesman Problem", Operations Research,
           2, 1954, pp.393-410.


Copyright Erwin Kalvelagen, 2007

$offtext

option optcr=0;
option iterlim=100000;
option reslim=10000;
option limrow=0;
option limcol=0;
option solprint=off;

$eolcom //

$if not set db $error db is not set

$onecho >  mdb2gms.txt
i=%db%
x=db.gdx
q1=select distinct code from cities
s1=i
q2=select city1,city2,distance from distances
p2=d
$offecho

$call =mdb2gms.exe @mdb2gms.txt

set i 'cities';
alias (i,j,k);
parameter d(i,j) 'distances';

$gdxin db.gdx
$load i,d

display i,d;


set lt(i,j) lower triangular;
lt(i,j)$(ord(i)>ord(j)) = yes;

free variable z;
binary variable x(i,j);
equations twomatch(i), obj;


obj.. z =e= sum(lt(i,j), d(i,j)*x(i,j));
twomatch(k).. sum(lt(i,k), x(i,k)) + sum(lt(k,j),x(k,j)) =e= 2;

model match /obj,twomatch/;

solve match minimizing z using mip;
file log /''/;
if (match.modelstat<>1,
    put log,  "**** ERROR: Model not solved to optimality"/;
    abort$1  "Model not solved to optimality";
);


display x.l;



*
* dynamic number of cuts
*

set cycle /cycle1*cycle50/;
set t tours  /t1*t100/;
set cutindex(cycle,t)              dynamic set for addressing cuts;
equation cut(cycle,t)              actual cuts;
parameter cutcoeff(cycle,t,i,j)    coefficients in cuts;
parameter cutlength(cycle,t)       rhs for cuts;

cut(cutindex).. sum((i,j), cutcoeff(cutindex,i,j)*x(i,j)) =l= cutlength(cutindex)-1;

model tsp /obj, twomatch, cut/;


*
* used to find and display tours
*
set s(i,j) current solution;
set tour(t,i,j) subtours;
set tt(t) current subtour;
set first(i,j) 'first (i,j) in S';
set reach(i,j) "(i,j)'s connected to tour(tt)";

scalar done /0/;
scalar continue;

parameter results;  // for reporting
results('relaxed','obj')=z.l;
results('relaxed','total constraints') = match.numequ;


loop(cycle,

*
* initialization
*
   tour(t,i,j) = no;
   s(i,j) = no;
   first(i,j) = no;
   reach(i,j) = no;
   done = 0;

   tt(t) = no;
   tt(t)$(ord(t)=1) = 1;         // initialize tt(1) = yes

   s(i,j)=yes$(x.l(i,j) > 0.5);  // initialize to current solution

   while(not done,

*
* pick any (i,j) from remaining solutions
*

      first(i,j)=no;
      loop((i,j)$(card(first)=0),
         first(i,j) = s(i,j);
      );
*   display first;

*
* this is he beginning of a new subtour
*
      tour(tt,first) = yes;
      s(first) = no;

      continue = 1;
      while(continue,
*
* find (i,j)'s connected to tour(tt)
* note that the loop over tt is just syntax: tt contains one element
*
         reach(i,j) = no;
         loop((tt,k),
            reach(s(i,j))$tour(tt,k,j) = yes;
            reach(s(i,j))$tour(tt,i,k) = yes;
            reach(s(i,j))$tour(tt,k,i) = yes;
            reach(s(i,j))$tour(tt,j,k) = yes;
         );
*       display tour,s,reach;


*
* if reach = empty then we can stop the while loop
*
         if (card(reach)=0, continue = 0);

*
* add them to the current subtour
* note that the loop over tt is just syntax: tt contains one element
*
         loop(tt,
            tour(tt,reach) = yes;
         );
         s(reach) = no;
      );

*   display tour;

*
* if no remaining solutions, we are done
*
      if (card(s)=0, done=1);

*
* new subtour
*
      tt(t) = tt(t-1);   // t := t + 1

   );


   cutindex(cycle,t)$(sum(tour(t,i,j),1)>0.5) = yes;
*
* check
*
   option tour:0:1:1;
   display tour;
   loop(cutindex(cycle,t),
      abort$(sum(tour(t,i,j),1) < 2.5) "subtour with 1 or 2 points";
   );

   if (sum(cutindex(cycle,t),1) > 1.5,

      cutcoeff(cycle,tour(t,i,j))$cutindex(cycle,t) = 1;
      cutlength(cutindex(cycle,t)) = sum(tour(t,i,j),1);

      option cutindex:0:1:1;
      option cutcoeff:0:1:1;
*     display cutindex,cutcoeff,cutlength;

      solve tsp minimizing z using mip;
      if (tsp.modelstat<>1,
          put log,  "**** ERROR: Model not solved to optimality"/;
          abort$1  "Model not solved to optimality";
      );
      display x.l;

      results(cycle,'obj') = z.l;
      results(cycle,'cuts added') = sum(cutindex(cycle,t), 1);
      results(cycle,'total constraints') = tsp.numequ;
*     display results;

   else

      display "Optimal solution found!";
      done = 1;

   );
);


display results;


$if exist solution.csv $call del solution.csv
file sol /solution.csv/;
sol.pc = 5;
put sol;
loop((i,j)$(x.l(i,j)>0.5),
   put i.tl,j.tl/;
);


largetext:
The blob contains the gamsinterrupt.dll.


