$ontext

   Flow formulation of minimum spanning tree problem.
   This is MIP but automatically integer so we solve as RMIP.

   Erwin Kalvelagen, april, 2002, revised 2009

   References:
     Magnanti, Thomas L.; Wolsey, Laurence A., "Optimal Trees",
     in M.O. Ball, T.L. Magnanti, C.L. Monma, and G.L. Nemhauser, editors,
     Network Models, volume 7 of Handbooks in Operations Research and Management Science,
     Chapter 9, pages 503--616. North-Holland, 1995.


$offtext

set i /i1*i51/;



alias (i,j,k);


table xy(i,*)
          x     y
    i1   37    52
    i2   49    49
    i3   52    64
    i4   20    26
    i5   40    30
    i6   21    47
    i7   17    63
    i8   31    62
    i9   52    33
   i10   51    21
   i11   42    41
   i12   31    32
   i13    5    25
   i14   12    42
   i15   36    16
   i16   52    41
   i17   27    23
   i18   17    33
   i19   13    13
   i20   57    58
   i21   62    42
   i22   42    57
   i23   16    57
   i24    8    52
   i25    7    38
   i26   27    68
   i27   30    48
   i28   43    67
   i29   58    48
   i30   58    27
   i31   37    69
   i32   38    46
   i33   46    10
   i34   61    33
   i35   62    63
   i36   63    69
   i37   32    22
   i38   45    35
   i39   59    15
   i40    5     6
   i41   10    17
   i42   21    10
   i43    5    64
   i44   30    15
   i45   39    10
   i46   32    39
   i47   25    32
   i48   25    55
   i49   48    28
   i50   56    37
   i51   30    40

;

parameter c(i,j);


set arcs(i,j);
arcs(i,j)$(not sameas(i,j)) = yes;

c(arcs(i,j)) = round(sqrt(sqr(xy(i,'x')-xy(j,'x'))+sqr(xy(i,'y')-xy(j,'y'))));



scalar n 'number of nodes';
n = card(i);


display xy;

alias(i,t);

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

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


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/;
solve m minimizing z using rmip;

display y.l, x.l;

parameter tree(i,*);
loop((i,j)$(x.l(i,j)>0.5),
   tree(i,'x0') = xy(i,'x');
   tree(i,'y0') = xy(i,'y');
   tree(i,'x1') = xy(j,'x');
   tree(i,'y1') = xy(j,'y');
);

display tree;