# ---------------------------------------- # DANTZIG-WOLFE DECOMPOSITION FOR # MULTI-COMMODITY TRANSPORTATION # ---------------------------------------- model multi2.mod; data multi2.dat; param nIter default 0; let {p in PROD} nPROP[p] := 0; let {p in PROD} price_convex[p] := 1; let {i in ORIG, j in DEST} price[i,j] := 0; option solver minos; option omit_zero_rows 1; option display_1col 0; option display_eps .000001; option presolve_eps 2e-13; # ---------------------------------------------------------- problem MasterI: Artificial, Weight, Excess, Multi, Convex; problem SubI {p in PROD}: Artif_Reduced_Cost[p], {i in ORIG, j in DEST} Trans[i,j,p], {i in ORIG} Supply[i,p], {j in DEST} Demand[j,p]; repeat { let nIter := nIter + 1; printf "\nPHASE I -- ITERATION %d\n", nIter; for {p in PROD} { printf "\nPRODUCT %s\n\n", p; solve SubI[p]; display Artif_Reduced_Cost[p]; printf "\n"; display {i in ORIG, j in DEST} Trans[i,j,p]; if Artif_Reduced_Cost[p] < - 0.00001 then { let nPROP[p] := nPROP[p] + 1; let {i in ORIG, j in DEST} prop_ship[i,j,p,nPROP[p]] := Trans[i,j,p]; let prop_cost[p,nPROP[p]] := sum {i in ORIG, j in DEST} cost[i,j,p] * Trans[i,j,p]; }; }; if min {p in PROD} Artif_Reduced_Cost[p] >= - 0.00001 then { printf "\n*** NO FEASIBLE SOLUTION ***\n"; break; }; solve MasterI; printf "\n"; display Weight; display Multi.dual; display {i in ORIG, j in DEST} limit[i,j] - sum {p in PROD, k in 1..nPROP[p]} prop_ship[i,j,p,k] * Weight[p,k]; if Excess <= 0.00001 then break; else { let {i in ORIG, j in DEST} price[i,j] := Multi[i,j].dual; let {p in PROD} price_convex[p] := Convex[p].dual; }; }; # ---------------------------------------------------------- printf "\nSETTING UP FOR PHASE II\n"; problem MasterII: Total_Cost, Weight, Multi, Convex; problem SubII {p in PROD}: Reduced_Cost[p], {i in ORIG, j in DEST} Trans[i,j,p], {i in ORIG} Supply[i,p], {j in DEST} Demand[j,p]; solve MasterII; printf "\n"; display Weight; display Multi.dual; display Multi.slack; let {i in ORIG, j in DEST} price[i,j] := Multi[i,j].dual; let {p in PROD} price_convex[p] := Convex[p].dual; repeat { let nIter := nIter + 1; printf "\nPHASE II -- ITERATION %d\n\n", nIter; for {p in PROD} { printf "\nPRODUCT %s\n\n", p; solve SubII[p]; display Reduced_Cost[p]; printf "\n"; display {i in ORIG, j in DEST} Trans[i,j,p]; if Reduced_Cost[p] < - 0.00001 then { let nPROP[p] := nPROP[p] + 1; let {i in ORIG, j in DEST} prop_ship[i,j,p,nPROP[p]] := Trans[i,j,p]; let prop_cost[p,nPROP[p]] := sum {i in ORIG, j in DEST} cost[i,j,p] * Trans[i,j,p]; }; }; if min {p in PROD} Reduced_Cost[p] >= - 0.00001 then break; solve MasterII; printf "\n"; display Weight; let {i in ORIG, j in DEST} price[i,j] := Multi[i,j].dual; let {p in PROD} price_convex[p] := Convex[p].dual; }; # ---------------------------------------------------------- printf "\nPHASE III\n"; let {i in ORIG, j in DEST, p in PROD} Trans[i,j,p] := sum {k in 1..nPROP[p]} prop_ship[i,j,p,k] * Weight[p,k]; param true_Total_Cost := sum {i in ORIG, j in DEST, p in PROD} cost[i,j,p] * Trans[i,j,p].val; printf "\n"; display true_Total_Cost; display Trans;