{ "cells": [ { "cell_type": "markdown", "id": "07234160", "metadata": {}, "source": [ "\n", "*This notebook contains material from [CBE60499](https://ndcbe.github.io/CBE60499);\n", "content is available [on Github](git@github.com:ndcbe/CBE60499.git).*\n" ] }, { "cell_type": "markdown", "id": "f2c995ee", "metadata": {}, "source": [ "\n", "< [2.1 Continuous Optimization](https://ndcbe.github.io/CBE60499/02.01-LP-NLP.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [2.3 Logical Modeling and Generalized Disjunctive Programs](https://ndcbe.github.io/CBE60499/02.03-GDP.html) >

\"Open

\"Download\"" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 1, "link": "[2.2 Integer Programs](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2-Integer-Programs)", "section": "2.2 Integer Programs" } }, "source": [ "# 2.2 Integer Programs" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 1, "link": "[2.2 Integer Programs](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2-Integer-Programs)", "section": "2.2 Integer Programs" } }, "outputs": [], "source": [ "# This code cell installs packages on Colab\n", "\n", "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()\n", " helper.download_figures(['feasible.png'])" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[2.2.1 Optimizing Across Process Alternatives](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1-Optimizing-Across-Process-Alternatives)", "section": "2.2.1 Optimizing Across Process Alternatives" } }, "source": [ "## 2.2.1 Optimizing Across Process Alternatives" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[2.2.1 Optimizing Across Process Alternatives](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1-Optimizing-Across-Process-Alternatives)", "section": "2.2.1 Optimizing Across Process Alternatives" } }, "source": [ "Reference: Example 15.3 from Biegler, Grossmann, Westerberg (1997)." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[2.2.1 Optimizing Across Process Alternatives](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1-Optimizing-Across-Process-Alternatives)", "section": "2.2.1 Optimizing Across Process Alternatives" } }, "source": [ "Assume that we have the choice of selecting two reactors (shown below) for the reaction $A \\rightarrow B$. Reactor I has a higher conversation (80%) but it is more expensive; reactor II has a lower conversion (66.7%) but is cheaper. The cost of feed $A$ is \\$5/kmol. Which process alternative (reactor I, reactor II, or both) has the minimum costs to make 10 kmol/hr of product B?" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[2.2.1 Optimizing Across Process Alternatives](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1-Optimizing-Across-Process-Alternatives)", "section": "2.2.1 Optimizing Across Process Alternatives" } }, "source": [ "Let $x_i$ be the size of reactor $i$.\n", "\n", "Continuous cost model: $c_i (x_i)^{0.6}$" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.1 Develop the optimization model](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.1-Develop-the-optimization-model)", "section": "2.2.1.1 Develop the optimization model" } }, "source": [ "### 2.2.1.1 Develop the optimization model\n", "* Draw a picture\n", "* Sets\n", "* Parameters\n", "* Variables\n", "* Objective\n", "* Constraints\n", "* Degree of freedom analysis" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.2 Solve with Continuous Cost Model in Pyomo](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.2-Solve-with-Continuous-Cost-Model-in-Pyomo)", "section": "2.2.1.2 Solve with Continuous Cost Model in Pyomo" } }, "source": [ "### 2.2.1.2 Solve with Continuous Cost Model in Pyomo" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.2 Solve with Continuous Cost Model in Pyomo](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.2-Solve-with-Continuous-Cost-Model-in-Pyomo)", "section": "2.2.1.2 Solve with Continuous Cost Model in Pyomo" } }, "source": [ "We start my defining the model in Pyomo." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.2 Solve with Continuous Cost Model in Pyomo](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.2-Solve-with-Continuous-Cost-Model-in-Pyomo)", "section": "2.2.1.2 Solve with Continuous Cost Model in Pyomo" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 Set Declarations\n", " REACTORS : Size=1, Index=None, Ordered=Insertion\n", " Key : Dimen : Domain : Size : Members\n", " None : 1 : Any : 2 : {1, 2}\n", "\n", "4 Param Declarations\n", " conversion : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " 1 : 0.8\n", " 2 : 0.67\n", " feed_cost : Size=1, Index=None, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " None : 5.0\n", " product_flowrate : Size=1, Index=None, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " None : 10.0\n", " reactor_cost : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " 1 : 5.5\n", " 2 : 4.0\n", "\n", "3 Var Declarations\n", " feed_flowrate : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : 0 : None : None : False : True : NonNegativeReals\n", " reactor_effluent : Size=2, Index=REACTORS\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 1 : 0 : None : None : False : True : NonNegativeReals\n", " 2 : 0 : None : None : False : True : NonNegativeReals\n", " reactor_feed : Size=2, Index=REACTORS\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 1 : 0 : None : None : False : True : NonNegativeReals\n", " 2 : 0 : None : None : False : True : NonNegativeReals\n", "\n", "1 Objective Declarations\n", " cost : Size=1, Index=None, Active=True\n", " Key : Active : Sense : Expression\n", " None : True : minimize : 5.5*reactor_feed[1]**0.6 + 4.0*reactor_feed[2]**0.6 + 5.0*feed_flowrate\n", "\n", "3 Constraint Declarations\n", " inlet_split : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 0.0 : feed_flowrate - (reactor_feed[1] + reactor_feed[2]) : 0.0 : True\n", " mixer : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 10.0 : reactor_effluent[1] + reactor_effluent[2] : 10.0 : True\n", " reactor_performance : Size=2, Index=REACTORS, Active=True\n", " Key : Lower : Body : Upper : Active\n", " 1 : 0.0 : reactor_effluent[1] - 0.8*reactor_feed[1] : 0.0 : True\n", " 2 : 0.0 : reactor_effluent[2] - 0.67*reactor_feed[2] : 0.0 : True\n", "\n", "12 Declarations: REACTORS reactor_cost product_flowrate conversion feed_cost feed_flowrate reactor_feed reactor_effluent inlet_split reactor_performance mixer cost\n" ] } ], "source": [ "import pyomo.environ as pyo\n", "\n", "nlp = pyo.ConcreteModel()\n", "\n", "## Define sets\n", "nlp.REACTORS = pyo.Set(initialize=range(1,3))\n", "\n", "## Define parameters (data)\n", "\n", "# $ / hr\n", "cost_coefficient = {1:5.5, 2:4.0}\n", "nlp.reactor_cost = pyo.Param(nlp.REACTORS, initialize=cost_coefficient)\n", "\n", "# kmol/hr B\n", "nlp.product_flowrate = pyo.Param(initialize=10.0)\n", "\n", "# conversion fraction\n", "reactor_conversion = {1:0.8, 2:0.67}\n", "nlp.conversion = pyo.Param(nlp.REACTORS, initialize=reactor_conversion)\n", "\n", "# feed cost, $/kmol of A\n", "nlp.feed_cost = pyo.Param(initialize=5.0)\n", "\n", "\n", "## Define variables\n", "\n", "# Feed flowrate into reactor, x0 in handout illustration\n", "nlp.feed_flowrate = pyo.Var(domain=pyo.NonNegativeReals)\n", "\n", "# Reactor feed, x1 and x2 in handout illustration\n", "nlp.reactor_feed = pyo.Var(nlp.REACTORS, domain=pyo.NonNegativeReals)\n", "\n", "# Reactor effluent (outlet), z1 and z2 in handout illustration\n", "nlp.reactor_effluent = pyo.Var(nlp.REACTORS, domain=pyo.NonNegativeReals)\n", "\n", "## Define constraints\n", "\n", "# YOUR SOLUTION HERE\n", "\n", "## Define objective\n", "nlp.cost = pyo.Objective(expr=sum(nlp.reactor_cost[r] * (nlp.reactor_feed[r])**(0.6) for r in nlp.REACTORS) +\n", " nlp.feed_cost * nlp.feed_flowrate)\n", "\n", "nlp.pprint()" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.3 Initialize to Favor Reaction 1 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.3-Initialize-to-Favor-Reaction-1-and-Solve)", "section": "2.2.1.3 Initialize to Favor Reaction 1 and Solve" } }, "source": [ "### 2.2.1.3 Initialize to Favor Reaction 1 and Solve" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.3 Initialize to Favor Reaction 1 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.3-Initialize-to-Favor-Reaction-1-and-Solve)", "section": "2.2.1.3 Initialize to Favor Reaction 1 and Solve" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 Set Declarations\n", " REACTORS : Size=1, Index=None, Ordered=Insertion\n", " Key : Dimen : Domain : Size : Members\n", " None : 1 : Any : 2 : {1, 2}\n", "\n", "4 Param Declarations\n", " conversion : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " 1 : 0.8\n", " 2 : 0.67\n", " feed_cost : Size=1, Index=None, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " None : 5.0\n", " product_flowrate : Size=1, Index=None, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " None : 10.0\n", " reactor_cost : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " 1 : 5.5\n", " 2 : 4.0\n", "\n", "3 Var Declarations\n", " feed_flowrate : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : 0 : 20.0 : None : False : False : NonNegativeReals\n", " reactor_effluent : Size=2, Index=REACTORS\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 1 : 0 : 16.0 : None : False : False : NonNegativeReals\n", " 2 : 0 : 0.0 : None : False : False : NonNegativeReals\n", " reactor_feed : Size=2, Index=REACTORS\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 1 : 0 : 20.0 : None : False : False : NonNegativeReals\n", " 2 : 0 : 0 : None : False : False : NonNegativeReals\n", "\n", "1 Objective Declarations\n", " cost : Size=1, Index=None, Active=True\n", " Key : Active : Sense : Expression\n", " None : True : minimize : 5.5*reactor_feed[1]**0.6 + 4.0*reactor_feed[2]**0.6 + 5.0*feed_flowrate\n", "\n", "3 Constraint Declarations\n", " inlet_split : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 0.0 : feed_flowrate - (reactor_feed[1] + reactor_feed[2]) : 0.0 : True\n", " mixer : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 10.0 : reactor_effluent[1] + reactor_effluent[2] : 10.0 : True\n", " reactor_performance : Size=2, Index=REACTORS, Active=True\n", " Key : Lower : Body : Upper : Active\n", " 1 : 0.0 : reactor_effluent[1] - 0.8*reactor_feed[1] : 0.0 : True\n", " 2 : 0.0 : reactor_effluent[2] - 0.67*reactor_feed[2] : 0.0 : True\n", "\n", "12 Declarations: REACTORS reactor_cost product_flowrate conversion feed_cost feed_flowrate reactor_feed reactor_effluent inlet_split reactor_performance mixer cost\n" ] } ], "source": [ "def initialize(model, reactor_choice=1):\n", " ''' Initialize all of the variables in the model to demonstrate local solutions\n", " \n", " Arguments:\n", " model: Pyomo model\n", " reactor_choice: 1 or 2\n", " \n", " Returns:\n", " nothing\n", " \n", " Action:\n", " initializes model\n", " \n", " '''\n", " \n", " # Guess 20 kmol/hr feed of A\n", " model.feed_flowrate = 20.0\n", " \n", " # Either assign all of the feed to reactor 1 or 2\n", " if reactor_choice == 1:\n", " model.reactor_feed[1] = 20.0\n", " model.reactor_feed[2] = 0\n", " elif reactor_choice == 2:\n", " model.reactor_feed[1] = 0\n", " model.reactor_feed[2] = 20.0\n", " else:\n", " raise ValueError(\"Argument reactor_choice needs value 1 or 2.\")\n", " \n", " # Based on the feed assignments, calculate effluent flowrate\n", " for r in model.REACTORS:\n", " model.reactor_effluent[r] = model.reactor_feed[r]() * model.conversion[r]\n", " \n", "initialize(nlp, reactor_choice=1)\n", "nlp.pprint()" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.3 Initialize to Favor Reaction 1 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.3-Initialize-to-Favor-Reaction-1-and-Solve)", "section": "2.2.1.3 Initialize to Favor Reaction 1 and Solve" } }, "source": [ "Now let's solve the model." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.3 Initialize to Favor Reaction 1 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.3-Initialize-to-Favor-Reaction-1-and-Solve)", "section": "2.2.1.3 Initialize to Favor Reaction 1 and Solve" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Error evaluating \"var =\" definition -1: can't evaluate pow'(0,0.6).\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Ipopt 3.13.2: \n", "\n", "******************************************************************************\n", "This program contains Ipopt, a library for large-scale nonlinear optimization.\n", " Ipopt is released as open source code under the Eclipse Public License (EPL).\n", " For more information visit http://projects.coin-or.org/Ipopt\n", "******************************************************************************\n", "\n", "This is Ipopt version 3.13.2, running with linear solver ma27.\n", "\n", "Number of nonzeros in equality constraint Jacobian...: 9\n", "Number of nonzeros in inequality constraint Jacobian.: 0\n", "Number of nonzeros in Lagrangian Hessian.............: 2\n", "\n", "ERROR: Solver (ipopt) returned non-zero return code (1)\n", "ERROR: See the solver log above for diagnostic information.\n" ] }, { "ename": "ApplicationError", "evalue": "Solver (ipopt) did not exit normally", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mApplicationError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0msolver\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpyo\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSolverFactory\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'ipopt'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mresults\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnlp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtee\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m/anaconda3/envs/spring2021/lib/python3.8/site-packages/pyomo/opt/base/solvers.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(self, *args, **kwds)\u001b[0m\n\u001b[1;32m 595\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_status\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'log'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0m_status\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 596\u001b[0m \u001b[0mlogger\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0merror\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Solver log:\\n\"\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_status\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 597\u001b[0;31m raise ApplicationError(\n\u001b[0m\u001b[1;32m 598\u001b[0m \"Solver (%s) did not exit normally\" % self.name)\n\u001b[1;32m 599\u001b[0m \u001b[0msolve_completion_time\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mApplicationError\u001b[0m: Solver (ipopt) did not exit normally" ] } ], "source": [ "solver = pyo.SolverFactory('ipopt')\n", "results = solver.solve(nlp, tee=True)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.3 Initialize to Favor Reaction 1 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.3-Initialize-to-Favor-Reaction-1-and-Solve)", "section": "2.2.1.3 Initialize to Favor Reaction 1 and Solve" } }, "source": [ "What happened? $0^{0.6}$ is not well defined. Work around? Let's set the lower bound to something really small:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.3 Initialize to Favor Reaction 1 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.3-Initialize-to-Favor-Reaction-1-and-Solve)", "section": "2.2.1.3 Initialize to Favor Reaction 1 and Solve" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 Set Declarations\n", " REACTORS : Size=1, Index=None, Ordered=Insertion\n", " Key : Dimen : Domain : Size : Members\n", " None : 1 : Any : 2 : {1, 2}\n", "\n", "4 Param Declarations\n", " conversion : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " 1 : 0.8\n", " 2 : 0.67\n", " feed_cost : Size=1, Index=None, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " None : 5.0\n", " product_flowrate : Size=1, Index=None, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " None : 10.0\n", " reactor_cost : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " 1 : 5.5\n", " 2 : 4.0\n", "\n", "3 Var Declarations\n", " feed_flowrate : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : 0 : 20.0 : None : False : False : NonNegativeReals\n", " reactor_effluent : Size=2, Index=REACTORS\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 1 : 0 : 16.0 : None : False : False : NonNegativeReals\n", " 2 : 0 : 0.0 : None : False : False : NonNegativeReals\n", " reactor_feed : Size=2, Index=REACTORS\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 1 : 1e-06 : 20.0 : None : False : False : NonNegativeReals\n", " 2 : 1e-06 : 1e-06 : None : False : False : NonNegativeReals\n", "\n", "1 Objective Declarations\n", " cost : Size=1, Index=None, Active=True\n", " Key : Active : Sense : Expression\n", " None : True : minimize : 5.5*reactor_feed[1]**0.6 + 4.0*reactor_feed[2]**0.6 + 5.0*feed_flowrate\n", "\n", "3 Constraint Declarations\n", " inlet_split : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 0.0 : feed_flowrate - (reactor_feed[1] + reactor_feed[2]) : 0.0 : True\n", " mixer : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 10.0 : reactor_effluent[1] + reactor_effluent[2] : 10.0 : True\n", " reactor_performance : Size=2, Index=REACTORS, Active=True\n", " Key : Lower : Body : Upper : Active\n", " 1 : 0.0 : reactor_effluent[1] - 0.8*reactor_feed[1] : 0.0 : True\n", " 2 : 0.0 : reactor_effluent[2] - 0.67*reactor_feed[2] : 0.0 : True\n", "\n", "12 Declarations: REACTORS reactor_cost product_flowrate conversion feed_cost feed_flowrate reactor_feed reactor_effluent inlet_split reactor_performance mixer cost\n" ] } ], "source": [ "small_number = 1E-6\n", "for r in nlp.REACTORS:\n", " \n", " # Set lower bound\n", " nlp.reactor_feed[r].setlb(small_number)\n", " \n", " # Adjust initial point if needed\n", " nlp.reactor_feed[r] = max(nlp.reactor_feed[r](), small_number)\n", " \n", "nlp.pprint()" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.3 Initialize to Favor Reaction 1 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.3-Initialize-to-Favor-Reaction-1-and-Solve)", "section": "2.2.1.3 Initialize to Favor Reaction 1 and Solve" } }, "source": [ "Now let's resolve:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.3 Initialize to Favor Reaction 1 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.3-Initialize-to-Favor-Reaction-1-and-Solve)", "section": "2.2.1.3 Initialize to Favor Reaction 1 and Solve" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ipopt 3.13.2: \n", "\n", "******************************************************************************\n", "This program contains Ipopt, a library for large-scale nonlinear optimization.\n", " Ipopt is released as open source code under the Eclipse Public License (EPL).\n", " For more information visit http://projects.coin-or.org/Ipopt\n", "******************************************************************************\n", "\n", "This is Ipopt version 3.13.2, running with linear solver ma27.\n", "\n", "Number of nonzeros in equality constraint Jacobian...: 9\n", "Number of nonzeros in inequality constraint Jacobian.: 0\n", "Number of nonzeros in Lagrangian Hessian.............: 2\n", "\n", "Total number of variables............................: 5\n", " variables with only lower bounds: 5\n", " variables with lower and upper bounds: 0\n", " variables with only upper bounds: 0\n", "Total number of equality constraints.................: 4\n", "Total number of inequality constraints...............: 0\n", " inequality constraints with only lower bounds: 0\n", " inequality constraints with lower and upper bounds: 0\n", " inequality constraints with only upper bounds: 0\n", "\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 0 1.3344037e+02 6.01e+00 8.32e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", " 1 8.9498771e+01 1.72e-15 1.11e+01 -1.0 7.77e+00 - 4.44e-02 1.00e+00f 1\n", " 2 8.9572341e+01 6.11e-16 5.52e+00 -1.0 4.18e-02 - 9.98e-01 5.00e-01f 2\n", " 3 8.9546293e+01 8.88e-16 8.09e-05 -1.0 7.46e-03 - 1.00e+00 1.00e+00f 1\n", " 4 8.7594136e+01 3.89e-16 5.50e+00 -2.5 5.48e-01 - 9.88e-01 6.12e-01f 1\n", " 5 8.7593412e+01 1.95e-16 1.87e-01 -2.5 1.87e-05 4.0 1.00e+00 1.00e+00f 1\n", " 6 8.7589094e+01 1.08e-15 2.70e-02 -2.5 1.08e-04 - 1.00e+00 1.00e+00f 1\n", " 7 8.7533830e+01 1.00e-15 8.54e+01 -3.8 1.31e-03 - 1.00e+00 6.28e-01f 1\n", " 8 8.7542171e+01 4.36e-16 2.02e+02 -3.8 8.12e-05 - 2.29e-03 5.00e-01f 2\n", " 9 8.7536494e+01 1.78e-15 1.91e+01 -3.8 3.28e-05 5.3 1.00e+00 1.00e+00f 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 10 8.7535635e+01 3.64e-17 7.35e+00 -3.8 6.31e-06 5.8 1.00e+00 5.00e-01f 2\n", " 11 8.7536085e+01 1.78e-15 8.56e-01 -3.8 1.58e-06 - 1.00e+00 1.00e+00f 1\n", " 12 8.7536020e+01 1.78e-15 1.36e-02 -3.8 2.39e-07 - 1.00e+00 1.00e+00f 1\n", " 13 8.7536022e+01 1.06e-15 1.73e-05 -3.8 8.30e-09 - 1.00e+00 1.00e+00h 1\n", " 14 8.7533756e+01 1.78e-15 3.24e+01 -5.7 1.03e-05 - 1.00e+00 5.98e-01f 1\n", " 15 8.7533783e+01 9.34e-16 5.52e-02 -5.7 4.47e-08 - 1.00e+00 1.00e+00f 1\n", " 16 8.7533768e+01 1.78e-15 1.74e-02 -5.7 2.57e-08 - 1.00e+00 1.00e+00f 1\n", " 17 8.7533768e+01 2.11e-17 1.13e-07 -5.7 6.42e-11 - 1.00e+00 1.00e+00h 1\n", " 18 8.7533756e+01 1.53e-16 6.30e-02 -8.6 1.93e-08 - 1.00e+00 9.76e-01f 1\n", " 19 8.7533756e+01 1.78e-15 1.79e-08 -8.6 2.50e-11 - 1.00e+00 1.00e+00f 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 20 8.7533756e+01 1.78e-15 7.27e-09 -9.0 1.59e-11 - 1.00e+00 1.00e+00h 1\n", "\n", "Number of Iterations....: 20\n", "\n", " (scaled) (unscaled)\n", "Objective...............: 1.4519923356722005e+01 8.7533756319141716e+01\n", "Dual infeasibility......: 7.2660668593016453e-09 4.3803683410370973e-08\n", "Constraint violation....: 1.7763568394002505e-15 1.7763568394002505e-15\n", "Complementarity.........: 9.0911618532142856e-10 5.4806318653791501e-09\n", "Overall NLP error.......: 7.2660668593016453e-09 4.3803683410370973e-08\n", "\n", "\n", "Number of objective function evaluations = 26\n", "Number of objective gradient evaluations = 21\n", "Number of equality constraint evaluations = 26\n", "Number of inequality constraint evaluations = 0\n", "Number of equality constraint Jacobian evaluations = 21\n", "Number of inequality constraint Jacobian evaluations = 0\n", "Number of Lagrangian Hessian evaluations = 20\n", "Total CPU secs in IPOPT (w/o function evaluations) = 0.005\n", "Total CPU secs in NLP function evaluations = 0.000\n", "\n", "EXIT: Optimal Solution Found.\n" ] } ], "source": [ "solver = pyo.SolverFactory('ipopt')\n", "results = solver.solve(nlp, tee=True)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.3 Initialize to Favor Reaction 1 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.3-Initialize-to-Favor-Reaction-1-and-Solve)", "section": "2.2.1.3 Initialize to Favor Reaction 1 and Solve" } }, "source": [ "Now we can print the variable names and values:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.3 Initialize to Favor Reaction 1 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.3-Initialize-to-Favor-Reaction-1-and-Solve)", "section": "2.2.1.3 Initialize to Favor Reaction 1 and Solve" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Variable Names\t\tValue\n", "feed_flowrate \t\t 12.50000016087647\n", "reactor_feed[1] \t\t 12.499999170867413\n", "reactor_feed[2] \t\t 1e-06\n", "reactor_effluent[1] \t\t 9.99999933669393\n", "reactor_effluent[2] \t\t 6.633060682547189e-07\n", "\n", "Objective Name\t\tValue\n", "cost \t\t 87.53376235430073\n" ] } ], "source": [ "def print_solution(model):\n", " '''Print variable names and values\n", " \n", " Arguments:\n", " model: Pyomo model\n", " \n", " '''\n", "\n", " print(\"Variable Names\\t\\tValue\")\n", " for c in model.component_data_objects(pyo.Var):\n", " print(c.name,\"\\t\\t\", pyo.value(c))\n", " \n", " print(\"\\nObjective Name\\t\\tValue\")\n", " for c in model.component_data_objects(pyo.Objective):\n", " print(c.name,\"\\t\\t\", pyo.value(c))\n", " \n", "print_solution(nlp)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.4 Initialize to Favor Reaction 2 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.4-Initialize-to-Favor-Reaction-2-and-Solve)", "section": "2.2.1.4 Initialize to Favor Reaction 2 and Solve" } }, "source": [ "### 2.2.1.4 Initialize to Favor Reaction 2 and Solve" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.4 Initialize to Favor Reaction 2 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.4-Initialize-to-Favor-Reaction-2-and-Solve)", "section": "2.2.1.4 Initialize to Favor Reaction 2 and Solve" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ipopt 3.13.2: \n", "\n", "******************************************************************************\n", "This program contains Ipopt, a library for large-scale nonlinear optimization.\n", " Ipopt is released as open source code under the Eclipse Public License (EPL).\n", " For more information visit http://projects.coin-or.org/Ipopt\n", "******************************************************************************\n", "\n", "This is Ipopt version 3.13.2, running with linear solver ma27.\n", "\n", "Number of nonzeros in equality constraint Jacobian...: 9\n", "Number of nonzeros in inequality constraint Jacobian.: 0\n", "Number of nonzeros in Lagrangian Hessian.............: 2\n", "\n", "Total number of variables............................: 5\n", " variables with only lower bounds: 5\n", " variables with lower and upper bounds: 0\n", " variables with only upper bounds: 0\n", "Total number of equality constraints.................: 4\n", "Total number of inequality constraints...............: 0\n", " inequality constraints with only lower bounds: 0\n", " inequality constraints with lower and upper bounds: 0\n", " inequality constraints with only upper bounds: 0\n", "\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 0 1.2448375e+02 3.41e+00 8.54e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", " 1 9.6785883e+01 0.00e+00 9.38e+00 -1.0 5.37e+00 - 6.58e-02 1.00e+00f 1\n", " 2 9.9320862e+01 2.22e-16 2.05e+00 -1.0 2.25e+00 - 1.41e-01 1.00e+00f 1\n", " 3 9.9386245e+01 1.78e-15 5.06e-04 -1.0 2.06e-01 - 1.00e+00 1.00e+00f 1\n", " 4 9.8712942e+01 1.78e-15 3.67e-02 -1.7 1.26e+00 -2.0 9.62e-01 1.00e+00f 1\n", " 5 9.5055675e+01 1.78e-15 3.30e+00 -2.5 5.41e+00 -1.6 1.00e+00 2.76e-01f 1\n", " 6 9.4951632e+01 3.25e-19 2.68e+00 -2.5 3.26e-03 2.5 1.00e+00 1.00e+00f 1\n", " 7 9.4957287e+01 1.78e-15 1.23e+00 -2.5 2.52e-04 2.9 1.00e+00 5.00e-01f 2\n", " 8 9.4955819e+01 1.78e-15 1.75e-03 -2.5 3.33e-05 - 1.00e+00 1.00e+00f 1\n", " 9 9.4877882e+01 1.78e-15 8.54e+01 -3.8 1.71e-03 - 1.00e+00 6.18e-01f 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 10 9.4886607e+01 1.78e-15 1.19e+02 -3.8 1.30e-04 - 1.66e-03 2.50e-01f 3\n", " 11 9.4883285e+01 1.78e-15 9.54e+01 -3.8 6.25e-05 5.1 1.00e+00 2.60e-01f 2\n", " 12 9.4881130e+01 9.32e-21 9.72e+01 -3.8 3.54e-05 - 1.00e+00 2.30e-01f 2\n", " 13 9.4880984e+01 1.78e-15 3.37e-02 -3.8 4.78e-07 - 1.00e+00 1.00e+00f 1\n", " 14 9.4880897e+01 1.78e-15 1.28e-02 -3.8 2.79e-07 - 1.00e+00 1.00e+00f 1\n", " 15 9.4880898e+01 1.78e-15 1.69e-06 -3.8 3.12e-09 - 1.00e+00 1.00e+00h 1\n", " 16 9.4877775e+01 1.48e-21 3.25e+01 -5.7 1.24e-05 - 1.00e+00 5.96e-01f 1\n", " 17 9.4877812e+01 1.06e-22 5.58e-02 -5.7 5.36e-08 - 1.00e+00 1.00e+00f 1\n", " 18 9.4877790e+01 0.00e+00 1.77e-02 -5.7 3.10e-08 - 1.00e+00 1.00e+00f 1\n", " 19 9.4877790e+01 1.06e-22 1.14e-07 -5.7 7.69e-11 - 1.00e+00 1.00e+00f 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 20 9.4877775e+01 1.06e-22 5.31e-02 -8.6 2.31e-08 - 1.00e+00 9.76e-01f 1\n", " 21 9.4877775e+01 1.78e-15 1.80e-08 -8.6 2.99e-11 - 1.00e+00 1.00e+00f 1\n", " 22 9.4877775e+01 1.78e-15 7.30e-09 -9.0 1.90e-11 - 1.00e+00 1.00e+00h 1\n", "\n", "Number of Iterations....: 22\n", "\n", " (scaled) (unscaled)\n", "Objective...............: 1.1445915871146484e+01 9.4877774550678950e+01\n", "Dual infeasibility......: 7.2954691177073983e-09 6.0473786631279065e-08\n", "Constraint violation....: 1.7763568394002505e-15 1.7763568394002505e-15\n", "Complementarity.........: 9.0911629756661984e-10 7.5358697453220019e-09\n", "Overall NLP error.......: 7.2954691177073983e-09 6.0473786631279065e-08\n", "\n", "\n", "Number of objective function evaluations = 31\n", "Number of objective gradient evaluations = 23\n", "Number of equality constraint evaluations = 31\n", "Number of inequality constraint evaluations = 0\n", "Number of equality constraint Jacobian evaluations = 23\n", "Number of inequality constraint Jacobian evaluations = 0\n", "Number of Lagrangian Hessian evaluations = 22\n", "Total CPU secs in IPOPT (w/o function evaluations) = 0.006\n", "Total CPU secs in NLP function evaluations = 0.000\n", "\n", "EXIT: Optimal Solution Found.\n", "Variable Names\t\tValue\n", "feed_flowrate \t\t 14.925372942237043\n", "reactor_feed[1] \t\t 1e-06\n", "reactor_feed[2] \t\t 14.925371952227968\n", "reactor_effluent[1] \t\t 7.920072602986078e-07\n", "reactor_effluent[2] \t\t 9.99999920799274\n", "\n", "Objective Name\t\tValue\n", "cost \t\t 94.87778284900739\n" ] } ], "source": [ "# Initialize\n", "initialize(nlp, reactor_choice=2)\n", "\n", "# Correct for bound\n", "# Note: I would have put this in the initialize function but I wanted to show\n", "# the error in class\n", "for r in nlp.REACTORS:\n", " # Adjust initial point if needed\n", " nlp.reactor_feed[r] = max(nlp.reactor_feed[r](), small_number)\n", "\n", "results = solver.solve(nlp, tee=True)\n", "print_solution(nlp)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.4 Initialize to Favor Reaction 2 and Solve](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.4-Initialize-to-Favor-Reaction-2-and-Solve)", "section": "2.2.1.4 Initialize to Favor Reaction 2 and Solve" } }, "source": [ "Which solution is better? Why are there multiple solutions?" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.5 Discrete Cost Model](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.5-Discrete-Cost-Model)", "section": "2.2.1.5 Discrete Cost Model" } }, "source": [ "### 2.2.1.5 Discrete Cost Model" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.5 Discrete Cost Model](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.5-Discrete-Cost-Model)", "section": "2.2.1.5 Discrete Cost Model" } }, "source": [ "We want to modify the model to:\n", "* Easily find the best global solution (reduce impacts of initialization)\n", "* Account for the fact there is a minimum reactor size we can purchase" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.6 Enumerate the solutions](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.6-Enumerate-the-solutions)", "section": "2.2.1.6 Enumerate the solutions" } }, "source": [ "### 2.2.1.6 Enumerate the solutions" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.6 Enumerate the solutions](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.6-Enumerate-the-solutions)", "section": "2.2.1.6 Enumerate the solutions" } }, "source": [ "As an illustration, enumerate through the following four options:\n", "1. No reactor\n", "2. Reactor I only\n", "3. Reactor II only\n", "4. Reactor I and II only\n", "\n", "For each option, ask:\n", "* Are the constraints feasible?\n", "* What is the objective?" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.7 Solve with Pyomo](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.7-Solve-with-Pyomo)", "section": "2.2.1.7 Solve with Pyomo" } }, "source": [ "### 2.2.1.7 Solve with Pyomo" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.7 Solve with Pyomo](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.7-Solve-with-Pyomo)", "section": "2.2.1.7 Solve with Pyomo" } }, "source": [ "Create and inspect the model." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.7 Solve with Pyomo](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.7-Solve-with-Pyomo)", "section": "2.2.1.7 Solve with Pyomo" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 Set Declarations\n", " REACTORS : Size=1, Index=None, Ordered=Insertion\n", " Key : Dimen : Domain : Size : Members\n", " None : 1 : Any : 2 : {1, 2}\n", "\n", "6 Param Declarations\n", " conversion : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " 1 : 0.8\n", " 2 : 0.67\n", " feed_cost : Size=1, Index=None, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " None : 5.0\n", " max_flowrate : Size=1, Index=None, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " None : 20.0\n", " product_flowrate : Size=1, Index=None, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " None : 10.0\n", " reactor_cost_fixed : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " 1 : 7.5\n", " 2 : 5.5\n", " reactor_cost_linear : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False\n", " Key : Value\n", " 1 : 6.4\n", " 2 : 6.0\n", "\n", "4 Var Declarations\n", " feed_flowrate : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : 0 : None : 20.0 : False : True : NonNegativeReals\n", " reactor_boolean : Size=2, Index=REACTORS\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 1 : 0 : None : 1 : False : True : Boolean\n", " 2 : 0 : None : 1 : False : True : Boolean\n", " reactor_effluent : Size=2, Index=REACTORS\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 1 : 0 : None : None : False : True : NonNegativeReals\n", " 2 : 0 : None : None : False : True : NonNegativeReals\n", " reactor_feed : Size=2, Index=REACTORS\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 1 : 0 : None : 20.0 : False : True : NonNegativeReals\n", " 2 : 0 : None : 20.0 : False : True : NonNegativeReals\n", "\n", "1 Objective Declarations\n", " cost : Size=1, Index=None, Active=True\n", " Key : Active : Sense : Expression\n", " None : True : minimize : 6.4*reactor_feed[1] + 7.5*reactor_boolean[1] + 6.0*reactor_feed[2] + 5.5*reactor_boolean[2]\n", "\n", "4 Constraint Declarations\n", " inlet_split : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 0.0 : feed_flowrate - (reactor_feed[1] + reactor_feed[2]) : 0.0 : True\n", " mixer : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 10.0 : reactor_effluent[1] + reactor_effluent[2] : 10.0 : True\n", " reactor_performance : Size=2, Index=REACTORS, Active=True\n", " Key : Lower : Body : Upper : Active\n", " 1 : 0.0 : reactor_effluent[1] - 0.8*reactor_feed[1] : 0.0 : True\n", " 2 : 0.0 : reactor_effluent[2] - 0.67*reactor_feed[2] : 0.0 : True\n", " toggle_reactor : Size=2, Index=REACTORS, Active=True\n", " Key : Lower : Body : Upper : Active\n", " 1 : -Inf : reactor_feed[1] - 20.0*reactor_boolean[1] : 0.0 : True\n", " 2 : -Inf : reactor_feed[2] - 20.0*reactor_boolean[2] : 0.0 : True\n", "\n", "16 Declarations: REACTORS max_flowrate reactor_cost_linear reactor_cost_fixed product_flowrate conversion feed_cost feed_flowrate reactor_feed reactor_effluent reactor_boolean inlet_split reactor_performance mixer toggle_reactor cost\n" ] } ], "source": [ "milp = pyo.ConcreteModel()\n", "\n", "## Define sets\n", "milp.REACTORS = pyo.Set(initialize=range(1,3))\n", "\n", "## Define parameters (data)\n", "\n", "# kmol/hour\n", "milp.max_flowrate = pyo.Param(initialize=20.0)\n", "\n", "# $ / hr\n", "cost_coefficient1 = {1:6.4, 2:6.0}\n", "milp.reactor_cost_linear = pyo.Param(milp.REACTORS, initialize=cost_coefficient1)\n", "\n", "# $\n", "cost_coefficient2 = {1:7.5, 2:5.5}\n", "milp.reactor_cost_fixed = pyo.Param(milp.REACTORS, initialize=cost_coefficient2)\n", "\n", "# kmol/hr B\n", "milp.product_flowrate = pyo.Param(initialize=10.0)\n", "\n", "# conversion fraction\n", "reactor_conversion = {1:0.8, 2:0.67}\n", "milp.conversion = pyo.Param(milp.REACTORS, initialize=reactor_conversion)\n", "\n", "# feed cost, $/kmol of A\n", "milp.feed_cost = pyo.Param(initialize=5.0)\n", "\n", "\n", "## Define variables\n", "\n", "# Feed flowrate into reactor, x0 in handout illustration\n", "milp.feed_flowrate = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0, milp.max_flowrate))\n", "\n", "# Reactor feed, x1 and x2 in handout illustration\n", "milp.reactor_feed = pyo.Var(milp.REACTORS, domain=pyo.NonNegativeReals, bounds=(0, milp.max_flowrate))\n", "\n", "# Reactor effluent (outlet), z1 and z2 in handout illustration\n", "milp.reactor_effluent = pyo.Var(milp.REACTORS, domain=pyo.NonNegativeReals)\n", "\n", "# Boolean variables\n", "# YOUR SOLUTION HERE\n", "milp.pprint()" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.7 Solve with Pyomo](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.7-Solve-with-Pyomo)", "section": "2.2.1.7 Solve with Pyomo" } }, "source": [ "Solve the model using `cbc`, `bonmin`, `gurobi` (need license), or `cplex` (need license)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "nbpages": { "level": 3, "link": "[2.2.1.7 Solve with Pyomo](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.1.7-Solve-with-Pyomo)", "section": "2.2.1.7 Solve with Pyomo" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using license file /Users/adowling/gurobi.lic\n", "Academic license - for non-commercial use only\n", "Read LP format model from file /var/folders/xy/24xvnyss36v3d8mw68tygxdw0000gp/T/tmpqhpzlzbd.pyomo.lp\n", "Reading time = 0.00 seconds\n", "x8: 7 rows, 8 columns, 14 nonzeros\n", "Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)\n", "Optimize a model with 7 rows, 8 columns and 14 nonzeros\n", "Model fingerprint: 0x49a58b0f\n", "Variable types: 6 continuous, 2 integer (2 binary)\n", "Coefficient statistics:\n", " Matrix range [7e-01, 2e+01]\n", " Objective range [6e+00, 8e+00]\n", " Bounds range [1e+00, 2e+01]\n", " RHS range [1e+00, 1e+01]\n", "Presolve removed 7 rows and 8 columns\n", "Presolve time: 0.01s\n", "Presolve: All rows and columns removed\n", "\n", "Explored 0 nodes (0 simplex iterations) in 0.01 seconds\n", "Thread count was 1 (of 6 available processors)\n", "\n", "Solution count 1: 87.5 \n", "\n", "Optimal solution found (tolerance 1.00e-04)\n", "Best objective 8.750000000000e+01, best bound 8.750000000000e+01, gap 0.0000%\n", "Variable Names\t\tValue\n", "feed_flowrate \t\t 12.5\n", "reactor_feed[1] \t\t 12.5\n", "reactor_feed[2] \t\t 0.0\n", "reactor_effluent[1] \t\t 10.0\n", "reactor_effluent[2] \t\t 0.0\n", "reactor_boolean[1] \t\t 1.0\n", "reactor_boolean[2] \t\t 0.0\n", "\n", "Objective Name\t\tValue\n", "cost \t\t 87.5\n" ] } ], "source": [ "\n", "# Set solver\n", "solver = pyo.SolverFactory('gurobi')\n", "# solver = pyo.SolverFactory('glpk')\n", "# solver = pyo.SolverFactory('cbc')\n", "# solver = pyo.SolverFactory('bonmin')\n", "\n", "# YOUR SOLUTION HERE\n", "\n", "print_solution(milp)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[2.2.2 Is rounding good enough?](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2-Is-rounding-good-enough?)", "section": "2.2.2 Is rounding good enough?" } }, "source": [ "## 2.2.2 Is rounding good enough?" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.1 Linear Program (Relaxation)](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.1-Linear-Program-(Relaxation))", "section": "2.2.2.1 Linear Program (Relaxation)" } }, "source": [ "### 2.2.2.1 Linear Program (Relaxation)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.1 Linear Program (Relaxation)](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.1-Linear-Program-(Relaxation))", "section": "2.2.2.1 Linear Program (Relaxation)" } }, "source": [ "$$\\begin{align}\\min_{x_1,x_2} \\quad & x_2 \\\\\n", "\\mathrm{s.t.} \\quad & 2 x_1 + x_2 \\geq 13 \\\\\n", "& 5 x_1 + 2 x_2 \\leq 30 \\\\\n", "& -x_1 + x_2 \\geq 5 \\\\\n", "& x_1, x_2 \\geq 0\n", "\\end{align}$$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.1 Linear Program (Relaxation)](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.1-Linear-Program-(Relaxation))", "section": "2.2.2.1 Linear Program (Relaxation)" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2 Var Declarations\n", " x1 : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : 0 : None : None : False : True : NonNegativeReals\n", " x2 : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : 0 : None : None : False : True : NonNegativeReals\n", "\n", "1 Objective Declarations\n", " obj : Size=1, Index=None, Active=True\n", " Key : Active : Sense : Expression\n", " None : True : minimize : x2\n", "\n", "3 Constraint Declarations\n", " con1 : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 13.0 : 2*x1 + x2 : +Inf : True\n", " con2 : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : -Inf : 5*x1 + 2*x2 : 30.0 : True\n", " con3 : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 5.0 : - x1 + x2 : +Inf : True\n", "\n", "6 Declarations: x1 x2 con1 con2 con3 obj\n" ] } ], "source": [ "import pyomo.environ as pyo\n", "\n", "m = pyo.ConcreteModel()\n", "\n", "# Declare variables with bounds\n", "m.x1 = pyo.Var(domain=pyo.NonNegativeReals)\n", "m.x2 = pyo.Var(domain=pyo.NonNegativeReals)\n", "\n", "# Constraint 1\n", "m.con1 = pyo.Constraint(expr=2*m.x1 + m.x2 >= 13)\n", "\n", "# Constraint 2\n", "m.con2 = pyo.Constraint(expr=5*m.x1 + 2*m.x2 <= 30)\n", "\n", "# Constraint 3\n", "m.con3 = pyo.Constraint(expr=-m.x1 + m.x2 >= 5)\n", "\n", "# Objective\n", "m.obj = pyo.Objective(expr=m.x2)\n", "\n", "# Print model\n", "m.pprint()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.1 Linear Program (Relaxation)](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.1-Linear-Program-(Relaxation))", "section": "2.2.2.1 Linear Program (Relaxation)" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GLPSOL: GLPK LP/MIP Solver, v4.65\n", "Parameter(s) specified in the command line:\n", " --write /var/folders/xy/24xvnyss36v3d8mw68tygxdw0000gp/T/tmpmp918eg6.glpk.raw\n", " --wglp /var/folders/xy/24xvnyss36v3d8mw68tygxdw0000gp/T/tmpjwou7gyx.glpk.glp\n", " --cpxlp /var/folders/xy/24xvnyss36v3d8mw68tygxdw0000gp/T/tmponqo37ck.pyomo.lp\n", "Reading problem data from '/var/folders/xy/24xvnyss36v3d8mw68tygxdw0000gp/T/tmponqo37ck.pyomo.lp'...\n", "4 rows, 3 columns, 7 non-zeros\n", "30 lines were read\n", "Writing problem data to '/var/folders/xy/24xvnyss36v3d8mw68tygxdw0000gp/T/tmpjwou7gyx.glpk.glp'...\n", "22 lines were written\n", "GLPK Simplex Optimizer, v4.65\n", "4 rows, 3 columns, 7 non-zeros\n", "Preprocessing...\n", "3 rows, 2 columns, 6 non-zeros\n", "Scaling...\n", " A: min|aij| = 1.000e+00 max|aij| = 5.000e+00 ratio = 5.000e+00\n", "Problem data seem to be well scaled\n", "Constructing initial basis...\n", "Size of triangular part is 3\n", " 0: obj = 0.000000000e+00 inf = 1.800e+01 (2)\n", " 2: obj = 7.666666667e+00 inf = 0.000e+00 (0)\n", "OPTIMAL LP SOLUTION FOUND\n", "Time used: 0.0 secs\n", "Memory used: 0.0 Mb (40424 bytes)\n", "Writing basic solution to '/var/folders/xy/24xvnyss36v3d8mw68tygxdw0000gp/T/tmpmp918eg6.glpk.raw'...\n", "16 lines were written\n", " \n", "x1 = 2.66666666666667\n", "x2 = 7.66666666666667\n" ] } ], "source": [ "# Set solver\n", "# solver = pyo.SolverFactory('gurobi')\n", "solver = pyo.SolverFactory('glpk')\n", "# solver = pyo.SolverFactory('ipopt')\n", "\n", "\n", "# Solve\n", "solver.solve(m,tee=True)\n", "\n", "# Print solution\n", "print(\" \")\n", "print(\"x1 = \",pyo.value(m.x1))\n", "print(\"x2 = \",pyo.value(m.x2))" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.2 Rounding](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.2-Rounding)", "section": "2.2.2.2 Rounding" } }, "source": [ "### 2.2.2.2 Rounding" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.2 Rounding](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.2-Rounding)", "section": "2.2.2.2 Rounding" } }, "source": [ "With your neighbor, discuss:\n", "* If you round $x_1$ and $x_2$ to an integer, are all of the constraints feasible?\n", "* How would you go about checking the optimality of a feasible integer solution?" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.3 Integer Program](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.3-Integer-Program)", "section": "2.2.2.3 Integer Program" } }, "source": [ "### 2.2.2.3 Integer Program" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.3 Integer Program](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.3-Integer-Program)", "section": "2.2.2.3 Integer Program" } }, "source": [ "Consider the following integer program:\n", "\n", "$$\\begin{align}\\min_{x_1,x_2} \\quad & x_2 \\\\\n", "\\mathrm{s.t.} \\quad & 2 x_1 + x_2 \\geq 13 \\\\\n", "& 5 x_1 + 2 x_2 \\leq 30 \\\\\n", "& -x_1 + x_2 \\geq 5 \\\\\n", "& x_1, x_2 \\in \\mathcal{Z} := \\{0,1,2,...\\}\n", "\\end{align}$$" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.3 Integer Program](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.3-Integer-Program)", "section": "2.2.2.3 Integer Program" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2 Var Declarations\n", " x1 : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : 1 : None : None : False : True : PositiveIntegers\n", " x2 : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : 1 : None : None : False : True : PositiveIntegers\n", "\n", "1 Objective Declarations\n", " obj : Size=1, Index=None, Active=True\n", " Key : Active : Sense : Expression\n", " None : True : minimize : x2\n", "\n", "3 Constraint Declarations\n", " con1 : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 13.0 : 2*x1 + x2 : +Inf : True\n", " con2 : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : -Inf : 5*x1 + 2*x2 : 30.0 : True\n", " con3 : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 5.0 : - x1 + x2 : +Inf : True\n", "\n", "6 Declarations: x1 x2 con1 con2 con3 obj\n" ] } ], "source": [ "m2 = pyo.ConcreteModel()\n", "\n", "# Declare variables as positive integers\n", "m2.x1 = pyo.Var(domain=pyo.PositiveIntegers)\n", "m2.x2 = pyo.Var(domain=pyo.PositiveIntegers)\n", "\n", "# Constraint 1\n", "m2.con1 = pyo.Constraint(expr=2*m2.x1 + m2.x2 >= 13)\n", "\n", "# Constraint 2\n", "m2.con2 = pyo.Constraint(expr=5*m2.x1 + 2*m2.x2 <= 30)\n", "\n", "# Constraint 3\n", "m2.con3 = pyo.Constraint(expr=-m2.x1 + m2.x2 >= 5)\n", "\n", "# Objective\n", "m2.obj = pyo.Objective(expr=m2.x2)\n", "\n", "m2.pprint()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.3 Integer Program](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.3-Integer-Program)", "section": "2.2.2.3 Integer Program" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using license file /Users/adowling/gurobi.lic\n", "Academic license - for non-commercial use only\n", "Read LP format model from file /var/folders/xy/24xvnyss36v3d8mw68tygxdw0000gp/T/tmpmsmdp9dm.pyomo.lp\n", "Reading time = 0.00 seconds\n", "x3: 4 rows, 3 columns, 7 nonzeros\n", "Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)\n", "Optimize a model with 4 rows, 3 columns and 7 nonzeros\n", "Model fingerprint: 0x3a0f60a9\n", "Variable types: 1 continuous, 2 integer (0 binary)\n", "Coefficient statistics:\n", " Matrix range [1e+00, 5e+00]\n", " Objective range [1e+00, 1e+00]\n", " Bounds range [1e+00, 1e+00]\n", " RHS range [1e+00, 3e+01]\n", "Presolve removed 1 rows and 1 columns\n", "Presolve time: 0.00s\n", "Presolved: 3 rows, 2 columns, 6 nonzeros\n", "Variable types: 0 continuous, 2 integer (0 binary)\n", "Found heuristic solution: objective 10.0000000\n", "Found heuristic solution: objective 9.0000000\n", "\n", "Explored 0 nodes (0 simplex iterations) in 0.01 seconds\n", "Thread count was 6 (of 6 available processors)\n", "\n", "Solution count 2: 9 10 \n", "\n", "Optimal solution found (tolerance 1.00e-04)\n", "Best objective 9.000000000000e+00, best bound 9.000000000000e+00, gap 0.0000%\n", " \n", "x1 = 2.0\n", "x2 = 9.0\n" ] } ], "source": [ "# Set solver\n", "solver = pyo.SolverFactory('gurobi')\n", "# solver = pyo.SolverFactory('glpk')\n", "# solver = pyo.SolverFactory('cbc')\n", "# solver = pyo.SolverFactory('bonmin')\n", "\n", "# Solve\n", "solver.solve(m2,tee=True)\n", "\n", "# Print solution\n", "print(\" \")\n", "print(\"x1 = \",m2.x1())\n", "print(\"x2 = \",m2.x2())" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.4 Why rounding does not always work](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.4-Why-rounding-does-not-always-work)", "section": "2.2.2.4 Why rounding does not always work" } }, "source": [ "### 2.2.2.4 Why rounding does not always work" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.4 Why rounding does not always work](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.4-Why-rounding-does-not-always-work)", "section": "2.2.2.4 Why rounding does not always work" } }, "source": [ "Reference: http://personal.lse.ac.uk/williahp/publications/wp%2010-118.pdf" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.4 Why rounding does not always work](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.4-Why-rounding-does-not-always-work)", "section": "2.2.2.4 Why rounding does not always work" } }, "source": [ "![feasible](./figures/feasible.png)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[2.2.2.4 Why rounding does not always work](https://ndcbe.github.io/CBE60499/02.02-IP.html#2.2.2.4-Why-rounding-does-not-always-work)", "section": "2.2.2.4 Why rounding does not always work" } }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "eef0501b", "metadata": {}, "source": [ "\n", "< [2.1 Continuous Optimization](https://ndcbe.github.io/CBE60499/02.01-LP-NLP.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [2.3 Logical Modeling and Generalized Disjunctive Programs](https://ndcbe.github.io/CBE60499/02.03-GDP.html) >

\"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.8.6" } }, "nbformat": 4, "nbformat_minor": 2 }