Additional AMPL Scripts: Enhance Your Modeling

Additional Scripts: Looping and Testing – 2

Implementing algorithms through AMPL scripts

Implements

Script

Uses

Gilmore-Gomory column generation procedure for the cutting-stock (roll trim) problem

cut1.run

# ----------------------------------------
# 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;

Open file in a new tab

cut1.mod

# ----------------------------------------
# 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;

Open file in a new tab

cut.dat

param roll_width := 110 ;

param: WIDTHS: orders :=
          20     48
          45     35
          50     24
          55     10
          75      8  ;

Open file in a new tab

Same as cut1.run, but using an alternative arrangement wherein problems are defined immediately before before their members are declared

cut2.run

# ----------------------------------------
# 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;

Open file in a new tab

cut2.mod

  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;

Open file in a new tab

cut.dat

param roll_width := 110 ;

param: WIDTHS: orders :=
          20     48
          45     35
          50     24
          55     10
          75      8  ;

Open file in a new tab

Same as cut1.run, but with better formatting of output

cut3.run

# ----------------------------------------
# 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));

Open file in a new tab

cut1.mod

# ----------------------------------------
# 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;

Open file in a new tab

cut.dat

param roll_width := 110 ;

param: WIDTHS: orders :=
          20     48
          45     35
          50     24
          55     10
          75      8  ;

Open file in a new tab

Dantzig-Wolfe decomposition for a multi-commodity transportation problem, using a single subproblem

multi1.run

# ----------------------------------------
# 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;

Open file in a new tab

multi1.mod

# ----------------------------------------
# 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];

Open file in a new tab

multi1.dat

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 ;

Open file in a new tab

Same as multi1.run, but using the same repeat loop for both phase I (infeasible) and phase II (feasible).

multi1a.run

# ----------------------------------------
# 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;

Open file in a new tab

multi1.mod

# ----------------------------------------
# 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];

Open file in a new tab

multi1.dat

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 ;

Open file in a new tab

Same as multi1.run, but using a separate subproblem for each product; subproblems are represented in AMPL by an indexed collection of named problems

multi2.run

# ----------------------------------------
# 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;

Open file in a new tab

multi2.mod

# ----------------------------------------
# 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;

Open file in a new tab

multi2.dat

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 ;

Open file in a new tab

Same as multi2.run, except that the separate subproblems are realized by changing the data to a single AMPL named problem

multi3.run

# ----------------------------------------
# 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;

Open file in a new tab

multi3.mod

# ----------------------------------------
# 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;

Open file in a new tab

multi3.dat

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 ;

Open file in a new tab

Benders decomposition for a stochastic programming variant of a multi-period production problem (see Exercise 4-5)

stoch1.run

# ----------------------------------------
# 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]);
   }

Open file in a new tab

stoch1.mod

# ----------------------------------------
# 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];

Open file in a new tab

stoch.dat

# ----------------------------------------
# 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]);
   }

Open file in a new tab

Same as stoch1.run, but using a separate subproblem for each scenario; subproblems are represented in AMPL by an indexed collection of named problems

stoch2.run

# ----------------------------------------
# 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]);
   }

Open file in a new tab

stoch2.mod

# ----------------------------------------
# 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];

Open file in a new tab

stoch.dat

# ----------------------------------------
# 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]);
   }

Open file in a new tab

Same as stoch2.run, except that the separate subproblems are realized by changing the data to a single AMPL named problem

stoch3.run

# ----------------------------------------
# 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]);
   }

Open file in a new tab

stoch3.mod

# ----------------------------------------
# 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];

Open file in a new tab

stoch.dat

# ----------------------------------------
# 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]);
   }

Open file in a new tab

Benders decomposition for a location-transportation problem (original model in trnloc.mod)

trnloc1.run

# ----------------------------------------
# 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;

Open file in a new tab

trnloc1.mod

# ----------------------------------------
# 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];

Open file in a new tab

trnloc.dat

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

trnloc2b.run

# ----------------------------------------
# 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;

Open file in a new tab

trnloc2b.mod

# ----------------------------------------
# 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;

Open file in a new tab

trnloc2.dat

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 ;

Open file in a new tab

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

trnloc2c.run

# ----------------------------------------
# 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;

Open file in a new tab

trnloc2c.mod

# ----------------------------------------
# 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;

Open file in a new tab

trnloc2.dat

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 ;

Open file in a new tab

Lagrangian relaxation for a location-transportation problem: LP relaxation bound is poor, and subproblem has the integrality property so no improvement can be made

trnloc2a.run

# ----------------------------------------
# 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;

Open file in a new tab

trnloc2a.mod

# ----------------------------------------
# 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;

Open file in a new tab

trnloc2.dat

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 ;

DOCUMENTATION & RESOURCES

OPEN SOURCE PROJECTS

Linear

Nonlinear

Global

Docs

Models and Code

Our Company

AMPL in Practice

Pricing

Not quite sure what license is right for you?

Docs

Models & Code

Our Company

AMPL in Practice

Free Licenses

Paid Licenses

Get Support

Linear

Nonlinear

Global

Docs

Models and Code

Our Company

AMPL in Practice

Pricing

Not quite sure what license is right for you?

Docs

Models & Code

Our Company

AMPL in Practice

Free Licenses

Paid Licenses

Get Support

Additional Scripts: Looping and Testing – 2

Implements

Script

Uses

Free Licenses

Paid Licenses

LEARN

Models & Code

RESOURCES

COMPANY