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;