diff --git a/idaes_examples/notebooks/docs/unit_models/reactors/cstr_doc.ipynb b/idaes_examples/notebooks/docs/unit_models/reactors/cstr_doc.ipynb index 79bd9890..4e52cd16 100644 --- a/idaes_examples/notebooks/docs/unit_models/reactors/cstr_doc.ipynb +++ b/idaes_examples/notebooks/docs/unit_models/reactors/cstr_doc.ipynb @@ -187,7 +187,7 @@ "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:" + "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.ipynb), 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:" ] }, { @@ -237,9 +237,7 @@ { "cell_type": "code", "execution_count": 7, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "m.fs.R101 = CSTR(\n", @@ -369,26 +367,803 @@ "15\n" ] } - ], - "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" + ], + "source": [ + "print(degrees_of_freedom(m))" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "# Check the degrees of freedom\n", + "assert degrees_of_freedom(m) == 15" + ] + }, + { + "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": 14, + "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": 15, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.H101.outlet.temperature.fix(328.15 * pyunits.K)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We'll add constraints defining the reactor volume and conversion in relation to the stream properties. Particularly, we want to use our CSTR performance relation: \n", + "\n", + "$V = \\frac{v_0 X} {k(1-X)}$, where the `CSTR` reaction volume $V$ will be specified, the inlet volumetric flow $v_0$ is determined by stream properties, $k$ is calculated by the reaction package, and $X$ will be calculated. Reactor volume is commonly selected as a specification in simulation problems, and choosing conversion is often to perform reactor design.\n", + "\n", + "For the `CSTR`, we have to define the conversion in terms of ethylene oxide as well as the `CSTR` reaction volume. This requires us to create new variables and constraints relating reactor properties to stream properties. Note that the `CSTR` reaction volume variable (m.fs.R101.volume) does not need to be defined here since it is internally defined by the `CSTR` model. Additionally, the heat duty is not fixed, since the heat of reaction depends on the reactor conversion (through the extent of reaction and heat of reaction). We'll estimate 80% conversion for our initial flowsheet:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "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.volume.fix(5.538 * pyunits.m**3)" + ] + }, + { + "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": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0\n" + ] } + ], + "source": [ + "print(degrees_of_freedom(m))" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "# Check the degrees of freedom\n", + "assert degrees_of_freedom(m) == 0" + ] + }, + { + "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": 19, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2024-05-22 14:30:32 [INFO] idaes.init.fs.OXIDE.properties: Starting initialization\n", + "2024-05-22 14:30:32 [INFO] idaes.init.fs.OXIDE.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:32 [INFO] idaes.init.fs.OXIDE.properties: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:32 [INFO] idaes.init.fs.OXIDE: Initialization Complete.\n", + "2024-05-22 14:30:32 [INFO] idaes.init.fs.ACID.properties: Starting initialization\n", + "2024-05-22 14:30:32 [INFO] idaes.init.fs.ACID.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:32 [INFO] idaes.init.fs.ACID.properties: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:32 [INFO] idaes.init.fs.ACID: Initialization Complete.\n", + "2024-05-22 14:30:32 [INFO] idaes.init.fs.M101.reagent_feed_state: Starting initialization\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.M101.reagent_feed_state: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.M101.catalyst_feed_state: Starting initialization\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.M101.catalyst_feed_state: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.M101.mixed_state: Starting initialization\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.M101.mixed_state: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.M101.mixed_state: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.M101: Initialization Complete: optimal - Optimal Solution Found\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.H101.control_volume.properties_in: Starting initialization\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.H101.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.H101.control_volume.properties_out: Starting initialization\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.H101.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.H101.control_volume: Initialization Complete\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.H101: Initialization Complete: optimal - Optimal Solution Found\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.R101.control_volume.properties_in: Starting initialization\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.R101.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:33 [INFO] idaes.init.fs.R101.control_volume.properties_out: Starting initialization\n", + "2024-05-22 14:30:34 [INFO] idaes.init.fs.R101.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:34 [INFO] idaes.init.fs.R101.control_volume.reactions: Initialization Complete.\n", + "2024-05-22 14:30:34 [INFO] idaes.init.fs.R101.control_volume: Initialization Complete\n", + "2024-05-22 14:30:34 [INFO] idaes.init.fs.R101: Initialization Complete: optimal - Optimal Solution Found\n", + "2024-05-22 14:30:34 [INFO] idaes.init.fs.PROD.properties: Starting initialization\n", + "2024-05-22 14:30:34 [INFO] idaes.init.fs.PROD.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:34 [INFO] idaes.init.fs.PROD.properties: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-22 14:30:34 [INFO] idaes.init.fs.PROD: Initialization Complete.\n" + ] + } + ], + "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": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt 3.13.2: nlp_scaling_method=gradient-based\n", + "tol=1e-06\n", + "max_iter=200\n", + "\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", + "This version of Ipopt was compiled from source code available at\n", + " https://github.com/IDAES/Ipopt as part of the Institute for the Design of\n", + " Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE\n", + " Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.\n", + "\n", + "This version of Ipopt was compiled using HSL, a collection of Fortran codes\n", + " for large-scale scientific computation. All technical papers, sales and\n", + " publicity material resulting from use of the HSL codes within IPOPT must\n", + " contain the following acknowledgement:\n", + " HSL, a collection of Fortran codes for large-scale scientific\n", + " computation. See http://www.hsl.rl.ac.uk.\n", + "******************************************************************************\n", + "\n", + "This is Ipopt version 3.13.2, running with linear solver ma27.\n", + "\n", + "Number of nonzeros in equality constraint Jacobian...: 345\n", + "Number of nonzeros in inequality constraint Jacobian.: 0\n", + "Number of nonzeros in Lagrangian Hessian.............: 393\n", + "\n", + "Total number of variables............................: 96\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 87\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 96\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 0.0000000e+00 1.24e+06 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 0.0000000e+00 2.99e+06 1.69e+01 -1.0 1.09e+07 - 6.77e-02 9.90e-01h 1\n", + " 2 0.0000000e+00 2.51e+04 2.90e+02 -1.0 6.71e+04 - 7.00e-01 9.90e-01h 1\n", + " 3 0.0000000e+00 2.56e+02 1.41e+03 -1.0 2.21e+04 - 9.75e-01 9.90e-01h 1\n", + " 4 0.0000000e+00 1.93e+00 3.22e+03 -1.0 1.36e+03 - 9.90e-01 9.92e-01h 1\n", + " 5 0.0000000e+00 2.08e-06 4.65e+03 -1.0 1.28e+01 - 9.92e-01 1.00e+00h 1\n", + "Cannot recompute multipliers for feasibility problem. Error in eq_mult_calculator\n", + "\n", + "Number of Iterations....: 5\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Dual infeasibility......: 1.6686895402311948e+06 1.6686895402311948e+06\n", + "Constraint violation....: 1.4285171770845562e-09 2.0787119865417480e-06\n", + "Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Overall NLP error.......: 1.4285171770845562e-09 1.6686895402311948e+06\n", + "\n", + "\n", + "Number of objective function evaluations = 6\n", + "Number of objective gradient evaluations = 6\n", + "Number of equality constraint evaluations = 6\n", + "Number of inequality constraint evaluations = 0\n", + "Number of equality constraint Jacobian evaluations = 6\n", + "Number of inequality constraint Jacobian evaluations = 0\n", + "Number of Lagrangian Hessian evaluations = 5\n", + "Total CPU secs in IPOPT (w/o function evaluations) = 0.003\n", + "Total CPU secs in NLP function evaluations = 0.000\n", + "\n", + "EXIT: Optimal Solution Found.\n" + ] + } + ], + "source": [ + "# Solve the model\n", + "results = solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "# Check solver solve status\n", + "from pyomo.environ import TerminationCondition\n", + "\n", + "assert results.solver.termination_condition == TerminationCondition.optimal" + ] + }, + { + "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": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "operating cost = $8.004 million per year\n" + ] + } + ], + "source": [ + "print(f\"operating cost = ${value(m.fs.operating_cost)/1e6:0.3f} million per year\")" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "import pytest\n", + "\n", + "assert value(m.fs.operating_cost) / 1e6 == pytest.approx(8.003020, rel=1e-3)" + ] + }, + { + "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": 24, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "====================================================================================\n", + "Unit : fs.R101 Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Unit Performance\n", + "\n", + " Variables: \n", + "\n", + " Key : Value : Units : Fixed : Bounds\n", + " Heat Duty : -5.8667e+06 : watt : False : (None, None)\n", + " Volume : 5.5380 : meter ** 3 : True : (None, None)\n", + "\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Inlet Outlet \n", + " Molar Flowrate ('Liq', 'ethylene_oxide') mole / second 58.000 11.600\n", + " Molar Flowrate ('Liq', 'water') mole / second 239.60 193.20\n", + " Molar Flowrate ('Liq', 'sulfuric_acid') mole / second 0.33401 0.33401\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 2.0000e-05 46.400\n", + " Temperature kelvin 328.15 329.30\n", + " Pressure pascal 1.0000e+05 1.0000e+05\n", + "====================================================================================\n", + "\n", + "Conversion achieved = 80.0%\n", + "\n", + "Assuming a 20% design factor for reactor volume,total CSTR volume required = 6.646 m^3 = 1755.582 gal\n" + ] + } + ], + "source": [ + "m.fs.R101.report()\n", + "\n", + "print()\n", + "print(f\"Conversion achieved = {value(m.fs.R101.conversion):.1%}\")\n", + "print()\n", + "print(\n", + " f\"Assuming a 20% design factor for reactor volume,\"\n", + " f\"total CSTR volume required = {value(1.2*m.fs.R101.volume[0]):0.3f}\"\n", + " f\" m^3 = {value(pyunits.convert(1.2*m.fs.R101.volume[0], to_units=pyunits.gal)):0.3f} gal\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "assert value(m.fs.R101.conversion) == pytest.approx(0.8000, rel=1e-3)\n", + "assert value(m.fs.R101.volume[0]) == pytest.approx(5.5380, rel=1e-3)\n", + "assert value(m.fs.R101.heat_duty[0]) / 1e6 == pytest.approx(-5.8659, rel=1e-3)\n", + "assert value(m.fs.R101.outlet.temperature[0]) / 1e2 == pytest.approx(3.2933, rel=1e-3)" + ] + }, + { + "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 reactor volume (considering operating/non-capital costs only) and reactor temperature (heater outlet)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us declare our objective function for this problem. " + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "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:" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "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.R101.volume.unfix()\n", + "m.fs.R101.volume.setlb(0 * pyunits.m**3)\n", + "m.fs.R101.volume.setub(pyunits.convert(5000 * pyunits.gal, to_units=pyunits.m**3))\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)" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "assert degrees_of_freedom(m) == 2" + ] + }, + { + "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": 29, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt 3.13.2: nlp_scaling_method=gradient-based\n", + "tol=1e-06\n", + "max_iter=200\n", + "\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", + "This version of Ipopt was compiled from source code available at\n", + " https://github.com/IDAES/Ipopt as part of the Institute for the Design of\n", + " Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE\n", + " Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.\n", + "\n", + "This version of Ipopt was compiled using HSL, a collection of Fortran codes\n", + " for large-scale scientific computation. All technical papers, sales and\n", + " publicity material resulting from use of the HSL codes within IPOPT must\n", + " contain the following acknowledgement:\n", + " HSL, a collection of Fortran codes for large-scale scientific\n", + " computation. See http://www.hsl.rl.ac.uk.\n", + "******************************************************************************\n", + "\n", + "This is Ipopt version 3.13.2, running with linear solver ma27.\n", + "\n", + "Number of nonzeros in equality constraint Jacobian...: 348\n", + "Number of nonzeros in inequality constraint Jacobian.: 1\n", + "Number of nonzeros in Lagrangian Hessian.............: 408\n", + "\n", + "Total number of variables............................: 98\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 89\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 96\n", + "Total number of inequality constraints...............: 1\n", + " inequality constraints with only lower bounds: 1\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 8.0035262e+06 1.74e+06 6.34e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 7.7522414e+06 1.74e+06 1.58e+01 -1.0 1.01e+10 - 1.05e-04 6.26e-06f 1\n", + " 2 7.3310760e+06 1.74e+06 8.39e+01 -1.0 6.50e+09 - 7.69e-05 1.06e-04f 1\n", + " 3 7.1953539e+06 1.74e+06 1.13e+02 -1.0 2.15e+09 - 2.87e-04 1.04e-04f 1\n", + " 4 7.1947130e+06 1.74e+06 1.30e+04 -1.0 3.23e+06 - 7.84e-02 3.61e-04f 1\n", + " 5 7.2867360e+06 1.61e+06 6.37e+04 -1.0 1.56e+06 - 1.30e-01 7.88e-02h 1\n", + " 6 8.3361464e+06 3.53e+04 1.70e+07 -1.0 1.41e+06 - 7.28e-01 9.90e-01h 1\n", + " 7 8.3189323e+06 3.60e+02 7.39e+05 -1.0 3.17e+04 - 9.86e-01 9.91e-01f 1\n", + " 8 8.3176450e+06 8.64e-02 4.26e+03 -1.0 2.14e+03 - 9.90e-01 1.00e+00f 1\n", + " 9 8.3176416e+06 5.00e-07 2.10e+03 -1.7 5.44e+00 - 9.91e-01 1.00e+00f 1\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 8.3176416e+06 1.49e-08 2.00e-07 -1.7 4.22e-05 - 1.00e+00 1.00e+00h 1\n", + " 11 8.3176415e+06 1.12e-08 3.69e-07 -5.7 3.56e-02 - 1.00e+00 1.00e+00f 1\n", + "\n", + "Number of Iterations....: 11\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 8.3176415170773193e+06 8.3176415170773193e+06\n", + "Dual infeasibility......: 3.6871028551502261e-07 3.6871028551502261e-07\n", + "Constraint violation....: 2.0987168619865631e-14 1.1175870895385742e-08\n", + "Complementarity.........: 1.8475671339363312e-06 1.8475671339363312e-06\n", + "Overall NLP error.......: 1.8402097991736475e-07 1.8475671339363312e-06\n", + "\n", + "\n", + "Number of objective function evaluations = 12\n", + "Number of objective gradient evaluations = 12\n", + "Number of equality constraint evaluations = 12\n", + "Number of inequality constraint evaluations = 12\n", + "Number of equality constraint Jacobian evaluations = 12\n", + "Number of inequality constraint Jacobian evaluations = 12\n", + "Number of Lagrangian Hessian evaluations = 11\n", + "Total CPU secs in IPOPT (w/o function evaluations) = 0.004\n", + "Total CPU secs in NLP function evaluations = 0.000\n", + "\n", + "EXIT: Optimal Solution Found.\n" + ] + } + ], + "source": [ + "results = solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "# Check for solver solve status\n", + "from pyomo.environ import TerminationCondition\n", + "\n", + "assert results.solver.termination_condition == TerminationCondition.optimal" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "operating cost = $8.318 million per year\n", + "\n", + "Heater results\n", + "\n", + "====================================================================================\n", + "Unit : fs.H101 Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Unit Performance\n", + "\n", + " Variables: \n", + "\n", + " Key : Value : Units : Fixed : Bounds\n", + " Heat Duty : 6.9784e+05 : watt : False : (None, None)\n", + "\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Inlet Outlet \n", + " Molar Flowrate ('Liq', 'ethylene_oxide') mole / second 58.000 58.000\n", + " Molar Flowrate ('Liq', 'water') mole / second 239.60 239.60\n", + " Molar Flowrate ('Liq', 'sulfuric_acid') mole / second 0.33401 0.33401\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 2.0000e-05 2.0000e-05\n", + " Temperature kelvin 298.15 328.15\n", + " Pressure pascal 1.0000e+05 1.0000e+05\n", + "====================================================================================\n", + "\n", + "CSTR reactor results\n", + "\n", + "====================================================================================\n", + "Unit : fs.R101 Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Unit Performance\n", + "\n", + " Variables: \n", + "\n", + " Key : Value : Units : Fixed : Bounds\n", + " Heat Duty : -6.3812e+06 : watt : False : (None, None)\n", + " Volume : 18.927 : meter ** 3 : False : (0.0, 18.927058919999997)\n", + "\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Inlet Outlet \n", + " Molar Flowrate ('Liq', 'ethylene_oxide') mole / second 58.000 5.8000\n", + " Molar Flowrate ('Liq', 'water') mole / second 239.60 187.40\n", + " Molar Flowrate ('Liq', 'sulfuric_acid') mole / second 0.33401 0.33401\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 2.0000e-05 52.200\n", + " Temperature kelvin 328.15 338.94\n", + " Pressure pascal 1.0000e+05 1.0000e+05\n", + "====================================================================================\n" + ] + } + ], + "source": [ + "print(f\"operating cost = ${value(m.fs.operating_cost)/1e6:0.3f} million per year\")\n", + "\n", + "print()\n", + "print(\"Heater results\")\n", + "\n", + "m.fs.H101.report()\n", + "\n", + "print()\n", + "print(\"CSTR reactor results\")\n", + "\n", + "m.fs.R101.report()" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "assert value(m.fs.operating_cost) / 1e6 == pytest.approx(8.317406, rel=1e-3)\n", + "assert value(m.fs.R101.volume[0]) == pytest.approx(18.927, rel=1e-3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Display optimal values for the decision variables and design variables:" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Optimal Values\n", + "\n", + "H101 outlet temperature = 328.150 K\n", + "\n", + "Assuming a 20% design factor for reactor volume,total CSTR volume required = 22.712 m^3 = 6000.000 gal\n", + "\n", + "Ethylene glycol produced = 225.415 MM lb/year\n", + "\n", + "Conversion achieved = 90.0%\n" + ] + } + ], + "source": [ + "print(\"Optimal Values\")\n", + "print()\n", + "\n", + "print(f\"H101 outlet temperature = {value(m.fs.H101.outlet.temperature[0]):0.3f} K\")\n", + "\n", + "print()\n", + "print(\n", + " f\"Assuming a 20% design factor for reactor volume,\"\n", + " f\"total CSTR volume required = {value(1.2*m.fs.R101.volume[0]):0.3f}\"\n", + " f\" m^3 = {value(pyunits.convert(1.2*m.fs.R101.volume[0], to_units=pyunits.gal)):0.3f} gal\"\n", + ")\n", + "\n", + "print()\n", + "print(f\"Ethylene glycol produced = {value(m.fs.eg_prod):0.3f} MM lb/year\")\n", + "\n", + "print()\n", + "print(f\"Conversion achieved = {value(m.fs.R101.conversion):.1%}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "assert value(m.fs.H101.outlet.temperature[0]) / 100 == pytest.approx(3.2815, rel=1e-3)\n", + "assert value(m.fs.R101.volume[0] * 1.2) == pytest.approx(22.712, rel=1e-3)\n", + "assert value(m.fs.eg_prod) == pytest.approx(225.415, rel=1e-3)\n", + "assert value(m.fs.R101.conversion) * 100 == pytest.approx(90.0, rel=1e-3)" + ] + }, + { + "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": { @@ -400,9 +1175,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.9.18" } }, "nbformat": 4, - "nbformat_minor": 3 + "nbformat_minor": 4 }