Refinery Optimization: How to Run an Oil Refinery

The NOMISS option in the first READ DATA statement ensures that only nonmissing values of the variable multiplier in the arc_data data set populate the arc_mult parameter. Because arc_mult is declared with an initial value of 1, parameters with no value default to 1.

proc optmodel;
   set <str,str> ARCS;
   num arc_mult {ARCS} init 1;
   read data arc_data nomiss into ARCS=[i j] arc_mult=multiplier;
   var Flow {ARCS} >= 0;

   set <str> FINAL_PRODUCTS;
   num profit {FINAL_PRODUCTS};
   read data final_product_data into FINAL_PRODUCTS=[product] profit;

The following FOR loop converts profit from pence per barrel to pounds per barrel, without altering the input data set:

   for {product in FINAL_PRODUCTS} profit[product] = profit[product] / 100;

Most arcs appear in the arc_data data set, but the following assignment statement uses the set operators UNION and CROSS to augment the ARCS set with an arc from each final product to the sink node:

   ARCS = ARCS union (FINAL_PRODUCTS cross {'sink'});

   set NODES = union {<i,j> in ARCS} {i,j};

   max TotalProfit = sum {i in FINAL_PRODUCTS} profit[i] * Flow[i,'sink'];

   con Flow_balance {i in NODES diff {'source','sink'}}:
      sum {<(i),j> in ARCS} Flow[i,j]
         = sum {<j,(i)> in ARCS} arc_mult[j,i] * Flow[j,i];

   set <str> CRUDES;
   var CrudeDistilled {CRUDES} >= 0;

Because the decision variable CrudeDistilled automatically contains the .ub suffix, you can populate CrudeDistilled.ub directly by using the following READ DATA statement, without having to declare an additional parameter to store this upper bound:

   read data crude_data into CRUDES=[crude] CrudeDistilled.ub=crude_ub;
   con Distillation {<i,j> in ARCS: i in CRUDES}:
      Flow[i,j] = CrudeDistilled[i];

The SETOF operator, used together with the concatenation operator () in the following statements, enables you to create one index set from another:

   set OILS = {'light_oil','heavy_oil'};
   set CRACKED_OILS = setof {i in OILS} i||'_cracked';
   var OilCracked {CRACKED_OILS} >= 0;
   con Cracking {<i,j> in ARCS: i in CRACKED_OILS}:
      Flow[i,j] = OilCracked[i];

   num octane {NODES} init .;
   read data octane_data nomiss into [i] octane;

   set <str> PETROLS;
   num octane_lb {PETROLS};
   read data petrol_data into PETROLS=[petrol] octane_lb;

   num vapour_pressure {NODES} init .;
   read data vapour_pressure_data nomiss into [oil] vapour_pressure;

As expressed on , both octane numbers and vapour pressures of the blended fuels are ratios of linear functions of the decision variables. To increase algorithmic performance and reliability, the following two CON statements linearize the nonlinear ratio constraints by clearing the denominators:

   con Blending_petrol {petrol in PETROLS}:
      sum {<i,(petrol)> in ARCS}
         octane[i] * arc_mult[i,petrol] * Flow[i,petrol]
   >= octane_lb[petrol] *
      sum {<i,(petrol)> in ARCS} arc_mult[i,petrol] * Flow[i,petrol];

   con Blending_jet_fuel:
      sum {<i,'jet_fuel'> in ARCS}
         vapour_pressure[i] * arc_mult[i,'jet_fuel'] * Flow[i,'jet_fuel']
   <= &vapour_pressure_ub *
      sum {<i,'jet_fuel'> in ARCS} arc_mult[i,'jet_fuel'] * Flow[i,'jet_fuel'];

Similarly, the following CON statement linearizes the fuel oil ratio constraints by clearing the denominators:

   num fuel_oil_coefficient {NODES} init 0;
   read data fuel_oil_ratio_data nomiss into [oil]
      fuel_oil_coefficient=coefficient;
   num sum_fuel_oil_coefficient
      = sum {<i,'fuel_oil'> in ARCS} fuel_oil_coefficient[i];
   con Blending_fuel_oil {<i,'fuel_oil'> in ARCS}:
      sum_fuel_oil_coefficient * Flow[i,'fuel_oil']
   = fuel_oil_coefficient[i] * sum {<j,'fuel_oil'> in ARCS} Flow[j,'fuel_oil'];

   con Crude_total_ub_con:
      sum {i in CRUDES} CrudeDistilled[i] <= &crude_total_ub;

The following CON statement uses the SAS function INDEX together with the colon operator (:) to select the subset of arcs whose tail node contains ‘naphtha’ in its name:

   con Naphtha_ub_con:
      sum {<i,'reformed_gasoline'> in ARCS: index(i,'naphtha') > 0}
         Flow[i,'reformed_gasoline']
   <= &naphtha_ub;

   con Cracked_oil_ub_con:
      sum {<i,'cracked_oil'> in ARCS} Flow[i,'cracked_oil'] <= &cracked_oil_ub;

   con Lube_oil_range_con:
      &lube_oil_lb <= Flow['lube_oil','sink'] <= &lube_oil_ub;

As expressed on , the premium ratio constraint involves a ratio of linear functions of the decision variables. The following CON statement linearizes the nonlinear ratio constraint by clearing the denominator:

   con Premium_ratio_con:
      sum {<'premium_petrol',j> in ARCS} Flow['premium_petrol',j]
   >= &premium_ratio *
      sum {<'regular_petrol',j> in ARCS} Flow['regular_petrol',j];

   num octane_sol {petrol in PETROLS} =
      (sum {<i,(petrol)> in ARCS}
         octane[i] * arc_mult[i,petrol] * Flow[i,petrol].sol) /
      (sum {<i,(petrol)> in ARCS} arc_mult[i,petrol] * Flow[i,petrol].sol);

   num vapour_pressure_sol =
      (sum {<i,'jet_fuel'> in ARCS} vapour_pressure[i] *
         arc_mult[i,'jet_fuel'] * Flow[i,'jet_fuel'].sol)
      / (sum {<i,'jet_fuel'> in ARCS} arc_mult[i,'jet_fuel'] *
         Flow[i,'jet_fuel'].sol);

   num fuel_oil_ratio_sol {<i,'fuel_oil'> in ARCS} =
      (arc_mult[i,'fuel_oil'] * Flow[i,'fuel_oil'].sol) /
      (sum {<j,'fuel_oil'> in ARCS} arc_mult[j,'fuel_oil'] *
         Flow[j,'fuel_oil'].sol);

   solve;
   print CrudeDistilled;
   print OilCracked Flow;
   print octane_sol octane_lb;

Although the previous PRINT statements print all values of the given parameters, the following two PRINT statements use an index set to print a specified subset of the values:

   print {<i,'jet_fuel'> in ARCS} vapour_pressure vapour_pressure_sol;
   print {<i,'fuel_oil'> in ARCS} fuel_oil_coefficient fuel_oil_ratio_sol;
   create data sol_data1 from [i j] Flow;
quit;

Figure 6.2 shows the output from the linear programming solver. For this test instance, it turns out that the optimal solution contains no fuel oil.