The first several PROC OPTMODEL statements declare index sets and parameters and then read the input data:

proc optmodel;
   set PERIODS;
   num length {PERIODS};
   num demand {PERIODS};
   read data period_data into PERIODS=[_N_] length demand;

   set TYPES;
   num num_avail {TYPES};
   num min_level {TYPES};
   num max_level {TYPES};
   num unit_cost {TYPES};
   num excess_cost {TYPES};
   num startup_cost {TYPES};
   read data type_data into TYPES=[_N_]
      num_avail min_level max_level unit_cost excess_cost startup_cost;

The following VAR statements declare decision variables, with bounds that depend on type:

   var NumWorking {PERIODS, type in TYPES} >= 0 <= num_avail[type] integer;
   var Excess {PERIODS, TYPES} >= 0;
   var NumStartup {PERIODS, type in TYPES} >= 0 <= num_avail[type] integer;

The following IMPVAR statement declares Output as an implicit variable:

   impvar Output {period in PERIODS, type in TYPES} = 
      min_level[type] * NumWorking[period,type] + Excess[period,type];

The following MIN statement declares the objective:

   min TotalCost =
      sum {period in PERIODS, type in TYPES} (
         unit_cost[type] * length[period] * NumWorking[period,type]
         + excess_cost[type] * length[period] * Excess[period,type]
         + startup_cost[type] * NumStartup[period,type]);

The following CON statements declare the constraints, with an IF-THEN/ELSE expression in the last constraint to account for a boundary condition in the first period:

   con Demand_con {period in PERIODS}:
      sum {type in TYPES} Output[period,type]
   >= demand[period];

   con Reserve_con {period in PERIODS}:
      sum {type in TYPES} max_level[type] * NumWorking[period,type]
   >= (1 + &reserve) * demand[period];

   con Excess_ub {period in PERIODS, type in TYPES}:
      Excess[period,type] 
   <= (max_level[type] - min_level[type]) * NumWorking[period,type];

   con Startup_con {period in PERIODS, type in TYPES}:
      NumStartup[period,type] 
   >= NumWorking[period,type] 
    - (if period - 1 in PERIODS then NumWorking[period-1,type] 
       else NumWorking[card(PERIODS),type]);

The following statements call the mixed integer linear programming solver, print the optimal solution, and create a data set that contains the optimal solution, with one observation per period-type pair:

   solve;
   print NumWorking NumStartup Excess Output;
   create data sol_data from [period type] NumWorking NumStartup Excess Output;

Figure 16.1 shows the output from the mixed integer linear programming solver.

The following statements fix the integer variables to their optimal values, call the linear programming solver, and use the .dual constraint suffix to compute marginal costs. The RELAXINT option in the SOLVE statement relaxes the integer variables to be continuous.

   fix NumWorking;
   fix NumStartup;
   solve with LP relaxint;
   print NumWorking NumStartup Excess Output;
   print {period in PERIODS} (demand_con[period].dual / length[period]);

Figure 16.2 shows the output from the linear programming solver when the integer variables are fixed to their optimal values. As expected, the same optimal solution is obtained; the point of calling the LP solver is to obtain the dual variables used to compute the marginal costs, as shown in Figure 16.3.

The following statements unfix the integer variables, call the linear programming solver, and use the .dual constraint suffix to compute marginal costs:

   unfix NumWorking;
   unfix NumStartup;
   solve with LP relaxint;
   print NumWorking NumStartup Excess Output;
   print {period in PERIODS} (demand_con[period].dual / length[period]);
   print {period in PERIODS} (reserve_con[period].dual / length[period]);
quit;

Figure 16.4 shows the output from the linear programming solver when the integer variables are unfixed and relaxed to be continuous. As expected, some relaxed integer variables now take fractional values, and the total cost is slightly less than before. Figure 16.5 and Figure 16.6 show the resulting marginal costs.