Refinery Optimization¶
Reference¶
SAS/OR example: http://go.documentation.sas.com/?docsetId=ormpex&docsetTarget=ormpex_ex6_toc.htm&docsetVersion=15.1&locale=en
SAS/OR code for example: http://support.sas.com/documentation/onlinedoc/or/ex_code/151/mpex06.html
Model¶
import sasoptpy as so
import pandas as pd
import numpy as np
def test(cas_conn, **kwargs):
m = so.Model(name='refinery_optimization', session=cas_conn)
crude_data = pd.DataFrame([
['crude1', 20000],
['crude2', 30000]
], columns=['crude', 'crude_ub']).set_index(['crude'])
arc_data = pd.DataFrame([
['source', 'crude1', 6],
['source', 'crude2', 6],
['crude1', 'light_naphtha', 0.1],
['crude1', 'medium_naphtha', 0.2],
['crude1', 'heavy_naphtha', 0.2],
['crude1', 'light_oil', 0.12],
['crude1', 'heavy_oil', 0.2],
['crude1', 'residuum', 0.13],
['crude2', 'light_naphtha', 0.15],
['crude2', 'medium_naphtha', 0.25],
['crude2', 'heavy_naphtha', 0.18],
['crude2', 'light_oil', 0.08],
['crude2', 'heavy_oil', 0.19],
['crude2', 'residuum', 0.12],
['light_naphtha', 'regular_petrol', np.nan],
['light_naphtha', 'premium_petrol', np.nan],
['medium_naphtha', 'regular_petrol', np.nan],
['medium_naphtha', 'premium_petrol', np.nan],
['heavy_naphtha', 'regular_petrol', np.nan],
['heavy_naphtha', 'premium_petrol', np.nan],
['light_naphtha', 'reformed_gasoline', 0.6],
['medium_naphtha', 'reformed_gasoline', 0.52],
['heavy_naphtha', 'reformed_gasoline', 0.45],
['light_oil', 'jet_fuel', np.nan],
['light_oil', 'fuel_oil', np.nan],
['heavy_oil', 'jet_fuel', np.nan],
['heavy_oil', 'fuel_oil', np.nan],
['light_oil', 'light_oil_cracked', 2],
['light_oil_cracked', 'cracked_oil', 0.68],
['light_oil_cracked', 'cracked_gasoline', 0.28],
['heavy_oil', 'heavy_oil_cracked', 2],
['heavy_oil_cracked', 'cracked_oil', 0.75],
['heavy_oil_cracked', 'cracked_gasoline', 0.2],
['cracked_oil', 'jet_fuel', np.nan],
['cracked_oil', 'fuel_oil', np.nan],
['reformed_gasoline', 'regular_petrol', np.nan],
['reformed_gasoline', 'premium_petrol', np.nan],
['cracked_gasoline', 'regular_petrol', np.nan],
['cracked_gasoline', 'premium_petrol', np.nan],
['residuum', 'lube_oil', 0.5],
['residuum', 'jet_fuel', np.nan],
['residuum', 'fuel_oil', np.nan],
], columns=['i', 'j', 'multiplier']).set_index(['i', 'j'])
octane_data = pd.DataFrame([
['light_naphtha', 90],
['medium_naphtha', 80],
['heavy_naphtha', 70],
['reformed_gasoline', 115],
['cracked_gasoline', 105],
], columns=['i', 'octane']).set_index(['i'])
petrol_data = pd.DataFrame([
['regular_petrol', 84],
['premium_petrol', 94],
], columns=['petrol', 'octane_lb']).set_index(['petrol'])
vapour_pressure_data = pd.DataFrame([
['light_oil', 1.0],
['heavy_oil', 0.6],
['cracked_oil', 1.5],
['residuum', 0.05],
], columns=['oil', 'vapour_pressure']).set_index(['oil'])
fuel_oil_ratio_data = pd.DataFrame([
['light_oil', 10],
['cracked_oil', 4],
['heavy_oil', 3],
['residuum', 1],
], columns=['oil', 'coefficient']).set_index(['oil'])
final_product_data = pd.DataFrame([
['premium_petrol', 700],
['regular_petrol', 600],
['jet_fuel', 400],
['fuel_oil', 350],
['lube_oil', 150],
], columns=['product', 'profit']).set_index(['product'])
vapour_pressure_ub = 1
crude_total_ub = 45000
naphtha_ub = 10000
cracked_oil_ub = 8000
lube_oil_lb = 500
lube_oil_ub = 1000
premium_ratio = 0.40
ARCS = arc_data.index.tolist()
arc_mult = arc_data['multiplier'].fillna(1)
FINAL_PRODUCTS = final_product_data.index.tolist()
final_product_data['profit'] = final_product_data['profit'] / 100
profit = final_product_data['profit']
ARCS = ARCS + [(i, 'sink') for i in FINAL_PRODUCTS]
flow = m.add_variables(ARCS, name='flow', lb=0)
NODES = np.unique([i for j in ARCS for i in j])
m.set_objective(so.expr_sum(profit[i] * flow[i, 'sink']
for i in FINAL_PRODUCTS
if (i, 'sink') in ARCS),
name='totalProfit', sense=so.MAX)
m.add_constraints((so.expr_sum(flow[a] for a in ARCS if a[0] == n) ==
so.expr_sum(arc_mult[a] * flow[a]
for a in ARCS if a[1] == n)
for n in NODES if n not in ['source', 'sink']),
name='flow_balance')
CRUDES = crude_data.index.tolist()
crudeDistilled = m.add_variables(CRUDES, name='crudesDistilled', lb=0)
crudeDistilled.set_bounds(ub=crude_data['crude_ub'])
m.add_constraints((flow[i, j] == crudeDistilled[i]
for (i, j) in ARCS if i in CRUDES), name='distillation')
OILS = ['light_oil', 'heavy_oil']
CRACKED_OILS = [i+'_cracked' for i in OILS]
oilCracked = m.add_variables(CRACKED_OILS, name='oilCracked', lb=0)
m.add_constraints((flow[i, j] == oilCracked[i] for (i, j) in ARCS
if i in CRACKED_OILS), name='cracking')
octane = octane_data['octane']
PETROLS = petrol_data.index.tolist()
octane_lb = petrol_data['octane_lb']
vapour_pressure = vapour_pressure_data['vapour_pressure']
m.add_constraints((so.expr_sum(octane[a[0]] * arc_mult[a] * flow[a]
for a in ARCS if a[1] == p)
>= octane_lb[p] *
so.expr_sum(arc_mult[a] * flow[a]
for a in ARCS if a[1] == p)
for p in PETROLS), name='blending_petrol')
m.add_constraint(so.expr_sum(vapour_pressure[a[0]] * arc_mult[a] * flow[a]
for a in ARCS if a[1] == 'jet_fuel') <=
vapour_pressure_ub *
so.expr_sum(arc_mult[a] * flow[a]
for a in ARCS if a[1] == 'jet_fuel'),
name='blending_jet_fuel')
fuel_oil_coefficient = fuel_oil_ratio_data['coefficient']
sum_fuel_oil_coefficient = sum(fuel_oil_coefficient)
m.add_constraints((sum_fuel_oil_coefficient * flow[a] ==
fuel_oil_coefficient[a[0]] * flow.sum('*', ['fuel_oil'])
for a in ARCS if a[1] == 'fuel_oil'),
name='blending_fuel_oil')
m.add_constraint(crudeDistilled.sum('*') <= crude_total_ub,
name='crude_total_ub')
m.add_constraint(so.expr_sum(flow[a] for a in ARCS
if a[0].find('naphtha') > -1 and
a[1] == 'reformed_gasoline')
<= naphtha_ub, name='naphtba_ub')
m.add_constraint(so.expr_sum(flow[a] for a in ARCS if a[1] ==
'cracked_oil') <=
cracked_oil_ub, name='cracked_oil_ub')
m.add_constraint(flow['lube_oil', 'sink'] == [lube_oil_lb, lube_oil_ub],
name='lube_oil_range')
m.add_constraint(flow.sum('premium_petrol', '*') >= premium_ratio *
flow.sum('regular_petrol', '*'), name='premium_ratio')
res = m.solve(**kwargs)
if res is not None:
print(so.get_solution_table(crudeDistilled))
print(so.get_solution_table(oilCracked))
print(so.get_solution_table(flow))
octane_sol = []
for p in PETROLS:
octane_sol.append(so.expr_sum(octane[a[0]] * arc_mult[a] *
flow[a].get_value() for a in ARCS
if a[1] == p) /
sum(arc_mult[a] * flow[a].get_value()
for a in ARCS if a[1] == p))
octane_sol = pd.Series(octane_sol, name='octane_sol', index=PETROLS)
print(so.get_solution_table(octane_sol, octane_lb))
print(so.get_solution_table(vapour_pressure))
vapour_pressure_sol = sum(vapour_pressure[a[0]] *
arc_mult[a] *
flow[a].get_value() for a in ARCS
if a[1] == 'jet_fuel') /\
sum(arc_mult[a] * flow[a].get_value() for a in ARCS
if a[1] == 'jet_fuel')
print('Vapour_pressure_sol: {:.4f}'.format(vapour_pressure_sol))
num_fuel_oil_ratio_sol = [arc_mult[a] * flow[a].get_value() /
sum(arc_mult[b] *
flow[b].get_value()
for b in ARCS if b[1] == 'fuel_oil')
for a in ARCS if a[1] == 'fuel_oil']
num_fuel_oil_ratio_sol = pd.Series(num_fuel_oil_ratio_sol,
name='num_fuel_oil_ratio_sol',
index=[a[0] for a in ARCS
if a[1] == 'fuel_oil'])
print(so.get_solution_table(fuel_oil_coefficient,
num_fuel_oil_ratio_sol))
return m.get_objective_value()
Output¶
In [1]: import os
In [2]: hostname = os.getenv('CASHOST')
In [3]: port = os.getenv('CASPORT')
In [4]: from swat import CAS
In [5]: cas_conn = CAS(hostname, port)
In [6]: import sasoptpy
In [7]: from examples.client_side.refinery_optimization import test
In [8]: test(cas_conn)
NOTE: Initialized model refinery_optimization.
NOTE: Added action set 'optimization'.
NOTE: Converting model refinery_optimization to OPTMODEL.
NOTE: Submitting OPTMODEL code to CAS server.
NOTE: Problem generation will use 8 threads.
NOTE: The problem has 51 variables (0 free, 0 fixed).
NOTE: The problem has 46 linear constraints (4 LE, 38 EQ, 3 GE, 1 range).
NOTE: The problem has 158 linear constraint coefficients.
NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).
NOTE: The OPTMODEL presolver is disabled for linear problems.
NOTE: The LP presolver value AUTOMATIC is applied.
NOTE: The LP presolver time is 0.01 seconds.
NOTE: The LP presolver removed 29 variables and 30 constraints.
NOTE: The LP presolver removed 86 constraint coefficients.
NOTE: The presolved problem has 22 variables, 16 constraints, and 72 constraint coefficients.
NOTE: The LP solver is called.
NOTE: The Dual Simplex algorithm is used.
Objective
Phase Iteration Value Time
D 2 1 9.878656E+05 0
P 2 18 2.113651E+05 0
NOTE: Optimal.
NOTE: Objective = 211365.13477.
NOTE: The Dual Simplex solve time is 0.02 seconds.
NOTE: The output table 'SOLUTION' in caslib 'CASUSER(casuser)' has 51 rows and 6 columns.
NOTE: The output table 'DUAL' in caslib 'CASUSER(casuser)' has 46 rows and 4 columns.
NOTE: The CAS table 'solutionSummary' in caslib 'CASUSER(casuser)' has 13 rows and 4 columns.
NOTE: The CAS table 'problemSummary' in caslib 'CASUSER(casuser)' has 18 rows and 4 columns.
crudesDistilled
crude1 15000.0
crude2 30000.0
oilCracked
light_oil_cracked 4200.0
heavy_oil_cracked 3800.0
flow
(source, crude1) 15000.000000
(source, crude2) 30000.000000
(crude1, light_naphtha) 15000.000000
(crude1, medium_naphtha) 15000.000000
(crude1, heavy_naphtha) 15000.000000
(crude1, light_oil) 15000.000000
(crude1, heavy_oil) 15000.000000
(crude1, residuum) 15000.000000
(crude2, light_naphtha) 30000.000000
(crude2, medium_naphtha) 30000.000000
(crude2, heavy_naphtha) 30000.000000
(crude2, light_oil) 30000.000000
(crude2, heavy_oil) 30000.000000
(crude2, residuum) 30000.000000
(light_naphtha, regular_petrol) 3293.112993
(light_naphtha, premium_petrol) 2706.887007
(medium_naphtha, regular_petrol) 10500.000000
(medium_naphtha, premium_petrol) 0.000000
(heavy_naphtha, regular_petrol) 1315.334140
(heavy_naphtha, premium_petrol) 1677.804016
(light_naphtha, reformed_gasoline) 0.000000
(medium_naphtha, reformed_gasoline) 0.000000
(heavy_naphtha, reformed_gasoline) 5406.861844
(light_oil, jet_fuel) 0.000000
(light_oil, fuel_oil) 0.000000
(heavy_oil, jet_fuel) 4900.000000
(heavy_oil, fuel_oil) 0.000000
(light_oil, light_oil_cracked) 4200.000000
(light_oil_cracked, cracked_oil) 4200.000000
(light_oil_cracked, cracked_gasoline) 4200.000000
(heavy_oil, heavy_oil_cracked) 3800.000000
(heavy_oil_cracked, cracked_oil) 3800.000000
(heavy_oil_cracked, cracked_gasoline) 3800.000000
(cracked_oil, jet_fuel) 5706.000000
(cracked_oil, fuel_oil) 0.000000
(reformed_gasoline, regular_petrol) 0.000000
(reformed_gasoline, premium_petrol) 2433.087830
(cracked_gasoline, regular_petrol) 1936.000000
(cracked_gasoline, premium_petrol) 0.000000
(residuum, lube_oil) 1000.000000
(residuum, jet_fuel) 4550.000000
(residuum, fuel_oil) 0.000000
(premium_petrol, sink) 6817.778853
(regular_petrol, sink) 17044.447133
(jet_fuel, sink) 15156.000000
(fuel_oil, sink) 0.000000
(lube_oil, sink) 500.000000
octane_sol octane_lb
petrol
regular_petrol 84.0 84
premium_petrol 94.0 94
vapour_pressure
oil
light_oil 1.0
heavy_oil 0.6
cracked_oil 1.5
residuum 0.05
Vapour_pressure_sol: 0.7737
coefficient num_fuel_oil_ratio_sol
light_oil 10 nan
cracked_oil 4 nan
heavy_oil 3 nan
residuum 1 nan
Out[8]: 211365.13476893297