danwolfe.gms : Dantzig Wolfe Decomposition and Grid Computing

Description

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.


Large Model of Type : LP


Category : GAMS Model library


Main file : danwolfe.gms

$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;