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

AMPL SUITE

Python API - amplpy

Stage

Profession

Industry

AMPL IN ACTION→

AMPL SUITE

Python API - amplpy

Stage

Profession

Industry

AMPL IN ACTION→

Additional Scripts: Looping and Testing – 2

Implements

Script

Uses