None
The task is to determine the most profitable blend of gasoline products from given set of refinery streams.
%matplotlib inline
import shutil
import sys
import os.path
if not shutil.which("pyomo"):
!pip install -q pyomo
assert(shutil.which("pyomo"))
if not (shutil.which("cbc") or os.path.isfile("cbc")):
if "google.colab" in sys.modules:
!apt-get install -y -qq coinor-cbc
else:
try:
!conda install -c conda-forge coincbc
except:
pass
assert(shutil.which("cbc") or os.path.isfile("cbc"))
import pandas as pd
import pyomo.environ as pyomo
|████████████████████████████████| 9.1 MB 2.6 MB/s |████████████████████████████████| 49 kB 3.0 MB/s Selecting previously unselected package coinor-libcoinutils3v5. (Reading database ... 148492 files and directories currently installed.) Preparing to unpack .../0-coinor-libcoinutils3v5_2.10.14+repack1-1_amd64.deb ... Unpacking coinor-libcoinutils3v5 (2.10.14+repack1-1) ... Selecting previously unselected package coinor-libosi1v5. Preparing to unpack .../1-coinor-libosi1v5_0.107.9+repack1-1_amd64.deb ... Unpacking coinor-libosi1v5 (0.107.9+repack1-1) ... Selecting previously unselected package coinor-libclp1. Preparing to unpack .../2-coinor-libclp1_1.16.11+repack1-1_amd64.deb ... Unpacking coinor-libclp1 (1.16.11+repack1-1) ... Selecting previously unselected package coinor-libcgl1. Preparing to unpack .../3-coinor-libcgl1_0.59.10+repack1-1_amd64.deb ... Unpacking coinor-libcgl1 (0.59.10+repack1-1) ... Selecting previously unselected package coinor-libcbc3. Preparing to unpack .../4-coinor-libcbc3_2.9.9+repack1-1_amd64.deb ... Unpacking coinor-libcbc3 (2.9.9+repack1-1) ... Selecting previously unselected package coinor-cbc. Preparing to unpack .../5-coinor-cbc_2.9.9+repack1-1_amd64.deb ... Unpacking coinor-cbc (2.9.9+repack1-1) ... Setting up coinor-libcoinutils3v5 (2.10.14+repack1-1) ... Setting up coinor-libosi1v5 (0.107.9+repack1-1) ... Setting up coinor-libclp1 (1.16.11+repack1-1) ... Setting up coinor-libcgl1 (0.59.10+repack1-1) ... Setting up coinor-libcbc3 (2.9.9+repack1-1) ... Setting up coinor-cbc (2.9.9+repack1-1) ... Processing triggers for man-db (2.8.3-2ubuntu0.1) ... Processing triggers for libc-bin (2.27-3ubuntu1.2) ... /sbin/ldconfig.real: /usr/local/lib/python3.7/dist-packages/ideep4py/lib/libmkldnn.so.0 is not a symbolic link
The gasoline products include regular and premium gasoline. In addition to the current price, the specifications include
products = {
'Regular' : {'price': 2.75, 'octane': 87, 'RVPmin': 0.0, 'RVPmax': 15.0, 'benzene': 1.1},
'Premium' : {'price': 2.85, 'octane': 91, 'RVPmin': 0.0, 'RVPmax': 15.0, 'benzene': 1.1},
}
print(pd.DataFrame.from_dict(products).T)
price octane RVPmin RVPmax benzene Regular 2.75 87.0 0.0 15.0 1.1 Premium 2.85 91.0 0.0 15.0 1.1
A typical refinery produces many intermediate streams that can be incorporated in a blended gasoline product. Here we provide data on seven streams that include:
The stream specifications include research octane and motor octane numbers for each blending component, the Reid vapor pressure, the benzene content, cost, and availability (in gallons per day). The road octane number is computed as the average of the RON and MON.
streams = {
'Butane' : {'RON': 93.0, 'MON': 92.0, 'RVP': 54.0, 'benzene': 0.00, 'cost': 0.85, 'avail': 30000},
'LSR' : {'RON': 78.0, 'MON': 76.0, 'RVP': 11.2, 'benzene': 0.73, 'cost': 2.05, 'avail': 35000},
'Isomerate' : {'RON': 83.0, 'MON': 81.1, 'RVP': 13.5, 'benzene': 0.00, 'cost': 2.20, 'avail': 0},
'Reformate' : {'RON':100.0, 'MON': 88.2, 'RVP': 3.2, 'benzene': 1.85, 'cost': 2.80, 'avail': 60000},
'Reformate LB' : {'RON': 93.7, 'MON': 84.0, 'RVP': 2.8, 'benzene': 0.12, 'cost': 2.75, 'avail': 0},
'FCC Naphtha' : {'RON': 92.1, 'MON': 77.1, 'RVP': 1.4, 'benzene': 1.06, 'cost': 2.60, 'avail': 70000},
'Alkylate' : {'RON': 97.3, 'MON': 95.9, 'RVP': 4.6, 'benzene': 0.00, 'cost': 2.75, 'avail': 40000},
}
# calculate road octane as (R+M)/2
for s in streams.keys():
streams[s]['octane'] = (streams[s]['RON'] + streams[s]['MON'])/2
# display feed information
print(pd.DataFrame.from_dict(streams).T)
RON MON RVP benzene cost avail octane Butane 93.0 92.0 54.0 0.00 0.85 30000.0 92.50 LSR 78.0 76.0 11.2 0.73 2.05 35000.0 77.00 Isomerate 83.0 81.1 13.5 0.00 2.20 0.0 82.05 Reformate 100.0 88.2 3.2 1.85 2.80 60000.0 94.10 Reformate LB 93.7 84.0 2.8 0.12 2.75 0.0 88.85 FCC Naphtha 92.1 77.1 1.4 1.06 2.60 70000.0 84.60 Alkylate 97.3 95.9 4.6 0.00 2.75 40000.0 96.60
This simplified blending model assumes the product attributes can be computed as linear volume weighted averages of the component properties. Let the decision variable $x_{s,p} \geq 0$ be the volume, in gallons, of blending component $s \in S$ used in the final product $p \in P$.
The objective is maximize profit, which is the difference between product revenue and stream costs.
\begin{align} \mbox{profit} & = \max_{x_{s,p}}\left( \sum_{p\in P} \mbox{Price}_p \sum_{s\in S} x_{s,p} - \sum_{s\in S} \mbox{Cost}_s \sum_{p\in P} x_{s,p}\right) \end{align}or \begin{align} \mbox{profit} & = \max_{x_{s,p}} \sum_{p\in P} \sum_{s\in S} x_{s,p}\mbox{Price}_p - \max_{x_{s,p}} \sum_{p\in P} \sum_{s\in S} x_{s,p}\mbox{Cost}_s \end{align}
The blending constraint for octane can be written as
\begin{align} \frac{\sum_{s \in S} x_{s,p} \mbox{Octane}_s}{\sum_{s \in S} x_{s,p}} & \geq \mbox{Octane}_p & \forall p \in P \end{align}where $\mbox{Octane}_s$ refers to the octane rating of stream $s$, whereas $\mbox{Octane}_p$ refers to the octane rating of product $p$. Multiplying through by the denominator, and consolidating terms gives
\begin{align} \sum_{s \in S} x_{s,p}\left(\mbox{Octane}_s - \mbox{Octane}_p\right) & \geq 0 & \forall p \in P \end{align}The same assumptions and development apply to the benzene constraint
\begin{align} \sum_{s \in S} x_{s,p}\left(\mbox{Benzene}_s - \mbox{Benzene}_p\right) & \leq 0 & \forall p \in P \end{align}Reid vapor pressure, however, follows a somewhat different mixing rule. For the Reid vapor pressure we have
\begin{align} \sum_{s \in S} x_{s,p}\left(\mbox{RVP}_s^{1.25} - \mbox{RVP}_{min,p}^{1.25}\right) & \geq 0 & \forall p \in P \\ \sum_{s \in S} x_{s,p}\left(\mbox{RVP}_s^{1.25} - \mbox{RVP}_{max,p}^{1.25}\right) & \leq 0 & \forall p \in P \end{align}This model is implemented in the following cell.
import pyomo.environ as pyomo
# create model
m = pyomo.ConcreteModel()
# create decision variables
S = streams.keys()
P = products.keys()
m.x = pyomo.Var(S,P, domain=pyomo.NonNegativeReals)
# objective
revenue = sum(sum(m.x[s,p]*products[p]['price'] for s in S) for p in P)
cost = sum(sum(m.x[s,p]*streams[s]['cost'] for s in S) for p in P)
m.profit = pyomo.Objective(expr = revenue - cost, sense=pyomo.maximize)
# constraints
m.cons = pyomo.ConstraintList()
for s in S:
m.cons.add(sum(m.x[s,p] for p in P) <= streams[s]['avail'])
for p in P:
m.cons.add(sum(m.x[s,p]*(streams[s]['octane'] - products[p]['octane']) for s in S) >= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['RVP']**1.25 - products[p]['RVPmin']**1.25) for s in S) >= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['RVP']**1.25 - products[p]['RVPmax']**1.25) for s in S) <= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['benzene'] - products[p]['benzene']) for s in S) <= 0)
# solve
solver = pyomo.SolverFactory('cbc')
solver.solve(m)
# display results
vol = sum(m.x[s,p]() for s in S for p in P)
print("Total Volume =", round(vol, 1), "gallons.")
print("Total Profit =", round(m.profit(), 1), "dollars.")
print("Profit =", round(100*m.profit()/vol,1), "cents per gallon.")
Total Volume = 235000.0 gallons. Total Profit = 100425.0 dollars. Profit = 42.7 cents per gallon.
stream_results = pd.DataFrame()
for s in S:
for p in P:
stream_results.loc[s,p] = round(m.x[s,p](), 1)
stream_results.loc[s,'Total'] = round(sum(m.x[s,p]() for p in P), 1)
stream_results.loc[s,'Available'] = streams[s]['avail']
stream_results['Unused (Slack)'] = stream_results['Available'] - stream_results['Total']
print(stream_results)
Regular Premium Total Available Unused (Slack) Butane 21754.6 8245.4 30000.0 30000.0 0.0 LSR 9211.6 25788.4 35000.0 35000.0 0.0 Isomerate 0.0 0.0 0.0 0.0 0.0 Reformate 19783.9 40216.1 60000.0 60000.0 0.0 Reformate LB 0.0 0.0 0.0 0.0 0.0 FCC Naphtha 70000.0 0.0 70000.0 70000.0 0.0 Alkylate -0.0 40000.0 40000.0 40000.0 0.0
product_results = pd.DataFrame()
for p in P:
product_results.loc[p,'Volume'] = round(sum(m.x[s,p]() for s in S), 1)
product_results.loc[p,'octane'] = round(sum(m.x[s,p]()*streams[s]['octane'] for s in S)
/product_results.loc[p,'Volume'], 1)
product_results.loc[p,'RVP'] = round((sum(m.x[s,p]()*streams[s]['RVP']**1.25 for s in S)
/product_results.loc[p,'Volume'])**0.8, 1)
product_results.loc[p,'benzene'] = round(sum(m.x[s,p]()*streams[s]['benzene'] for s in S)
/product_results.loc[p,'Volume'], 1)
print(product_results)
Volume octane RVP benzene Regular 120750.0 87.0 15.0 1.0 Premium 114250.0 91.0 10.6 0.8
The marketing team says there is an opportunity to create a mid-grade gasoline product with a road octane number of 89 that would sell for $2.82/gallon, and with all other specifications the same. Could an additional profit be created?
Create a new cell (or cells) below to compute a solution to this exercise.
products = {
'Regular' : {'price': 2.75, 'octane': 87, 'RVPmin': 0.0, 'RVPmax': 15.0, 'benzene': 1.1},
'Midgrade' : {'price': 2.82, 'octane': 89, 'RVPmin': 0.0, 'RVPmax': 15.0, 'benzene': 1.1},
'Premium' : {'price': 2.85, 'octane': 91, 'RVPmin': 0.0, 'RVPmax': 15.0, 'benzene': 1.1},
}
print(pd.DataFrame.from_dict(products).T)
streams = {
'Butane' : {'RON': 93.0, 'MON': 92.0, 'RVP': 54.0, 'benzene': 0.00, 'cost': 0.85, 'avail': 30000},
'LSR' : {'RON': 78.0, 'MON': 76.0, 'RVP': 11.2, 'benzene': 0.73, 'cost': 2.05, 'avail': 35000},
'Isomerate' : {'RON': 83.0, 'MON': 81.1, 'RVP': 13.5, 'benzene': 0.00, 'cost': 2.20, 'avail': 0},
'Reformate' : {'RON':100.0, 'MON': 88.2, 'RVP': 3.2, 'benzene': 1.85, 'cost': 2.80, 'avail': 60000},
'Reformate LB' : {'RON': 93.7, 'MON': 84.0, 'RVP': 2.8, 'benzene': 0.12, 'cost': 2.75, 'avail': 0},
'FCC Naphtha' : {'RON': 92.1, 'MON': 77.1, 'RVP': 1.4, 'benzene': 1.06, 'cost': 2.60, 'avail': 70000},
'Alkylate' : {'RON': 97.3, 'MON': 95.9, 'RVP': 4.6, 'benzene': 0.00, 'cost': 2.75, 'avail': 40000},
}
# calculate road octane as (R+M)/2
for s in streams.keys():
streams[s]['octane'] = (streams[s]['RON'] + streams[s]['MON'])/2
# display feed information
print(pd.DataFrame.from_dict(streams).T)
# create model
m = pyomo.ConcreteModel()
# create decision variables
S = streams.keys()
P = products.keys()
m.x = pyomo.Var(S,P, domain=pyomo.NonNegativeReals)
# objective
revenue = sum(sum(m.x[s,p]*products[p]['price'] for s in S) for p in P)
cost = sum(sum(m.x[s,p]*streams[s]['cost'] for s in S) for p in P)
m.profit = pyomo.Objective(expr = revenue - cost, sense=pyomo.maximize)
# constraints
m.cons = pyomo.ConstraintList()
for s in S:
m.cons.add(sum(m.x[s,p] for p in P) <= streams[s]['avail'])
for p in P:
m.cons.add(sum(m.x[s,p]*(streams[s]['octane'] - products[p]['octane']) for s in S) >= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['RVP']**1.25 - products[p]['RVPmin']**1.25) for s in S) >= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['RVP']**1.25 - products[p]['RVPmax']**1.25) for s in S) <= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['benzene'] - products[p]['benzene']) for s in S) <= 0)
# solve
solver = pyomo.SolverFactory('cbc')
solver.solve(m)
print("Profit = $", round(m.profit(),2))
stream_results = pd.DataFrame()
for s in S:
for p in P:
stream_results.loc[s,p] = round(m.x[s,p](), 1)
stream_results.loc[s,'Total'] = round(sum(m.x[s,p]() for p in P), 1)
stream_results.loc[s,'Available'] = streams[s]['avail']
stream_results['Unused (Slack)'] = stream_results['Available'] - stream_results['Total']
print(stream_results)
product_results = pd.DataFrame()
for p in P:
product_results.loc[p,'Volume'] = round(sum(m.x[s,p]() for s in S), 1)
product_results.loc[p,'octane'] = round(sum(m.x[s,p]()*streams[s]['octane'] for s in S)
/product_results.loc[p,'Volume'], 1)
product_results.loc[p,'RVP'] = round((sum(m.x[s,p]()*streams[s]['RVP']**1.25 for s in S)
/product_results.loc[p,'Volume'])**0.8, 1)
product_results.loc[p,'benzene'] = round(sum(m.x[s,p]()*streams[s]['benzene'] for s in S)
/product_results.loc[p,'Volume'], 1)
print(product_results)
vol = sum(m.x[s,p]() for s in S for p in P)
print("Total Profit =", round(m.profit(), 1), "dollars.")
print("Total Volume =", round(vol, 1), "gallons.")
print("Profit =", round(100*m.profit()/vol,1), "cents per gallon.")
price octane RVPmin RVPmax benzene Regular 2.75 87.0 0.0 15.0 1.1 Midgrade 2.82 89.0 0.0 15.0 1.1 Premium 2.85 91.0 0.0 15.0 1.1 RON MON RVP benzene cost avail octane Butane 93.0 92.0 54.0 0.00 0.85 30000.0 92.50 LSR 78.0 76.0 11.2 0.73 2.05 35000.0 77.00 Isomerate 83.0 81.1 13.5 0.00 2.20 0.0 82.05 Reformate 100.0 88.2 3.2 1.85 2.80 60000.0 94.10 Reformate LB 93.7 84.0 2.8 0.12 2.75 0.0 88.85 FCC Naphtha 92.1 77.1 1.4 1.06 2.60 70000.0 84.60 Alkylate 97.3 95.9 4.6 0.00 2.75 40000.0 96.60 Profit = $ 104995.0 Regular Midgrade Premium Total Available Unused (Slack) Butane 849.5 29150.5 0.0 30000.0 30000.0 0.0 LSR 2646.2 32353.8 0.0 35000.0 35000.0 0.0 Isomerate 0.0 0.0 0.0 0.0 0.0 0.0 Reformate 2820.7 57179.3 0.0 60000.0 60000.0 0.0 Reformate LB 0.0 0.0 0.0 0.0 0.0 0.0 FCC Naphtha 0.0 70000.0 0.0 70000.0 70000.0 0.0 Alkylate 183.6 39816.4 0.0 40000.0 40000.0 0.0 Volume octane RVP benzene Regular 6500.0 87.0 15.0 1.1 Midgrade 228500.0 89.0 12.8 0.9 Premium 0.0 NaN NaN NaN Total Profit = 104995.0 dollars. Total Volume = 235000.0 gallons. Profit = 44.7 cents per gallon.
/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:68: RuntimeWarning: invalid value encountered in double_scalars /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:70: RuntimeWarning: invalid value encountered in double_scalars /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:72: RuntimeWarning: invalid value encountered in double_scalars
New environmental regulations have reduced the allowable benzene levels from 1.1 vol% to 0.62 vol%, and the maximum Reid vapor pressure from 15.0 to 9.0.
Compared to the base case (i.e., without the midgrade product), how does this change profitability?
products = {
'Regular' : {'price': 2.75, 'octane': 87, 'RVPmin': 0.0, 'RVPmax': 9.0, 'benzene': 0.62},
'Premium' : {'price': 2.85, 'octane': 91, 'RVPmin': 0.0, 'RVPmax': 9.0, 'benzene': 0.62},
}
print(pd.DataFrame.from_dict(products).T)
streams = {
'Butane' : {'RON': 93.0, 'MON': 92.0, 'RVP': 54.0, 'benzene': 0.00, 'cost': 0.85, 'avail': 30000},
'LSR' : {'RON': 78.0, 'MON': 76.0, 'RVP': 11.2, 'benzene': 0.73, 'cost': 2.05, 'avail': 35000},
'Isomerate' : {'RON': 83.0, 'MON': 81.1, 'RVP': 13.5, 'benzene': 0.00, 'cost': 2.20, 'avail': 0},
'Reformate' : {'RON':100.0, 'MON': 88.2, 'RVP': 3.2, 'benzene': 1.85, 'cost': 2.80, 'avail': 60000},
'Reformate LB' : {'RON': 93.7, 'MON': 84.0, 'RVP': 2.8, 'benzene': 0.12, 'cost': 2.75, 'avail': 0},
'FCC Naphtha' : {'RON': 92.1, 'MON': 77.1, 'RVP': 1.4, 'benzene': 1.06, 'cost': 2.60, 'avail': 70000},
'Alkylate' : {'RON': 97.3, 'MON': 95.9, 'RVP': 4.6, 'benzene': 0.00, 'cost': 2.75, 'avail': 40000},
}
# calculate road octane as (R+M)/2
for s in streams.keys():
streams[s]['octane'] = (streams[s]['RON'] + streams[s]['MON'])/2
# display feed information
print(pd.DataFrame.from_dict(streams).T)
# create model
m = pyomo.ConcreteModel()
# create decision variables
S = streams.keys()
P = products.keys()
m.x = pyomo.Var(S,P, domain=pyomo.NonNegativeReals)
# objective
revenue = sum(sum(m.x[s,p]*products[p]['price'] for s in S) for p in P)
cost = sum(sum(m.x[s,p]*streams[s]['cost'] for s in S) for p in P)
m.profit = pyomo.Objective(expr = revenue - cost, sense=pyomo.maximize)
# constraints
m.cons = pyomo.ConstraintList()
for s in S:
m.cons.add(sum(m.x[s,p] for p in P) <= streams[s]['avail'])
for p in P:
m.cons.add(sum(m.x[s,p]*(streams[s]['octane'] - products[p]['octane']) for s in S) >= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['RVP']**1.25 - products[p]['RVPmin']**1.25) for s in S) >= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['RVP']**1.25 - products[p]['RVPmax']**1.25) for s in S) <= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['benzene'] - products[p]['benzene']) for s in S) <= 0)
# solve
solver = pyomo.SolverFactory('cbc')
solver.solve(m)
print("Profit = $", round(m.profit(),2))
stream_results = pd.DataFrame()
for s in S:
for p in P:
stream_results.loc[s,p] = round(m.x[s,p](), 1)
stream_results.loc[s,'Total'] = round(sum(m.x[s,p]() for p in P), 1)
stream_results.loc[s,'Available'] = streams[s]['avail']
stream_results['Unused (Slack)'] = stream_results['Available'] - stream_results['Total']
print(stream_results)
product_results = pd.DataFrame()
for p in P:
product_results.loc[p,'Volume'] = round(sum(m.x[s,p]() for s in S), 1)
product_results.loc[p,'octane'] = round(sum(m.x[s,p]()*streams[s]['octane'] for s in S)
/product_results.loc[p,'Volume'], 1)
product_results.loc[p,'RVP'] = round((sum(m.x[s,p]()*streams[s]['RVP']**1.25 for s in S)
/product_results.loc[p,'Volume'])**0.8, 1)
product_results.loc[p,'benzene'] = round(sum(m.x[s,p]()*streams[s]['benzene'] for s in S)
/product_results.loc[p,'Volume'], 1)
print(product_results)
vol = sum(m.x[s,p]() for s in S for p in P)
print("Total Profit =", round(m.profit(), 1), "dollars.")
print("Total Volume =", round(vol, 1), "gallons.")
print("Profit =", round(100*m.profit()/vol,1), "cents per gallon.")
price octane RVPmin RVPmax benzene Regular 2.75 87.0 0.0 9.0 0.62 Premium 2.85 91.0 0.0 9.0 0.62 RON MON RVP benzene cost avail octane Butane 93.0 92.0 54.0 0.00 0.85 30000.0 92.50 LSR 78.0 76.0 11.2 0.73 2.05 35000.0 77.00 Isomerate 83.0 81.1 13.5 0.00 2.20 0.0 82.05 Reformate 100.0 88.2 3.2 1.85 2.80 60000.0 94.10 Reformate LB 93.7 84.0 2.8 0.12 2.75 0.0 88.85 FCC Naphtha 92.1 77.1 1.4 1.06 2.60 70000.0 84.60 Alkylate 97.3 95.9 4.6 0.00 2.75 40000.0 96.60 Profit = $ 44493.62 Regular Premium Total Available Unused (Slack) Butane 8187.5 0.0 8187.5 30000.0 21812.5 LSR 28305.3 0.0 28305.3 35000.0 6694.7 Isomerate 0.0 0.0 0.0 0.0 0.0 Reformate 0.0 0.0 0.0 60000.0 60000.0 Reformate LB 0.0 0.0 0.0 0.0 0.0 FCC Naphtha 60824.3 0.0 60824.3 70000.0 9175.7 Alkylate 40000.0 0.0 40000.0 40000.0 0.0 Volume octane RVP benzene Regular 137317.1 87.0 9.0 0.6 Premium 0.0 NaN NaN NaN Total Profit = 44493.6 dollars. Total Volume = 137317.1 gallons. Profit = 32.4 cents per gallon.
/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:67: RuntimeWarning: invalid value encountered in double_scalars /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:69: RuntimeWarning: invalid value encountered in double_scalars /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:71: RuntimeWarning: invalid value encountered in double_scalars
Given the new product specifications in Exercise 2, let's consider using different refinery streams. In place of Reformate, the refinery could produce Reformate LB. (That is, one or the other of the two streams could be 60000 gallons per day, but not both). Same for LSR and Reformate. How should the refinery be operated to maximize profitability?
products = {
'Regular' : {'price': 2.75, 'octane': 87, 'RVPmin': 0.0, 'RVPmax': 9.0, 'benzene': 0.62},
'Premium' : {'price': 2.85, 'octane': 91, 'RVPmin': 0.0, 'RVPmax': 9.0, 'benzene': 0.62},
}
print(pd.DataFrame.from_dict(products).T)
streams = {
'Butane' : {'RON': 93.0, 'MON': 92.0, 'RVP': 54.0, 'benzene': 0.00, 'cost': 0.85, 'avail': 30000},
'LSR' : {'RON': 78.0, 'MON': 76.0, 'RVP': 11.2, 'benzene': 0.73, 'cost': 2.05, 'avail': 35000},
'Isomerate' : {'RON': 83.0, 'MON': 81.1, 'RVP': 13.5, 'benzene': 0.00, 'cost': 2.20, 'avail': 0},
'Reformate' : {'RON':100.0, 'MON': 88.2, 'RVP': 3.2, 'benzene': 1.85, 'cost': 2.80, 'avail': 0},
'Reformate LB' : {'RON': 93.7, 'MON': 84.0, 'RVP': 2.8, 'benzene': 0.12, 'cost': 2.75, 'avail': 60000},
'FCC Naphtha' : {'RON': 92.1, 'MON': 77.1, 'RVP': 1.4, 'benzene': 1.06, 'cost': 2.60, 'avail': 70000},
'Alkylate' : {'RON': 97.3, 'MON': 95.9, 'RVP': 4.6, 'benzene': 0.00, 'cost': 2.75, 'avail': 40000},
}
# calculate road octane as (R+M)/2
for s in streams.keys():
streams[s]['octane'] = (streams[s]['RON'] + streams[s]['MON'])/2
# display feed information
print(pd.DataFrame.from_dict(streams).T)
# create model
m = pyomo.ConcreteModel()
# create decision variables
S = streams.keys()
P = products.keys()
m.x = pyomo.Var(S,P, domain=pyomo.NonNegativeReals)
# objective
revenue = sum(sum(m.x[s,p]*products[p]['price'] for s in S) for p in P)
cost = sum(sum(m.x[s,p]*streams[s]['cost'] for s in S) for p in P)
m.profit = pyomo.Objective(expr = revenue - cost, sense=pyomo.maximize)
# constraints
m.cons = pyomo.ConstraintList()
for s in S:
m.cons.add(sum(m.x[s,p] for p in P) <= streams[s]['avail'])
for p in P:
m.cons.add(sum(m.x[s,p]*(streams[s]['octane'] - products[p]['octane']) for s in S) >= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['RVP']**1.25 - products[p]['RVPmin']**1.25) for s in S) >= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['RVP']**1.25 - products[p]['RVPmax']**1.25) for s in S) <= 0)
m.cons.add(sum(m.x[s,p]*(streams[s]['benzene'] - products[p]['benzene']) for s in S) <= 0)
# solve
solver = pyomo.SolverFactory('cbc')
solver.solve(m)
print("Profit = $", round(m.profit(),2))
stream_results = pd.DataFrame()
for s in S:
for p in P:
stream_results.loc[s,p] = round(m.x[s,p](), 1)
stream_results.loc[s,'Total'] = round(sum(m.x[s,p]() for p in P), 1)
stream_results.loc[s,'Available'] = streams[s]['avail']
stream_results['Unused (Slack)'] = stream_results['Available'] - stream_results['Total']
print(stream_results)
product_results = pd.DataFrame()
for p in P:
product_results.loc[p,'Volume'] = round(sum(m.x[s,p]() for s in S), 1)
product_results.loc[p,'octane'] = round(sum(m.x[s,p]()*streams[s]['octane'] for s in S)
/product_results.loc[p,'Volume'], 1)
product_results.loc[p,'RVP'] = round((sum(m.x[s,p]()*streams[s]['RVP']**1.25 for s in S)
/product_results.loc[p,'Volume'])**0.8, 1)
product_results.loc[p,'benzene'] = round(sum(m.x[s,p]()*streams[s]['benzene'] for s in S)
/product_results.loc[p,'Volume'], 1)
print(product_results)
vol = sum(m.x[s,p]() for s in S for p in P)
print("Total Profit =", round(m.profit(), 1), "dollars.")
print("Total Volume =", round(vol, 1), "gallons.")
print("Profit =", round(100*m.profit()/vol,1), "cents per gallon.")
price octane RVPmin RVPmax benzene Regular 2.75 87.0 0.0 9.0 0.62 Premium 2.85 91.0 0.0 9.0 0.62 RON MON RVP benzene cost avail octane Butane 93.0 92.0 54.0 0.00 0.85 30000.0 92.50 LSR 78.0 76.0 11.2 0.73 2.05 35000.0 77.00 Isomerate 83.0 81.1 13.5 0.00 2.20 0.0 82.05 Reformate 100.0 88.2 3.2 1.85 2.80 0.0 94.10 Reformate LB 93.7 84.0 2.8 0.12 2.75 60000.0 88.85 FCC Naphtha 92.1 77.1 1.4 1.06 2.60 70000.0 84.60 Alkylate 97.3 95.9 4.6 0.00 2.75 40000.0 96.60 Profit = $ 63791.16 Regular Premium Total Available Unused (Slack) Butane 13906.5 506.3 14412.8 30000.0 15587.2 LSR 31086.6 3913.4 35000.0 35000.0 0.0 Isomerate 0.0 0.0 0.0 0.0 0.0 Reformate 0.0 0.0 0.0 0.0 0.0 Reformate LB 60000.0 0.0 60000.0 60000.0 0.0 FCC Naphtha 70000.0 0.0 70000.0 70000.0 0.0 Alkylate 30352.1 9647.9 40000.0 40000.0 0.0 Volume octane RVP benzene Regular 205345.2 87.0 9.0 0.5 Premium 14067.7 91.0 9.0 0.2 Total Profit = 63791.2 dollars. Total Volume = 219412.8 gallons. Profit = 29.1 cents per gallon.