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;
``````
GAMS Development Corp.
GAMS Software GmbH

General Information and Sales
U.S. (+1) 202 342-0180
Europe: (+49) 221 949-9170