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