Cutting Stock (omod5)
/***************************************************************/
/* */
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: omod5 */
/* TITLE: Cutting Stock (omod5) */
/* PRODUCT: OR */
/* SYSTEM: ALL */
/* KEYS: OR */
/* PROCS: OPTMODEL */
/* DATA: */
/* */
/* SUPPORT: UPDATE: */
/* REF: */
/* MISC: Example 5 from the OPTMODEL chapter of */
/* Mathematical Programming. */
/* */
/***************************************************************/
/* cutting-stock problem */
/* uses delayed column generation from
Chvatal's Linear Programming (1983), page 198 */
%macro csp(capacity);
proc optmodel printlevel=0;
/* declare parameters and sets */
num capacity = &capacity;
set ITEMS;
num demand {ITEMS};
num width {ITEMS};
num num_patterns init card(ITEMS);
set PATTERNS = 1..num_patterns;
num a {ITEMS, PATTERNS};
num c {ITEMS} init 0;
num epsilon = 1E-6;
/* read input data */
read data indata into ITEMS=[_N_] demand width;
/* generate trivial initial columns */
for {i in ITEMS, j in PATTERNS}
a[i,j] = (if (i = j) then floor(capacity/width[i]) else 0);
/* define master problem */
var x {PATTERNS} >= 0 integer;
minimize NumberOfRaws = sum {j in PATTERNS} x[j];
con demand_con {i in ITEMS}:
sum {j in PATTERNS} a[i,j] * x[j] >= demand[i];
problem Master include x NumberOfRaws demand_con;
/* define column generation subproblem */
var y {ITEMS} >= 0 integer;
maximize KnapsackObjective = sum {i in ITEMS} c[i] * y[i];
con knapsack_con:
sum {i in ITEMS} width[i] * y[i] <= capacity;
problem Knapsack include y KnapsackObjective knapsack_con;
/* main loop */
do while (1);
print _page_ a;
/* master problem */
/* minimize sum_j x[j]
subj. to sum_j a[i,j] * x[j] >= demand[i]
x[j] >= 0 and integer */
use problem Master;
put "solve master problem";
solve with lp relaxint /
presolver=none solver=ps basis=warmstart printfreq=1;
print x;
print demand_con.dual;
for {i in ITEMS} c[i] = demand_con[i].dual;
/* knapsack problem */
/* maximize sum_i c[i] * y[i]
subj. to sum_i width[i] * y[i] <= capacity
y[i] >= 0 and integer */
use problem Knapsack;
put "solve column generation subproblem";
solve with milp / printfreq=0;
for {i in ITEMS} y[i] = round(y[i]);
print y;
print KnapsackObjective;
if KnapsackObjective <= 1 + epsilon then leave;
/* include new pattern */
num_patterns = num_patterns + 1;
for {i in ITEMS} a[i,num_patterns] = y[i];
end;
/* solve IP, using rounded-up LP solution as warm start */
use problem Master;
for {j in PATTERNS} x[j] = ceil(x[j].sol);
put "solve (restricted) master problem as IP";
solve with milp / primalin;
/* cleanup solution and save to output data set */
for {j in PATTERNS} x[j] = round(x[j].sol);
create data solution from [pattern]={j in PATTERNS: x[j] > 0}
raws=x {i in ITEMS} <col('i'||i)=a[i,j]>;
quit;
%mend csp;
/* Chvatal, p.199 */
data indata;
input demand width;
datalines;
78 25.5
40 22.5
30 20
30 15
;
run;
%csp(91)
/* LP solution is integer */
/* Chvatal, p.195 */
data indata;
input demand width;
datalines;
97 45
610 36
395 31
211 14
;
run;
%csp(100)
/* LP solution is fractional */