diff --git a/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor.ipynb b/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor.ipynb index 2335f6a4..b65de628 100644 --- a/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor.ipynb +++ b/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor.ipynb @@ -306,7 +306,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now, let us add expressions to compute the reactor cooling cost (\\\\$/s) assuming a cost of 2.12E-5 \\\\$/kW, and the heating utility cost (\\\\$/s) assuming 2.2E-4 \\\\$/kW. Note that the heat duty is in units of watt (J/s). The total operating cost will be the sum of the two, expressed in \\\\$/year assuming 8000 operating hours per year (~10\\% downtime, which is fairly common for small scale chemical plants):" + "Now, let us add expressions to compute the reactor cooling cost in USD/s assuming a cost of 2.12E-5 USD/kW, and the heating utility cost in USD/s assuming 2.2E-4 USD/kW. Note that the heat duty is in units of watt (J/s). The total operating cost will be the sum of the two, expressed in USD/year assuming 8000 operating hours per year (~10 percent downtime, which is fairly common for small scale chemical plants):" ] }, { diff --git a/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_doc.ipynb b/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_doc.ipynb index d8a8866f..e171570f 100644 --- a/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_doc.ipynb +++ b/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_doc.ipynb @@ -1,693 +1,693 @@ { - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [ - "header", - "hide-cell" - ] - }, - "outputs": [], - "source": [ - "###############################################################################\n", - "# The Institute for the Design of Advanced Energy Systems Integrated Platform\n", - "# Framework (IDAES IP) was produced under the DOE Institute for the\n", - "# Design of Advanced Energy Systems (IDAES).\n", - "#\n", - "# Copyright (c) 2018-2023 by the software owners: The Regents of the\n", - "# University of California, through Lawrence Berkeley National Laboratory,\n", - "# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon\n", - "# University, West Virginia University Research Corporation, et al.\n", - "# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md\n", - "# for full copyright and license information.\n", - "###############################################################################" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "# Flowsheet Stoichiometric Reactor Simulation and Optimization of Ethylene Glycol Production\n", - "Author: Brandon Paul \n", - "Maintainer: Brandon Paul \n", - "Updated: 2023-06-01 \n", - "\n", - "## Learning Outcomes\n", - "\n", - "\n", - "- Call and implement the IDAES StochiometricReactor unit model\n", - "- Construct a steady-state flowsheet using the IDAES unit model library\n", - "- Connecting unit models in a flowsheet using Arcs\n", - "- Formulate and solve an optimization problem\n", - " - Defining an objective function\n", - " - Setting variable bounds\n", - " - Adding additional constraints \n", - "\n", - "\n", - "## Problem Statement\n", - "\n", - "Following the previous example implementing a [Continuous Stirred Tank Reactor (CSTR) unit model](http://localhost:8888/notebooks/GitHub/examples-pse/src/Examples/UnitModels/Reactors/cstr_testing_doc.md), we can alter the flowsheet to use a stochiometric (or yield) reactor. As before, this example is adapted from Fogler, H.S., Elements of Chemical Reaction Engineering 5th ed., 2016, Prentice Hall, p. 157-160 with the following chemical reaction, property packages and flowsheet. Unlike the previous two reactors which apply performance equations to calculate reaction extent, this simplified reactor model neglects all geometric properties and allows the user to specify a yield per reaction. The state variables chosen for the property package are **molar flows of each component by phase in each stream, temperature of each stream and pressure of each stream**. The components considered are: **ethylene oxide, water, sulfuric acid and ethylene glycol** and the process occurs in liquid phase only. Therefore, every stream has 4 flow variables, 1 temperature and 1 pressure variable.\n", - "\n", - "Chemical reaction:\n", - "\n", - "**C2H4O + H2O + H2SO4 \u2192 C2H6O2 + H2SO4**\n", - "\n", - "Property Packages:\n", - "\n", - "- egprod_ideal.py\n", - "- egprod_reaction.py\n", - "\n", - "Flowsheet:\n", - "\n", - "![](egprod_flowsheet.png)\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Importing Required Pyomo and IDAES components\n", - "\n", - "\n", - "To construct a flowsheet, we will need several components from the Pyomo and IDAES packages. Let us first import the following components from Pyomo:\n", - "- Constraint (to write constraints)\n", - "- Var (to declare variables)\n", - "- ConcreteModel (to create the concrete model object)\n", - "- Expression (to evaluate values as a function of variables defined in the model)\n", - "- Objective (to define an objective function for optimization)\n", - "- TransformationFactory (to apply certain transformations)\n", - "- Arc (to connect two unit models)\n", - "\n", - "For further details on these components, please refer to the pyomo documentation: https://pyomo.readthedocs.io/en/stable/\n", - "\n", - "From idaes, we will be needing the `FlowsheetBlock` and the following unit models:\n", - "- Mixer\n", - "- Heater\n", - "- StoichiometricReactor\n", - "\n", - "We will also be needing some utility tools to put together the flowsheet and calculate the degrees of freedom, tools for model expressions and calling variable values, and built-in functions to define property packages, add unit containers to objects and define our initialization scheme.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from pyomo.environ import (\n", - " Constraint,\n", - " Var,\n", - " ConcreteModel,\n", - " Expression,\n", - " Objective,\n", - " TransformationFactory,\n", - " value,\n", - " units as pyunits,\n", - ")\n", - "from pyomo.network import Arc\n", - "\n", - "from idaes.core import FlowsheetBlock\n", - "from idaes.models.properties.modular_properties.base.generic_property import (\n", - " GenericParameterBlock,\n", - ")\n", - "from idaes.models.properties.modular_properties.base.generic_reaction import (\n", - " GenericReactionParameterBlock,\n", - ")\n", - "from idaes.models.unit_models import Feed, Mixer, Heater, StoichiometricReactor, Product\n", - "\n", - "from idaes.core.solvers import get_solver\n", - "from idaes.core.util.model_statistics import degrees_of_freedom\n", - "from idaes.core.util.initialization import propagate_state" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Importing Required Thermophysical and Reaction Packages\n", - "\n", - "The final step is to import the thermophysical and reaction packages. We have created a custom thermophysical package that support ideal vapor and liquid behavior for this system, and in this case we will restrict it to ideal liquid behavior only. \n", - "\n", - "Let us import the following modules from the same directory as this Jupyter notebook:\n", - "- egprod_ideal as thermo_props\n", - "- egprod_reaction as reaction_props" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import egprod_ideal as thermo_props\n", - "import egprod_reaction as reaction_props" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Constructing the Flowsheet\n", - "\n", - "We have now imported all the components, unit models, and property modules we need to construct a flowsheet. Let us create a ConcreteModel and add the flowsheet block. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m = ConcreteModel()\n", - "m.fs = FlowsheetBlock(dynamic=False)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We now need to add the property packages to the flowsheet. Unlike the basic [Flash unit model example](http://localhost:8888/notebooks/GitHub/examples-pse/src/Tutorials/Basics/flash_unit_solution_testing_doc.md), where we only had a thermophysical property package, for this flowsheet we will also need to add a reaction property package. We will use the [Modular Property Framework](https://idaes-pse.readthedocs.io/en/stable/explanations/components/property_package/index.html#generic-property-package-framework) and [Modular Reaction Framework](https://idaes-pse.readthedocs.io/en/stable/explanations/components/property_package/index.html#generic-reaction-package-framework). The get_prop method for the natural gas property module automatically returns the correct dictionary using a component list argument. The GenericParameterBlock and GenericReactionParameterBlock methods build states blocks from passed parameter data; the reaction block unpacks using **reaction_props.config_dict to allow for optional or empty keyword arguments:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "m.fs.thermo_params = GenericParameterBlock(**thermo_props.config_dict)\n", - "m.fs.reaction_params = GenericReactionParameterBlock(\n", - " property_package=m.fs.thermo_params, **reaction_props.config_dict\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Adding Unit Models\n", - "\n", - "Let us start adding the unit models we have imported to the flowsheet. Here, we are adding a `Mixer`, a `Heater` and a `StoichiometricReactor`. Note that all unit models need to be given a property package argument. In addition to that, there are several arguments depending on the unit model, please refer to the documentation for more details on [IDAES Unit Models](https://idaes-pse.readthedocs.io/en/stable/reference_guides/model_libraries/index.html). For example, the `Mixer` is given a `list` consisting of names to the two inlets." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "m.fs.OXIDE = Feed(property_package=m.fs.thermo_params)\n", - "m.fs.ACID = Feed(property_package=m.fs.thermo_params)\n", - "m.fs.PROD = Product(property_package=m.fs.thermo_params)\n", - "m.fs.M101 = Mixer(\n", - " property_package=m.fs.thermo_params, inlet_list=[\"reagent_feed\", \"catalyst_feed\"]\n", - ")\n", - "m.fs.H101 = Heater(\n", - " property_package=m.fs.thermo_params,\n", - " has_pressure_change=False,\n", - " has_phase_equilibrium=False,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.R101 = StoichiometricReactor(\n", - " property_package=m.fs.thermo_params,\n", - " reaction_package=m.fs.reaction_params,\n", - " has_heat_of_reaction=True,\n", - " has_heat_transfer=True,\n", - " has_pressure_change=False,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Connecting Unit Models Using Arcs\n", - "\n", - "We have now added all the unit models we need to the flowsheet. However, we have not yet specified how the units are to be connected. To do this, we will be using the `Arc` which is a pyomo component that takes in two arguments: `source` and `destination`. Let us connect the outlet of the `Mixer` to the inlet of the `Heater`, and the outlet of the `Heater` to the inlet of the `StoichiometricReactor`. Additionally, we will connect the `Feed` and `Product` blocks to the flowsheet:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.s01 = Arc(source=m.fs.OXIDE.outlet, destination=m.fs.M101.reagent_feed)\n", - "m.fs.s02 = Arc(source=m.fs.ACID.outlet, destination=m.fs.M101.catalyst_feed)\n", - "m.fs.s03 = Arc(source=m.fs.M101.outlet, destination=m.fs.H101.inlet)\n", - "m.fs.s04 = Arc(source=m.fs.H101.outlet, destination=m.fs.R101.inlet)\n", - "m.fs.s05 = Arc(source=m.fs.R101.outlet, destination=m.fs.PROD.inlet)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We have now connected the unit model block using the arcs. However, we also need to link the state variables on connected ports. Pyomo provides a convenient method `TransformationFactory` to write these equality constraints for us between two ports:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "TransformationFactory(\"network.expand_arcs\").apply_to(m)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Adding Expressions to Compute Operating Costs\n", - "\n", - "In this section, we will add a few Expressions that allows us to evaluate the performance. `Expressions` provide a convenient way of calculating certain values that are a function of the variables defined in the model. For more details on `Expressions`, please refer to the [Pyomo Expression documentation]( https://pyomo.readthedocs.io/en/stable/pyomo_modeling_components/Expressions.html).\n", - "\n", - "For this flowsheet, we are interested in computing ethylene glycol production in millions of pounds per year, as well as the total costs due to cooling and heating utilities." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let us first add an `Expression` to convert the product flow from mol/s to MM lb/year of ethylene glycol. We see that the molecular weight exists in the thermophysical property package, so we may use that value for our calculations." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.eg_prod = Expression(\n", - " expr=pyunits.convert(\n", - " m.fs.PROD.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", - " * m.fs.thermo_params.ethylene_glycol.mw, # MW defined in properties as kg/mol\n", - " to_units=pyunits.Mlb / pyunits.yr,\n", - " )\n", - ") # converting kg/s to MM lb/year" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now, let us add expressions to compute the reactor cooling cost (\\\\$/s) assuming a cost of 2.12E-5 \\\\$/kW, and the heating utility cost (\\\\$/s) assuming 2.2E-4 \\\\$/kW. Note that the heat duty is in units of watt (J/s). The total operating cost will be the sum of the two, expressed in \\\\$/year assuming 8000 operating hours per year (~10\\% downtime, which is fairly common for small scale chemical plants):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.cooling_cost = Expression(\n", - " expr=2.12e-8 * (-m.fs.R101.heat_duty[0])\n", - ") # the reaction is exothermic, so R101 duty is negative\n", - "m.fs.heating_cost = Expression(\n", - " expr=2.2e-7 * m.fs.H101.heat_duty[0]\n", - ") # the stream must be heated to T_rxn, so H101 duty is positive\n", - "m.fs.operating_cost = Expression(\n", - " expr=(3600 * 8000 * (m.fs.heating_cost + m.fs.cooling_cost))\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Fixing Feed Conditions\n", - "\n", - "Let us first check how many degrees of freedom exist for this flowsheet using the `degrees_of_freedom` tool we imported earlier. We expect each stream to have 6 degrees of freedom, the mixer to have 0 (after both streams are accounted for), the heater to have 1 (just the duty, since the inlet is also the outlet of M101), and the reactor to have 1 (duty or overall conversion, since the inlet is also the outlet of H101). In this case, the reactor has an extra degree of freedom since we have not yet defined the yield of the sole rate-kinetics reaction. Therefore, we have 15 degrees of freedom to specify: temperature, pressure and flow of all four components on both streams; outlet heater temperature; reactor conversion and duty." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "print(degrees_of_freedom(m))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We will now be fixing the feed stream to the conditions shown in the flowsheet above. As mentioned in other tutorials, the IDAES framework expects a time index value for every referenced internal stream or unit variable, even in steady-state systems with a single time point $ t = 0 $ (`t = [0]` is the default when creating a `FlowsheetBlock` without passing a `time_set` argument). The non-present components in each stream are assigned a very small non-zero value to help with convergence and initializing. Based on stoichiometric ratios for the reaction, 80% conversion and 200 MM lb/year (46.4 mol/s) of ethylene glycol, we will initialize our simulation with the following calculated values:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.OXIDE.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_oxide\"].fix(\n", - " 58.0 * pyunits.mol / pyunits.s\n", - ")\n", - "m.fs.OXIDE.outlet.flow_mol_phase_comp[0, \"Liq\", \"water\"].fix(\n", - " 39.6 * pyunits.mol / pyunits.s\n", - ") # calculated from 16.1 mol EO / cudm in stream\n", - "m.fs.OXIDE.outlet.flow_mol_phase_comp[0, \"Liq\", \"sulfuric_acid\"].fix(\n", - " 1e-5 * pyunits.mol / pyunits.s\n", - ")\n", - "m.fs.OXIDE.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"].fix(\n", - " 1e-5 * pyunits.mol / pyunits.s\n", - ")\n", - "m.fs.OXIDE.outlet.temperature.fix(298.15 * pyunits.K)\n", - "m.fs.OXIDE.outlet.pressure.fix(1e5 * pyunits.Pa)\n", - "\n", - "m.fs.ACID.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_oxide\"].fix(\n", - " 1e-5 * pyunits.mol / pyunits.s\n", - ")\n", - "m.fs.ACID.outlet.flow_mol_phase_comp[0, \"Liq\", \"water\"].fix(\n", - " 200 * pyunits.mol / pyunits.s\n", - ")\n", - "m.fs.ACID.outlet.flow_mol_phase_comp[0, \"Liq\", \"sulfuric_acid\"].fix(\n", - " 0.334 * pyunits.mol / pyunits.s\n", - ") # calculated from 0.9 wt% SA in stream\n", - "m.fs.ACID.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"].fix(\n", - " 1e-5 * pyunits.mol / pyunits.s\n", - ")\n", - "m.fs.ACID.outlet.temperature.fix(298.15 * pyunits.K)\n", - "m.fs.ACID.outlet.pressure.fix(1e5 * pyunits.Pa)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Fixing Unit Model Specifications\n", - "\n", - "Now that we have fixed our inlet feed conditions, we will now be fixing the operating conditions for the unit models in the flowsheet. Let us fix the outlet temperature of H101 to 328.15 K. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.H101.outlet.temperature.fix(328.15 * pyunits.K)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We will need to specify both initial reactant extent (conversion or yield) and heat duty values (these are the only two free variables to choose from). The reaction extent can be specified directly, as a molar or mass yield ratio of product to a particular reactant, or fractional conversion of a particular reactant. Here, we choose fractional conversion in terms of ethylene oxide. Since heat duty and the outlet reactor temperature are interdependent, we can choose to specify this quantity instead. While the reaction kinetic parameters exist in the property package, we also do not need to add a rate constant expression since generation is explicitly defined through the conversion/yield. Note that our initial problem will solve with zero *temperature change* but will be infeasible with zero *heat duty*; this is due to the heat of reaction enforced by allowing heat transfer and mandating a non-zero conversion." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.R101.conversion = Var(\n", - " initialize=0.80, bounds=(0, 1), units=pyunits.dimensionless\n", - ") # fraction\n", - "\n", - "m.fs.R101.conv_constraint = Constraint(\n", - " expr=m.fs.R101.conversion\n", - " * m.fs.R101.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_oxide\"]\n", - " == (\n", - " m.fs.R101.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_oxide\"]\n", - " - m.fs.R101.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_oxide\"]\n", - " )\n", - ")\n", - "\n", - "m.fs.R101.conversion.fix(0.80)\n", - "\n", - "m.fs.R101.outlet.temperature.fix(328.15 * pyunits.K) # equal inlet reactor temperature" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For initialization, we solve a square problem (degrees of freedom = 0). Let's check the degrees of freedom below:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(degrees_of_freedom(m))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, we need to initialize the each unit operation in sequence to solve the flowsheet. As in best practice, unit operations are initialized or solved, and outlet properties are propagated to connected inlet streams via arc definitions as follows:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Initialize and solve each unit operation\n", - "m.fs.OXIDE.initialize()\n", - "propagate_state(arc=m.fs.s01)\n", - "\n", - "m.fs.ACID.initialize()\n", - "propagate_state(arc=m.fs.s01)\n", - "\n", - "m.fs.M101.initialize()\n", - "propagate_state(arc=m.fs.s03)\n", - "\n", - "m.fs.H101.initialize()\n", - "propagate_state(arc=m.fs.s04)\n", - "\n", - "m.fs.R101.initialize()\n", - "propagate_state(arc=m.fs.s05)\n", - "\n", - "m.fs.PROD.initialize()\n", - "\n", - "# set solver\n", - "solver = get_solver()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "# Solve the model\n", - "results = solver.solve(m, tee=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Analyze the Results of the Square Problem\n", - "\n", - "\n", - "What is the total operating cost? " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(f\"operating cost = ${value(m.fs.operating_cost)/1e6:0.6f} million per year\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For this operating cost, what conversion did we achieve of ethylene oxide to ethylene glycol? " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.R101.report()\n", - "\n", - "print()\n", - "print(f\"Conversion achieved = {value(m.fs.R101.conversion):.1%}\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Optimizing Ethylene Glycol Production\n", - "\n", - "Now that the flowsheet has been squared and solved, we can run a small optimization problem to minimize our production costs. Suppose we require at least 200 million pounds/year of ethylene glycol produced and 90% conversion of ethylene oxide, allowing for variable and reactor temperature (heater outlet)." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let us declare our objective function for this problem. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.objective = Objective(expr=m.fs.operating_cost)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now, we need to add the design constraints and unfix the decision variables as we had solved a square problem (degrees of freedom = 0) until now, as well as set bounds for the design variables (reactor outlet temperature is set by state variable bounds in property package):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.eg_prod_con = Constraint(\n", - " expr=m.fs.eg_prod >= 200 * pyunits.Mlb / pyunits.yr\n", - ") # MM lb/year\n", - "m.fs.R101.conversion.fix(0.90)\n", - "\n", - "m.fs.H101.outlet.temperature.unfix()\n", - "m.fs.H101.outlet.temperature[0].setlb(328.15 * pyunits.K)\n", - "m.fs.H101.outlet.temperature[0].setub(\n", - " 470.45 * pyunits.K\n", - ") # highest component boiling point (ethylene glycol)\n", - "\n", - "m.fs.R101.outlet.temperature.unfix()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "We have now defined the optimization problem and we are now ready to solve this problem. \n", - "\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "results = solver.solve(m, tee=True)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(f\"operating cost = ${value(m.fs.operating_cost)/1e6:0.6f} million per year\")\n", - "\n", - "print()\n", - "print(\"Heater results\")\n", - "\n", - "m.fs.H101.report()\n", - "\n", - "print()\n", - "print(\"Stoichiometric reactor results\")\n", - "\n", - "m.fs.R101.report()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Display optimal values for the decision variables and design variables:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(\"Optimal Values\")\n", - "print()\n", - "\n", - "print(f\"H101 outlet temperature = {value(m.fs.H101.outlet.temperature[0]):0.6f} K\")\n", - "\n", - "print()\n", - "print(f\"R101 outlet temperature = {value(m.fs.R101.outlet.temperature[0]):0.6f} K\")\n", - "\n", - "print()\n", - "print(f\"Ethylene glycol produced = {value(m.fs.eg_prod):0.6f} MM lb/year\")\n", - "\n", - "print()\n", - "print(f\"Conversion achieved = {value(m.fs.R101.conversion):.1%}\")" - ] - }, - { - "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.18" - } - }, - "nbformat": 4, - "nbformat_minor": 3 + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "header", + "hide-cell" + ] + }, + "outputs": [], + "source": [ + "###############################################################################\n", + "# The Institute for the Design of Advanced Energy Systems Integrated Platform\n", + "# Framework (IDAES IP) was produced under the DOE Institute for the\n", + "# Design of Advanced Energy Systems (IDAES).\n", + "#\n", + "# Copyright (c) 2018-2023 by the software owners: The Regents of the\n", + "# University of California, through Lawrence Berkeley National Laboratory,\n", + "# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon\n", + "# University, West Virginia University Research Corporation, et al.\n", + "# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md\n", + "# for full copyright and license information.\n", + "###############################################################################" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "# Flowsheet Stoichiometric Reactor Simulation and Optimization of Ethylene Glycol Production\n", + "Author: Brandon Paul \n", + "Maintainer: Brandon Paul \n", + "Updated: 2023-06-01 \n", + "\n", + "## Learning Outcomes\n", + "\n", + "\n", + "- Call and implement the IDAES StochiometricReactor unit model\n", + "- Construct a steady-state flowsheet using the IDAES unit model library\n", + "- Connecting unit models in a flowsheet using Arcs\n", + "- Formulate and solve an optimization problem\n", + " - Defining an objective function\n", + " - Setting variable bounds\n", + " - Adding additional constraints \n", + "\n", + "\n", + "## Problem Statement\n", + "\n", + "Following the previous example implementing a [Continuous Stirred Tank Reactor (CSTR) unit model](http://localhost:8888/notebooks/GitHub/examples-pse/src/Examples/UnitModels/Reactors/cstr_testing_doc.md), we can alter the flowsheet to use a stochiometric (or yield) reactor. As before, this example is adapted from Fogler, H.S., Elements of Chemical Reaction Engineering 5th ed., 2016, Prentice Hall, p. 157-160 with the following chemical reaction, property packages and flowsheet. Unlike the previous two reactors which apply performance equations to calculate reaction extent, this simplified reactor model neglects all geometric properties and allows the user to specify a yield per reaction. The state variables chosen for the property package are **molar flows of each component by phase in each stream, temperature of each stream and pressure of each stream**. The components considered are: **ethylene oxide, water, sulfuric acid and ethylene glycol** and the process occurs in liquid phase only. Therefore, every stream has 4 flow variables, 1 temperature and 1 pressure variable.\n", + "\n", + "Chemical reaction:\n", + "\n", + "**C2H4O + H2O + H2SO4 → C2H6O2 + H2SO4**\n", + "\n", + "Property Packages:\n", + "\n", + "- egprod_ideal.py\n", + "- egprod_reaction.py\n", + "\n", + "Flowsheet:\n", + "\n", + "![](egprod_flowsheet.png)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Importing Required Pyomo and IDAES components\n", + "\n", + "\n", + "To construct a flowsheet, we will need several components from the Pyomo and IDAES packages. Let us first import the following components from Pyomo:\n", + "- Constraint (to write constraints)\n", + "- Var (to declare variables)\n", + "- ConcreteModel (to create the concrete model object)\n", + "- Expression (to evaluate values as a function of variables defined in the model)\n", + "- Objective (to define an objective function for optimization)\n", + "- TransformationFactory (to apply certain transformations)\n", + "- Arc (to connect two unit models)\n", + "\n", + "For further details on these components, please refer to the pyomo documentation: https://pyomo.readthedocs.io/en/stable/\n", + "\n", + "From idaes, we will be needing the `FlowsheetBlock` and the following unit models:\n", + "- Mixer\n", + "- Heater\n", + "- StoichiometricReactor\n", + "\n", + "We will also be needing some utility tools to put together the flowsheet and calculate the degrees of freedom, tools for model expressions and calling variable values, and built-in functions to define property packages, add unit containers to objects and define our initialization scheme.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyomo.environ import (\n", + " Constraint,\n", + " Var,\n", + " ConcreteModel,\n", + " Expression,\n", + " Objective,\n", + " TransformationFactory,\n", + " value,\n", + " units as pyunits,\n", + ")\n", + "from pyomo.network import Arc\n", + "\n", + "from idaes.core import FlowsheetBlock\n", + "from idaes.models.properties.modular_properties.base.generic_property import (\n", + " GenericParameterBlock,\n", + ")\n", + "from idaes.models.properties.modular_properties.base.generic_reaction import (\n", + " GenericReactionParameterBlock,\n", + ")\n", + "from idaes.models.unit_models import Feed, Mixer, Heater, StoichiometricReactor, Product\n", + "\n", + "from idaes.core.solvers import get_solver\n", + "from idaes.core.util.model_statistics import degrees_of_freedom\n", + "from idaes.core.util.initialization import propagate_state" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Importing Required Thermophysical and Reaction Packages\n", + "\n", + "The final step is to import the thermophysical and reaction packages. We have created a custom thermophysical package that support ideal vapor and liquid behavior for this system, and in this case we will restrict it to ideal liquid behavior only. \n", + "\n", + "Let us import the following modules from the same directory as this Jupyter notebook:\n", + "- egprod_ideal as thermo_props\n", + "- egprod_reaction as reaction_props" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import egprod_ideal as thermo_props\n", + "import egprod_reaction as reaction_props" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Constructing the Flowsheet\n", + "\n", + "We have now imported all the components, unit models, and property modules we need to construct a flowsheet. Let us create a ConcreteModel and add the flowsheet block. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m = ConcreteModel()\n", + "m.fs = FlowsheetBlock(dynamic=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now need to add the property packages to the flowsheet. Unlike the basic [Flash unit model example](http://localhost:8888/notebooks/GitHub/examples-pse/src/Tutorials/Basics/flash_unit_solution_testing_doc.md), where we only had a thermophysical property package, for this flowsheet we will also need to add a reaction property package. We will use the [Modular Property Framework](https://idaes-pse.readthedocs.io/en/stable/explanations/components/property_package/index.html#generic-property-package-framework) and [Modular Reaction Framework](https://idaes-pse.readthedocs.io/en/stable/explanations/components/property_package/index.html#generic-reaction-package-framework). The get_prop method for the natural gas property module automatically returns the correct dictionary using a component list argument. The GenericParameterBlock and GenericReactionParameterBlock methods build states blocks from passed parameter data; the reaction block unpacks using **reaction_props.config_dict to allow for optional or empty keyword arguments:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "m.fs.thermo_params = GenericParameterBlock(**thermo_props.config_dict)\n", + "m.fs.reaction_params = GenericReactionParameterBlock(\n", + " property_package=m.fs.thermo_params, **reaction_props.config_dict\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Adding Unit Models\n", + "\n", + "Let us start adding the unit models we have imported to the flowsheet. Here, we are adding a `Mixer`, a `Heater` and a `StoichiometricReactor`. Note that all unit models need to be given a property package argument. In addition to that, there are several arguments depending on the unit model, please refer to the documentation for more details on [IDAES Unit Models](https://idaes-pse.readthedocs.io/en/stable/reference_guides/model_libraries/index.html). For example, the `Mixer` is given a `list` consisting of names to the two inlets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "m.fs.OXIDE = Feed(property_package=m.fs.thermo_params)\n", + "m.fs.ACID = Feed(property_package=m.fs.thermo_params)\n", + "m.fs.PROD = Product(property_package=m.fs.thermo_params)\n", + "m.fs.M101 = Mixer(\n", + " property_package=m.fs.thermo_params, inlet_list=[\"reagent_feed\", \"catalyst_feed\"]\n", + ")\n", + "m.fs.H101 = Heater(\n", + " property_package=m.fs.thermo_params,\n", + " has_pressure_change=False,\n", + " has_phase_equilibrium=False,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.R101 = StoichiometricReactor(\n", + " property_package=m.fs.thermo_params,\n", + " reaction_package=m.fs.reaction_params,\n", + " has_heat_of_reaction=True,\n", + " has_heat_transfer=True,\n", + " has_pressure_change=False,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Connecting Unit Models Using Arcs\n", + "\n", + "We have now added all the unit models we need to the flowsheet. However, we have not yet specified how the units are to be connected. To do this, we will be using the `Arc` which is a pyomo component that takes in two arguments: `source` and `destination`. Let us connect the outlet of the `Mixer` to the inlet of the `Heater`, and the outlet of the `Heater` to the inlet of the `StoichiometricReactor`. Additionally, we will connect the `Feed` and `Product` blocks to the flowsheet:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.s01 = Arc(source=m.fs.OXIDE.outlet, destination=m.fs.M101.reagent_feed)\n", + "m.fs.s02 = Arc(source=m.fs.ACID.outlet, destination=m.fs.M101.catalyst_feed)\n", + "m.fs.s03 = Arc(source=m.fs.M101.outlet, destination=m.fs.H101.inlet)\n", + "m.fs.s04 = Arc(source=m.fs.H101.outlet, destination=m.fs.R101.inlet)\n", + "m.fs.s05 = Arc(source=m.fs.R101.outlet, destination=m.fs.PROD.inlet)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have now connected the unit model block using the arcs. However, we also need to link the state variables on connected ports. Pyomo provides a convenient method `TransformationFactory` to write these equality constraints for us between two ports:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TransformationFactory(\"network.expand_arcs\").apply_to(m)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Adding Expressions to Compute Operating Costs\n", + "\n", + "In this section, we will add a few Expressions that allows us to evaluate the performance. `Expressions` provide a convenient way of calculating certain values that are a function of the variables defined in the model. For more details on `Expressions`, please refer to the [Pyomo Expression documentation]( https://pyomo.readthedocs.io/en/stable/pyomo_modeling_components/Expressions.html).\n", + "\n", + "For this flowsheet, we are interested in computing ethylene glycol production in millions of pounds per year, as well as the total costs due to cooling and heating utilities." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us first add an `Expression` to convert the product flow from mol/s to MM lb/year of ethylene glycol. We see that the molecular weight exists in the thermophysical property package, so we may use that value for our calculations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.eg_prod = Expression(\n", + " expr=pyunits.convert(\n", + " m.fs.PROD.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", + " * m.fs.thermo_params.ethylene_glycol.mw, # MW defined in properties as kg/mol\n", + " to_units=pyunits.Mlb / pyunits.yr,\n", + " )\n", + ") # converting kg/s to MM lb/year" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, let us add expressions to compute the reactor cooling cost in USD/s assuming a cost of 2.12E-5 USD/kW, and the heating utility cost in USD/s assuming 2.2E-4 USD/kW. Note that the heat duty is in units of watt (J/s). The total operating cost will be the sum of the two, expressed in USD/year assuming 8000 operating hours per year (~10 percent downtime, which is fairly common for small scale chemical plants):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.cooling_cost = Expression(\n", + " expr=2.12e-8 * (-m.fs.R101.heat_duty[0])\n", + ") # the reaction is exothermic, so R101 duty is negative\n", + "m.fs.heating_cost = Expression(\n", + " expr=2.2e-7 * m.fs.H101.heat_duty[0]\n", + ") # the stream must be heated to T_rxn, so H101 duty is positive\n", + "m.fs.operating_cost = Expression(\n", + " expr=(3600 * 8000 * (m.fs.heating_cost + m.fs.cooling_cost))\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Fixing Feed Conditions\n", + "\n", + "Let us first check how many degrees of freedom exist for this flowsheet using the `degrees_of_freedom` tool we imported earlier. We expect each stream to have 6 degrees of freedom, the mixer to have 0 (after both streams are accounted for), the heater to have 1 (just the duty, since the inlet is also the outlet of M101), and the reactor to have 1 (duty or overall conversion, since the inlet is also the outlet of H101). In this case, the reactor has an extra degree of freedom since we have not yet defined the yield of the sole rate-kinetics reaction. Therefore, we have 15 degrees of freedom to specify: temperature, pressure and flow of all four components on both streams; outlet heater temperature; reactor conversion and duty." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "print(degrees_of_freedom(m))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will now be fixing the feed stream to the conditions shown in the flowsheet above. As mentioned in other tutorials, the IDAES framework expects a time index value for every referenced internal stream or unit variable, even in steady-state systems with a single time point $ t = 0 $ (`t = [0]` is the default when creating a `FlowsheetBlock` without passing a `time_set` argument). The non-present components in each stream are assigned a very small non-zero value to help with convergence and initializing. Based on stoichiometric ratios for the reaction, 80% conversion and 200 MM lb/year (46.4 mol/s) of ethylene glycol, we will initialize our simulation with the following calculated values:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.OXIDE.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_oxide\"].fix(\n", + " 58.0 * pyunits.mol / pyunits.s\n", + ")\n", + "m.fs.OXIDE.outlet.flow_mol_phase_comp[0, \"Liq\", \"water\"].fix(\n", + " 39.6 * pyunits.mol / pyunits.s\n", + ") # calculated from 16.1 mol EO / cudm in stream\n", + "m.fs.OXIDE.outlet.flow_mol_phase_comp[0, \"Liq\", \"sulfuric_acid\"].fix(\n", + " 1e-5 * pyunits.mol / pyunits.s\n", + ")\n", + "m.fs.OXIDE.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"].fix(\n", + " 1e-5 * pyunits.mol / pyunits.s\n", + ")\n", + "m.fs.OXIDE.outlet.temperature.fix(298.15 * pyunits.K)\n", + "m.fs.OXIDE.outlet.pressure.fix(1e5 * pyunits.Pa)\n", + "\n", + "m.fs.ACID.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_oxide\"].fix(\n", + " 1e-5 * pyunits.mol / pyunits.s\n", + ")\n", + "m.fs.ACID.outlet.flow_mol_phase_comp[0, \"Liq\", \"water\"].fix(\n", + " 200 * pyunits.mol / pyunits.s\n", + ")\n", + "m.fs.ACID.outlet.flow_mol_phase_comp[0, \"Liq\", \"sulfuric_acid\"].fix(\n", + " 0.334 * pyunits.mol / pyunits.s\n", + ") # calculated from 0.9 wt% SA in stream\n", + "m.fs.ACID.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"].fix(\n", + " 1e-5 * pyunits.mol / pyunits.s\n", + ")\n", + "m.fs.ACID.outlet.temperature.fix(298.15 * pyunits.K)\n", + "m.fs.ACID.outlet.pressure.fix(1e5 * pyunits.Pa)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Fixing Unit Model Specifications\n", + "\n", + "Now that we have fixed our inlet feed conditions, we will now be fixing the operating conditions for the unit models in the flowsheet. Let us fix the outlet temperature of H101 to 328.15 K. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.H101.outlet.temperature.fix(328.15 * pyunits.K)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will need to specify both initial reactant extent (conversion or yield) and heat duty values (these are the only two free variables to choose from). The reaction extent can be specified directly, as a molar or mass yield ratio of product to a particular reactant, or fractional conversion of a particular reactant. Here, we choose fractional conversion in terms of ethylene oxide. Since heat duty and the outlet reactor temperature are interdependent, we can choose to specify this quantity instead. While the reaction kinetic parameters exist in the property package, we also do not need to add a rate constant expression since generation is explicitly defined through the conversion/yield. Note that our initial problem will solve with zero *temperature change* but will be infeasible with zero *heat duty*; this is due to the heat of reaction enforced by allowing heat transfer and mandating a non-zero conversion." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.R101.conversion = Var(\n", + " initialize=0.80, bounds=(0, 1), units=pyunits.dimensionless\n", + ") # fraction\n", + "\n", + "m.fs.R101.conv_constraint = Constraint(\n", + " expr=m.fs.R101.conversion\n", + " * m.fs.R101.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_oxide\"]\n", + " == (\n", + " m.fs.R101.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_oxide\"]\n", + " - m.fs.R101.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_oxide\"]\n", + " )\n", + ")\n", + "\n", + "m.fs.R101.conversion.fix(0.80)\n", + "\n", + "m.fs.R101.outlet.temperature.fix(328.15 * pyunits.K) # equal inlet reactor temperature" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For initialization, we solve a square problem (degrees of freedom = 0). Let's check the degrees of freedom below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(degrees_of_freedom(m))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we need to initialize the each unit operation in sequence to solve the flowsheet. As in best practice, unit operations are initialized or solved, and outlet properties are propagated to connected inlet streams via arc definitions as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Initialize and solve each unit operation\n", + "m.fs.OXIDE.initialize()\n", + "propagate_state(arc=m.fs.s01)\n", + "\n", + "m.fs.ACID.initialize()\n", + "propagate_state(arc=m.fs.s01)\n", + "\n", + "m.fs.M101.initialize()\n", + "propagate_state(arc=m.fs.s03)\n", + "\n", + "m.fs.H101.initialize()\n", + "propagate_state(arc=m.fs.s04)\n", + "\n", + "m.fs.R101.initialize()\n", + "propagate_state(arc=m.fs.s05)\n", + "\n", + "m.fs.PROD.initialize()\n", + "\n", + "# set solver\n", + "solver = get_solver()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Solve the model\n", + "results = solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Analyze the Results of the Square Problem\n", + "\n", + "\n", + "What is the total operating cost? " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"operating cost = ${value(m.fs.operating_cost)/1e6:0.6f} million per year\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For this operating cost, what conversion did we achieve of ethylene oxide to ethylene glycol? " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.R101.report()\n", + "\n", + "print()\n", + "print(f\"Conversion achieved = {value(m.fs.R101.conversion):.1%}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Optimizing Ethylene Glycol Production\n", + "\n", + "Now that the flowsheet has been squared and solved, we can run a small optimization problem to minimize our production costs. Suppose we require at least 200 million pounds/year of ethylene glycol produced and 90% conversion of ethylene oxide, allowing for variable and reactor temperature (heater outlet)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us declare our objective function for this problem. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.objective = Objective(expr=m.fs.operating_cost)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, we need to add the design constraints and unfix the decision variables as we had solved a square problem (degrees of freedom = 0) until now, as well as set bounds for the design variables (reactor outlet temperature is set by state variable bounds in property package):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.eg_prod_con = Constraint(\n", + " expr=m.fs.eg_prod >= 200 * pyunits.Mlb / pyunits.yr\n", + ") # MM lb/year\n", + "m.fs.R101.conversion.fix(0.90)\n", + "\n", + "m.fs.H101.outlet.temperature.unfix()\n", + "m.fs.H101.outlet.temperature[0].setlb(328.15 * pyunits.K)\n", + "m.fs.H101.outlet.temperature[0].setub(\n", + " 470.45 * pyunits.K\n", + ") # highest component boiling point (ethylene glycol)\n", + "\n", + "m.fs.R101.outlet.temperature.unfix()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "We have now defined the optimization problem and we are now ready to solve this problem. \n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "results = solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"operating cost = ${value(m.fs.operating_cost)/1e6:0.6f} million per year\")\n", + "\n", + "print()\n", + "print(\"Heater results\")\n", + "\n", + "m.fs.H101.report()\n", + "\n", + "print()\n", + "print(\"Stoichiometric reactor results\")\n", + "\n", + "m.fs.R101.report()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Display optimal values for the decision variables and design variables:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Optimal Values\")\n", + "print()\n", + "\n", + "print(f\"H101 outlet temperature = {value(m.fs.H101.outlet.temperature[0]):0.6f} K\")\n", + "\n", + "print()\n", + "print(f\"R101 outlet temperature = {value(m.fs.R101.outlet.temperature[0]):0.6f} K\")\n", + "\n", + "print()\n", + "print(f\"Ethylene glycol produced = {value(m.fs.eg_prod):0.6f} MM lb/year\")\n", + "\n", + "print()\n", + "print(f\"Conversion achieved = {value(m.fs.R101.conversion):.1%}\")" + ] + }, + { + "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.18" + } + }, + "nbformat": 4, + "nbformat_minor": 3 } \ No newline at end of file diff --git a/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_test.ipynb b/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_test.ipynb index 4135cb99..9ef0bc1c 100644 --- a/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_test.ipynb +++ b/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_test.ipynb @@ -306,7 +306,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now, let us add expressions to compute the reactor cooling cost (\\\\$/s) assuming a cost of 2.12E-5 \\\\$/kW, and the heating utility cost (\\\\$/s) assuming 2.2E-4 \\\\$/kW. Note that the heat duty is in units of watt (J/s). The total operating cost will be the sum of the two, expressed in \\\\$/year assuming 8000 operating hours per year (~10\\% downtime, which is fairly common for small scale chemical plants):" + "Now, let us add expressions to compute the reactor cooling cost in USD/s assuming a cost of 2.12E-5 USD/kW, and the heating utility cost in USD/s assuming 2.2E-4 USD/kW. Note that the heat duty is in units of watt (J/s). The total operating cost will be the sum of the two, expressed in USD/year assuming 8000 operating hours per year (~10 percent downtime, which is fairly common for small scale chemical plants):" ] }, { diff --git a/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_usr.ipynb b/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_usr.ipynb index 3ce8eb06..8a0e88b0 100644 --- a/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_usr.ipynb +++ b/idaes_examples/notebooks/docs/unit_models/reactors/stoichiometric_reactor_usr.ipynb @@ -306,7 +306,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now, let us add expressions to compute the reactor cooling cost (\\\\$/s) assuming a cost of 2.12E-5 \\\\$/kW, and the heating utility cost (\\\\$/s) assuming 2.2E-4 \\\\$/kW. Note that the heat duty is in units of watt (J/s). The total operating cost will be the sum of the two, expressed in \\\\$/year assuming 8000 operating hours per year (~10\\% downtime, which is fairly common for small scale chemical plants):" + "Now, let us add expressions to compute the reactor cooling cost in USD/s assuming a cost of 2.12E-5 USD/kW, and the heating utility cost in USD/s assuming 2.2E-4 USD/kW. Note that the heat duty is in units of watt (J/s). The total operating cost will be the sum of the two, expressed in USD/year assuming 8000 operating hours per year (~10 percent downtime, which is fairly common for small scale chemical plants):" ] }, {