Implementing algorithms through AMPL scripts
Gilmore-Gomory column generation procedure for the cutting-stock (roll trim) problem
# ----------------------------------------
# GILMORE-GOMORY METHOD FOR
# CUTTING STOCK PROBLEM
# ----------------------------------------
option solver cplex;
option solution_round 6;
model cut1.mod;
data cut.dat;
problem Cutting_Opt: Cut, Number, Fill;
option relax_integrality 1;
option presolve 0;
problem Pattern_Gen: Use, Reduced_Cost, Width_Limit;
option relax_integrality 0;
option presolve 1;
let nPAT := 0;
for {i in WIDTHS} {
let nPAT := nPAT + 1;
let nbr[i,nPAT] := floor (roll_width/i);
let {i2 in WIDTHS: i2 <> i} nbr[i2,nPAT] := 0;
};
repeat {
solve Cutting_Opt;
let {i in WIDTHS} price[i] := Fill[i].dual;
solve Pattern_Gen;
if Reduced_Cost < -0.00001 then {
let nPAT := nPAT + 1;
let {i in WIDTHS} nbr[i,nPAT] := Use[i];
}
else break;
};
display nbr;
display Cut;
option Cutting_Opt.relax_integrality 0;
option Cutting_Opt.presolve 10;
solve Cutting_Opt;
display Cut;# ----------------------------------------
# CUTTING STOCK USING PATTERNS
# ----------------------------------------
param roll_width > 0; # width of raw rolls
set WIDTHS; # set of widths to be cut
param orders {WIDTHS} > 0; # number of each width to be cut
param nPAT integer >= 0; # number of patterns
set PATTERNS := 1..nPAT; # set of patterns
param nbr {WIDTHS,PATTERNS} integer >= 0;
check {j in PATTERNS}:
sum {i in WIDTHS} i * nbr[i,j] <= roll_width;
# defn of patterns: nbr[i,j] = number
# of rolls of width i in pattern j
var Cut {PATTERNS} integer >= 0; # rolls cut using each pattern
minimize Number: # minimize total raw rolls cut
sum {j in PATTERNS} Cut[j];
subj to Fill {i in WIDTHS}:
sum {j in PATTERNS} nbr[i,j] * Cut[j] >= orders[i];
# for each width, total
# rolls cut meets total orders
# ----------------------------------------
# KNAPSACK SUBPROBLEM FOR CUTTING STOCK
# ----------------------------------------
param price {WIDTHS} default 0.0;
var Use {WIDTHS} integer >= 0;
minimize Reduced_Cost:
1 - sum {i in WIDTHS} price[i] * Use[i];
subj to Width_Limit:
sum {i in WIDTHS} i * Use[i] <= roll_width;param roll_width := 110 ;
param: WIDTHS: orders :=
20 48
45 35
50 24
55 10
75 8 ;Same as cut1.run, but using an alternative arrangement wherein problems are defined immediately before before their members are declared
# ----------------------------------------
# GILMORE-GOMORY METHOD FOR
# CUTTING STOCK PROBLEM
# ----------------------------------------
option solver cplex;
option solution_round 6;
model cut2.mod;
data cut.dat;
problem Cutting_Opt;
option relax_integrality 1;
option presolve 0;
problem Pattern_Gen;
option relax_integrality 0;
option presolve 1;
let nPAT := 0;
for {i in WIDTHS} {
let nPAT := nPAT + 1;
let nbr[i,nPAT] := floor (roll_width/i);
let {i2 in WIDTHS: i2 <> i} nbr[i2,nPAT] := 0;
};
repeat {
solve Cutting_Opt;
let {i in WIDTHS} price[i] := Fill[i].dual;
solve Pattern_Gen;
if Reduced_Cost < -0.00001 then {
let nPAT := nPAT + 1;
let {i in WIDTHS} nbr[i,nPAT] := Use[i];
}
else break;
};
display nbr;
display Cut;
option Cutting_Opt.relax_integrality 0;
option Cutting_Opt.presolve 10;
solve Cutting_Opt;
display Cut; problem Cutting_Opt;
# ----------------------------------------
param nPAT integer >= 0, default 0;
param roll_width;
set PATTERNS := 1..nPAT;
set WIDTHS;
param orders {WIDTHS} > 0;
param nbr {WIDTHS,PATTERNS} integer >= 0;
check {j in PATTERNS}: sum {i in WIDTHS} i * nbr[i,j] <= roll_width;
var Cut {PATTERNS} integer >= 0;
minimize Number: sum {j in PATTERNS} Cut[j];
subj to Fill {i in WIDTHS}:
sum {j in PATTERNS} nbr[i,j] * Cut[j] >= orders[i];
problem Pattern_Gen;
# ----------------------------------------
param price {WIDTHS} default 0;
var Use {WIDTHS} integer >= 0;
minimize Reduced_Cost:
1 - sum {i in WIDTHS} price[i] * Use[i];
subj to Width_Limit:
sum {i in WIDTHS} i * Use[i] <= roll_width;param roll_width := 110 ;
param: WIDTHS: orders :=
20 48
45 35
50 24
55 10
75 8 ;Same as cut1.run, but with better formatting of output
# ----------------------------------------
# GILMORE-GOMORY METHOD FOR
# CUTTING STOCK PROBLEM
# ----------------------------------------
option solver cplex;
option solution_round 6;
option solver_msg 0;
model cut1.mod;
data cut.dat;
problem Cutting_Opt: Cut, Number, Fill;
option relax_integrality 1;
option presolve 0;
problem Pattern_Gen: Use, Reduced_Cost, Width_Limit;
option relax_integrality 0;
option presolve 1;
let nPAT := 0;
for {i in WIDTHS} {
let nPAT := nPAT + 1;
let nbr[i,nPAT] := floor (roll_width/i);
let {i2 in WIDTHS: i2 <> i} nbr[i2,nPAT] := 0;
};
repeat {
solve Cutting_Opt >/dev/null;
let {i in WIDTHS} price[i] := Fill[i].dual;
solve Pattern_Gen >/dev/null;
printf "\n%7.2f%11.2e ", Number, Reduced_Cost;
if Reduced_Cost < -0.00001 then {
let nPAT := nPAT + 1;
let {i in WIDTHS} nbr[i,nPAT] := Use[i];
}
else break;
for {i in WIDTHS} printf "%3i", Use[i];
};
let {j in PATTERNS} Cut[j] := ceil(Cut[j]);
printf "\n\n\nRounded up to integer: %3i rolls", sum {j in PATTERNS} Cut[j];
printf "\n\nCut ";
printf {j in PATTERNS}: "%4i", Cut[j];
printf "\n\n";
for {i in WIDTHS} {
printf "%3i ", i;
printf {j in PATTERNS}: "%4i", nbr[i,j];
printf "\n";
}
printf "\nWASTE = %5.2f%%\n\n\n",
100 * (1 - (sum {i in WIDTHS} i * orders[i]) / (roll_width * Number));
option Cutting_Opt.relax_integrality 0;
option Cutting_Opt.presolve 10;
solve Cutting_Opt >/dev/null;
printf "Best integer: %3i rolls", sum {j in PATTERNS} Cut[j];
printf "\n\nCut ";
printf {j in PATTERNS}: "%4i", Cut[j];
printf "\n\n";
for {i in WIDTHS} {
printf "%3i ", i;
printf {j in PATTERNS}: "%4i", nbr[i,j];
printf "\n";
}
printf "\nWASTE = %5.2f%%\n\n",
100 * (1 - (sum {i in WIDTHS} i * orders[i]) / (roll_width * Number));# ----------------------------------------
# CUTTING STOCK USING PATTERNS
# ----------------------------------------
param roll_width > 0; # width of raw rolls
set WIDTHS; # set of widths to be cut
param orders {WIDTHS} > 0; # number of each width to be cut
param nPAT integer >= 0; # number of patterns
set PATTERNS := 1..nPAT; # set of patterns
param nbr {WIDTHS,PATTERNS} integer >= 0;
check {j in PATTERNS}:
sum {i in WIDTHS} i * nbr[i,j] <= roll_width;
# defn of patterns: nbr[i,j] = number
# of rolls of width i in pattern j
var Cut {PATTERNS} integer >= 0; # rolls cut using each pattern
minimize Number: # minimize total raw rolls cut
sum {j in PATTERNS} Cut[j];
subj to Fill {i in WIDTHS}:
sum {j in PATTERNS} nbr[i,j] * Cut[j] >= orders[i];
# for each width, total
# rolls cut meets total orders
# ----------------------------------------
# KNAPSACK SUBPROBLEM FOR CUTTING STOCK
# ----------------------------------------
param price {WIDTHS} default 0.0;
var Use {WIDTHS} integer >= 0;
minimize Reduced_Cost:
1 - sum {i in WIDTHS} price[i] * Use[i];
subj to Width_Limit:
sum {i in WIDTHS} i * Use[i] <= roll_width;param roll_width := 110 ;
param: WIDTHS: orders :=
20 48
45 35
50 24
55 10
75 8 ;Dantzig-Wolfe decomposition for a multi-commodity transportation problem, using a single subproblem
# ----------------------------------------
# DANTZIG-WOLFE DECOMPOSITION FOR
# MULTI-COMMODITY TRANSPORTATION
# ----------------------------------------
model multi1.mod;
data multi1.dat;
let nPROP := 0;
let price_convex := 1;
let {i in ORIG, j in DEST} price[i,j] := 0;
option solver minos;
option omit_zero_rows 1;
option display_1col 10;
option display_eps .000001;
# ----------------------------------------------------------
problem MasterI: Artificial, Weight, Excess, Multi, Convex;
problem SubI: Artif_Reduced_Cost, Trans, Supply, Demand;
repeat { printf "\nPHASE I -- ITERATION %d\n\n", nPROP+1;
solve SubI;
printf "\n";
display Trans;
if Artif_Reduced_Cost >= - 0.00001 then {
printf "\n*** NO FEASIBLE SOLUTION ***\n";
break;
}
else {
let nPROP := nPROP + 1;
let {i in ORIG, j in DEST}
prop_ship[i,j,nPROP] := sum {p in PROD} Trans[i,j,p];
let prop_cost[nPROP] :=
sum {i in ORIG, j in DEST, p in PROD} cost[i,j,p] * Trans[i,j,p];
};
solve MasterI;
printf "\n";
display Weight; display Multi.dual;
display {i in ORIG, j in DEST}
limit[i,j] - sum {k in 1..nPROP} prop_ship[i,j,k] * Weight[k];
if Excess <= 0.00001 then break;
else {
let {i in ORIG, j in DEST} price[i,j] := Multi[i,j].dual;
let price_convex := Convex.dual;
};
};
# ----------------------------------------------------------
printf "\nSETTING UP FOR PHASE II\n\n";
problem MasterII: Total_Cost, Weight, Multi, Convex;
problem SubII: Reduced_Cost, Trans, Supply, Demand;
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 price_convex := Convex.dual;
repeat { printf "\nPHASE II -- ITERATION %d\n\n", nPROP+1;
solve SubII;
printf "\n";
display Trans;
if Reduced_Cost >= - 0.00001 then {
printf "\n*** OPTIMAL SOLUTION ***\n";
break;
}
else {
let nPROP := nPROP + 1;
let {i in ORIG, j in DEST}
prop_ship[i,j,nPROP] := sum {p in PROD} Trans[i,j,p];
let prop_cost[nPROP] :=
sum {i in ORIG, j in DEST, p in PROD} cost[i,j,p] * Trans[i,j,p];
};
solve MasterII;
printf "\n";
display Weight;
let {i in ORIG, j in DEST} price[i,j] := Multi[i,j].dual;
let price_convex := Convex.dual;
};
# ----------------------------------------------------------
printf "\nPHASE III\n\n";
problem MasterIII: Opt_Cost, Trans, Supply, Demand, Opt_Multi;
let {i in ORIG, j in DEST}
opt_ship[i,j] := sum {k in 1..nPROP} prop_ship[i,j,k] * Weight[k];
solve MasterIII;
printf "\n";
display Trans;# ----------------------------------------
# MULTI-COMMODITY FLOW USING
# DANTZIG-WOLFE DECOMPOSITION
# ----------------------------------------
### SUBPROBLEM ###
set ORIG; # origins
set DEST; # destinations
set PROD; # products
param supply {ORIG,PROD} >= 0; # amounts available at origins
param demand {DEST,PROD} >= 0; # amounts required at destinations
check {p in PROD}:
sum {i in ORIG} supply[i,p] = sum {j in DEST} demand[j,p];
param price_convex; # dual price on convexity constr
param price {ORIG,DEST} <= 0.000001; # dual price on shipment limit
param cost {ORIG,DEST,PROD} >= 0; # shipment costs per unit
var Trans {ORIG,DEST,PROD} >= 0; # units to be shipped
minimize Artif_Reduced_Cost:
sum {i in ORIG, j in DEST, p in PROD}
(- price[i,j]) * Trans[i,j,p] - price_convex;
minimize Reduced_Cost:
sum {i in ORIG, j in DEST, p in PROD}
(cost[i,j,p] - price[i,j]) * Trans[i,j,p] - price_convex;
subject to Supply {i in ORIG, p in PROD}:
sum {j in DEST} Trans[i,j,p] = supply[i,p];
subject to Demand {j in DEST, p in PROD}:
sum {i in ORIG} Trans[i,j,p] = demand[j,p];
### MASTER PROBLEM ###
param limit {ORIG,DEST} >= 0; # max shipped on each link
param nPROP integer >= 0;
param prop_ship {ORIG,DEST,1..nPROP} >= -0.000001;
param prop_cost {1..nPROP} >= 0;
# For each proposal from the subproblem:
# amount it ships over each link, and its cost
var Weight {1..nPROP} >= 0;
var Excess >= 0;
minimize Artificial: Excess;
minimize Total_Cost:
sum {k in 1..nPROP} prop_cost[k] * Weight[k];
subject to Multi {i in ORIG, j in DEST}:
sum {k in 1..nPROP} prop_ship[i,j,k] * Weight[k] - Excess <= limit[i,j];
subject to Convex: sum {k in 1..nPROP} Weight[k] = 1;
### PHASE III PROBLEM ###
param opt_ship {ORIG,DEST} >= -0.000001;
minimize Opt_Cost:
sum {i in ORIG, j in DEST, p in PROD} cost[i,j,p] * Trans[i,j,p];
subject to Opt_Multi {i in ORIG, j in DEST}:
sum {p in PROD} Trans[i,j,p] = opt_ship[i,j];set ORIG := GARY CLEV PITT ;
set DEST := FRA DET LAN WIN STL FRE LAF ;
set PROD := bands coils plate ;
param supply (tr): GARY CLEV PITT :=
bands 400 700 800
coils 800 1600 1800
plate 200 300 300 ;
param demand (tr):
FRA DET LAN WIN STL FRE LAF :=
bands 300 300 100 75 650 225 250
coils 500 750 400 250 950 850 500
plate 100 100 0 50 200 100 250 ;
param limit default 625 ;
param cost :=
[*,*,bands]: FRA DET LAN WIN STL FRE LAF :=
GARY 30 10 8 10 11 71 6
CLEV 22 7 10 7 21 82 13
PITT 19 11 12 10 25 83 15
[*,*,coils]: FRA DET LAN WIN STL FRE LAF :=
GARY 39 14 11 14 16 82 8
CLEV 27 9 12 9 26 95 17
PITT 24 14 17 13 28 99 20
[*,*,plate]: FRA DET LAN WIN STL FRE LAF :=
GARY 41 15 12 16 17 86 8
CLEV 29 9 13 9 28 99 18
PITT 26 14 17 13 31 104 20 ;Same as multi1.run, but using the same repeat loop for both phase I (infeasible) and phase II (feasible).
# ----------------------------------------
# DANTZIG-WOLFE DECOMPOSITION FOR
# MULTI-COMMODITY TRANSPORTATION
# ----------------------------------------
model multi1.mod;
data multi1.dat;
let nPROP := 0;
let price_convex := 1;
let {i in ORIG, j in DEST} price[i,j] := 0;
param phase symbolic, default "I";
option solver minos;
option omit_zero_rows 1;
option display_1col 10;
option display_eps .000001;
# ----------------------------------------------------------
problem Master: Artificial, Weight, Excess, Multi, Convex;
problem Sub: Artif_Reduced_Cost, Trans, Supply, Demand;
repeat { printf "\nPHASE %s -- ITERATION %d\n\n", phase, nPROP+1;
solve Sub;
printf "\n";
display Trans;
if phase = "I" then {
if Artif_Reduced_Cost >= - 0.00001 then {
printf "\n*** NO FEASIBLE SOLUTION ***\n";
break;
}
}
else {
if Reduced_Cost >= - 0.00001 then {
printf "\n*** OPTIMAL SOLUTION ***\n";
break;
}
};
let nPROP := nPROP + 1;
let {i in ORIG, j in DEST}
prop_ship[i,j,nPROP] := sum {p in PROD} Trans[i,j,p];
let prop_cost[nPROP] :=
sum {i in ORIG, j in DEST, p in PROD} cost[i,j,p] * Trans[i,j,p];
solve Master;
printf "\n";
display Weight;
if phase = "I" then {
display Multi.dual;
display {i in ORIG, j in DEST}
limit[i,j] - sum {k in 1..nPROP} prop_ship[i,j,k] * Weight[k];
if Excess <= 0.00001 then {
printf "\nSETTING UP FOR PHASE II\n\n";
let phase := "II";
problem Master;
drop Artificial; restore Total_Cost; fix Excess;
problem Sub;
drop Artif_Reduced_Cost; restore Reduced_Cost;
# problem Master; show obj; solve;
solve Master;
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 price_convex := Convex.dual;
};
# ----------------------------------------------------------
printf "\nPHASE III\n\n";
problem MasterIII: Opt_Cost, Trans, Supply, Demand, Opt_Multi;
let {i in ORIG, j in DEST}
opt_ship[i,j] := sum {k in 1..nPROP} prop_ship[i,j,k] * Weight[k];
solve MasterIII;
printf "\n";
display Trans;# ----------------------------------------
# MULTI-COMMODITY FLOW USING
# DANTZIG-WOLFE DECOMPOSITION
# ----------------------------------------
### SUBPROBLEM ###
set ORIG; # origins
set DEST; # destinations
set PROD; # products
param supply {ORIG,PROD} >= 0; # amounts available at origins
param demand {DEST,PROD} >= 0; # amounts required at destinations
check {p in PROD}:
sum {i in ORIG} supply[i,p] = sum {j in DEST} demand[j,p];
param price_convex; # dual price on convexity constr
param price {ORIG,DEST} <= 0.000001; # dual price on shipment limit
param cost {ORIG,DEST,PROD} >= 0; # shipment costs per unit
var Trans {ORIG,DEST,PROD} >= 0; # units to be shipped
minimize Artif_Reduced_Cost:
sum {i in ORIG, j in DEST, p in PROD}
(- price[i,j]) * Trans[i,j,p] - price_convex;
minimize Reduced_Cost:
sum {i in ORIG, j in DEST, p in PROD}
(cost[i,j,p] - price[i,j]) * Trans[i,j,p] - price_convex;
subject to Supply {i in ORIG, p in PROD}:
sum {j in DEST} Trans[i,j,p] = supply[i,p];
subject to Demand {j in DEST, p in PROD}:
sum {i in ORIG} Trans[i,j,p] = demand[j,p];
### MASTER PROBLEM ###
param limit {ORIG,DEST} >= 0; # max shipped on each link
param nPROP integer >= 0;
param prop_ship {ORIG,DEST,1..nPROP} >= -0.000001;
param prop_cost {1..nPROP} >= 0;
# For each proposal from the subproblem:
# amount it ships over each link, and its cost
var Weight {1..nPROP} >= 0;
var Excess >= 0;
minimize Artificial: Excess;
minimize Total_Cost:
sum {k in 1..nPROP} prop_cost[k] * Weight[k];
subject to Multi {i in ORIG, j in DEST}:
sum {k in 1..nPROP} prop_ship[i,j,k] * Weight[k] - Excess <= limit[i,j];
subject to Convex: sum {k in 1..nPROP} Weight[k] = 1;
### PHASE III PROBLEM ###
param opt_ship {ORIG,DEST} >= -0.000001;
minimize Opt_Cost:
sum {i in ORIG, j in DEST, p in PROD} cost[i,j,p] * Trans[i,j,p];
subject to Opt_Multi {i in ORIG, j in DEST}:
sum {p in PROD} Trans[i,j,p] = opt_ship[i,j];set ORIG := GARY CLEV PITT ;
set DEST := FRA DET LAN WIN STL FRE LAF ;
set PROD := bands coils plate ;
param supply (tr): GARY CLEV PITT :=
bands 400 700 800
coils 800 1600 1800
plate 200 300 300 ;
param demand (tr):
FRA DET LAN WIN STL FRE LAF :=
bands 300 300 100 75 650 225 250
coils 500 750 400 250 950 850 500
plate 100 100 0 50 200 100 250 ;
param limit default 625 ;
param cost :=
[*,*,bands]: FRA DET LAN WIN STL FRE LAF :=
GARY 30 10 8 10 11 71 6
CLEV 22 7 10 7 21 82 13
PITT 19 11 12 10 25 83 15
[*,*,coils]: FRA DET LAN WIN STL FRE LAF :=
GARY 39 14 11 14 16 82 8
CLEV 27 9 12 9 26 95 17
PITT 24 14 17 13 28 99 20
[*,*,plate]: FRA DET LAN WIN STL FRE LAF :=
GARY 41 15 12 16 17 86 8
CLEV 29 9 13 9 28 99 18
PITT 26 14 17 13 31 104 20 ;Same as multi1.run, but using a separate subproblem for each product; subproblems are represented in AMPL by an indexed collection of named problems
# ----------------------------------------
# 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 omit_zero_rows 1;
option display_1col 0;
option display_eps .000001;
# ----------------------------------------------------------
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];
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];
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;# ----------------------------------------
# MULTI-COMMODITY FLOW USING
# DANTZIG-WOLFE DECOMPOSITION
# (one subproblem for each product)
# ----------------------------------------
### SUBPROBLEM ###
set ORIG; # origins
set DEST; # destinations
set PROD; # products
param supply {ORIG,PROD} >= 0; # amounts available at origins
param demand {DEST,PROD} >= 0; # amounts required at destinations
check {p in PROD}:
sum {i in ORIG} supply[i,p] = sum {j in DEST} demand[j,p];
param price_convex {PROD}; # dual price on convexity constr
param price {ORIG,DEST} <= 0.000001; # dual price on shipment limit
param cost {ORIG,DEST,PROD} >= 0; # shipment costs per unit
var Trans {ORIG,DEST,PROD} >= 0; # units to be shipped
minimize Artif_Reduced_Cost {p in PROD}:
sum {i in ORIG, j in DEST}
(- price[i,j]) * Trans[i,j,p] - price_convex[p];
minimize Reduced_Cost {p in PROD}:
sum {i in ORIG, j in DEST}
(cost[i,j,p] - price[i,j]) * Trans[i,j,p] - price_convex[p];
subject to Supply {i in ORIG, p in PROD}:
sum {j in DEST} Trans[i,j,p] = supply[i,p];
subject to Demand {j in DEST, p in PROD}:
sum {i in ORIG} Trans[i,j,p] = demand[j,p];
### MASTER PROBLEM ###
param limit {ORIG,DEST} >= 0; # max shipped on each link
param nPROP {PROD} integer >= 0;
param prop_ship {ORIG, DEST, p in PROD, 1..nPROP[p]} >= -0.000001;
param prop_cost {p in PROD, 1..nPROP[p]} >= 0;
# For each proposal from each subproblem:
# amount it ships over each link, and its cost
var Weight {p in PROD, 1..nPROP[p]} >= 0;
var Excess >= 0;
minimize Artificial: Excess;
minimize Total_Cost:
sum {p in PROD, k in 1..nPROP[p]} prop_cost[p,k] * Weight[p,k];
subject to Multi {i in ORIG, j in DEST}:
sum {p in PROD, k in 1..nPROP[p]}
prop_ship[i,j,p,k] * Weight[p,k] - Excess <= limit[i,j];
subject to Convex {p in PROD}: sum {k in 1..nPROP[p]} Weight[p,k] = 1;set ORIG := GARY CLEV PITT ;
set DEST := FRA DET LAN WIN STL FRE LAF ;
set PROD := bands coils plate ;
param supply (tr): GARY CLEV PITT :=
bands 400 700 800
coils 800 1600 1800
plate 200 300 300 ;
param demand (tr):
FRA DET LAN WIN STL FRE LAF :=
bands 300 300 100 75 650 225 250
coils 500 750 400 250 950 850 500
plate 100 100 0 50 200 100 250 ;
param limit default 625 ;
param cost :=
[*,*,bands]: FRA DET LAN WIN STL FRE LAF :=
GARY 30 10 8 10 11 71 6
CLEV 22 7 10 7 21 82 13
PITT 19 11 12 10 25 83 15
[*,*,coils]: FRA DET LAN WIN STL FRE LAF :=
GARY 39 14 11 14 16 82 8
CLEV 27 9 12 9 26 95 17
PITT 24 14 17 13 28 99 20
[*,*,plate]: FRA DET LAN WIN STL FRE LAF :=
GARY 41 15 12 16 17 86 8
CLEV 29 9 13 9 28 99 18
PITT 26 14 17 13 31 104 20 ;Same as multi2.run, except that the separate subproblems are realized by changing the data to a single AMPL named problem
# ----------------------------------------
# DANTZIG-WOLFE DECOMPOSITION FOR
# MULTI-COMMODITY TRANSPORTATION
# ----------------------------------------
model multi3.mod;
data multi3.dat;
param nIter default 0;
param nNegRC;
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 omit_zero_rows 1;
option display_1col 0;
option display_eps .000001;
# ----------------------------------------------------------
problem MasterI: Artificial, Weight, Excess, Multi, Convex;
problem SubI: Artif_Reduced_Cost, Trans, Supply, Demand;
repeat {
let nIter := nIter + 1;
let nNegRC := 0;
printf "\nPHASE I -- ITERATION %d\n", nIter;
problem SubI;
for {p in PROD} { printf "\nPRODUCT %s\n\n", p;
let P := p;
fix Trans; unfix {i in ORIG, j in DEST} Trans[i,j,p];
solve SubI;
printf "\n";
display {i in ORIG, j in DEST} Trans[i,j,p];
if Artif_Reduced_Cost < - 0.00001 then {
let nNegRC := nNegRC + 1;
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 nNegRC = 0 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: Reduced_Cost, Trans, Supply, Demand;
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;
let nNegRC := 0;
printf "\nPHASE II -- ITERATION %d\n\n", nIter;
problem SubII;
for {p in PROD} { printf "\nPRODUCT %s\n\n", p;
let P := p;
fix Trans; unfix {i in ORIG, j in DEST} Trans[i,j,p];
solve SubII;
printf "\n";
display {i in ORIG, j in DEST} Trans[i,j,p];
if Reduced_Cost < - 0.00001 then {
let nNegRC := nNegRC + 1;
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 nNegRC = 0 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;# ----------------------------------------
# MULTI-COMMODITY FLOW USING
# DANTZIG-WOLFE DECOMPOSITION
# (one single-product subproblem model)
# ----------------------------------------
### SUBPROBLEM ###
set ORIG; # origins
set DEST; # destinations
set PROD; # products
param supply {ORIG,PROD} >= 0; # amounts available at origins
param demand {DEST,PROD} >= 0; # amounts required at destinations
check {p in PROD}:
sum {i in ORIG} supply[i,p] = sum {j in DEST} demand[j,p];
param price_convex {PROD}; # dual price on convexity constr
param price {ORIG,DEST} <= 0.000001; # dual price on shipment limit
param cost {ORIG,DEST,PROD} >= 0; # shipment costs per unit
var Trans {ORIG,DEST,PROD} >= 0; # units to be shipped
# ----------------------------------------
param P symbolic within PROD;
minimize Artif_Reduced_Cost:
sum {i in ORIG, j in DEST}
(- price[i,j]) * Trans[i,j,P] - price_convex[P];
minimize Reduced_Cost:
sum {i in ORIG, j in DEST}
(cost[i,j,P] - price[i,j]) * Trans[i,j,P] - price_convex[P];
subject to Supply {i in ORIG}:
sum {j in DEST} Trans[i,j,P] = supply[i,P];
subject to Demand {j in DEST}:
sum {i in ORIG} Trans[i,j,P] = demand[j,P];
### MASTER PROBLEM ###
param limit {ORIG,DEST} >= 0; # max shipped on each link
param nPROP {PROD} integer >= 0;
param prop_ship {ORIG, DEST, p in PROD, 1..nPROP[p]} >= -0.000001;
param prop_cost {p in PROD, 1..nPROP[p]} >= 0;
# For each proposal from each subproblem:
# amount it ships over each link, and its cost
var Weight {p in PROD, 1..nPROP[p]} >= 0;
var Excess >= 0;
minimize Artificial: Excess;
minimize Total_Cost:
sum {p in PROD, k in 1..nPROP[p]} prop_cost[p,k] * Weight[p,k];
subject to Multi {i in ORIG, j in DEST}:
sum {p in PROD, k in 1..nPROP[p]}
prop_ship[i,j,p,k] * Weight[p,k] - Excess <= limit[i,j];
subject to Convex {p in PROD}: sum {k in 1..nPROP[p]} Weight[p,k] = 1;set ORIG := GARY CLEV PITT ;
set DEST := FRA DET LAN WIN STL FRE LAF ;
set PROD := bands coils plate ;
param supply (tr): GARY CLEV PITT :=
bands 400 700 800
coils 800 1600 1800
plate 200 300 300 ;
param demand (tr):
FRA DET LAN WIN STL FRE LAF :=
bands 300 300 100 75 650 225 250
coils 500 750 400 250 950 850 500
plate 100 100 0 50 200 100 250 ;
param limit default 625 ;
param cost :=
[*,*,bands]: FRA DET LAN WIN STL FRE LAF :=
GARY 30 10 8 10 11 71 6
CLEV 22 7 10 7 21 82 13
PITT 19 11 12 10 25 83 15
[*,*,coils]: FRA DET LAN WIN STL FRE LAF :=
GARY 39 14 11 14 16 82 8
CLEV 27 9 12 9 26 95 17
PITT 24 14 17 13 28 99 20
[*,*,plate]: FRA DET LAN WIN STL FRE LAF :=
GARY 41 15 12 16 17 86 8
CLEV 29 9 13 9 28 99 18
PITT 26 14 17 13 31 104 20 ;Benders decomposition for a stochastic programming variant of a multi-period production problem (see Exercise 4-5)
# ----------------------------------------
# STOCHASTIC PROGRAMMING PROBLEM
# USING BENDERS DECOMPOSITION
# ----------------------------------------
model stoch1.mod;
data stoch.dat;
option solver cplex;
option omit_zero_rows 1;
option display_eps .000001;
problem Master: Make1, Inv1, Sell1, Min_Stage2_Profit,
Expected_Profit, Cut_Defn, Time1, Balance1;
problem Sub: Make, Inv, Sell, Stage2_Profit, Time, Balance2, Balance;
let nCUT := 0;
let Min_Stage2_Profit := Infinity;
let {p in PROD} inv1[p] := 0;
param GAP default Infinity;
for {1..50} { printf "\nITERATION %d\n\n", nCUT+1;
solve Sub;
printf "\n";
if Stage2_Profit < Min_Stage2_Profit - 0.00001 then {
let GAP := min (GAP, Min_Stage2_Profit - Stage2_Profit);
option display_1col 0;
display GAP, Make, Sell, Inv;
let nCUT := nCUT + 1;
let {t in 2..T, s in SCEN} time_price[t,s,nCUT] := Time[t,s].dual;
let {p in PROD, s in SCEN} bal2_price[p,s,nCUT] := Balance2[p,s].dual;
let {p in PROD, t in 2..T, s in SCEN}
sell_lim_price[p,t,s,nCUT] := Sell[p,t,s].urc;
}
else break;
printf "\nRE-SOLVING MASTER PROBLEM\n\n";
solve Master;
printf "\n";
option display_1col 20;
display Make1, Inv1, Sell1;
let {p in PROD} inv1[p] := Inv1[p];
};
printf "\nOPTIMAL SOLUTION FOUND\nExpected Profit = %f\n\n", Expected_Profit;
option display_1col 0;
param MAKE {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Make1[p].val else Make[p,t,s].val;
param SELL {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Sell1[p].val else Sell[p,t,s].val;
param INV {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Inv1[p].val else Inv[p,t,s].val;
for {s in SCEN} {
printf "SCENARIO %s\n", s;
display {p in PROD, t in 1..T}
(MAKE[p,t,s], SELL[p,t,s], INV[p,t,s]);
}# ----------------------------------------
# STOCHASTIC PROGRAMMING PROBLEM
# USING BENDERS DECOMPOSITION
# ----------------------------------------
### SUBPROBLEM (always feasible) ###
set PROD; # products
param T > 0; # number of weeks
set SCEN; # number of scenarios
param rate {PROD} > 0; # tons per hour produced
param avail {1..T} >= 0; # hours available in week
param market {PROD,1..T} >= 0; # limit on tons sold in week
param prodcost {PROD} >= 0; # cost per ton produced
param invcost {PROD} >= 0; # carrying cost/ton of inventory
param revenue {PROD,1..T,SCEN} >= 0; # projected revenue/ton
param prob {SCEN} >= 0, <= 1;
check: 0.99999 < sum {s in SCEN} prob[s] < 1.00001;
param inv1 {PROD} >= 0; # inventory at end of first period
var Make {PROD,2..T,SCEN} >= 0; # tons produced
var Inv {PROD,2..T,SCEN} >= 0; # tons inventoried
var Sell {p in PROD, t in 2..T, SCEN} # tons sold
>= 0, <= market[p,t];
maximize Stage2_Profit:
sum {s in SCEN} prob[s] *
sum {p in PROD, t in 2..T} (revenue[p,t,s]*Sell[p,t,s] -
prodcost[p]*Make[p,t,s] - invcost[p]*Inv[p,t,s]);
subject to Time {t in 2..T, s in SCEN}:
sum {p in PROD} (1/rate[p]) * Make[p,t,s] <= avail[t];
subject to Balance2 {p in PROD, s in SCEN}:
Make[p,2,s] + inv1[p] = Sell[p,2,s] + Inv[p,2,s];
subject to Balance {p in PROD, t in 3..T, s in SCEN}:
Make[p,t,s] + Inv[p,t-1,s] = Sell[p,t,s] + Inv[p,t,s];
### MASTER PROBLEM ###
param inv0 {PROD} >= 0; # initial inventory
param nCUT >= 0 integer;
param time_price {2..T,SCEN,1..nCUT} >= -0.000001;
param bal2_price {PROD,SCEN,1..nCUT};
param sell_lim_price {PROD,2..T,SCEN,1..nCUT} >= -0.000001;
var Make1 {PROD} >= 0;
var Inv1 {PROD} >= 0;
var Sell1 {p in PROD} >= 0, <= market[p,1];
var Min_Stage2_Profit >= 0;
maximize Expected_Profit:
sum {s in SCEN} prob[s] *
sum {p in PROD} (revenue[p,1,s]*Sell1[p] -
prodcost[p]*Make1[p] - invcost[p]*Inv1[p]) +
Min_Stage2_Profit;
subj to Cut_Defn {k in 1..nCUT}:
Min_Stage2_Profit <=
sum {t in 2..T, s in SCEN} time_price[t,s,k] * avail[t] +
sum {p in PROD, s in SCEN} bal2_price[p,s,k] * (-Inv1[p]) +
sum {p in PROD, t in 2..T, s in SCEN}
sell_lim_price[p,t,s,k] * market[p,t];
subject to Time1:
sum {p in PROD} (1/rate[p]) * Make1[p] <= avail[1];
subject to Balance1 {p in PROD}:
Make1[p] + inv0[p] = Sell1[p] + Inv1[p];
# ----------------------------------------
# STOCHASTIC PROGRAMMING PROBLEM
# USING BENDERS DECOMPOSITION
# ----------------------------------------
model stoch2.mod;
data stoch.dat;
option solver cplex;
option omit_zero_rows 1;
option display_eps .000001;
problem Master: Make1, Inv1, Sell1, Min_Stage2_Profit,
Expected_Profit, Cut_Defn, Time1, Balance1;
problem Sub {s in SCEN}:
{p in PROD, t in 2..T} Make[p,t,s],
{p in PROD, t in 2..T} Inv[p,t,s],
{p in PROD, t in 2..T} Sell[p,t,s],
Exp_Stage2_Profit[s],
{t in 2..T} Time[t,s],
{p in PROD} Balance2[p,s],
{p in PROD, t in 3..T} Balance[p,t,s];
let nCUT := 0;
let Min_Stage2_Profit := Infinity;
let {p in PROD} inv1[p] := 0;
param GAP default Infinity;
param newGAP;
for {1..50} { printf "\nITERATION %d\n\n", nCUT+1;
let nCUT := nCUT + 1;
for {s in SCEN} {
solve Sub[s];
let {t in 2..T} time_price[t,s,nCUT] := Time[t,s].dual;
let {p in PROD} bal2_price[p,s,nCUT] := Balance2[p,s].dual;
let {p in PROD, t in 2..T}
sell_lim_price[p,t,s,nCUT] := Sell[p,t,s].urc;
}
let newGAP := Min_Stage2_Profit - sum {s in SCEN} Exp_Stage2_Profit[s];
printf "\n";
if newGAP > 0.00001 then {
let GAP := min (GAP, newGAP);
option display_1col 0;
display GAP, Make, Sell, Inv;
}
else break;
printf "\nRE-SOLVING MASTER PROBLEM\n\n";
solve Master;
printf "\n";
option display_1col 20;
display Make1, Inv1, Sell1;
let {p in PROD} inv1[p] := Inv1[p];
};
printf "\nOPTIMAL SOLUTION FOUND\nExpected Profit = %f\n\n", Expected_Profit;
option display_1col 0;
param MAKE {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Make1[p].val else Make[p,t,s].val;
param SELL {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Sell1[p].val else Sell[p,t,s].val;
param INV {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Inv1[p].val else Inv[p,t,s].val;
for {s in SCEN} {
printf "SCENARIO %s\n", s;
display {p in PROD, t in 1..T}
(MAKE[p,t,s], SELL[p,t,s], INV[p,t,s]);
}Same as stoch1.run, but using a separate subproblem for each scenario; subproblems are represented in AMPL by an indexed collection of named problems
# ----------------------------------------
# STOCHASTIC PROGRAMMING PROBLEM
# USING BENDERS DECOMPOSITION
# ----------------------------------------
model stoch2.mod;
data stoch.dat;
option solver cplex;
option omit_zero_rows 1;
option display_eps .000001;
problem Master: Make1, Inv1, Sell1, Min_Stage2_Profit,
Expected_Profit, Cut_Defn, Time1, Balance1;
problem Sub {s in SCEN}:
{p in PROD, t in 2..T} Make[p,t,s],
{p in PROD, t in 2..T} Inv[p,t,s],
{p in PROD, t in 2..T} Sell[p,t,s],
Exp_Stage2_Profit[s],
{t in 2..T} Time[t,s],
{p in PROD} Balance2[p,s],
{p in PROD, t in 3..T} Balance[p,t,s];
let nCUT := 0;
let Min_Stage2_Profit := Infinity;
let {p in PROD} inv1[p] := 0;
param GAP default Infinity;
param newGAP;
for {1..50} { printf "\nITERATION %d\n\n", nCUT+1;
let nCUT := nCUT + 1;
for {s in SCEN} {
solve Sub[s];
let {t in 2..T} time_price[t,s,nCUT] := Time[t,s].dual;
let {p in PROD} bal2_price[p,s,nCUT] := Balance2[p,s].dual;
let {p in PROD, t in 2..T}
sell_lim_price[p,t,s,nCUT] := Sell[p,t,s].urc;
}
let newGAP := Min_Stage2_Profit - sum {s in SCEN} Exp_Stage2_Profit[s];
printf "\n";
if newGAP > 0.00001 then {
let GAP := min (GAP, newGAP);
option display_1col 0;
display GAP, Make, Sell, Inv;
}
else break;
printf "\nRE-SOLVING MASTER PROBLEM\n\n";
solve Master;
printf "\n";
option display_1col 20;
display Make1, Inv1, Sell1;
let {p in PROD} inv1[p] := Inv1[p];
};
printf "\nOPTIMAL SOLUTION FOUND\nExpected Profit = %f\n\n", Expected_Profit;
option display_1col 0;
param MAKE {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Make1[p].val else Make[p,t,s].val;
param SELL {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Sell1[p].val else Sell[p,t,s].val;
param INV {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Inv1[p].val else Inv[p,t,s].val;
for {s in SCEN} {
printf "SCENARIO %s\n", s;
display {p in PROD, t in 1..T}
(MAKE[p,t,s], SELL[p,t,s], INV[p,t,s]);
}# ----------------------------------------
# STOCHASTIC PROGRAMMING PROBLEM
# USING BENDERS DECOMPOSITION
# ----------------------------------------
### SUBPROBLEM ###
set PROD; # products
param T > 0; # number of weeks
set SCEN; # number of scenarios
param rate {PROD} > 0; # tons per hour produced
param avail {1..T} >= 0; # hours available in week
param market {PROD,1..T} >= 0; # limit on tons sold in week
param prodcost {PROD} >= 0; # cost per ton produced
param invcost {PROD} >= 0; # carrying cost/ton of inventory
param revenue {PROD,1..T,SCEN} >= 0; # projected revenue/ton
param prob {SCEN} >= 0, <= 1;
check: 0.99999 < sum {s in SCEN} prob[s] < 1.00001;
param inv1 {PROD} >= 0; # inventory at end of first period
var Make {PROD,2..T,SCEN} >= 0; # tons produced
var Inv {PROD,2..T,SCEN} >= 0; # tons inventoried
var Sell {p in PROD, t in 2..T, SCEN} # tons sold
>= 0, <= market[p,t];
maximize Exp_Stage2_Profit {s in SCEN}: prob[s] *
sum {p in PROD, t in 2..T} (revenue[p,t,s]*Sell[p,t,s] -
prodcost[p]*Make[p,t,s] - invcost[p]*Inv[p,t,s]);
subject to Time {t in 2..T, s in SCEN}:
sum {p in PROD} (1/rate[p]) * Make[p,t,s] <= avail[t];
subject to Balance2 {p in PROD, s in SCEN}:
Make[p,2,s] + inv1[p] = Sell[p,2,s] + Inv[p,2,s];
subject to Balance {p in PROD, t in 3..T, s in SCEN}:
Make[p,t,s] + Inv[p,t-1,s] = Sell[p,t,s] + Inv[p,t,s];
### MASTER PROBLEM ###
param inv0 {PROD} >= 0; # initial inventory
param nCUT >= 0 integer;
param time_price {2..T,SCEN,1..nCUT} >= -0.000001 default 0;
param bal2_price {PROD,SCEN,1..nCUT} default 0;
param sell_lim_price {PROD,2..T,SCEN,1..nCUT} >= -0.000001 default 0;
var Make1 {PROD} >= 0;
var Inv1 {PROD} >= 0;
var Sell1 {p in PROD} >= 0, <= market[p,1];
var Min_Stage2_Profit >= 0;
maximize Expected_Profit:
sum {s in SCEN} prob[s] *
sum {p in PROD} (revenue[p,1,s]*Sell1[p] -
prodcost[p]*Make1[p] - invcost[p]*Inv1[p]) +
Min_Stage2_Profit;
subj to Cut_Defn {k in 1..nCUT}:
Min_Stage2_Profit <=
sum {t in 2..T, s in SCEN} time_price[t,s,k] * avail[t] +
sum {p in PROD, s in SCEN} bal2_price[p,s,k] * (-Inv1[p]) +
sum {p in PROD, t in 2..T, s in SCEN}
sell_lim_price[p,t,s,k] * market[p,t];
subject to Time1:
sum {p in PROD} (1/rate[p]) * Make1[p] <= avail[1];
subject to Balance1 {p in PROD}:
Make1[p] + inv0[p] = Sell1[p] + Inv1[p];
# ----------------------------------------
# STOCHASTIC PROGRAMMING PROBLEM
# USING BENDERS DECOMPOSITION
# ----------------------------------------
model stoch2.mod;
data stoch.dat;
option solver cplex;
option omit_zero_rows 1;
option display_eps .000001;
problem Master: Make1, Inv1, Sell1, Min_Stage2_Profit,
Expected_Profit, Cut_Defn, Time1, Balance1;
problem Sub {s in SCEN}:
{p in PROD, t in 2..T} Make[p,t,s],
{p in PROD, t in 2..T} Inv[p,t,s],
{p in PROD, t in 2..T} Sell[p,t,s],
Exp_Stage2_Profit[s],
{t in 2..T} Time[t,s],
{p in PROD} Balance2[p,s],
{p in PROD, t in 3..T} Balance[p,t,s];
let nCUT := 0;
let Min_Stage2_Profit := Infinity;
let {p in PROD} inv1[p] := 0;
param GAP default Infinity;
param newGAP;
for {1..50} { printf "\nITERATION %d\n\n", nCUT+1;
let nCUT := nCUT + 1;
for {s in SCEN} {
solve Sub[s];
let {t in 2..T} time_price[t,s,nCUT] := Time[t,s].dual;
let {p in PROD} bal2_price[p,s,nCUT] := Balance2[p,s].dual;
let {p in PROD, t in 2..T}
sell_lim_price[p,t,s,nCUT] := Sell[p,t,s].urc;
}
let newGAP := Min_Stage2_Profit - sum {s in SCEN} Exp_Stage2_Profit[s];
printf "\n";
if newGAP > 0.00001 then {
let GAP := min (GAP, newGAP);
option display_1col 0;
display GAP, Make, Sell, Inv;
}
else break;
printf "\nRE-SOLVING MASTER PROBLEM\n\n";
solve Master;
printf "\n";
option display_1col 20;
display Make1, Inv1, Sell1;
let {p in PROD} inv1[p] := Inv1[p];
};
printf "\nOPTIMAL SOLUTION FOUND\nExpected Profit = %f\n\n", Expected_Profit;
option display_1col 0;
param MAKE {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Make1[p].val else Make[p,t,s].val;
param SELL {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Sell1[p].val else Sell[p,t,s].val;
param INV {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Inv1[p].val else Inv[p,t,s].val;
for {s in SCEN} {
printf "SCENARIO %s\n", s;
display {p in PROD, t in 1..T}
(MAKE[p,t,s], SELL[p,t,s], INV[p,t,s]);
}Same as stoch2.run, except that the separate subproblems are realized by changing the data to a single AMPL named problem
# ----------------------------------------
# STOCHASTIC PROGRAMMING PROBLEM
# USING BENDERS DECOMPOSITION
# ----------------------------------------
model stoch3.mod;
data stoch.dat;
option solver cplex;
option omit_zero_rows 1;
option display_eps .000001;
problem Master: Make1, Inv1, Sell1, Min_Stage2_Profit,
Expected_Profit, Cut_Defn, Time1, Balance1;
problem Sub: Make, Inv, Sell, Exp_Stage2_Profit, Time, Balance2, Balance;
let nCUT := 0;
let Min_Stage2_Profit := Infinity;
let {p in PROD} inv1[p] := 0;
param GAP default Infinity;
param newGAP;
for {1..50} { printf "\nITERATION %d\n\n", nCUT+1;
let nCUT := nCUT + 1;
let newGAP := Min_Stage2_Profit;
for {s in SCEN} {
let S := s;
solve Sub;
let newGAP := newGAP - Exp_Stage2_Profit;
let {t in 2..T} time_price[t,S,nCUT] := Time[t].dual;
let {p in PROD} bal2_price[p,S,nCUT] := Balance2[p].dual;
let {p in PROD, t in 2..T}
sell_lim_price[p,t,S,nCUT] := Sell[p,t,S].urc;
}
printf "\n";
if newGAP > 0.00001 then {
let GAP := min (GAP, newGAP);
option display_1col 0;
display GAP, Make, Sell, Inv;
}
else break;
printf "\nRE-SOLVING MASTER PROBLEM\n\n";
solve Master;
printf "\n";
option display_1col 20;
display Make1, Inv1, Sell1;
let {p in PROD} inv1[p] := Inv1[p];
};
printf "\nOPTIMAL SOLUTION FOUND\nExpected Profit = %f\n\n", Expected_Profit;
option display_1col 0;
param MAKE {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Make1[p].val else Make[p,t,s].val;
param SELL {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Sell1[p].val else Sell[p,t,s].val;
param INV {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Inv1[p].val else Inv[p,t,s].val;
for {s in SCEN} {
printf "SCENARIO %s\n", s;
display {p in PROD, t in 1..T}
(MAKE[p,t,s], SELL[p,t,s], INV[p,t,s]);
}# ----------------------------------------
# STOCHASTIC PROGRAMMING PROBLEM
# USING BENDERS DECOMPOSITION
# ----------------------------------------
### SUBPROBLEMS ###
set PROD; # products
param T > 0; # number of weeks
set SCEN; # number of scenarios
param rate {PROD} > 0; # tons per hour produced
param avail {1..T} >= 0; # hours available in week
param market {PROD,1..T} >= 0; # limit on tons sold in week
param prodcost {PROD} >= 0; # cost per ton produced
param invcost {PROD} >= 0; # carrying cost/ton of inventory
param revenue {PROD,1..T,SCEN} >= 0; # projected revenue/ton
param prob {SCEN} >= 0, <= 1;
check: 0.99999 < sum {s in SCEN} prob[s] < 1.00001;
param inv1 {PROD} >= 0; # inventory at end of first period
var Make {PROD,2..T,SCEN} >= 0; # tons produced
var Inv {PROD,2..T,SCEN} >= 0; # tons inventoried
var Sell {p in PROD, t in 2..T,SCEN} # tons sold
>= 0, <= market[p,t];
# ----------------------------------------
param S symbolic in SCEN; # scenario of subproblem
maximize Exp_Stage2_Profit: prob[S] *
sum {p in PROD, t in 2..T} (revenue[p,t,S]*Sell[p,t,S] -
prodcost[p]*Make[p,t,S] - invcost[p]*Inv[p,t,S]);
subject to Time {t in 2..T}:
sum {p in PROD} (1/rate[p]) * Make[p,t,S] <= avail[t];
subject to Balance2 {p in PROD}:
Make[p,2,S] + inv1[p] = Sell[p,2,S] + Inv[p,2,S];
subject to Balance {p in PROD, t in 3..T}:
Make[p,t,S] + Inv[p,t-1,S] = Sell[p,t,S] + Inv[p,t,S];
### MASTER PROBLEM ###
param inv0 {PROD} >= 0; # initial inventory
param nCUT >= 0 integer;
param time_price {2..T,SCEN,1..nCUT} >= -0.000001 default 0;
param bal2_price {PROD,SCEN,1..nCUT} default 0;
param sell_lim_price {PROD,2..T,SCEN,1..nCUT} >= -0.000001 default 0;
var Make1 {PROD} >= 0;
var Inv1 {PROD} >= 0;
var Sell1 {p in PROD} >= 0, <= market[p,1];
var Min_Stage2_Profit >= 0;
maximize Expected_Profit:
sum {s in SCEN} prob[s] *
sum {p in PROD} (revenue[p,1,s]*Sell1[p] -
prodcost[p]*Make1[p] - invcost[p]*Inv1[p]) +
Min_Stage2_Profit;
subj to Cut_Defn {k in 1..nCUT}:
Min_Stage2_Profit <=
sum {t in 2..T, s in SCEN} time_price[t,s,k] * avail[t] +
sum {p in PROD, s in SCEN} bal2_price[p,s,k] * (-Inv1[p]) +
sum {p in PROD, t in 2..T, s in SCEN}
sell_lim_price[p,t,s,k] * market[p,t];
subject to Time1:
sum {p in PROD} (1/rate[p]) * Make1[p] <= avail[1];
subject to Balance1 {p in PROD}:
Make1[p] + inv0[p] = Sell1[p] + Inv1[p];
# ----------------------------------------
# STOCHASTIC PROGRAMMING PROBLEM
# USING BENDERS DECOMPOSITION
# ----------------------------------------
model stoch2.mod;
data stoch.dat;
option solver cplex;
option omit_zero_rows 1;
option display_eps .000001;
problem Master: Make1, Inv1, Sell1, Min_Stage2_Profit,
Expected_Profit, Cut_Defn, Time1, Balance1;
problem Sub {s in SCEN}:
{p in PROD, t in 2..T} Make[p,t,s],
{p in PROD, t in 2..T} Inv[p,t,s],
{p in PROD, t in 2..T} Sell[p,t,s],
Exp_Stage2_Profit[s],
{t in 2..T} Time[t,s],
{p in PROD} Balance2[p,s],
{p in PROD, t in 3..T} Balance[p,t,s];
let nCUT := 0;
let Min_Stage2_Profit := Infinity;
let {p in PROD} inv1[p] := 0;
param GAP default Infinity;
param newGAP;
for {1..50} { printf "\nITERATION %d\n\n", nCUT+1;
let nCUT := nCUT + 1;
for {s in SCEN} {
solve Sub[s];
let {t in 2..T} time_price[t,s,nCUT] := Time[t,s].dual;
let {p in PROD} bal2_price[p,s,nCUT] := Balance2[p,s].dual;
let {p in PROD, t in 2..T}
sell_lim_price[p,t,s,nCUT] := Sell[p,t,s].urc;
}
let newGAP := Min_Stage2_Profit - sum {s in SCEN} Exp_Stage2_Profit[s];
printf "\n";
if newGAP > 0.00001 then {
let GAP := min (GAP, newGAP);
option display_1col 0;
display GAP, Make, Sell, Inv;
}
else break;
printf "\nRE-SOLVING MASTER PROBLEM\n\n";
solve Master;
printf "\n";
option display_1col 20;
display Make1, Inv1, Sell1;
let {p in PROD} inv1[p] := Inv1[p];
};
printf "\nOPTIMAL SOLUTION FOUND\nExpected Profit = %f\n\n", Expected_Profit;
option display_1col 0;
param MAKE {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Make1[p].val else Make[p,t,s].val;
param SELL {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Sell1[p].val else Sell[p,t,s].val;
param INV {p in PROD, t in 1..T, s in SCEN}
:= if t = 1 then Inv1[p].val else Inv[p,t,s].val;
for {s in SCEN} {
printf "SCENARIO %s\n", s;
display {p in PROD, t in 1..T}
(MAKE[p,t,s], SELL[p,t,s], INV[p,t,s]);
}Benders decomposition for a location-transportation problem (original model in trnloc.mod)
# ----------------------------------------
# BENDERS DECOMPOSITION FOR
# THE LOCATION-TRANSPORTATION PROBLEM
# ----------------------------------------
model trnloc1.mod;
data trnloc.dat;
option omit_zero_rows 1;
option display_eps .000001;
problem Master: Build, Max_Ship_Cost, Total_Cost, Cut_Defn;
option display_1col 20;
option solver cplex, cplex_options '';
problem Sub: Ship, Ship_Cost, Supply, Demand;
option presolve 0;
option display_1col 0;
option solver cplex, cplex_options 'netopt 2 netfind 2';
let nCUT := 0;
let Max_Ship_Cost := 0;
let {i in ORIG} build[i] := 1;
param GAP default Infinity;
suffix dunbdd;
repeat { printf "\nITERATION %d\n\n", nCUT+1;
solve Sub;
printf "\n";
if Sub.result = "infeasible" then {
let nCUT := nCUT + 1;
let cut_type[nCUT] := "ray";
let {i in ORIG} supply_price[i,nCUT] := Supply[i].dunbdd;
let {j in DEST} demand_price[j,nCUT] := Demand[j].dunbdd;
expand Cut_Defn[nCUT];
}
else if Ship_Cost > Max_Ship_Cost + 0.00001 then {
let GAP := min (GAP, Ship_Cost - Max_Ship_Cost);
display GAP, Ship;
let nCUT := nCUT + 1;
let cut_type[nCUT] := "point";
let {i in ORIG} supply_price[i,nCUT] := Supply[i].dual;
let {j in DEST} demand_price[j,nCUT] := Demand[j].dual;
}
else break;
printf "\nRE-SOLVING MASTER PROBLEM\n\n";
solve Master;
printf "\n";
display Build;
let {i in ORIG} build[i] := Build[i];
};
option display_1col 0;
display Ship;# ----------------------------------------
# LOCATION-TRANSPORTATION PROBLEM
# USING BENDERS DECOMPOSITION
# (using primal formulation of subproblem)
# ----------------------------------------
### SUBPROBLEM ###
set ORIG; # shipment origins (warehouses)
set DEST; # shipment destinations (stores)
param supply {ORIG} > 0;
param demand {DEST} > 0;
param fix_cost {ORIG} > 0;
param var_cost {ORIG,DEST} > 0;
param build {ORIG} binary; # = 1 iff warehouse built at i
var Ship {ORIG,DEST} >= 0; # amounts shipped
minimize Ship_Cost:
sum {i in ORIG, j in DEST} var_cost[i,j] * Ship[i,j];
subj to Supply {i in ORIG}:
sum {j in DEST} Ship[i,j] <= supply[i] * build[i];
subj to Demand {j in DEST}:
sum {i in ORIG} Ship[i,j] = demand[j];
### MASTER PROBLEM ###
param nCUT >= 0 integer;
param cut_type {1..nCUT} symbolic within {"point","ray"};
param supply_price {ORIG,1..nCUT} <= 0.000001;
param demand_price {DEST,1..nCUT};
var Build {ORIG} binary; # = 1 iff warehouse built at i
var Max_Ship_Cost >= 0;
minimize Total_Cost:
sum {i in ORIG} fix_cost[i] * Build[i] + Max_Ship_Cost;
subj to Cut_Defn {k in 1..nCUT}:
if cut_type[k] = "point" then Max_Ship_Cost >=
sum {i in ORIG} supply_price[i,k] * supply[i] * Build[i] +
sum {j in DEST} demand_price[j,k] * demand[j];param: ORIG: supply :=
1 23070 6 21090 11 22690 16 17110 21 16720
2 18290 7 16650 12 19360 17 15380 22 21540
3 20010 8 18420 13 22330 18 18690 23 16500
4 15080 9 19160 14 15440 19 20720 24 15310
5 17540 10 18860 15 19330 20 21220 25 18970 ;
param: DEST: demand :=
A3 12000
A6 12000
A8 14000
A9 13500
B2 25000
B4 29000 ;
param fix_cost default 500000 ;
param var_cost:
A3 A6 A8 A9 B2 B4 :=
1 73.78 14.76 86.82 91.19 51.03 76.49
2 60.28 20.92 76.43 83.99 58.84 68.86
3 58.18 21.64 69.84 72.39 61.64 58.39
4 50.37 21.74 61.49 65.72 60.48 56.68
5 42.73 35.19 44.11 58.08 65.76 55.51
6 44.62 39.21 44.44 48.32 76.12 51.17
7 49.31 51.72 36.27 42.96 84.52 49.61
8 50.79 59.25 22.53 33.22 94.30 49.66
9 51.93 72.13 21.66 29.39 93.52 49.63
10 65.90 13.07 79.59 86.07 46.83 69.55
11 50.79 9.99 67.83 78.81 49.34 60.79
12 47.51 12.95 59.57 67.71 51.13 54.65
13 39.36 19.01 56.39 62.37 57.25 47.91
14 33.55 30.16 40.66 48.50 60.83 42.51
15 34.17 40.46 40.23 47.10 66.22 38.94
16 41.68 53.03 22.56 30.89 77.22 35.88
17 42.75 62.94 18.58 27.02 80.36 40.11
18 46.46 71.17 17.17 21.16 91.65 41.56
19 56.83 8.84 83.99 91.88 41.38 67.79
20 46.21 2.92 68.94 76.86 38.89 60.38
21 41.67 11.69 61.05 70.06 43.24 48.48
22 25.57 17.59 54.93 57.07 44.93 43.97
23 28.16 29.39 38.64 46.48 50.16 34.20
24 26.97 41.62 29.72 40.61 59.56 31.21
25 34.24 54.09 22.13 28.43 69.68 24.09 ;Same as trnloc2a.run, but model has upper limits on the Ship variables: LP relaxation bound is still poor, but subproblem does not have the integrality property and considerable improvement is made
# ----------------------------------------
# LAGRANGIAN RELAXATION FOR
# THE LOCATION-TRANSPORTATION PROBLEM
# ----------------------------------------
printf "\nLP RELAXATION\n\n";
model trnloc2b.mod;
data trnloc2.dat;
option omit_zero_rows 1;
option display_eps .000001;
option solution_round 8;
option solver cplex;
option solver_msg 0;
option relax_integrality 1;
objective Shipping_Cost;
solve;
param LB; param UB;
let LB := Shipping_Cost.val;
let UB := sum {j in CITY} max {i in CITY} ship_cost[i,j];
option relax_integrality 0;
problem LowerBound: Build, Ship, Supply, Limit, Lagrangian;
problem UpperBound: Ship, Supply, Demand, Limit, Shipping_Cost;
let {j in CITY} mult[j] := 0;
param slack {CITY};
param scale default 1;
param norm;
param step;
param same default 0;
param same_limit := 3;
param iter_limit := 20;
param LBlog {0..iter_limit}; let LBlog[0] := LB;
param UBlog {0..iter_limit}; let UBlog[0] := UB;
param scalelog {1..iter_limit};
param steplog {1..iter_limit};
for {k in 1..iter_limit} { printf "\nITERATION %d\n\n", k;
solve LowerBound;
let {j in CITY} slack[j] := sum {i in CITY} Ship[i,j] - demand[j];
if Lagrangian > LB + 0.000001 then {
let LB := Lagrangian;
let same := 0; }
else let same := same + 1;
if same = same_limit then {
let scale := scale / 2;
let same := 0;
};
let norm := sum {j in CITY} slack[j]^2;
let step := scale * (UB - Lagrangian) / norm;
let {j in CITY} mult[j] := mult[j] less step * slack[j];
if sum {i in CITY} supply[i] * Build[i]
>= sum {j in CITY} demand[j] - 1e-8 then {
solve UpperBound;
let UB := min (UB, Shipping_Cost);
}
let LBlog[k] := LB;
let UBlog[k] := UB;
let scalelog[k] := scale;
let steplog[k] := step;
if UB - LB < 0.9999 then {break};
if step < 1.0e-8 then {break};
}
printf "\n\n";
display LBlog, UBlog, scalelog, steplog;# ----------------------------------------
# LOCATION-TRANSPORTATION PROBLEM
# USING LAGRANGIAN RELAXATION
# (with upper bounds on Ship variables)
# ----------------------------------------
set CITY;
param build_limit integer;
param demand {i in CITY} integer > 0;
param supply {CITY} integer > 0;
param ship_cost {i in CITY, j in CITY} >= 0;
param mult {CITY} >= 0; # Lagrange multipliers for Demand constr
var Build {CITY} integer >= 0 <= 1; # = 1 iff warehouse built at i
var Ship {i in CITY, j in CITY} # amounts shipped
>= 0, <= demand[j];
minimize Lagrangian:
sum {i in CITY, j in CITY} ship_cost[i,j] * Ship[i,j] +
sum {j in CITY} mult[j] * (demand[j] - sum {i in CITY} Ship[i,j]);
minimize Shipping_Cost:
sum {i in CITY, j in CITY} ship_cost[i,j] * Ship[i,j];
subj to Supply {i in CITY}:
sum {j in CITY} Ship[i,j] <= supply[i] * Build[i];
subj to Demand {j in CITY}:
sum {i in CITY} Ship[i,j] >= demand[j];
subj to Limit: sum {i in CITY} Build[i] <= build_limit;param build_limit := 8;
param: CITY: supply demand :=
1 19 4
2 13 1
3 20 2
4 18 7
5 23 9
6 25 9
7 22 3
8 28 9
9 18 5
10 16 2
11 17 4
12 14 1
13 22 9
14 13 3
15 18 1
16 20 7
17 15 5
18 23 9
19 25 2
20 19 4
21 18 9
22 19 6
23 22 5
24 22 4
25 21 5
26 19 6
27 17 8
28 21 9
29 15 5
30 25 6 ;
param ship_cost
: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 :=
1 0 61 35 35 53 16 11 58 28 34 43 18 48 30 68
2 61 0 96 88 29 77 50 92 57 75 48 47 76 32 77
3 35 96 0 18 88 21 45 72 58 53 65 49 51 66 93
4 35 88 18 0 87 27 42 84 62 62 51 42 33 62 100
5 53 29 88 87 0 67 45 66 37 54 62 47 85 26 49
6 16 77 21 27 67 0 27 57 36 34 55 32 51 46 73
7 11 50 45 42 45 27 0 64 29 40 35 8 46 20 68
8 58 92 72 84 66 57 64 0 37 24 99 72 106 65 31
9 28 57 58 62 37 36 29 37 0 18 62 36 74 28 40
10 34 75 53 62 54 34 40 24 18 0 75 48 81 45 41
11 43 48 65 51 62 55 35 99 62 75 0 27 29 39 99
12 18 47 49 42 47 32 8 72 36 48 27 0 40 21 74
13 48 76 51 33 85 51 46 106 74 81 29 40 0 60 114
14 30 32 66 62 26 46 20 65 28 45 39 21 60 0 60
15 68 77 93 100 49 73 68 31 40 41 99 74 114 60 0
16 15 74 22 20 68 10 24 67 41 43 46 27 41 44 80
17 27 34 61 55 35 43 16 72 34 50 30 13 51 10 69
18 53 75 76 85 48 57 55 19 26 24 88 62 100 51 17
19 26 87 14 26 77 10 37 60 45 39 63 42 55 56 79
20 32 68 40 25 73 37 31 90 59 66 27 27 15 48 99
21 13 70 33 39 57 12 23 47 24 23 57 31 59 38 61
22 41 38 76 76 14 54 35 54 23 40 59 39 79 19 42
23 54 8 88 81 27 69 43 87 51 69 41 39 69 25 75
24 36 87 43 55 67 29 44 30 30 15 79 52 80 55 53
25 51 71 76 83 44 57 52 23 24 24 85 60 98 48 17
26 68 82 91 100 54 73 69 26 41 39 101 77 115 63 6
27 29 39 64 63 24 43 22 56 19 37 48 26 66 10 51
28 20 78 17 16 73 12 29 69 46 46 49 31 40 49 85
29 49 104 43 59 85 38 60 38 48 32 92 67 89 72 67
30 48 106 33 51 90 34 59 49 54 39 89 65 82 74 77
: 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 :=
1 15 27 53 26 32 13 41 54 36 51 68 29 20 49 48
2 74 34 75 87 68 70 38 8 87 71 82 39 78 104 106
3 22 61 76 14 40 33 76 88 43 76 91 64 17 43 33
4 20 55 85 26 25 39 76 81 55 83 100 63 16 59 51
5 68 35 48 77 73 57 14 27 67 44 54 24 73 85 90
6 10 43 57 10 37 12 54 69 29 57 73 43 12 38 34
7 24 16 55 37 31 23 35 43 44 52 69 22 29 60 59
8 67 72 19 60 90 47 54 87 30 23 26 56 69 38 49
9 41 34 26 45 59 24 23 51 30 24 41 19 46 48 54
10 43 50 24 39 66 23 40 69 15 24 39 37 46 32 39
11 46 30 88 63 27 57 59 41 79 85 101 48 49 92 89
12 27 13 62 42 27 31 39 39 52 60 77 26 31 67 65
13 41 51 100 55 15 59 79 69 80 98 115 66 40 89 82
14 44 10 51 56 48 38 19 25 55 48 63 10 49 72 74
15 80 69 17 79 99 61 42 75 53 17 6 51 85 67 77
16 0 39 65 17 26 20 56 66 39 63 80 44 5 48 43
17 39 0 59 53 39 38 29 27 58 56 72 18 44 74 74
18 65 59 0 63 85 45 37 70 37 4 15 42 69 51 60
19 17 53 63 0 42 21 64 80 30 63 78 53 15 34 27
20 26 39 85 42 0 44 65 61 65 83 100 52 27 75 69
21 20 38 45 21 44 0 44 63 22 44 61 34 24 36 36
22 56 29 37 64 65 44 0 34 53 33 46 13 61 71 76
23 66 27 70 80 61 63 34 0 80 67 79 32 71 97 99
24 39 58 37 30 65 22 53 80 0 38 51 48 41 18 25
25 63 56 4 63 83 44 33 67 38 0 17 38 68 53 62
26 80 72 15 78 100 61 46 79 51 17 0 54 84 63 73
27 44 18 42 53 52 34 13 32 48 38 54 0 49 66 69
28 5 44 69 15 27 24 61 71 41 68 84 49 0 48 42
29 48 74 51 34 75 36 71 97 18 53 63 66 48 0 12
30 43 74 60 27 69 36 76 99 25 62 73 69 42 12 0 ;Same as trnloc2b.run, but model has 0-1 constraints disaggregated: LP relaxation bound is good, but subproblem has the integrality property and no further improvement can be made
# ----------------------------------------
# LAGRANGIAN RELAXATION FOR
# THE LOCATION-TRANSPORTATION PROBLEM
# ----------------------------------------
printf "\nLP RELAXATION\n\n";
model trnloc2c.mod;
data trnloc2.dat;
option omit_zero_rows 1;
option display_eps .000001;
option solution_round 8;
option solver cplex;
option solver_msg 0;
option relax_integrality 1;
objective Shipping_Cost;
solve;
param LB; param UB;
let LB := Shipping_Cost.val;
let UB := sum {j in CITY} max {i in CITY} ship_cost[i,j];
option relax_integrality 0;
problem LowerBound: Build, Ship, Supply, Build_Def, Limit, Lagrangian;
problem UpperBound: Ship, Supply, Build_Def, Demand, Limit, Shipping_Cost;
let {j in CITY} mult[j] := 0;
param slack {CITY};
param scale default 1;
param norm;
param step;
param same default 0;
param same_limit := 3;
param iter_limit := 20;
param LBlog {0..iter_limit}; let LBlog[0] := LB;
param UBlog {0..iter_limit}; let UBlog[0] := UB;
param scalelog {1..iter_limit};
param steplog {1..iter_limit};
for {k in 1..iter_limit} { printf "\nITERATION %d\n\n", k;
solve LowerBound;
let {j in CITY} slack[j] := sum {i in CITY} Ship[i,j] - demand[j];
if Lagrangian > LB + 0.000001 then {
let LB := Lagrangian;
let same := 0; }
else let same := same + 1;
if same = same_limit then {
let scale := scale / 2;
let same := 0;
};
let norm := sum {j in CITY} slack[j]^2;
let step := scale * (UB - Lagrangian) / norm;
let {j in CITY} mult[j] := mult[j] less step * slack[j];
if sum {i in CITY} supply[i] * Build[i]
>= sum {j in CITY} demand[j] - 1e-8 then {
solve UpperBound;
let UB := min (UB, Shipping_Cost);
}
let LBlog[k] := LB;
let UBlog[k] := UB;
let scalelog[k] := scale;
let steplog[k] := step;
}
printf "\n\n";
display LBlog, UBlog, scalelog, steplog;# ----------------------------------------
# LOCATION-TRANSPORTATION PROBLEM
# USING LAGRANGIAN RELAXATION
# (with separate Supply, Build_Def constr)
# ----------------------------------------
set CITY;
param build_limit integer;
param demand {i in CITY} integer > 0;
param supply {CITY} integer > 0;
param ship_cost {i in CITY, j in CITY} >= 0;
param mult {CITY} >= 0; # Lagrange multipliers for Demand constr
var Build {CITY} integer >= 0 <= 1; # = 1 iff warehouse built at i
var Ship {i in CITY, j in CITY} # amounts shipped
>= 0, <= demand[j];
minimize Lagrangian:
sum {i in CITY, j in CITY} ship_cost[i,j] * Ship[i,j] +
sum {j in CITY} mult[j] * (demand[j] - sum {i in CITY} Ship[i,j]);
minimize Shipping_Cost:
sum {i in CITY, j in CITY} ship_cost[i,j] * Ship[i,j];
subj to Supply {i in CITY}:
sum {j in CITY} Ship[i,j] <= supply[i];
subj to Demand {j in CITY}:
sum {i in CITY} Ship[i,j] >= demand[j];
subj to Build_Def {i in CITY, j in CITY}:
Ship[i,j] <= demand[j] * Build[i];
subj to Limit: sum {i in CITY} Build[i] <= build_limit;param build_limit := 8;
param: CITY: supply demand :=
1 19 4
2 13 1
3 20 2
4 18 7
5 23 9
6 25 9
7 22 3
8 28 9
9 18 5
10 16 2
11 17 4
12 14 1
13 22 9
14 13 3
15 18 1
16 20 7
17 15 5
18 23 9
19 25 2
20 19 4
21 18 9
22 19 6
23 22 5
24 22 4
25 21 5
26 19 6
27 17 8
28 21 9
29 15 5
30 25 6 ;
param ship_cost
: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 :=
1 0 61 35 35 53 16 11 58 28 34 43 18 48 30 68
2 61 0 96 88 29 77 50 92 57 75 48 47 76 32 77
3 35 96 0 18 88 21 45 72 58 53 65 49 51 66 93
4 35 88 18 0 87 27 42 84 62 62 51 42 33 62 100
5 53 29 88 87 0 67 45 66 37 54 62 47 85 26 49
6 16 77 21 27 67 0 27 57 36 34 55 32 51 46 73
7 11 50 45 42 45 27 0 64 29 40 35 8 46 20 68
8 58 92 72 84 66 57 64 0 37 24 99 72 106 65 31
9 28 57 58 62 37 36 29 37 0 18 62 36 74 28 40
10 34 75 53 62 54 34 40 24 18 0 75 48 81 45 41
11 43 48 65 51 62 55 35 99 62 75 0 27 29 39 99
12 18 47 49 42 47 32 8 72 36 48 27 0 40 21 74
13 48 76 51 33 85 51 46 106 74 81 29 40 0 60 114
14 30 32 66 62 26 46 20 65 28 45 39 21 60 0 60
15 68 77 93 100 49 73 68 31 40 41 99 74 114 60 0
16 15 74 22 20 68 10 24 67 41 43 46 27 41 44 80
17 27 34 61 55 35 43 16 72 34 50 30 13 51 10 69
18 53 75 76 85 48 57 55 19 26 24 88 62 100 51 17
19 26 87 14 26 77 10 37 60 45 39 63 42 55 56 79
20 32 68 40 25 73 37 31 90 59 66 27 27 15 48 99
21 13 70 33 39 57 12 23 47 24 23 57 31 59 38 61
22 41 38 76 76 14 54 35 54 23 40 59 39 79 19 42
23 54 8 88 81 27 69 43 87 51 69 41 39 69 25 75
24 36 87 43 55 67 29 44 30 30 15 79 52 80 55 53
25 51 71 76 83 44 57 52 23 24 24 85 60 98 48 17
26 68 82 91 100 54 73 69 26 41 39 101 77 115 63 6
27 29 39 64 63 24 43 22 56 19 37 48 26 66 10 51
28 20 78 17 16 73 12 29 69 46 46 49 31 40 49 85
29 49 104 43 59 85 38 60 38 48 32 92 67 89 72 67
30 48 106 33 51 90 34 59 49 54 39 89 65 82 74 77
: 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 :=
1 15 27 53 26 32 13 41 54 36 51 68 29 20 49 48
2 74 34 75 87 68 70 38 8 87 71 82 39 78 104 106
3 22 61 76 14 40 33 76 88 43 76 91 64 17 43 33
4 20 55 85 26 25 39 76 81 55 83 100 63 16 59 51
5 68 35 48 77 73 57 14 27 67 44 54 24 73 85 90
6 10 43 57 10 37 12 54 69 29 57 73 43 12 38 34
7 24 16 55 37 31 23 35 43 44 52 69 22 29 60 59
8 67 72 19 60 90 47 54 87 30 23 26 56 69 38 49
9 41 34 26 45 59 24 23 51 30 24 41 19 46 48 54
10 43 50 24 39 66 23 40 69 15 24 39 37 46 32 39
11 46 30 88 63 27 57 59 41 79 85 101 48 49 92 89
12 27 13 62 42 27 31 39 39 52 60 77 26 31 67 65
13 41 51 100 55 15 59 79 69 80 98 115 66 40 89 82
14 44 10 51 56 48 38 19 25 55 48 63 10 49 72 74
15 80 69 17 79 99 61 42 75 53 17 6 51 85 67 77
16 0 39 65 17 26 20 56 66 39 63 80 44 5 48 43
17 39 0 59 53 39 38 29 27 58 56 72 18 44 74 74
18 65 59 0 63 85 45 37 70 37 4 15 42 69 51 60
19 17 53 63 0 42 21 64 80 30 63 78 53 15 34 27
20 26 39 85 42 0 44 65 61 65 83 100 52 27 75 69
21 20 38 45 21 44 0 44 63 22 44 61 34 24 36 36
22 56 29 37 64 65 44 0 34 53 33 46 13 61 71 76
23 66 27 70 80 61 63 34 0 80 67 79 32 71 97 99
24 39 58 37 30 65 22 53 80 0 38 51 48 41 18 25
25 63 56 4 63 83 44 33 67 38 0 17 38 68 53 62
26 80 72 15 78 100 61 46 79 51 17 0 54 84 63 73
27 44 18 42 53 52 34 13 32 48 38 54 0 49 66 69
28 5 44 69 15 27 24 61 71 41 68 84 49 0 48 42
29 48 74 51 34 75 36 71 97 18 53 63 66 48 0 12
30 43 74 60 27 69 36 76 99 25 62 73 69 42 12 0 ;Lagrangian relaxation for a location-transportation problem: LP relaxation bound is poor, and subproblem has the integrality property so no improvement can be made
# ----------------------------------------
# LAGRANGIAN RELAXATION FOR
# THE LOCATION-TRANSPORTATION PROBLEM
# ----------------------------------------
printf "\nLP RELAXATION\n\n";
model trnloc2a.mod;
data trnloc2.dat;
option omit_zero_rows 1;
option display_eps .000001;
option solution_round 8;
option solver cplex;
option solver_msg 0;
option relax_integrality 1;
objective Shipping_Cost;
solve;
param LB; param UB;
let LB := Shipping_Cost.val;
let UB := sum {j in CITY} max {i in CITY} ship_cost[i,j];
option relax_integrality 0;
problem LowerBound: Build, Ship, Supply, Limit, Lagrangian;
problem UpperBound: Ship, Supply, Demand, Limit, Shipping_Cost;
let {j in CITY} mult[j] := 0;
param slack {CITY};
param scale default 1;
param norm;
param step;
param same default 0;
param same_limit := 3;
param iter_limit := 20;
param LBlog {0..iter_limit}; let LBlog[0] := LB;
param UBlog {0..iter_limit}; let UBlog[0] := UB;
param scalelog {1..iter_limit};
param steplog {1..iter_limit};
for {k in 1..iter_limit} { printf "\nITERATION %d\n\n", k;
solve LowerBound;
let {j in CITY} slack[j] := sum {i in CITY} Ship[i,j] - demand[j];
if Lagrangian > LB + 0.000001 then {
let LB := Lagrangian;
let same := 0; }
else let same := same + 1;
if same = same_limit then {
let scale := scale / 2;
let same := 0;
};
let norm := sum {j in CITY} slack[j]^2;
let step := scale * (UB - Lagrangian) / norm;
let {j in CITY} mult[j] := mult[j] less step * slack[j];
if sum {i in CITY} supply[i] * Build[i]
>= sum {j in CITY} demand[j] - 1e-8 then {
solve UpperBound;
let UB := min (UB, Shipping_Cost);
}
let LBlog[k] := LB;
let UBlog[k] := UB;
let scalelog[k] := scale;
let steplog[k] := step;
}
printf "\n\n";
display LBlog, UBlog, scalelog, steplog;# ----------------------------------------
# LOCATION-TRANSPORTATION PROBLEM
# USING LAGRANGIAN RELAXATION
# ----------------------------------------
set CITY;
param build_limit integer;
param demand {i in CITY} integer > 0;
param supply {CITY} integer > 0;
param ship_cost {i in CITY, j in CITY} >= 0;
param mult {CITY} >= 0; # Lagrange multipliers for Demand constr
var Build {CITY} integer >= 0 <= 1; # = 1 iff warehouse built at i
var Ship {i in CITY, j in CITY} >= 0; # amounts shipped
minimize Lagrangian:
sum {i in CITY, j in CITY} ship_cost[i,j] * Ship[i,j] +
sum {j in CITY} mult[j] * (demand[j] - sum {i in CITY} Ship[i,j]);
minimize Shipping_Cost:
sum {i in CITY, j in CITY} ship_cost[i,j] * Ship[i,j];
subj to Supply {i in CITY}:
sum {j in CITY} Ship[i,j] <= supply[i] * Build[i];
subj to Demand {j in CITY}:
sum {i in CITY} Ship[i,j] >= demand[j];
subj to Limit: sum {i in CITY} Build[i] <= build_limit;param build_limit := 8;
param: CITY: supply demand :=
1 19 4
2 13 1
3 20 2
4 18 7
5 23 9
6 25 9
7 22 3
8 28 9
9 18 5
10 16 2
11 17 4
12 14 1
13 22 9
14 13 3
15 18 1
16 20 7
17 15 5
18 23 9
19 25 2
20 19 4
21 18 9
22 19 6
23 22 5
24 22 4
25 21 5
26 19 6
27 17 8
28 21 9
29 15 5
30 25 6 ;
param ship_cost
: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 :=
1 0 61 35 35 53 16 11 58 28 34 43 18 48 30 68
2 61 0 96 88 29 77 50 92 57 75 48 47 76 32 77
3 35 96 0 18 88 21 45 72 58 53 65 49 51 66 93
4 35 88 18 0 87 27 42 84 62 62 51 42 33 62 100
5 53 29 88 87 0 67 45 66 37 54 62 47 85 26 49
6 16 77 21 27 67 0 27 57 36 34 55 32 51 46 73
7 11 50 45 42 45 27 0 64 29 40 35 8 46 20 68
8 58 92 72 84 66 57 64 0 37 24 99 72 106 65 31
9 28 57 58 62 37 36 29 37 0 18 62 36 74 28 40
10 34 75 53 62 54 34 40 24 18 0 75 48 81 45 41
11 43 48 65 51 62 55 35 99 62 75 0 27 29 39 99
12 18 47 49 42 47 32 8 72 36 48 27 0 40 21 74
13 48 76 51 33 85 51 46 106 74 81 29 40 0 60 114
14 30 32 66 62 26 46 20 65 28 45 39 21 60 0 60
15 68 77 93 100 49 73 68 31 40 41 99 74 114 60 0
16 15 74 22 20 68 10 24 67 41 43 46 27 41 44 80
17 27 34 61 55 35 43 16 72 34 50 30 13 51 10 69
18 53 75 76 85 48 57 55 19 26 24 88 62 100 51 17
19 26 87 14 26 77 10 37 60 45 39 63 42 55 56 79
20 32 68 40 25 73 37 31 90 59 66 27 27 15 48 99
21 13 70 33 39 57 12 23 47 24 23 57 31 59 38 61
22 41 38 76 76 14 54 35 54 23 40 59 39 79 19 42
23 54 8 88 81 27 69 43 87 51 69 41 39 69 25 75
24 36 87 43 55 67 29 44 30 30 15 79 52 80 55 53
25 51 71 76 83 44 57 52 23 24 24 85 60 98 48 17
26 68 82 91 100 54 73 69 26 41 39 101 77 115 63 6
27 29 39 64 63 24 43 22 56 19 37 48 26 66 10 51
28 20 78 17 16 73 12 29 69 46 46 49 31 40 49 85
29 49 104 43 59 85 38 60 38 48 32 92 67 89 72 67
30 48 106 33 51 90 34 59 49 54 39 89 65 82 74 77
: 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 :=
1 15 27 53 26 32 13 41 54 36 51 68 29 20 49 48
2 74 34 75 87 68 70 38 8 87 71 82 39 78 104 106
3 22 61 76 14 40 33 76 88 43 76 91 64 17 43 33
4 20 55 85 26 25 39 76 81 55 83 100 63 16 59 51
5 68 35 48 77 73 57 14 27 67 44 54 24 73 85 90
6 10 43 57 10 37 12 54 69 29 57 73 43 12 38 34
7 24 16 55 37 31 23 35 43 44 52 69 22 29 60 59
8 67 72 19 60 90 47 54 87 30 23 26 56 69 38 49
9 41 34 26 45 59 24 23 51 30 24 41 19 46 48 54
10 43 50 24 39 66 23 40 69 15 24 39 37 46 32 39
11 46 30 88 63 27 57 59 41 79 85 101 48 49 92 89
12 27 13 62 42 27 31 39 39 52 60 77 26 31 67 65
13 41 51 100 55 15 59 79 69 80 98 115 66 40 89 82
14 44 10 51 56 48 38 19 25 55 48 63 10 49 72 74
15 80 69 17 79 99 61 42 75 53 17 6 51 85 67 77
16 0 39 65 17 26 20 56 66 39 63 80 44 5 48 43
17 39 0 59 53 39 38 29 27 58 56 72 18 44 74 74
18 65 59 0 63 85 45 37 70 37 4 15 42 69 51 60
19 17 53 63 0 42 21 64 80 30 63 78 53 15 34 27
20 26 39 85 42 0 44 65 61 65 83 100 52 27 75 69
21 20 38 45 21 44 0 44 63 22 44 61 34 24 36 36
22 56 29 37 64 65 44 0 34 53 33 46 13 61 71 76
23 66 27 70 80 61 63 34 0 80 67 79 32 71 97 99
24 39 58 37 30 65 22 53 80 0 38 51 48 41 18 25
25 63 56 4 63 83 44 33 67 38 0 17 38 68 53 62
26 80 72 15 78 100 61 46 79 51 17 0 54 84 63 73
27 44 18 42 53 52 34 13 32 48 38 54 0 49 66 69
28 5 44 69 15 27 24 61 71 41 68 84 49 0 48 42
29 48 74 51 34 75 36 71 97 18 53 63 66 48 0 12
30 43 74 60 27 69 36 76 99 25 62 73 69 42 12 0 ;