$title Dantzig Wolfe Decomposition and Grid Computing (DANWOLFE,SEQ=328)

$onText
This example implements the well known Dantzig Wolfe decomposition method. This
method exploits the structure in linear programs where the coefficient matrix
has the special form:

     B0  B1  B2 ... Bk
         A1
             A2
                .
                 .  Ak

The specific model is a multi-commodity network flow problem where Ak
corresponds to a commodity flow and Bk represents arc capacities.


Dantzig, G B, and Wolfe, P, Decomposition principle for linear
programs. Operations Research 8 (1960), 101-111.

Keywords: linear programming, Dantzig Wolfe decomposition, multi-commodity network flow problem,
          network optimization, grid computing
$offText

$eolCom //

$setDDList nodes comm maxtime
$if not set nodes   $set nodes    20
$if not set comm    $set comm      5
$if not set maxtime $set maxtime 100
$if not errorfree $abort wrong double dash parameters: --nodes=n --comm=n --maxtime=secs

Set
   i      'nodes'       / n1*n%nodes% /
   k      'commodities' / k1*k%comm%  /
   e(i,i) 'edges';

Alias (i,j);

Parameter
   cost(i,j) 'cost for edge use'
   bal(k,i)  'balance'
   kdem(k)   'demand'
   cap(i,j)  'bundle capacity';

Variable
   x(k,i,j) 'multi commodity flow'
   z        'objective';

Positive Variable x;

Equation
   defbal(k,i) 'balancing constraint'
   defcap(i,j) 'bundling capacity'
   defobj;

defobj..      z =e= sum((k,e), cost(e)*x(k,e));

defbal(k,i).. sum(e(i,j), x(k,e)) - sum(e(j,i),x(k,e)) =e= bal(k,i);

defcap(e)..   sum(k, x(k,e)) =l= cap(e);

Model mcf 'multi-commodity flow problem' / all /;

**** make a random instance
Scalar
   inum
   edgedensity / 0.3 /;

e(i,j)  = uniform(0,1) < edgedensity;
e(i,i)  = no;
cost(e) = uniform(1,10);
cap(e)  = uniform(50,100)*log(card(k));

loop(k,
   kdem(k) = uniform(50,150);
   inum    = uniformInt(1,card(i)); bal(k,i)$(ord(i) = inum) = kdem(k);
   inum    = uniformInt(1,card(i)); bal(k,i)$(ord(i) = inum) = bal(k,i) - kdem(k);
   kdem(k) = sum(i$(bal(k,i) > 0), bal(k,i));
);

**** see if the random model is feasible
option limRow = 0, limCol = 0, solPrint = off, solveLink = %solveLink.callModule%;

solve mcf min z using lp;

abort$(mcf.modelStat <> %modelStat.optimal%) 'problem not feasible. Increase edge density.'

Parameter xsingle(k,i,j) 'single solve';
xsingle(k,e) = x.l(k,e)$[x.l(k,e) > 1e-6];

display$(card(i) < 30) xsingle;

**** define master model
Set
   p           'paths idents' / p1*p100 /
   ap(k,p)     'active path'
   pe(k,p,i,j) 'edge path incidence vector';

Parameter pcost(k,p) 'path cost';

Positive Variable
   xp(k,p)
   slack(k);

Equation
   mdefcap(i,j) 'bundle constraint'
   mdefbal(k)   'balance constraint'
   mdefobj      'objective';

mdefobj..    z =e= sum(ap, pcost(ap)*xp(ap)) + sum(k, 999*slack(k));

mdefbal(k).. sum(ap(k,p), xp(ap)) + slack(k) =e= kdem(k);

mdefcap(e).. sum(pe(ap,e), xp(ap)) =l= cap(e);

Model master / mdefobj, mdefbal, mdefcap /;

**** define pricing model: shortest path
Parameter ebal(i);

Positive Variable xe(i,j);

Equation
   pdefbal(i) 'balance constraint'
   pdefobj    'objective';

pdefobj..    z =e= sum(e, (cost(e) - mdefcap.m(e))*xe(e));

pdefbal(i).. sum(e(i,j), xe(e)) - sum(e(j,i),xe(e)) =e= ebal(i);

Model pricing / pdefobj, pdefbal /;

**** serial decomposition method
Scalar
   done 'loop indicator'
   iter 'iteration counter';

Set nextp(k,p) 'next path to be added';

* clear path data
done          =   0;
iter          =   0;
ap(k,p)       =  no;
pe(k,p,e)     =  no;
pcost(k,p)    =   0;
nextp(k,p)    =  no;
nextp(k,'p1') = yes;

while(not done,
   iter = iter + 1;
   solve master using lp minimizing z;
   done = 1;
   loop(k$kdem(k),
      ebal(i) = bal(k,i)/kdem(k);
      solve pricing using lp minimizing z;
      pricing.solPrint = %solPrint.quiet%;  // turn off all outputs fpr pricing model
      if(mdefbal.m(k) - z.l > 1e-6,         // add new path
         ap(nextp(k,p))    = yes;
         pe(nextp(k,p),e)  = round(xe.l(e));
         pcost(nextp(k,p)) = sum(pe(nextp,e), cost(e));
         nextp(k,p)        = nextp(k,p-1);  // bump the path to the next free one
         abort$(sum(nextp(k,p),1) = 0) 'set p too small';
         done = 0;
      );
   );
);

abort$(abs(master.objVal - mcf.objVal) > 1e-3) 'different objective values', master.objVal, mcf.objVal;

Parameter xserial(k,i,j);
xserial(k,e) = sum(pe(ap(k,p),e), xp.l(ap));

display$(card(i) < 30) xserial;

**** Now the same with the Grid option
Parameter h(k) 'model handles';
pricing.solveLink = %solveLink.asyncGrid%;
pricing.solPrint  = %solPrint.summary%;

* clear path data
done          =   0;
iter          =   0;
ap(k,p)       =  no;
pe(k,p,e)     =  no;
pcost(k,p)    =   0;
nextp(k,p)    =  no;
nextp(k,'p1') = yes;

while(not done,
   iter = iter + 1;
   solve master using lp minimizing z;
   done = 1;
   pricing.number = 0;
   loop(k$kdem(k),
      ebal(i) = bal(k,i)/kdem(k);
      solve pricing using lp minimizing z;
      pricing.solPrint = %solPrint.quiet%; // turn off all outputs for pricing model
      h(k) = pricing.handle;
   );
   repeat
      loop(k$h(k),
         if(handlecollect(h(k)),
            if(mdefbal.m(k) - z.l > 1e-6,
               ap(nextp(k,p))    = yes;
               pe(nextp(k,p),e)  = round(xe.l(e));
               pcost(nextp(k,p)) = sum(pe(nextp,e), cost(e));
               nextp(k,p)        = nextp(k,p-1);
               abort$(sum(nextp(k,p),1) = 0) 'set p too small';
               done = 0;
            );
            display$handledelete(h(k)) 'trouble deleting handles' ;
            h(k) = 0;
         );
      );
      display$sleep(card(h)*0.1) 'sleep a bit';
   until card(h) = 0 or timeelapsed > %maxtime%;
   abort$(card(h) and timeelapsed <= %maxtime%) 'not all solves returned', h;
);

if(timeelapsed > %maxtime%,
   display 'Algorithm interrupted because maxtime=%maxtime% has been exceeded!';
);

abort$(abs(master.objVal - mcf.objVal) > 1e-3 and timeelapsed <= %maxtime%) 'different objective values', master.objVal, mcf.objVal;

Parameter xgrid(k,i,j);
xgrid(k,e) = sum(pe(ap(k,p),e), xp.l(ap));

display$(card(i) < 30) xgrid;

Parameter xall 'summary of flows';
xall(k,e,'single') = xsingle(k,e);
xall(k,e,'serial') = xserial(k,e);
xall(k,e,'grid')   = xgrid(k,e);

option  xall:3:3:1;
display xall;