{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Gasoline Blending\n", "\n", "** THIS WILL BE REVISED IN CLASS 3/22 (spill over to 3/24) **\n", "\n", "The task is to determine the most profitable blend of gasoline products from given set of refinery streams.\n", "\n", "![](https://www.researchgate.net/profile/Jeffrey-Kelly-2/publication/260356707/figure/fig4/AS:267835926642706@1440868469487/Typical-gasoline-blending-system-Color-figure-can-be-viewed-in-the-online-issue-which_W640.jpg)\n" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "import sys\n", "if \"google.colab\" in sys.modules:\n", " !wget \"https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py\"\n", " import helper\n", " helper.install_idaes()\n", " helper.install_ipopt()\n", " helper.install_glpk()\n", " helper.install_cbc()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import pyomo.environ as pyo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gasoline Product Specifications\n", "\n", "The gasoline products include regular and premium gasoline. In addition to the current price, the specifications include\n", "\n", "* **octane** the minimum road octane number. Road octane is the computed as the average of the Research Octane Number (RON) and Motor Octane Number (MON).\n", "* **Reid Vapor Pressure** Upper and lower limits are specified for the Reid vapor pressure. The Reid vapor pressure is the absolute pressure exerted by the liquid at 100°F.\n", "* **benzene** the maximum volume percentage of benzene allowed in the final product. Benzene helps to increase octane rating, but is also a treacherous environmental contaminant.\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
priceoctaneRVPminRVPmaxbenzene
Regular2.7587.00.015.01.1
Premium2.8591.00.015.01.1
\n", "
" ], "text/plain": [ " price octane RVPmin RVPmax benzene\n", "Regular 2.75 87.0 0.0 15.0 1.1\n", "Premium 2.85 91.0 0.0 15.0 1.1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "products = pd.DataFrame({\n", " 'Regular' : {'price': 2.75, 'octane': 87, 'RVPmin': 0.0, 'RVPmax': 15.0, 'benzene': 1.1},\n", " 'Premium' : {'price': 2.85, 'octane': 91, 'RVPmin': 0.0, 'RVPmax': 15.0, 'benzene': 1.1},\n", "}).T\n", "\n", "display(products)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stream Specifications\n", "\n", "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:\n", "\n", "* **Butane** n-butane is a C4 product stream produced from the light components of the crude being processed by the refinery. Butane is a highly volatile of gasoline.\n", "* **LSR** Light straight run naptha is a 90°F to 190°F cut from the crude distillation column primarily consisting of straight chain C5-C6 hydrocarbons.\n", "* **Isomerate** is the result of isomerizing LSR to produce branched molecules that results in higher octane number.\n", "* **Reformate** is result of catalytic reforming heavy straight run napthenes to produce a high octane blending component, as well by-product hydrogen used elsewhere in the refinery for hydro-treating.\n", "* **Reformate LB** is a is a low benzene variant of reformate.\n", "* **FCC Naphta** is the product of a fluidized catalytic cracking unit designed to produce gasoline blending components from long chain hydrocarbons present in the crude oil being processed by the refinery.\n", "* **Alkylate** The alkylation unit reacts iso-butane with low-molecular weight alkenes to produce a high octane blending component for gasoline.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
RONMONRVPbenzenecostavailoctane
Butane93.092.054.00.000.8530000.092.50
LSR78.076.011.20.732.0535001.077.00
Isomerate83.081.113.50.002.200.082.05
Reformate100.088.23.21.852.8060000.094.10
Reformate LB93.784.02.80.122.750.088.85
FCC Naphtha92.177.11.41.062.6070000.084.60
Alkylate97.395.94.60.002.7540000.096.60
\n", "
" ], "text/plain": [ " RON MON RVP benzene cost avail octane\n", "Butane 93.0 92.0 54.0 0.00 0.85 30000.0 92.50\n", "LSR 78.0 76.0 11.2 0.73 2.05 35001.0 77.00\n", "Isomerate 83.0 81.1 13.5 0.00 2.20 0.0 82.05\n", "Reformate 100.0 88.2 3.2 1.85 2.80 60000.0 94.10\n", "Reformate LB 93.7 84.0 2.8 0.12 2.75 0.0 88.85\n", "FCC Naphtha 92.1 77.1 1.4 1.06 2.60 70000.0 84.60\n", "Alkylate 97.3 95.9 4.6 0.00 2.75 40000.0 96.60" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "streams = pd.DataFrame({\n", " 'Butane' : {'RON': 93.0, 'MON': 92.0, 'RVP': 54.0, 'benzene': 0.00, 'cost': 0.85, 'avail': 30000},\n", " 'LSR' : {'RON': 78.0, 'MON': 76.0, 'RVP': 11.2, 'benzene': 0.73, 'cost': 2.05, 'avail': 35001},\n", " 'Isomerate' : {'RON': 83.0, 'MON': 81.1, 'RVP': 13.5, 'benzene': 0.00, 'cost': 2.20, 'avail': 0},\n", " 'Reformate' : {'RON':100.0, 'MON': 88.2, 'RVP': 3.2, 'benzene': 1.85, 'cost': 2.80, 'avail': 60000},\n", " 'Reformate LB' : {'RON': 93.7, 'MON': 84.0, 'RVP': 2.8, 'benzene': 0.12, 'cost': 2.75, 'avail': 0},\n", " 'FCC Naphtha' : {'RON': 92.1, 'MON': 77.1, 'RVP': 1.4, 'benzene': 1.06, 'cost': 2.60, 'avail': 70000},\n", " 'Alkylate' : {'RON': 97.3, 'MON': 95.9, 'RVP': 4.6, 'benzene': 0.00, 'cost': 2.75, 'avail': 40000},\n", "}).T\n", "\n", "streams['octane'] = (streams['RON'] + streams['MON'])/2\n", "\n", "display(streams)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Questions we want to Answer\n", "\n", "1. What is the maximum profit possible using the current product specifications and available streams?\n", "\n", "2. What are the marginal values of each blending stream? That is, how much would you be willing to pay for each additional gallon of the blending streams?\n", "\n", "3. A 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. Would an additional profit be created? What at what price point does the mid-grade product enhance profits?\n", "\n", "4. 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. What is the impact on profits?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Blending Model\n", "\n", "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$.\n", "\n", "### Objective\n", "\n", "The objective is maximize profit, which is the difference between product revenue and stream costs. \n", "\n", "\\begin{align}\n", "\\mbox{profit} & = \\max_{x_{s,p}}\\left( \\sum_{p\\in P} \\text{Price}_p \\sum_{s\\in S} x_{s,p}\n", "- \\sum_{s\\in S} \\text{Cost}_s \\sum_{p\\in P} x_{s,p}\\right) \\\\\n", "& = \\max_{x_{s,p}}\\left(\\sum_{p\\in P}\\sum_{s\\in S}x_{s,p}(\\text{Price}_p - \\text{Cost}_s)\\right)\n", "\\end{align}\n", "\n", "### Raw Materials\n", "\n", "The first constraints in any blending problem are normally the limits on available raw materials. Letting $\\text{\n", "\n", "The blending constraint for octane can be written as \n", "\n", "\\begin{align}\n", "\\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\n", "\\end{align}\n", "\n", "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\n", "\n", "\\begin{align}\n", "\\sum_{s \\in S} x_{s,p}\\left(\\mbox{Octane}_s - \\mbox{Octane}_p\\right) & \\geq 0 & \\forall p \\in P\n", "\\end{align}\n", "\n", "The same assumptions and development apply to the benzene constraint\n", "\n", "\\begin{align}\n", "\\sum_{s \\in S} x_{s,p}\\left(\\mbox{Benzene}_s - \\mbox{Benzene}_p\\right) & \\leq 0 & \\forall p \\in P\n", "\\end{align}\n", "\n", "Reid vapor pressure, however, follows a somewhat different mixing rule. For the Reid vapor pressure we have\n", "\n", "\\begin{align}\n", "\\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 \\\\\n", "\\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\n", "\\end{align}\n", "\n", "This model is implemented in the following cell." ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "profit = 100425.45002\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
RegularPremium
Butane2.175503e+048244.9672
LSR9.212971e+0325788.0290
Isomerate0.000000e+000.0000
Reformate1.978550e+0440214.5030
Reformate LB0.000000e+000.0000
FCC Naphtha7.000000e+040.0000
Alkylate-7.275958e-1240000.0000
\n", "
" ], "text/plain": [ " Regular Premium\n", "Butane 2.175503e+04 8244.9672\n", "LSR 9.212971e+03 25788.0290\n", "Isomerate 0.000000e+00 0.0000\n", "Reformate 1.978550e+04 40214.5030\n", "Reformate LB 0.000000e+00 0.0000\n", "FCC Naphtha 7.000000e+04 0.0000\n", "Alkylate -7.275958e-12 40000.0000" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import pyomo.environ as pyo\n", "\n", "def gas_blending(products, streams):\n", "\n", " m = pyo.ConcreteModel(\"Gasoline Blending\")\n", " \n", " m.PRODUCTS = pyo.Set(initialize=products.index)\n", " m.STREAMS = pyo.Set(initialize=streams.index)\n", " \n", " m.x = pyo.Var(m.STREAMS, m.PRODUCTS, domain=pyo.NonNegativeReals)\n", " \n", " @m.Objective(sense=pyo.maximize)\n", " def profit(m):\n", " return sum(sum(m.x[s, p]*(products.loc[p, 'price'] - streams.loc[s, 'cost']) for s in m.STREAMS) for p in m.PRODUCTS)\n", " \n", " @m.Constraint(m.STREAMS)\n", " def raw_material_available(m, s):\n", " return sum(m.x[s, p] for p in m.PRODUCTS) <= streams.loc[s, 'avail']\n", " \n", " @m.Constraint(m.PRODUCTS)\n", " def octane(m, p):\n", " return sum(m.x[s, p]*(streams.loc[s, 'octane'] - products.loc[p, 'octane']) for s in m.STREAMS) >= 0\n", " \n", " @m.Constraint(m.PRODUCTS)\n", " def benzene(m, p):\n", " return sum(m.x[s, p]*(streams.loc[s, 'benzene'] - products.loc[p, 'benzene']) for s in m.STREAMS) <= 0\n", " \n", " @m.Constraint(m.PRODUCTS)\n", " def min_reid_vapor_pressure(m, p):\n", " return sum(m.x[s, p]*(streams.loc[s, 'RVP']**1.25 - products.loc[p, 'RVPmin']**1.25) for s in m.STREAMS) >= 0\n", " \n", " @m.Constraint(m.PRODUCTS)\n", " def max_reid_vapor_pressure(m, p):\n", " return sum(m.x[s, p]*(streams.loc[s, 'RVP']**1.25 - products.loc[p, 'RVPmax']**1.25) for s in m.STREAMS) <= 0\n", "\n", " solver = pyo.SolverFactory('cbc')\n", " solver.solve(m)\n", " \n", " return m\n", "\n", "m = gas_blending(products, streams)\n", "\n", " \n", "soln = pd.DataFrame({s: {p: m.x[s,p]() for p in m.PRODUCTS} for s in m.STREAMS}).T\n", "print(\"profit = \", m.profit())\n", "display(soln)" ] }, { "cell_type": "code", "execution_count": 162, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total Volume = 217971.0 gallons.\n", "Total Profit = 60022.2 dollars.\n", "Profit = 27.5 cents per gallon.\n" ] } ], "source": [ "# create decision variables\n", "S = streams.keys()\n", "P = products.keys()\n", "m.x = pyomo.Var(S,P, domain=pyomo.NonNegativeReals)\n", " \n", "# objective\n", "revenue = sum(sum(m.x[s,p]*products[p]['price'] for s in S) for p in P)\n", "cost = sum(sum(m.x[s,p]*streams[s]['cost'] for s in S) for p in P)\n", "m.profit = pyomo.Objective(expr = revenue - cost, sense=pyomo.maximize)\n", "\n", "# constraints\n", "m.cons = pyomo.ConstraintList()\n", "for s in S:\n", " m.cons.add(sum(m.x[s,p] for p in P) <= streams[s]['avail'])\n", "for p in P:\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['octane'] - products[p]['octane']) for s in S) >= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['RVP']**1.25 - products[p]['RVPmin']**1.25) for s in S) >= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['RVP']**1.25 - products[p]['RVPmax']**1.25) for s in S) <= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['benzene'] - products[p]['benzene']) for s in S) <= 0)\n", "\n", "# solve\n", "solver = pyomo.SolverFactory('glpk')\n", "solver.solve(m)\n", "\n", "# display results\n", "vol = sum(m.x[s,p]() for s in S for p in P)\n", "print(\"Total Volume =\", round(vol, 1), \"gallons.\")\n", "print(\"Total Profit =\", round(m.profit(), 1), \"dollars.\")\n", "print(\"Profit =\", round(100*m.profit()/vol,1), \"cents per gallon.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Display Results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Results for each Stream" ] }, { "cell_type": "code", "execution_count": 152, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Regular Premium Total Available Unused (Slack)\n", "Butane 17924.9 12075.1 30000.0 30000.0 0.0\n", "LSR 35000.0 0.0 35000.0 35000.0 0.0\n", "Isomerate 0.0 0.0 0.0 0.0 0.0\n", "Reformate 43599.3 16400.7 60000.0 60000.0 0.0\n", "Reformate LB 0.0 0.0 0.0 0.0 0.0\n", "FCC Naphtha 24225.8 45774.2 70000.0 70000.0 0.0\n", "Alkylate 0.0 40000.0 40000.0 40000.0 0.0\n" ] } ], "source": [ "stream_results = pd.DataFrame()\n", "for s in S:\n", " for p in P:\n", " stream_results.loc[s,p] = round(m.x[s,p](), 1)\n", " stream_results.loc[s,'Total'] = round(sum(m.x[s,p]() for p in P), 1)\n", " stream_results.loc[s,'Available'] = streams[s]['avail']\n", " \n", "stream_results['Unused (Slack)'] = stream_results['Available'] - stream_results['Total']\n", "print(stream_results)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Results for each Product" ] }, { "cell_type": "code", "execution_count": 153, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Volume octane RVP benzene\n", "Regular 120750.0 87.0 15.0 1.1\n", "Premium 114250.0 91.0 10.6 0.7\n" ] } ], "source": [ "product_results = pd.DataFrame()\n", "for p in P:\n", " product_results.loc[p,'Volume'] = round(sum(m.x[s,p]() for s in S), 1)\n", " product_results.loc[p,'octane'] = round(sum(m.x[s,p]()*streams[s]['octane'] for s in S)\n", " /product_results.loc[p,'Volume'], 1)\n", " product_results.loc[p,'RVP'] = round((sum(m.x[s,p]()*streams[s]['RVP']**1.25 for s in S)\n", " /product_results.loc[p,'Volume'])**0.8, 1)\n", " product_results.loc[p,'benzene'] = round(sum(m.x[s,p]()*streams[s]['benzene'] for s in S)\n", " /product_results.loc[p,'Volume'], 1)\n", "print(product_results)" ] }, { "cell_type": "markdown", "metadata": { "tags": [ "exercise" ] }, "source": [ "## Exercise 1.\n", "\n", "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?\n", "\n", "Create a new cell (or cells) below to compute a solution to this exercise." ] }, { "cell_type": "code", "execution_count": 181, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " RVPmax RVPmin benzene octane price\n", "Midgrade 15.0 0.0 1.1 89.0 2.82\n", "Premium 15.0 0.0 1.1 91.0 2.85\n", "Regular 15.0 0.0 1.1 87.0 2.75\n", " MON RON RVP avail benzene cost octane\n", "Alkylate 95.9 97.3 4.6 40000.0 0.00 2.75 96.60\n", "Butane 92.0 93.0 54.0 30000.0 0.00 0.85 92.50\n", "FCC Naphtha 77.1 92.1 1.4 70000.0 1.06 2.60 84.60\n", "Isomerate 81.1 83.0 13.5 0.0 0.00 2.20 82.05\n", "LSR 76.0 78.0 11.2 35000.0 0.73 2.05 77.00\n", "Reformate 88.2 100.0 3.2 60000.0 1.85 2.80 94.10\n", "Reformate LB 84.0 93.7 2.8 0.0 0.12 2.75 88.85\n", "Profit = $ 104995.0\n", " Regular Midgrade Premium Total Available Unused (Slack)\n", "Butane 0.0 30000.0 0.0 30000.0 30000.0 0.0\n", "LSR 2833.3 32166.7 0.0 35000.0 35000.0 0.0\n", "Isomerate 0.0 0.0 0.0 0.0 0.0 0.0\n", "Reformate 2746.9 57253.1 0.0 60000.0 60000.0 0.0\n", "Reformate LB 0.0 0.0 0.0 0.0 0.0 0.0\n", "FCC Naphtha 0.0 70000.0 0.0 70000.0 70000.0 0.0\n", "Alkylate 919.8 39080.2 0.0 40000.0 40000.0 0.0\n", " Volume octane RVP benzene\n", "Regular 6500.0 87.0 7.2 1.1\n", "Midgrade 228500.0 89.0 13.0 0.9\n", "Premium 0.0 NaN NaN NaN\n", "Total Profit = 104995.0 dollars.\n", "Total Volume = 235000.0 gallons.\n", "Profit = 44.7 cents per gallon.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/jeff/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:68: RuntimeWarning: invalid value encountered in double_scalars\n", "/Users/jeff/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:70: RuntimeWarning: invalid value encountered in double_scalars\n", "/Users/jeff/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:72: RuntimeWarning: invalid value encountered in double_scalars\n" ] } ], "source": [ "products = {\n", " 'Regular' : {'price': 2.75, 'octane': 87, 'RVPmin': 0.0, 'RVPmax': 15.0, 'benzene': 1.1},\n", " 'Midgrade' : {'price': 2.82, 'octane': 89, 'RVPmin': 0.0, 'RVPmax': 15.0, 'benzene': 1.1},\n", " 'Premium' : {'price': 2.85, 'octane': 91, 'RVPmin': 0.0, 'RVPmax': 15.0, 'benzene': 1.1},\n", "}\n", "\n", "print(pd.DataFrame.from_dict(products).T)\n", "\n", "streams = {\n", " 'Butane' : {'RON': 93.0, 'MON': 92.0, 'RVP': 54.0, 'benzene': 0.00, 'cost': 0.85, 'avail': 30000},\n", " 'LSR' : {'RON': 78.0, 'MON': 76.0, 'RVP': 11.2, 'benzene': 0.73, 'cost': 2.05, 'avail': 35000},\n", " 'Isomerate' : {'RON': 83.0, 'MON': 81.1, 'RVP': 13.5, 'benzene': 0.00, 'cost': 2.20, 'avail': 0},\n", " 'Reformate' : {'RON':100.0, 'MON': 88.2, 'RVP': 3.2, 'benzene': 1.85, 'cost': 2.80, 'avail': 60000},\n", " 'Reformate LB' : {'RON': 93.7, 'MON': 84.0, 'RVP': 2.8, 'benzene': 0.12, 'cost': 2.75, 'avail': 0},\n", " 'FCC Naphtha' : {'RON': 92.1, 'MON': 77.1, 'RVP': 1.4, 'benzene': 1.06, 'cost': 2.60, 'avail': 70000},\n", " 'Alkylate' : {'RON': 97.3, 'MON': 95.9, 'RVP': 4.6, 'benzene': 0.00, 'cost': 2.75, 'avail': 40000},\n", "}\n", "\n", "# calculate road octane as (R+M)/2\n", "for s in streams.keys():\n", " streams[s]['octane'] = (streams[s]['RON'] + streams[s]['MON'])/2\n", " \n", "# display feed information\n", "print(pd.DataFrame.from_dict(streams).T)\n", "\n", "# create model\n", "m = pyomo.ConcreteModel()\n", "\n", "# create decision variables\n", "S = streams.keys()\n", "P = products.keys()\n", "m.x = pyomo.Var(S,P, domain=pyomo.NonNegativeReals)\n", " \n", "# objective\n", "revenue = sum(sum(m.x[s,p]*products[p]['price'] for s in S) for p in P)\n", "cost = sum(sum(m.x[s,p]*streams[s]['cost'] for s in S) for p in P)\n", "m.profit = pyomo.Objective(expr = revenue - cost, sense=pyomo.maximize)\n", "\n", "# constraints\n", "m.cons = pyomo.ConstraintList()\n", "for s in S:\n", " m.cons.add(sum(m.x[s,p] for p in P) <= streams[s]['avail'])\n", "for p in P:\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['octane'] - products[p]['octane']) for s in S) >= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['RVP']**1.25 - products[p]['RVPmin']**1.25) for s in S) >= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['RVP']**1.25 - products[p]['RVPmax']**1.25) for s in S) <= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['benzene'] - products[p]['benzene']) for s in S) <= 0)\n", "\n", "# solve\n", "solver = pyomo.SolverFactory('glpk')\n", "solver.solve(m)\n", "print(\"Profit = $\", round(m.profit(),2))\n", "\n", "stream_results = pd.DataFrame()\n", "for s in S:\n", " for p in P:\n", " stream_results.loc[s,p] = round(m.x[s,p](), 1)\n", " stream_results.loc[s,'Total'] = round(sum(m.x[s,p]() for p in P), 1)\n", " stream_results.loc[s,'Available'] = streams[s]['avail']\n", " \n", "stream_results['Unused (Slack)'] = stream_results['Available'] - stream_results['Total']\n", "print(stream_results)\n", "\n", "product_results = pd.DataFrame()\n", "for p in P:\n", " product_results.loc[p,'Volume'] = round(sum(m.x[s,p]() for s in S), 1)\n", " product_results.loc[p,'octane'] = round(sum(m.x[s,p]()*streams[s]['octane'] for s in S)\n", " /product_results.loc[p,'Volume'], 1)\n", " product_results.loc[p,'RVP'] = round((sum(m.x[s,p]()*streams[s]['RVP']**1.25 for s in S)\n", " /product_results.loc[p,'Volume'])**0.8, 1)\n", " product_results.loc[p,'benzene'] = round(sum(m.x[s,p]()*streams[s]['benzene'] for s in S)\n", " /product_results.loc[p,'Volume'], 1)\n", "print(product_results)\n", "\n", "vol = sum(m.x[s,p]() for s in S for p in P)\n", "print(\"Total Profit =\", round(m.profit(), 1), \"dollars.\")\n", "print(\"Total Volume =\", round(vol, 1), \"gallons.\")\n", "print(\"Profit =\", round(100*m.profit()/vol,1), \"cents per gallon.\")" ] }, { "cell_type": "markdown", "metadata": { "tags": [ "exercise" ] }, "source": [ "## Exercise 2.\n", "\n", "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.\n", "\n", "Compared to the base case (i.e., without the midgrade product), how does this change profitability? " ] }, { "cell_type": "code", "execution_count": 180, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " RVPmax RVPmin benzene octane price\n", "Premium 9.0 0.0 0.62 91.0 2.85\n", "Regular 9.0 0.0 0.62 87.0 2.75\n", " MON RON RVP avail benzene cost octane\n", "Alkylate 95.9 97.3 4.6 40000.0 0.00 2.75 96.60\n", "Butane 92.0 93.0 54.0 30000.0 0.00 0.85 92.50\n", "FCC Naphtha 77.1 92.1 1.4 70000.0 1.06 2.60 84.60\n", "Isomerate 81.1 83.0 13.5 0.0 0.00 2.20 82.05\n", "LSR 76.0 78.0 11.2 35000.0 0.73 2.05 77.00\n", "Reformate 88.2 100.0 3.2 60000.0 1.85 2.80 94.10\n", "Reformate LB 84.0 93.7 2.8 0.0 0.12 2.75 88.85\n", "Profit = $ 44493.62\n", " Regular Premium Total Available Unused (Slack)\n", "Butane 8187.5 0.0 8187.5 30000.0 21812.5\n", "LSR 28305.3 0.0 28305.3 35000.0 6694.7\n", "Isomerate 0.0 0.0 0.0 0.0 0.0\n", "Reformate 0.0 0.0 0.0 60000.0 60000.0\n", "Reformate LB 0.0 0.0 0.0 0.0 0.0\n", "FCC Naphtha 60824.3 0.0 60824.3 70000.0 9175.7\n", "Alkylate 40000.0 0.0 40000.0 40000.0 0.0\n", " Volume octane RVP benzene\n", "Regular 137317.1 87.0 9.0 0.6\n", "Premium 0.0 NaN NaN NaN\n", "Total Profit = 44493.6 dollars.\n", "Total Volume = 137317.1 gallons.\n", "Profit = 32.4 cents per gallon.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/jeff/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:67: RuntimeWarning: invalid value encountered in double_scalars\n", "/Users/jeff/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:69: RuntimeWarning: invalid value encountered in double_scalars\n", "/Users/jeff/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:71: RuntimeWarning: invalid value encountered in double_scalars\n" ] } ], "source": [ "products = {\n", " 'Regular' : {'price': 2.75, 'octane': 87, 'RVPmin': 0.0, 'RVPmax': 9.0, 'benzene': 0.62},\n", " 'Premium' : {'price': 2.85, 'octane': 91, 'RVPmin': 0.0, 'RVPmax': 9.0, 'benzene': 0.62},\n", "}\n", "\n", "print(pd.DataFrame.from_dict(products).T)\n", "\n", "streams = {\n", " 'Butane' : {'RON': 93.0, 'MON': 92.0, 'RVP': 54.0, 'benzene': 0.00, 'cost': 0.85, 'avail': 30000},\n", " 'LSR' : {'RON': 78.0, 'MON': 76.0, 'RVP': 11.2, 'benzene': 0.73, 'cost': 2.05, 'avail': 35000},\n", " 'Isomerate' : {'RON': 83.0, 'MON': 81.1, 'RVP': 13.5, 'benzene': 0.00, 'cost': 2.20, 'avail': 0},\n", " 'Reformate' : {'RON':100.0, 'MON': 88.2, 'RVP': 3.2, 'benzene': 1.85, 'cost': 2.80, 'avail': 60000},\n", " 'Reformate LB' : {'RON': 93.7, 'MON': 84.0, 'RVP': 2.8, 'benzene': 0.12, 'cost': 2.75, 'avail': 0},\n", " 'FCC Naphtha' : {'RON': 92.1, 'MON': 77.1, 'RVP': 1.4, 'benzene': 1.06, 'cost': 2.60, 'avail': 70000},\n", " 'Alkylate' : {'RON': 97.3, 'MON': 95.9, 'RVP': 4.6, 'benzene': 0.00, 'cost': 2.75, 'avail': 40000},\n", "}\n", "\n", "# calculate road octane as (R+M)/2\n", "for s in streams.keys():\n", " streams[s]['octane'] = (streams[s]['RON'] + streams[s]['MON'])/2\n", " \n", "# display feed information\n", "print(pd.DataFrame.from_dict(streams).T)\n", "\n", "# create model\n", "m = pyomo.ConcreteModel()\n", "\n", "# create decision variables\n", "S = streams.keys()\n", "P = products.keys()\n", "m.x = pyomo.Var(S,P, domain=pyomo.NonNegativeReals)\n", " \n", "# objective\n", "revenue = sum(sum(m.x[s,p]*products[p]['price'] for s in S) for p in P)\n", "cost = sum(sum(m.x[s,p]*streams[s]['cost'] for s in S) for p in P)\n", "m.profit = pyomo.Objective(expr = revenue - cost, sense=pyomo.maximize)\n", "\n", "# constraints\n", "m.cons = pyomo.ConstraintList()\n", "for s in S:\n", " m.cons.add(sum(m.x[s,p] for p in P) <= streams[s]['avail'])\n", "for p in P:\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['octane'] - products[p]['octane']) for s in S) >= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['RVP']**1.25 - products[p]['RVPmin']**1.25) for s in S) >= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['RVP']**1.25 - products[p]['RVPmax']**1.25) for s in S) <= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['benzene'] - products[p]['benzene']) for s in S) <= 0)\n", "\n", "# solve\n", "solver = pyomo.SolverFactory('glpk')\n", "solver.solve(m)\n", "print(\"Profit = $\", round(m.profit(),2))\n", "\n", "stream_results = pd.DataFrame()\n", "for s in S:\n", " for p in P:\n", " stream_results.loc[s,p] = round(m.x[s,p](), 1)\n", " stream_results.loc[s,'Total'] = round(sum(m.x[s,p]() for p in P), 1)\n", " stream_results.loc[s,'Available'] = streams[s]['avail']\n", " \n", "stream_results['Unused (Slack)'] = stream_results['Available'] - stream_results['Total']\n", "print(stream_results)\n", "\n", "product_results = pd.DataFrame()\n", "for p in P:\n", " product_results.loc[p,'Volume'] = round(sum(m.x[s,p]() for s in S), 1)\n", " product_results.loc[p,'octane'] = round(sum(m.x[s,p]()*streams[s]['octane'] for s in S)\n", " /product_results.loc[p,'Volume'], 1)\n", " product_results.loc[p,'RVP'] = round((sum(m.x[s,p]()*streams[s]['RVP']**1.25 for s in S)\n", " /product_results.loc[p,'Volume'])**0.8, 1)\n", " product_results.loc[p,'benzene'] = round(sum(m.x[s,p]()*streams[s]['benzene'] for s in S)\n", " /product_results.loc[p,'Volume'], 1)\n", "print(product_results)\n", "\n", "vol = sum(m.x[s,p]() for s in S for p in P)\n", "print(\"Total Profit =\", round(m.profit(), 1), \"dollars.\")\n", "print(\"Total Volume =\", round(vol, 1), \"gallons.\")\n", "print(\"Profit =\", round(100*m.profit()/vol,1), \"cents per gallon.\")" ] }, { "cell_type": "markdown", "metadata": { "tags": [ "exercise" ] }, "source": [ "## Exercise 3.\n", "\n", "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?" ] }, { "cell_type": "code", "execution_count": 184, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " RVPmax RVPmin benzene octane price\n", "Premium 9.0 0.0 0.62 91.0 2.85\n", "Regular 9.0 0.0 0.62 87.0 2.75\n", " MON RON RVP avail benzene cost octane\n", "Alkylate 95.9 97.3 4.6 40000.0 0.00 2.75 96.60\n", "Butane 92.0 93.0 54.0 30000.0 0.00 0.85 92.50\n", "FCC Naphtha 77.1 92.1 1.4 70000.0 1.06 2.60 84.60\n", "Isomerate 81.1 83.0 13.5 0.0 0.00 2.20 82.05\n", "LSR 76.0 78.0 11.2 35000.0 0.73 2.05 77.00\n", "Reformate 88.2 100.0 3.2 0.0 1.85 2.80 94.10\n", "Reformate LB 84.0 93.7 2.8 60000.0 0.12 2.75 88.85\n", "Profit = $ 63791.16\n", " Regular Premium Total Available Unused (Slack)\n", "Butane 13290.3 1122.5 14412.8 30000.0 15587.2\n", "LSR 35000.0 0.0 35000.0 35000.0 0.0\n", "Isomerate 0.0 0.0 0.0 0.0 0.0\n", "Reformate 0.0 0.0 0.0 0.0 0.0\n", "Reformate LB 60000.0 0.0 60000.0 60000.0 0.0\n", "FCC Naphtha 63818.6 6181.4 70000.0 70000.0 0.0\n", "Alkylate 33236.2 6763.8 40000.0 40000.0 0.0\n", " Volume octane RVP benzene\n", "Regular 205345.2 87.0 9.0 0.5\n", "Premium 14067.7 91.0 9.0 0.5\n", "Total Profit = 63791.2 dollars.\n", "Total Volume = 219412.8 gallons.\n", "Profit = 29.1 cents per gallon.\n" ] } ], "source": [ "products = {\n", " 'Regular' : {'price': 2.75, 'octane': 87, 'RVPmin': 0.0, 'RVPmax': 9.0, 'benzene': 0.62},\n", " 'Premium' : {'price': 2.85, 'octane': 91, 'RVPmin': 0.0, 'RVPmax': 9.0, 'benzene': 0.62},\n", "}\n", "\n", "print(pd.DataFrame.from_dict(products).T)\n", "\n", "streams = {\n", " 'Butane' : {'RON': 93.0, 'MON': 92.0, 'RVP': 54.0, 'benzene': 0.00, 'cost': 0.85, 'avail': 30000},\n", " 'LSR' : {'RON': 78.0, 'MON': 76.0, 'RVP': 11.2, 'benzene': 0.73, 'cost': 2.05, 'avail': 35000},\n", " 'Isomerate' : {'RON': 83.0, 'MON': 81.1, 'RVP': 13.5, 'benzene': 0.00, 'cost': 2.20, 'avail': 0},\n", " 'Reformate' : {'RON':100.0, 'MON': 88.2, 'RVP': 3.2, 'benzene': 1.85, 'cost': 2.80, 'avail': 0},\n", " 'Reformate LB' : {'RON': 93.7, 'MON': 84.0, 'RVP': 2.8, 'benzene': 0.12, 'cost': 2.75, 'avail': 60000},\n", " 'FCC Naphtha' : {'RON': 92.1, 'MON': 77.1, 'RVP': 1.4, 'benzene': 1.06, 'cost': 2.60, 'avail': 70000},\n", " 'Alkylate' : {'RON': 97.3, 'MON': 95.9, 'RVP': 4.6, 'benzene': 0.00, 'cost': 2.75, 'avail': 40000},\n", "}\n", "\n", "# calculate road octane as (R+M)/2\n", "for s in streams.keys():\n", " streams[s]['octane'] = (streams[s]['RON'] + streams[s]['MON'])/2\n", " \n", "# display feed information\n", "print(pd.DataFrame.from_dict(streams).T)\n", "\n", "# create model\n", "m = pyomo.ConcreteModel()\n", "\n", "# create decision variables\n", "S = streams.keys()\n", "P = products.keys()\n", "m.x = pyomo.Var(S,P, domain=pyomo.NonNegativeReals)\n", " \n", "# objective\n", "revenue = sum(sum(m.x[s,p]*products[p]['price'] for s in S) for p in P)\n", "cost = sum(sum(m.x[s,p]*streams[s]['cost'] for s in S) for p in P)\n", "m.profit = pyomo.Objective(expr = revenue - cost, sense=pyomo.maximize)\n", "\n", "# constraints\n", "m.cons = pyomo.ConstraintList()\n", "for s in S:\n", " m.cons.add(sum(m.x[s,p] for p in P) <= streams[s]['avail'])\n", "for p in P:\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['octane'] - products[p]['octane']) for s in S) >= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['RVP']**1.25 - products[p]['RVPmin']**1.25) for s in S) >= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['RVP']**1.25 - products[p]['RVPmax']**1.25) for s in S) <= 0)\n", " m.cons.add(sum(m.x[s,p]*(feeds[s]['benzene'] - products[p]['benzene']) for s in S) <= 0)\n", " \n", "# solve\n", "solver = pyomo.SolverFactory('glpk')\n", "solver.solve(m)\n", "print(\"Profit = $\", round(m.profit(),2))\n", "\n", "stream_results = pd.DataFrame()\n", "for s in S:\n", " for p in P:\n", " stream_results.loc[s,p] = round(m.x[s,p](), 1)\n", " stream_results.loc[s,'Total'] = round(sum(m.x[s,p]() for p in P), 1)\n", " stream_results.loc[s,'Available'] = streams[s]['avail']\n", " \n", "stream_results['Unused (Slack)'] = stream_results['Available'] - stream_results['Total']\n", "print(stream_results)\n", "\n", "product_results = pd.DataFrame()\n", "for p in P:\n", " product_results.loc[p,'Volume'] = round(sum(m.x[s,p]() for s in S), 1)\n", " product_results.loc[p,'octane'] = round(sum(m.x[s,p]()*streams[s]['octane'] for s in S)\n", " /product_results.loc[p,'Volume'], 1)\n", " product_results.loc[p,'RVP'] = round((sum(m.x[s,p]()*streams[s]['RVP']**1.25 for s in S)\n", " /product_results.loc[p,'Volume'])**0.8, 1)\n", " product_results.loc[p,'benzene'] = round(sum(m.x[s,p]()*streams[s]['benzene'] for s in S)\n", " /product_results.loc[p,'Volume'], 1)\n", "print(product_results)\n", "\n", "vol = sum(m.x[s,p]() for s in S for p in P)\n", "print(\"Total Profit =\", round(m.profit(), 1), \"dollars.\")\n", "print(\"Total Volume =\", round(vol, 1), \"gallons.\")\n", "print(\"Profit =\", round(100*m.profit()/vol,1), \"cents per gallon.\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Tags", "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 4 }