{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "*This notebook contains course material from [CBE30338](https://jckantor.github.io/CBE30338)\n", "by Jeffrey Kantor (jeff at nd.edu); the content is available [on Github](https://github.com/jckantor/CBE30338.git).\n", "The text is released under the [CC-BY-NC-ND-4.0 license](https://creativecommons.org/licenses/by-nc-nd/4.0/legalcode),\n", "and code is released under the [MIT license](https://opensource.org/licenses/MIT).*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [Gasoline Blending](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/06.08-Gasoline-Blending.ipynb) | [Contents](toc.ipynb) | [Simulation and Optimal Control](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/07.00-Simulation-and-Optimal-Control.ipynb) >

\"Open

\"Download\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Pyomo Examples" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example: Nonlinear Optimization of Series Reaction in a Continuous Stirred Tank Reactor" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 Var Declarations\n", " q : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : None : 8.944271909985442 : None : False : False : Reals\n", "\n", "1 Objective Declarations\n", " objective : Size=1, Index=None, Active=True\n", " Key : Active : Sense : Expression\n", " None : True : maximize : 40.0*q/(q + 4.0)/(q + 20.0)\n", "\n", "2 Declarations: q objective\n", "\n", "Flowrate at maximum CB = 8.944271909985442 liters per minute.\n", "\n", "Maximum CB = 0.954915028125263 moles per liter.\n", "\n", "Productivity = 8.541019662483748 moles per minute.\n" ] } ], "source": [ "from pyomo.environ import *\n", "\n", "V = 40 # liters\n", "kA = 0.5 # 1/min\n", "kB = 0.1 # l/min\n", "CAf = 2.0 # moles/liter\n", "\n", "# create a model instance\n", "model = ConcreteModel()\n", "\n", "# create x and y variables in the model\n", "model.q = Var()\n", "\n", "# add a model objective\n", "model.objective = Objective(expr = model.q*V*kA*CAf/(model.q + V*kB)/(model.q + V*kA), sense=maximize)\n", "\n", "# compute a solution using ipopt for nonlinear optimization\n", "results = SolverFactory('ipopt').solve(model)\n", "model.pprint()\n", "\n", "\n", "# print solutions\n", "qmax = model.q()\n", "CBmax = model.objective()\n", "print('\\nFlowrate at maximum CB = ', qmax, 'liters per minute.')\n", "print('\\nMaximum CB =', CBmax, 'moles per liter.')\n", "print('\\nProductivity = ', qmax*CBmax, 'moles per minute.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 19.3: Linear Programming Refinery" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " Price Max Production\n", "Gasoline 72 24000\n", "Kerosine 48 2000\n", "Fuel Oil 42 6000\n", "Residual 20 100000\n", "\n", " Processing Cost Feed Costs\n", "Crude 1 1.0 48\n", "Crude 2 2.0 30\n", "\n", " Crude 1 Crude 2\n", "Gasoline 0.80 0.44\n", "Kerosine 0.05 0.10\n", "Fuel Oil 0.10 0.36\n", "Residual 0.05 0.10\n" ] } ], "source": [ "import pandas as pd\n", "\n", "PRODUCTS = ['Gasoline', 'Kerosine', 'Fuel Oil', 'Residual']\n", "FEEDS = ['Crude 1', 'Crude 2']\n", "\n", "products = pd.DataFrame(index=PRODUCTS)\n", "products['Price'] = [72, 48, 42, 20]\n", "products['Max Production'] = [24000, 2000, 6000, 100000]\n", "\n", "crudes = pd.DataFrame(index=FEEDS)\n", "crudes['Processing Cost'] = [1.00, 2.00]\n", "crudes['Feed Costs'] = [48, 30]\n", "\n", "yields = pd.DataFrame(index=PRODUCTS)\n", "yields['Crude 1'] = [0.80, 0.05, 0.10, 0.05]\n", "yields['Crude 2'] = [0.44, 0.10, 0.36, 0.10]\n", "\n", "print('\\n', products)\n", "print('\\n', crudes)\n", "print('\\n', yields)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'Gasoline': 72, 'Kerosine': 48}" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "price = {}\n", "price['Gasoline'] = 72\n", "price['Kerosine'] = 48\n", "price\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "72" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "products.loc['Gasoline','Price']" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/jeff/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:9: FutureWarning: \n", ".ix is deprecated. Please use\n", ".loc for label based indexing or\n", ".iloc for positional indexing\n", "\n", "See the documentation here:\n", "http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#ix-indexer-is-deprecated\n", " if __name__ == '__main__':\n", "/Users/jeff/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: \n", ".ix is deprecated. Please use\n", ".loc for label based indexing or\n", ".iloc for positional indexing\n", "\n", "See the documentation here:\n", "http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#ix-indexer-is-deprecated\n", " # Remove the CWD from sys.path while we load stuff.\n" ] } ], "source": [ "# model formulation\n", "model = ConcreteModel()\n", "\n", "# variables\n", "model.x = Var(FEEDS, domain=NonNegativeReals)\n", "model.y = Var(PRODUCTS, domain=NonNegativeReals)\n", "\n", "# objective\n", "income = sum(products.ix[p, 'Price'] * model.y[p] for p in PRODUCTS)\n", "raw_materials_cost = sum(crudes.ix[f,'Feed Costs'] * model.x[f] for f in FEEDS)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'processing_costs' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprocessing_cost\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprocessing_costs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mf\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mFEEDS\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mprofit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mincome\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mraw_materials_cost\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mprocessing_cost\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mobjective\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mObjective\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprofit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msense\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaximize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;31m# constraints\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m(.0)\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprocessing_cost\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprocessing_costs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mf\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mFEEDS\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mprofit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mincome\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mraw_materials_cost\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mprocessing_cost\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mobjective\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mObjective\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprofit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msense\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaximize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;31m# constraints\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'processing_costs' is not defined" ] } ], "source": [ "processing_cost = sum(processing_costs[f] * model.x[f] for f in FEEDS)\n", "profit = income - raw_materials_cost - processing_cost\n", "model.objective = Objective(expr = profit, sense=maximize)\n", "\n", "# constraints\n", "model.constraints = ConstraintList()\n", "for p in PRODUCTS:\n", " model.constraints.add(0 <= model.y[p] <= max_production[p])\n", "for p in PRODUCTS:\n", " model.constraints.add(model.y[p] == sum(model.x[f] * product_yield[(f,p)] for f in FEEDS))\n", "\n", "solver = SolverFactory('glpk')\n", "solver.solve(model)\n", "model.pprint()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "ename": "SyntaxError", "evalue": "invalid syntax (, line 17)", "output_type": "error", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m17\u001b[0m\n\u001b[0;31m product_yield =\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" ] } ], "source": [ "from pyomo.environ import *\n", "import numpy as np\n", "\n", "# problem data\n", "FEEDS = ['Crude #1', 'Crude #2']\n", "PRODUCTS = ['Gasoline', 'Kerosine', 'Fuel Oil', 'Residual']\n", "\n", "# feed costs\n", "feed_costs = {'Crude #1': 48,\n", " 'Crude #2': 30}\n", "\n", "# processing costs\n", "processing_costs = {'Crude #1': 1.00,\n", " 'Crude #2': 2.00}\n", "\n", "# yield data\n", "product_yield = \n", "product_yield = {('Crude #1', 'Gasoline'): 0.80,\n", " ('Crude #1', 'Kerosine'): 0.05,\n", " ('Crude #1', 'Fuel Oil'): 0.10,\n", " ('Crude #1', 'Residual'): 0.05,\n", " ('Crude #2', 'Gasoline'): 0.44,\n", " ('Crude #2', 'Kerosine'): 0.10,\n", " ('Crude #2', 'Fuel Oil'): 0.36,\n", " ('Crude #2', 'Residual'): 0.10}\n", "\n", "# product sales prices\n", "sales_price = {'Gasoline': 72,\n", " 'Kerosine': 48,\n", " 'Fuel Oil': 42,\n", " 'Residual': 20}\n", "\n", "# production limits\n", "max_production = {'Gasoline': 24000,\n", " 'Kerosine': 2000,\n", " 'Fuel Oil': 6000,\n", " 'Residual': 100000}\n", "\n", "# model formulation\n", "model = ConcreteModel()\n", "\n", "# variables\n", "model.x = Var(FEEDS, domain=NonNegativeReals)\n", "model.y = Var(PRODUCTS, domain=NonNegativeReals)\n", "\n", "# objective\n", "income = sum(sales_price[p] * model.y[p] for p in PRODUCTS)\n", "raw_materials_cost = sum(feed_costs[f] * model.x[f] for f in FEEDS)\n", "processing_cost = sum(processing_costs[f] * model.x[f] for f in FEEDS)\n", "\n", "profit = income - raw_materials_cost - processing_cost\n", "model.objective = Objective(expr = profit, sense=maximize)\n", "\n", "# constraints\n", "model.constraints = ConstraintList()\n", "for p in PRODUCTS:\n", " model.constraints.add(0 <= model.y[p] <= max_production[p])\n", "for p in PRODUCTS:\n", " model.constraints.add(model.y[p] == sum(model.x[f] * product_yield[(f,p)] for f in FEEDS))\n", "\n", "solver = SolverFactory('glpk')\n", "solver.solve(model)\n", "model.pprint()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "573517.2413793121" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "profit()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2078344.8275862068" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "income()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1464827.5862068948" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "raw_materials_cost()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "39999.99999999996" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "processing_cost()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example: Making Change\n", "\n", "One of the important modeling features of Pyomo is the ability to index variables and constraints. The" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'dime': 0, 'nickel': 1, 'penny': 4, 'quarter': 41}" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from pyomo.environ import *\n", "\n", "def make_change(amount, coins):\n", " model = ConcreteModel()\n", " model.x = Var(coins.keys(), domain=NonNegativeIntegers)\n", " model.total = Objective(expr = sum(model.x[c] for c in coins.keys()), sense=minimize)\n", " model.amount = Constraint(expr = sum(model.x[c]*coins[c] for c in coins.keys()) == amount)\n", " SolverFactory('glpk').solve(model)\n", " return {c: int(model.x[c]()) for c in coins.keys()} \n", " \n", "# problem data\n", "coins = {\n", " 'penny': 1,\n", " 'nickel': 5,\n", " 'dime': 10,\n", " 'quarter': 25\n", "}\n", " \n", "make_change(1034, coins)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example: Linear Production Model with Constraints with Duals" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# ==========================================================\n", "# = Solver Results =\n", "# ==========================================================\n", "# ----------------------------------------------------------\n", "# Problem Information\n", "# ----------------------------------------------------------\n", "Problem: \n", "- Name: unknown\n", " Lower bound: 2600.0\n", " Upper bound: 2600.0\n", " Number of objectives: 1\n", " Number of constraints: 4\n", " Number of variables: 3\n", " Number of nonzeros: 6\n", " Sense: maximize\n", "# ----------------------------------------------------------\n", "# Solver Information\n", "# ----------------------------------------------------------\n", "Solver: \n", "- Status: ok\n", " Termination condition: optimal\n", " Statistics: \n", " Branch and bound: \n", " Number of bounded subproblems: 0\n", " Number of created subproblems: 0\n", " Error rc: 0\n", " Time: 0.01401972770690918\n", "# ----------------------------------------------------------\n", "# Solution Information\n", "# ----------------------------------------------------------\n", "Solution: \n", "- number of solutions: 0\n", " number of solutions displayed: 0\n" ] } ], "source": [ "from pyomo.environ import *\n", "model = ConcreteModel()\n", "\n", "# for access to dual solution for constraints\n", "model.dual = Suffix(direction=Suffix.IMPORT)\n", "\n", "# declare decision variables\n", "model.x = Var(domain=NonNegativeReals)\n", "model.y = Var(domain=NonNegativeReals)\n", "\n", "# declare objective\n", "model.profit = Objective(expr = 40*model.x + 30*model.y, sense = maximize)\n", "\n", "# declare constraints\n", "model.demand = Constraint(expr = model.x <= 40)\n", "model.laborA = Constraint(expr = model.x + model.y <= 80)\n", "model.laborB = Constraint(expr = 2*model.x + model.y <= 100)\n", "\n", "# solve\n", "SolverFactory('glpk').solve(model).write()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Constraint value lslack uslack dual\n", "demand 20.00 -inf 20.00 0.00\n", "laborA 80.00 -inf 0.00 20.00\n", "laborB 100.00 -inf 0.00 10.00\n" ] } ], "source": [ "str = \" {0:7.2f} {1:7.2f} {2:7.2f} {3:7.2f}\"\n", "\n", "print(\"Constraint value lslack uslack dual\")\n", "for c in [model.demand, model.laborA, model.laborB]:\n", " print(c, str.format(c(), c.lslack(), c.uslack(), model.dual[c]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [Gasoline Blending](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/06.08-Gasoline-Blending.ipynb) | [Contents](toc.ipynb) | [Simulation and Optimal Control](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/07.00-Simulation-and-Optimal-Control.ipynb) >

\"Open

\"Download\"" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.4" } }, "nbformat": 4, "nbformat_minor": 2 }