diff --git a/idaes_examples/notebooks/docs/unit_models/operations/skeleton_unit_doc.ipynb b/idaes_examples/notebooks/docs/unit_models/operations/skeleton_unit_doc.ipynb index 3885b9ea..4d2896d4 100644 --- a/idaes_examples/notebooks/docs/unit_models/operations/skeleton_unit_doc.ipynb +++ b/idaes_examples/notebooks/docs/unit_models/operations/skeleton_unit_doc.ipynb @@ -446,701 +446,647 @@ "metadata": {}, "outputs": [ { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# IDAES Skeleton Unit Model\n", - "Maintainer: Brandon Paul \n", - "Author: Brandon Paul \n", - "Updated: 2023-06-01 \n", - "\n", - "This notebook demonstrates usage of the IDAES Skeleton Unit Model, which provides a generic \"bare bones\" unit for user-defined models and custom variable and constraint sets. To allow maximum versatility, this unit may be defined as a surrogate model or a custom equation-oriented model. Users must add ports and variables that match connected models, and this is facilitated through a provided method to add port-variable sets.\n", - "\n", - "For users who wish to train surrogates with IDAES tools and insert obtained models into a flowsheet, see more detailed information on [IDAES Surrogate Tools](https://idaes-pse.readthedocs.io/en/stable/explanations/modeling_extensions/surrogate/index.html)." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# 1. Motivation\n", - "\n", - "In many cases, a specific application requires a unique unit operation that does not exist in the IDAES repository. Custom user models may source from external scripts, import surrogate equations or use first-principles calculations. However, IDAES flowsheets adhere to a standardized modeling hierarchy and simple Pyomo models do not always follow these conventions. Additionally, simple flowsheet submodels often require integration with other IDAES unit models which requires consistency between corresponding port variables, stream properties and physical unit sets, as well as proper usage of `ControlVolume` blocks.\n", - "\n", - "The IDAES `SkeletonUnitModel` allows custom creation of user models blocks that do not require `ControlVolume` blocks, and enabling connection with standard IDAES unit models that do contain `ControlVolume` blocks. To motivate the usefulness and versatility of this tool, we will consider a simple pervaporation unit. The custom model does not require rigourous thermodynamic calculations contained in adjacent unit models, and using a Skeleton model allows definition of only required variables and constraints. The new block does require state variable connections for the inlet and outlet streams. We will demonstrate this scenario below to highlight the usage and benefits of the Skeleton model." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# 2. Example - Pervaporation\n", - "\n", - "Pervaporation is a low-energy separation process, and is particularly advantageous over distillation for azeotropic solutions or aqueous mixtures of heavy alcohols. Ethylene glycol is more environmentally friendly than typical chloride- and bromide-based dessicants, and is a common choice for commericial recovery of water from flue gas via liquid spray columns. Due to ethylene glycol's high boiling point, diffusion-based water recovery is economically favorable compared to distillation-based processes. The following example and flux correlation are taken from the literature source below:\n", - "\n", - "Jennifer Runhong Du, Amit Chakma, X. Feng, Dehydration of ethylene glycol by pervaporation using poly(N,N-dimethylaminoethyl methacrylate)/polysulfone composite membranes, Separation and Purification Technology, Volume 64, Issue 1, 2008, Pages 63-70, ISSN 1383-5866, https://doi.org/10.1016/j.seppur.2008.08.004.\n", - "\n", - "The process is adapted from the literature, utilizing an inlet aqueous glycol feed circulated through a feed tank-membrane-feed tank recycle loop while permeate is continuously extracted by the membrane. To demonstrate the usefulness of the Skeleton model, we will model this system as a Mixer and custom Pervaporation unit per the diagram below and define the flux as an empirical custom mass balance term rather than requiring rigorous diffusion calculations. We will also circumvent the need for a vapor phase and VLE calculations by manually calculating the duty to condense and collect permeate vapor, and use correlations for steady-state fluxes to avoid a recycle requiring tear calculations.\n", - "\n", - "![](pervaporation_process.png)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 2.1 Pyomo and IDAES Imports\n", - "We will begin with relevant imports. We will need basic Pyomo and IDAES components:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import pytest\n", - "from pyomo.environ import (\n", - " check_optimal_termination,\n", - " ConcreteModel,\n", - " Constraint,\n", - " Expression,\n", - " Objective,\n", - " maximize,\n", - " Var,\n", - " Set,\n", - " TransformationFactory,\n", - " value,\n", - " exp,\n", - " units as pyunits,\n", - ")\n", - "from pyomo.network import Arc\n", - "from idaes.core import FlowsheetBlock\n", - "from idaes.models.unit_models import Feed, SkeletonUnitModel, Mixer, Product\n", - "from idaes.core.util.model_statistics import degrees_of_freedom\n", - "from idaes.core.util.initialization import propagate_state\n", - "from idaes.core.solvers import get_solver\n", - "from pyomo.util.check_units import assert_units_consistent\n", - "\n", - "# import thermophysical properties\n", - "import eg_h2o_ideal as thermo_props\n", - "from idaes.models.properties.modular_properties import GenericParameterBlock\n", - "from idaes.core.util.constants import Constants" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 2.2 Build Flowsheet\n", - "\n", - "We will build a simple model manually defining state variables relations entering and exiting the pervaporation unit. As shown below, we may define our pre-separation mixer as usual:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# build the flowsheet\n", - "m = ConcreteModel()\n", - "m.fs = FlowsheetBlock(dynamic=False)\n", - "\n", - "m.fs.thermo_params = GenericParameterBlock(**thermo_props.config_dict)\n", - "\n", - "m.fs.WATER = Feed(property_package=m.fs.thermo_params)\n", - "m.fs.GLYCOL = Feed(property_package=m.fs.thermo_params)\n", - "\n", - "m.fs.M101 = Mixer(\n", - " property_package=m.fs.thermo_params, inlet_list=[\"water_feed\", \"glycol_feed\"]\n", - ")\n", - "\n", - "m.fs.RETENTATE = Product(property_package=m.fs.thermo_params)\n", - "m.fs.PERMEATE = Product(property_package=m.fs.thermo_params)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 2.2 Defining Skeleton Model and Connections\n", - "\n", - "Now that our flowsheet exists, we can manually define variables, units, constraints and ports for our custom pervaporation unit model. By using a Skeleton model, we avoid rigorous mass and energy balances and phase equilibrium which impact model tractability. Instead, we define state variable relations as below - note that we include the fluxes as outlet flow terms. In this model, the variables specify an `FpcTP` system where molar flow of each component, temperature and pressure are selected as state variables:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# define Skeleton model for pervaporation unit\n", - "m.fs.pervap = SkeletonUnitModel(dynamic=False)\n", - "m.fs.pervap.comp_list = Set(initialize=[\"water\", \"ethylene_glycol\"])\n", - "m.fs.pervap.phase_list = Set(initialize=[\"Liq\"])\n", - "\n", - "# input vars for skeleton\n", - "# m.fs.time is a pre-initialized Set belonging to the FlowsheetBlock; for dynamic=False, time=[0]\n", - "m.fs.pervap.flow_in = Var(\n", - " m.fs.time,\n", - " m.fs.pervap.phase_list,\n", - " m.fs.pervap.comp_list,\n", - " initialize=1.0,\n", - " units=pyunits.mol / pyunits.s,\n", - ")\n", - "m.fs.pervap.temperature_in = Var(m.fs.time, initialize=298.15, units=pyunits.K)\n", - "m.fs.pervap.pressure_in = Var(m.fs.time, initialize=101e3, units=pyunits.Pa)\n", - "\n", - "# output vars for skeleton\n", - "m.fs.pervap.perm_flow = Var(\n", - " m.fs.time,\n", - " m.fs.pervap.phase_list,\n", - " m.fs.pervap.comp_list,\n", - " initialize=1.0,\n", - " units=pyunits.mol / pyunits.s,\n", - ")\n", - "m.fs.pervap.ret_flow = Var(\n", - " m.fs.time,\n", - " m.fs.pervap.phase_list,\n", - " m.fs.pervap.comp_list,\n", - " initialize=1.0,\n", - " units=pyunits.mol / pyunits.s,\n", - ")\n", - "m.fs.pervap.temperature_out = Var(m.fs.time, initialize=298.15, units=pyunits.K)\n", - "m.fs.pervap.pressure_out = Var(m.fs.time, initialize=101e3, units=pyunits.Pa)\n", - "m.fs.pervap.vacuum = Var(m.fs.time, initialize=1.3e3, units=pyunits.Pa)\n", - "\n", - "# dictionaries relating state properties to custom variables\n", - "inlet_dict = {\n", - " \"flow_mol_phase_comp\": m.fs.pervap.flow_in,\n", - " \"temperature\": m.fs.pervap.temperature_in,\n", - " \"pressure\": m.fs.pervap.pressure_in,\n", - "}\n", - "retentate_dict = {\n", - " \"flow_mol_phase_comp\": m.fs.pervap.ret_flow,\n", - " \"temperature\": m.fs.pervap.temperature_out,\n", - " \"pressure\": m.fs.pervap.pressure_out,\n", - "}\n", - "permeate_dict = {\n", - " \"flow_mol_phase_comp\": m.fs.pervap.perm_flow,\n", - " \"temperature\": m.fs.pervap.temperature_out,\n", - " \"pressure\": m.fs.pervap.vacuum,\n", - "}\n", - "\n", - "m.fs.pervap.add_ports(name=\"inlet\", member_dict=inlet_dict)\n", - "m.fs.pervap.add_ports(name=\"retentate\", member_dict=retentate_dict)\n", - "m.fs.pervap.add_ports(name=\"permeate\", member_dict=permeate_dict)\n", - "\n", - "# internal vars for skeleton\n", - "energy_activation_dict = {\n", - " (0, \"Liq\", \"water\"): 51e3,\n", - " (0, \"Liq\", \"ethylene_glycol\"): 53e3,\n", - "}\n", - "m.fs.pervap.energy_activation = Var(\n", - " m.fs.time,\n", - " m.fs.pervap.phase_list,\n", - " m.fs.pervap.comp_list,\n", - " initialize=energy_activation_dict,\n", - " units=pyunits.J / pyunits.mol,\n", - ")\n", - "m.fs.pervap.energy_activation.fix()\n", - "\n", - "permeance_dict = {\n", - " (0, \"Liq\", \"water\"): 5611320,\n", - " (0, \"Liq\", \"ethylene_glycol\"): 22358.88,\n", - "} # calculated from literature data\n", - "m.fs.pervap.permeance = Var(\n", - " m.fs.time,\n", - " m.fs.pervap.phase_list,\n", - " m.fs.pervap.comp_list,\n", - " initialize=permeance_dict,\n", - " units=pyunits.mol / pyunits.s / pyunits.m**2,\n", - ")\n", - "m.fs.pervap.permeance.fix()\n", - "\n", - "m.fs.pervap.area = Var(m.fs.time, initialize=6, units=pyunits.m**2)\n", - "m.fs.pervap.area.fix()\n", - "\n", - "latent_heat_dict = {\n", - " (0, \"Liq\", \"water\"): 40.660e3,\n", - " (0, \"Liq\", \"ethylene_glycol\"): 56.9e3,\n", - "}\n", - "m.fs.pervap.latent_heat_of_vaporization = Var(\n", - " m.fs.time,\n", - " m.fs.pervap.phase_list,\n", - " m.fs.pervap.comp_list,\n", - " initialize=latent_heat_dict,\n", - " units=pyunits.J / pyunits.mol,\n", - ")\n", - "m.fs.pervap.latent_heat_of_vaporization.fix()\n", - "m.fs.pervap.heat_duty = Var(\n", - " m.fs.time, initialize=1, units=pyunits.J / pyunits.s\n", - ") # we will calculate this later" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's define our surrogate equations for flux and permeance, and link them to the port variables. Users can use this structure to write custom relations between inlet and outlet streams; for example, here we define the outlet flow of the pervaporation unit as a sum of the inlet flow and calculated recovery fluxes. By defining model constraints in lieu of rigorous mass balances, we add the flux as a custom mass balance term via an empirical correlation and calculate only the condensation duty rather than implementing full energy balance calculations:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Surrogate and first principles model equations\n", - "\n", - "# flux equation (gas constant is defined as J/mol-K)\n", - "\n", - "\n", - "def rule_permeate_flux(pervap, t, p, i):\n", - " return pervap.permeate.flow_mol_phase_comp[t, p, i] / pervap.area[t] == (\n", - " pervap.permeance[t, p, i]\n", - " * exp(\n", - " -pervap.energy_activation[t, p, i]\n", - " / (Constants.gas_constant * pervap.inlet.temperature[t])\n", - " )\n", - " )\n", - "\n", - "\n", - "m.fs.pervap.eq_permeate_flux = Constraint(\n", - " m.fs.time, m.fs.pervap.phase_list, m.fs.pervap.comp_list, rule=rule_permeate_flux\n", - ")\n", - "\n", - "# permeate condensation equation\n", - "# heat duty based on condensing all of permeate product vapor\n", - "# avoids the need for a Heater or HeatExchanger unit model\n", - "\n", - "\n", - "def rule_duty(pervap, t):\n", - " return pervap.heat_duty[t] == sum(\n", - " pervap.latent_heat_of_vaporization[t, p, i]\n", - " * pervap.permeate.flow_mol_phase_comp[t, p, i]\n", - " for p in pervap.phase_list\n", - " for i in pervap.comp_list\n", - " )\n", - "\n", - "\n", - "m.fs.pervap.eq_duty = Constraint(m.fs.time, rule=rule_duty)\n", - "\n", - "# flow equation adding total recovery as a custom mass balance term\n", - "def rule_retentate_flow(pervap, t, p, i):\n", - " return pervap.retentate.flow_mol_phase_comp[t, p, i] == (\n", - " pervap.inlet.flow_mol_phase_comp[t, p, i]\n", - " - pervap.permeate.flow_mol_phase_comp[t, p, i]\n", - " )\n", - "\n", - "\n", - "m.fs.pervap.eq_retentate_flow = Constraint(\n", - " m.fs.time, m.fs.pervap.phase_list, m.fs.pervap.comp_list, rule=rule_retentate_flow\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, let's define the Arc connecting our two models (IDAES Mixer and custom Pervaporation) and build the flowsheet network:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.s01 = Arc(source=m.fs.WATER.outlet, destination=m.fs.M101.water_feed)\n", - "m.fs.s02 = Arc(source=m.fs.GLYCOL.outlet, destination=m.fs.M101.glycol_feed)\n", - "m.fs.s03 = Arc(source=m.fs.M101.outlet, destination=m.fs.pervap.inlet)\n", - "m.fs.s04 = Arc(source=m.fs.pervap.permeate, destination=m.fs.PERMEATE.inlet)\n", - "m.fs.s05 = Arc(source=m.fs.pervap.retentate, destination=m.fs.RETENTATE.inlet)\n", - "TransformationFactory(\"network.expand_arcs\").apply_to(m)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's see how many degrees of freedom the flowsheet has:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(degrees_of_freedom(m))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 2.3 Inlet Specifications\n", - "\n", - "To obtain a square problem with zero degrees of freedom, we specify the inlet water flow, ethylene glycol flow, temperature and pressure for each feed stream, as well as the permeate stream pressure:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.WATER.outlet.flow_mol_phase_comp[0, \"Liq\", \"water\"].fix(0.34) # mol/s\n", - "m.fs.WATER.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"].fix(1e-6) # mol/s\n", - "m.fs.WATER.outlet.temperature.fix(318.15) # K\n", - "m.fs.WATER.outlet.pressure.fix(101.325e3) # Pa\n", - "\n", - "m.fs.GLYCOL.outlet.flow_mol_phase_comp[0, \"Liq\", \"water\"].fix(1e-6) # mol/s\n", - "m.fs.GLYCOL.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"].fix(0.66) # mol/s\n", - "m.fs.GLYCOL.outlet.temperature.fix(318.15) # K\n", - "m.fs.GLYCOL.outlet.pressure.fix(101.325e3) # Pa" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Additionally, we need to pass rules defining the temperature and pressure outlets of the pervaporation unit:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Add a constraint to calculate the outlet temperature.\n", - "# Here, assume outlet temperature is the same as inlet temperature for illustration\n", - "# in reality, temperature change from latent heat loss through membrane is negligible\n", - "\n", - "\n", - "def rule_temp_out(pervap, t):\n", - " return pervap.inlet.temperature[t] == pervap.retentate.temperature[t]\n", - "\n", - "\n", - "m.fs.pervap.temperature_out_calculation = Constraint(m.fs.time, rule=rule_temp_out)\n", - "\n", - "# Add a constraint to calculate the retentate pressure\n", - "# Here, assume the retentate pressure is the same as the inlet pressure for illustration\n", - "# in reality, pressure change from mass loss through membrane is negligible\n", - "\n", - "\n", - "def rule_pres_out(pervap, t):\n", - " return pervap.inlet.pressure[t] == pervap.retentate.pressure[t]\n", - "\n", - "\n", - "m.fs.pervap.pressure_out_calculation = Constraint(m.fs.time, rule=rule_pres_out)\n", - "\n", - "# fix permeate vacuum pressure\n", - "m.fs.PERMEATE.inlet.pressure.fix(1.3e3)\n", - "\n", - "assert degrees_of_freedom(m) == 0" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 2.4 Custom Initialization\n", - "In addition to allowing custom variable and constraint defintions, the Skeleton model enables implementation of a custom initialization scheme. Complex unit operations may present unique tractability issues, and users have precise control over piecewise unit model solving." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Add this to the imports\n", - "from pyomo.util.calc_var_value import calculate_variable_from_constraint\n", - "\n", - "\n", - "def my_initialize(unit, **kwargs):\n", - " # Callback for user provided initialization sequence\n", - " # Fix the inlet state\n", - " unit.inlet.flow_mol_phase_comp.fix()\n", - " unit.inlet.pressure.fix()\n", - " unit.inlet.temperature.fix()\n", - "\n", - " # Calculate the values of the remaining variables\n", - " for t in m.fs.time:\n", - "\n", - " calculate_variable_from_constraint(\n", - " unit.permeate.flow_mol_phase_comp[t, \"Liq\", \"water\"],\n", - " unit.eq_permeate_flux[t, \"Liq\", \"water\"],\n", - " )\n", - "\n", - " calculate_variable_from_constraint(\n", - " unit.permeate.flow_mol_phase_comp[t, \"Liq\", \"ethylene_glycol\"],\n", - " unit.eq_permeate_flux[t, \"Liq\", \"ethylene_glycol\"],\n", - " )\n", - "\n", - " calculate_variable_from_constraint(unit.heat_duty[t], unit.eq_duty[t])\n", - "\n", - " calculate_variable_from_constraint(\n", - " unit.retentate.flow_mol_phase_comp[t, \"Liq\", \"water\"],\n", - " unit.eq_retentate_flow[t, \"Liq\", \"water\"],\n", - " )\n", - "\n", - " calculate_variable_from_constraint(\n", - " unit.retentate.flow_mol_phase_comp[t, \"Liq\", \"ethylene_glycol\"],\n", - " unit.eq_retentate_flow[t, \"Liq\", \"ethylene_glycol\"],\n", - " )\n", - "\n", - " calculate_variable_from_constraint(\n", - " unit.retentate.temperature[t], unit.temperature_out_calculation[t]\n", - " )\n", - "\n", - " calculate_variable_from_constraint(\n", - " unit.retentate.pressure[t], unit.pressure_out_calculation[t]\n", - " )\n", - "\n", - " assert degrees_of_freedom(unit) == 0\n", - " if degrees_of_freedom(unit) == 0:\n", - " res = solver.solve(unit, tee=True)\n", - " unit.inlet.flow_mol_phase_comp.unfix()\n", - " unit.inlet.temperature.unfix()\n", - " unit.inlet.pressure.unfix()\n", - " print(\"Custom initialization routine complete: \", res.solver.message)\n", - "\n", - "\n", - "solver = get_solver()\n", - "\n", - "m.fs.WATER.initialize()\n", - "propagate_state(m.fs.s01)\n", - "\n", - "m.fs.GLYCOL.initialize()\n", - "propagate_state(m.fs.s02)\n", - "\n", - "m.fs.pervap.config.initializer = my_initialize\n", - "my_initialize(m.fs.pervap)\n", - "propagate_state(m.fs.s03)\n", - "\n", - "m.fs.PERMEATE.initialize()\n", - "propagate_state(m.fs.s04)\n", - "\n", - "m.fs.RETENTATE.initialize()\n", - "\n", - "results = solver.solve(m, tee=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's check the results:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# print results\n", - "\n", - "m.fs.WATER.report()\n", - "m.fs.GLYCOL.report()\n", - "m.fs.PERMEATE.report()\n", - "m.fs.RETENTATE.report()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# separation factor for results analysis\n", - "m.fs.inlet_water_frac = Expression(\n", - " expr=(\n", - " m.fs.pervap.inlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", - " / sum(\n", - " m.fs.pervap.inlet.flow_mol_phase_comp[0, \"Liq\", i]\n", - " for i in m.fs.pervap.comp_list\n", - " )\n", - " )\n", - ")\n", - "m.fs.permeate_water_frac = Expression(\n", - " expr=(\n", - " m.fs.pervap.permeate.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", - " / sum(\n", - " m.fs.pervap.permeate.flow_mol_phase_comp[0, \"Liq\", i]\n", - " for i in m.fs.pervap.comp_list\n", - " )\n", - " )\n", - ")\n", - "m.fs.separation_factor = Expression(\n", - " expr=(m.fs.permeate_water_frac / (1 - m.fs.permeate_water_frac))\n", - " / (m.fs.inlet_water_frac / (1 - m.fs.inlet_water_frac))\n", - ")\n", - "\n", - "print(f\"Inlet water mole fraction: {value(m.fs.inlet_water_frac)}\")\n", - "print(f\"Permeate water mole fraction: {value(m.fs.permeate_water_frac)}\")\n", - "print(f\"Separation factor: {value(m.fs.separation_factor)}\")\n", - "print(f\"Condensation duty: {value(m.fs.pervap.heat_duty[0]/1000)} kW\")\n", - "print(\n", - " f\"Duty per mole water recovered: {value(m.fs.pervap.heat_duty[0]/(1000*m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, 'Liq', 'water']*3600))} kW-h / mol\"\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# check results\n", - "assert check_optimal_termination(results)\n", - "assert_units_consistent(m)\n", - "\n", - "assert value(\n", - " m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", - ") == pytest.approx(0.1426, rel=1e-3)\n", - "assert value(\n", - " m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", - ") == pytest.approx(0.0002667, rel=1e-3)\n", - "assert value(\n", - " m.fs.RETENTATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", - ") == pytest.approx(0.1974, rel=1e-3)\n", - "assert value(\n", - " m.fs.RETENTATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", - ") == pytest.approx(0.6597, rel=1e-3)\n", - "assert value(m.fs.separation_factor) == pytest.approx(1038, rel=1e-3)\n", - "assert value(m.fs.pervap.heat_duty[0]) == pytest.approx(5813, rel=1e-3)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# 3. Optimization\n", - "\n", - "Suppose we wish to characterize the membrane behavior by calculating the maximum inlet water mole fraction allowing a separation factor of at least 100 (typical value for high-efficiency separation processes such as gas separation of CO2/N2). We need to fix total inlet flow to ensure physically-sound solutions. We can quickly modify and resolve the model, and check some key results:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# unfix inlet flows but fix total to prevent divergence during solve\n", - "m.fs.WATER.outlet.flow_mol_phase_comp[0, \"Liq\", \"water\"].unfix()\n", - "m.fs.GLYCOL.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"].unfix()\n", - "m.fs.total_flow = Constraint(\n", - " expr=m.fs.WATER.outlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", - " + m.fs.GLYCOL.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", - " == 1 * pyunits.mol / pyunits.s\n", - ")\n", - "\n", - "# set criteria for separation factor\n", - "m.fs.sep_min = Constraint(expr=m.fs.separation_factor >= 100)\n", - "\n", - "# set objective - defaults to minimization\n", - "m.fs.obj = Objective(expr=m.fs.inlet_water_frac, sense=maximize)\n", - "\n", - "results = solver.solve(m, tee=True)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# print results\n", - "\n", - "m.fs.WATER.report()\n", - "m.fs.GLYCOL.report()\n", - "m.fs.PERMEATE.report()\n", - "m.fs.RETENTATE.report()" - ] - }, + "name": "stdout", + "output_type": "stream", + "text": [ + "2024-05-20 14:15:19 [INFO] idaes.init.fs.WATER.properties: Starting initialization\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.WATER.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.WATER.properties: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.WATER: Initialization Complete.\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.GLYCOL.properties: Starting initialization\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.GLYCOL.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.GLYCOL.properties: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.GLYCOL: Initialization Complete.\n", + "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...: 11\n", + "Number of nonzeros in inequality constraint Jacobian.: 0\n", + "Number of nonzeros in Lagrangian Hessian.............: 0\n", + "\n", + "Total number of variables............................: 7\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 0\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 7\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.13e-16 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + "\n", + "Number of Iterations....: 0\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Constraint violation....: 1.1275702593849246e-16 1.1275702593849246e-16\n", + "Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Overall NLP error.......: 1.1275702593849246e-16 1.1275702593849246e-16\n", + "\n", + "\n", + "Number of objective function evaluations = 1\n", + "Number of objective gradient evaluations = 1\n", + "Number of equality constraint evaluations = 1\n", + "Number of inequality constraint evaluations = 0\n", + "Number of equality constraint Jacobian evaluations = 1\n", + "Number of inequality constraint Jacobian evaluations = 0\n", + "Number of Lagrangian Hessian evaluations = 0\n", + "Total CPU secs in IPOPT (w/o function evaluations) = 0.000\n", + "Total CPU secs in NLP function evaluations = 0.000\n", + "\n", + "EXIT: Optimal Solution Found.\n", + "Custom initialization routine complete: Ipopt 3.13.2\\x3a Optimal Solution Found\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.PERMEATE.properties: Starting initialization\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.PERMEATE.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.PERMEATE.properties: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.PERMEATE: Initialization Complete.\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.RETENTATE.properties: Starting initialization\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.RETENTATE.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.RETENTATE.properties: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:19 [INFO] idaes.init.fs.RETENTATE: Initialization Complete.\n", + "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...: 113\n", + "Number of nonzeros in inequality constraint Jacobian.: 0\n", + "Number of nonzeros in Lagrangian Hessian.............: 74\n", + "\n", + "Total number of variables............................: 47\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 33\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 47\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 7.42e+07 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 0.0000000e+00 7.10e+05 2.45e+01 -1.0 1.01e+05 - 5.55e-01 9.90e-01h 1\n", + " 2 0.0000000e+00 2.30e+05 1.02e+02 -1.0 1.00e+03 - 7.67e-01 6.65e-01h 1\n", + " 3 0.0000000e+00 5.49e+03 6.43e+01 -1.0 1.57e+03 - 9.90e-01 9.90e-01h 1\n", + " 4 0.0000000e+00 7.38e+01 7.18e+02 -1.0 1.92e+03 - 9.90e-01 1.00e+00h 1\n", + " 5 0.0000000e+00 1.39e-03 1.95e+00 -1.0 1.76e+02 - 1.00e+00 1.00e+00h 1\n", + " 6 0.0000000e+00 3.38e-10 5.22e-03 -3.8 2.31e-01 - 1.00e+00 1.00e+00h 1\n", + "\n", + "Number of Iterations....: 6\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Constraint violation....: 1.6906692712481686e-10 3.3813385424963371e-10\n", + "Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Overall NLP error.......: 1.6906692712481686e-10 3.3813385424963371e-10\n", + "\n", + "\n", + "Number of objective function evaluations = 7\n", + "Number of objective gradient evaluations = 7\n", + "Number of equality constraint evaluations = 7\n", + "Number of inequality constraint evaluations = 0\n", + "Number of equality constraint Jacobian evaluations = 7\n", + "Number of inequality constraint Jacobian evaluations = 0\n", + "Number of Lagrangian Hessian evaluations = 6\n", + "Total CPU secs in IPOPT (w/o function evaluations) = 0.001\n", + "Total CPU secs in NLP function evaluations = 0.000\n", + "\n", + "EXIT: Optimal Solution Found.\n" + ] + } + ], + "source": [ + "# Add this to the imports\n", + "from pyomo.util.calc_var_value import calculate_variable_from_constraint\n", + "\n", + "\n", + "def my_initialize(unit, **kwargs):\n", + " # Callback for user provided initialization sequence\n", + " # Fix the inlet state\n", + " unit.inlet.flow_mol_phase_comp.fix()\n", + " unit.inlet.pressure.fix()\n", + " unit.inlet.temperature.fix()\n", + "\n", + " # Calculate the values of the remaining variables\n", + " for t in m.fs.time:\n", + "\n", + " calculate_variable_from_constraint(\n", + " unit.permeate.flow_mol_phase_comp[t, \"Liq\", \"water\"],\n", + " unit.eq_permeate_flux[t, \"Liq\", \"water\"],\n", + " )\n", + "\n", + " calculate_variable_from_constraint(\n", + " unit.permeate.flow_mol_phase_comp[t, \"Liq\", \"ethylene_glycol\"],\n", + " unit.eq_permeate_flux[t, \"Liq\", \"ethylene_glycol\"],\n", + " )\n", + "\n", + " calculate_variable_from_constraint(unit.heat_duty[t], unit.eq_duty[t])\n", + "\n", + " calculate_variable_from_constraint(\n", + " unit.retentate.flow_mol_phase_comp[t, \"Liq\", \"water\"],\n", + " unit.eq_retentate_flow[t, \"Liq\", \"water\"],\n", + " )\n", + "\n", + " calculate_variable_from_constraint(\n", + " unit.retentate.flow_mol_phase_comp[t, \"Liq\", \"ethylene_glycol\"],\n", + " unit.eq_retentate_flow[t, \"Liq\", \"ethylene_glycol\"],\n", + " )\n", + "\n", + " calculate_variable_from_constraint(\n", + " unit.retentate.temperature[t], unit.temperature_out_calculation[t]\n", + " )\n", + "\n", + " calculate_variable_from_constraint(\n", + " unit.retentate.pressure[t], unit.pressure_out_calculation[t]\n", + " )\n", + "\n", + " assert degrees_of_freedom(unit) == 0\n", + " if degrees_of_freedom(unit) == 0:\n", + " res = solver.solve(unit, tee=True)\n", + " unit.inlet.flow_mol_phase_comp.unfix()\n", + " unit.inlet.temperature.unfix()\n", + " unit.inlet.pressure.unfix()\n", + " print(\"Custom initialization routine complete: \", res.solver.message)\n", + "\n", + "\n", + "solver = get_solver()\n", + "\n", + "m.fs.WATER.initialize()\n", + "propagate_state(m.fs.s01)\n", + "\n", + "m.fs.GLYCOL.initialize()\n", + "propagate_state(m.fs.s02)\n", + "\n", + "m.fs.pervap.config.initializer = my_initialize\n", + "my_initialize(m.fs.pervap)\n", + "propagate_state(m.fs.s03)\n", + "\n", + "m.fs.PERMEATE.initialize()\n", + "propagate_state(m.fs.s04)\n", + "\n", + "m.fs.RETENTATE.initialize()\n", + "\n", + "results = solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's check the results:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(f\"Inlet water mole fraction: {value(m.fs.inlet_water_frac)}\")\n", - "print(f\"Permeate water mole fraction: {value(m.fs.permeate_water_frac)}\")\n", - "print(f\"Separation factor: {value(m.fs.separation_factor)}\")\n", - "print(f\"Condensation duty: {value(m.fs.pervap.heat_duty[0]/1000)} kW\")\n", - "print(\n", - " f\"Duty per mole water recovered: {value(m.fs.pervap.heat_duty[0]/(1000*m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, 'Liq', 'water']*3600))} kW-h / mol\"\n", - ")" - ] - }, + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "====================================================================================\n", + "Unit : fs.WATER Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Outlet \n", + " Molar Flowrate ('Liq', 'water') mole / second 0.34000\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 1.0000e-06\n", + " Temperature kelvin 318.15\n", + " Pressure pascal 1.0132e+05\n", + "====================================================================================\n", + "\n", + "====================================================================================\n", + "Unit : fs.GLYCOL Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Outlet \n", + " Molar Flowrate ('Liq', 'water') mole / second 1.0000e-06\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 0.66000\n", + " Temperature kelvin 318.15\n", + " Pressure pascal 1.0132e+05\n", + "====================================================================================\n", + "\n", + "====================================================================================\n", + "Unit : fs.PERMEATE Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Inlet \n", + " Molar Flowrate ('Liq', 'water') mole / second 0.14259\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 0.00026675\n", + " Temperature kelvin 318.15\n", + " Pressure pascal 1300.0\n", + "====================================================================================\n", + "\n", + "====================================================================================\n", + "Unit : fs.RETENTATE Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Inlet \n", + " Molar Flowrate ('Liq', 'water') mole / second 0.19742\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 0.65973\n", + " Temperature kelvin 318.15\n", + " Pressure pascal 1.0132e+05\n", + "====================================================================================\n" + ] + } + ], + "source": [ + "# print results\n", + "\n", + "m.fs.WATER.report()\n", + "m.fs.GLYCOL.report()\n", + "m.fs.PERMEATE.report()\n", + "m.fs.RETENTATE.report()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# check results\n", - "assert check_optimal_termination(results)\n", - "assert_units_consistent(m)\n", - "\n", - "assert value(\n", - " m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", - ") == pytest.approx(0.1426, rel=1e-3)\n", - "assert value(\n", - " m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", - ") == pytest.approx(0.0002667, rel=1e-3)\n", - "assert value(\n", - " m.fs.RETENTATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", - ") == pytest.approx(0.6998, rel=1e-3)\n", - "assert value(\n", - " m.fs.RETENTATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", - ") == pytest.approx(0.1573, rel=1e-3)\n", - "assert value(m.fs.separation_factor) == pytest.approx(100.0, rel=1e-3)\n", - "assert value(m.fs.pervap.heat_duty[0]) == pytest.approx(5813, rel=1e-3)" - ] - }, + "name": "stdout", + "output_type": "stream", + "text": [ + "Inlet water mole fraction: 0.34000031999936\n", + "Permeate water mole fraction: 0.998132696792035\n", + "Separation factor: 1037.6188153503224\n", + "Condensation duty: 5.812711092458081 kW\n", + "Duty per mole water recovered: 0.011324013423286048 kW-h / mol\n" + ] + } + ], + "source": [ + "# separation factor for results analysis\n", + "m.fs.inlet_water_frac = Expression(\n", + " expr=(\n", + " m.fs.pervap.inlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", + " / sum(\n", + " m.fs.pervap.inlet.flow_mol_phase_comp[0, \"Liq\", i]\n", + " for i in m.fs.pervap.comp_list\n", + " )\n", + " )\n", + ")\n", + "m.fs.permeate_water_frac = Expression(\n", + " expr=(\n", + " m.fs.pervap.permeate.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", + " / sum(\n", + " m.fs.pervap.permeate.flow_mol_phase_comp[0, \"Liq\", i]\n", + " for i in m.fs.pervap.comp_list\n", + " )\n", + " )\n", + ")\n", + "m.fs.separation_factor = Expression(\n", + " expr=(m.fs.permeate_water_frac / (1 - m.fs.permeate_water_frac))\n", + " / (m.fs.inlet_water_frac / (1 - m.fs.inlet_water_frac))\n", + ")\n", + "\n", + "print(f\"Inlet water mole fraction: {value(m.fs.inlet_water_frac)}\")\n", + "print(f\"Permeate water mole fraction: {value(m.fs.permeate_water_frac)}\")\n", + "print(f\"Separation factor: {value(m.fs.separation_factor)}\")\n", + "print(f\"Condensation duty: {value(m.fs.pervap.heat_duty[0]/1000)} kW\")\n", + "print(\n", + " f\"Duty per mole water recovered: {value(m.fs.pervap.heat_duty[0]/(1000*m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, 'Liq', 'water']*3600))} kW-h / mol\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "# check results\n", + "assert check_optimal_termination(results)\n", + "assert_units_consistent(m)\n", + "\n", + "assert value(\n", + " m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", + ") == pytest.approx(0.1426, rel=1e-3)\n", + "assert value(\n", + " m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", + ") == pytest.approx(0.0002667, rel=1e-3)\n", + "assert value(\n", + " m.fs.RETENTATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", + ") == pytest.approx(0.1974, rel=1e-3)\n", + "assert value(\n", + " m.fs.RETENTATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", + ") == pytest.approx(0.6597, rel=1e-3)\n", + "assert value(m.fs.separation_factor) == pytest.approx(1038, rel=1e-3)\n", + "assert value(m.fs.pervap.heat_duty[0]) == pytest.approx(5813, rel=1e-3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 3. Optimization\n", + "\n", + "Suppose we wish to characterize the membrane behavior by calculating the maximum inlet water mole fraction allowing a separation factor of at least 100 (typical value for high-efficiency separation processes such as gas separation of CO2/N2). We need to fix total inlet flow to ensure physically-sound solutions. We can quickly modify and resolve the model, and check some key results:" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# 4. Summary\n", - "\n", - "The IDAES Skeleton Unit Model is a powerful tool for implementing relatively simple first-princples, surrogate-based or empirical unit operations. More crucially, users can add their own custom models and integrate them into a larger IDAES flowsheet without adding control volumes or rigorous flow balance and equilibrium calculations when not required. The pervaporation example displays a case where all model equations are empirical correlations or simple manual calculations, with a small number of state variable and port connections, and the Skeleton model avoids complex calculations that impact model tractability. The example also demonstrates adding a custom initialization scheme to handle internally model degrees of freedom, a feature providing greater user control than with most IDAES unit models." - ] - }, + "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...: 121\n", + "Number of nonzeros in inequality constraint Jacobian.: 4\n", + "Number of nonzeros in Lagrangian Hessian.............: 88\n", + "\n", + "Total number of variables............................: 49\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 35\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 48\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 -3.4000032e-01 7.26e+03 6.14e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 -3.7417153e-01 7.06e+02 3.96e+01 -1.0 4.48e-02 - 9.90e-01 8.90e-01h 1\n", + " 2 -5.0289931e-01 2.22e+02 1.14e+02 -1.0 2.30e-01 - 9.90e-01 6.53e-01h 1\n", + " 3 -6.1963465e-01 9.31e+00 1.14e+04 -1.0 1.37e-01 - 9.91e-01 9.93e-01h 1\n", + " 4 -6.1883656e-01 1.48e-02 1.01e+02 -1.0 1.25e-03 - 1.00e+00 9.99e-01h 1\n", + " 5 -6.1153811e-01 2.14e-04 7.97e+02 -1.0 8.51e-03 - 1.00e+00 1.00e+00f 1\n", + " 6 -6.1019822e-01 1.82e-06 1.48e+01 -1.0 1.56e-03 - 1.00e+00 1.00e+00h 1\n", + " 7 -7.3642396e-01 1.07e-02 2.29e+04 -2.5 1.47e-01 - 7.70e-01 1.00e+00f 1\n", + " 8 -8.1765759e-01 2.06e-02 9.69e+02 -2.5 1.47e-01 - 1.00e+00 6.43e-01f 1\n", + " 9 -8.3576879e-01 3.92e-03 4.60e+01 -2.5 2.11e-02 - 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.3954010e-01 1.75e-04 1.05e+00 -2.5 4.40e-03 - 1.00e+00 1.00e+00h 1\n", + " 11 -8.3930724e-01 8.95e-08 2.93e-03 -2.5 2.72e-04 - 1.00e+00 1.00e+00h 1\n", + " 12 -8.4239161e-01 6.13e-05 9.69e+01 -3.8 3.73e-03 - 1.00e+00 9.63e-01f 1\n", + " 13 -8.4225198e-01 7.39e-08 1.37e-02 -3.8 1.63e-04 - 1.00e+00 1.00e+00f 1\n", + " 14 -8.4225232e-01 5.82e-11 3.86e-08 -3.8 3.92e-07 - 1.00e+00 1.00e+00h 1\n", + " 15 -8.4240230e-01 1.51e-07 2.22e-02 -5.7 1.75e-04 - 1.00e+00 1.00e+00f 1\n", + " 16 -8.4240161e-01 9.09e-13 1.67e-08 -5.7 8.12e-07 - 1.00e+00 1.00e+00h 1\n", + " 17 -8.4240336e-01 9.09e-13 3.07e-06 -7.0 2.05e-06 - 1.00e+00 1.00e+00f 1\n", + " 18 -8.4240336e-01 5.82e-11 1.08e-12 -7.0 1.12e-10 - 1.00e+00 1.00e+00h 1\n", + "\n", + "Number of Iterations....: 18\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: -8.4240336294854945e-01 -8.4240336294854945e-01\n", + "Dual infeasibility......: 1.0800249583553523e-12 1.0800249583553523e-12\n", + "Constraint violation....: 1.2738317640843951e-14 5.8207660913467407e-11\n", + "Complementarity.........: 9.0909090917798691e-08 9.0909090917798691e-08\n", + "Overall NLP error.......: 9.0909090917798691e-08 9.0909090917798691e-08\n", + "\n", + "\n", + "Number of objective function evaluations = 19\n", + "Number of objective gradient evaluations = 19\n", + "Number of equality constraint evaluations = 19\n", + "Number of inequality constraint evaluations = 19\n", + "Number of equality constraint Jacobian evaluations = 19\n", + "Number of inequality constraint Jacobian evaluations = 19\n", + "Number of Lagrangian Hessian evaluations = 18\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": [ + "# unfix inlet flows but fix total to prevent divergence during solve\n", + "m.fs.WATER.outlet.flow_mol_phase_comp[0, \"Liq\", \"water\"].unfix()\n", + "m.fs.GLYCOL.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"].unfix()\n", + "m.fs.total_flow = Constraint(\n", + " expr=m.fs.WATER.outlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", + " + m.fs.GLYCOL.outlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", + " == 1 * pyunits.mol / pyunits.s\n", + ")\n", + "\n", + "# set criteria for separation factor\n", + "m.fs.sep_min = Constraint(expr=m.fs.separation_factor >= 100)\n", + "\n", + "# set objective - defaults to minimization\n", + "m.fs.obj = Objective(expr=m.fs.inlet_water_frac, sense=maximize)\n", + "\n", + "results = solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "====================================================================================\n", + "Unit : fs.WATER Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Outlet \n", + " Molar Flowrate ('Liq', 'water') mole / second 0.84240\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 1.0000e-06\n", + " Temperature kelvin 318.15\n", + " Pressure pascal 1.0132e+05\n", + "====================================================================================\n", + "\n", + "====================================================================================\n", + "Unit : fs.GLYCOL Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Outlet \n", + " Molar Flowrate ('Liq', 'water') mole / second 1.0000e-06\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 0.15760\n", + " Temperature kelvin 318.15\n", + " Pressure pascal 1.0132e+05\n", + "====================================================================================\n", + "\n", + "====================================================================================\n", + "Unit : fs.PERMEATE Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Inlet \n", + " Molar Flowrate ('Liq', 'water') mole / second 0.14259\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 0.00026675\n", + " Temperature kelvin 318.15\n", + " Pressure pascal 1300.0\n", + "====================================================================================\n", + "\n", + "====================================================================================\n", + "Unit : fs.RETENTATE Time: 0.0\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Inlet \n", + " Molar Flowrate ('Liq', 'water') mole / second 0.69982\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 0.15733\n", + " Temperature kelvin 318.15\n", + " Pressure pascal 1.0132e+05\n", + "====================================================================================\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 results\n", + "\n", + "m.fs.WATER.report()\n", + "m.fs.GLYCOL.report()\n", + "m.fs.PERMEATE.report()\n", + "m.fs.RETENTATE.report()" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Inlet water mole fraction: 0.8424033629485495\n", + "Permeate water mole fraction: 0.9981326967920352\n", + "Separation factor: 100.00006747653238\n", + "Condensation duty: 5.812711092447902 kW\n", + "Duty per mole water recovered: 0.011324013423286044 kW-h / mol\n" + ] } + ], + "source": [ + "print(f\"Inlet water mole fraction: {value(m.fs.inlet_water_frac)}\")\n", + "print(f\"Permeate water mole fraction: {value(m.fs.permeate_water_frac)}\")\n", + "print(f\"Separation factor: {value(m.fs.separation_factor)}\")\n", + "print(f\"Condensation duty: {value(m.fs.pervap.heat_duty[0]/1000)} kW\")\n", + "print(\n", + " f\"Duty per mole water recovered: {value(m.fs.pervap.heat_duty[0]/(1000*m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, 'Liq', 'water']*3600))} kW-h / mol\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "# check results\n", + "assert check_optimal_termination(results)\n", + "assert_units_consistent(m)\n", + "\n", + "assert value(\n", + " m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", + ") == pytest.approx(0.1426, rel=1e-3)\n", + "assert value(\n", + " m.fs.PERMEATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", + ") == pytest.approx(0.0002667, rel=1e-3)\n", + "assert value(\n", + " m.fs.RETENTATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"water\"]\n", + ") == pytest.approx(0.6998, rel=1e-3)\n", + "assert value(\n", + " m.fs.RETENTATE.inlet.flow_mol_phase_comp[0, \"Liq\", \"ethylene_glycol\"]\n", + ") == pytest.approx(0.1573, rel=1e-3)\n", + "assert value(m.fs.separation_factor) == pytest.approx(100.0, rel=1e-3)\n", + "assert value(m.fs.pervap.heat_duty[0]) == pytest.approx(5813, rel=1e-3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 4. Summary\n", + "\n", + "The IDAES Skeleton Unit Model is a powerful tool for implementing relatively simple first-princples, surrogate-based or empirical unit operations. More crucially, users can add their own custom models and integrate them into a larger IDAES flowsheet without adding control volumes or rigorous flow balance and equilibrium calculations when not required. The pervaporation example displays a case where all model equations are empirical correlations or simple manual calculations, with a small number of state variable and port connections, and the Skeleton model avoids complex calculations that impact model tractability. The example also demonstrates adding a custom initialization scheme to handle internally model degrees of freedom, a feature providing greater user control than with most IDAES unit models." + ] + }, + { + "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": { @@ -1152,9 +1098,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 } diff --git a/idaes_examples/notebooks/docs/unit_models/reactors/plug_flow_reactor_doc.ipynb b/idaes_examples/notebooks/docs/unit_models/reactors/plug_flow_reactor_doc.ipynb index 9a8a70a0..139fbbe7 100644 --- a/idaes_examples/notebooks/docs/unit_models/reactors/plug_flow_reactor_doc.ipynb +++ b/idaes_examples/notebooks/docs/unit_models/reactors/plug_flow_reactor_doc.ipynb @@ -49,7 +49,7 @@ "\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 plug flow reactor (PFR). 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 a CSTR which assumes well-mixed liquid behavior, the concentration profiles will vary spatially in one dimension. In actuality, following start-up flow reactor exhibit dynamic behavior as they approach a steady-state equilibrium; we will assume our system has already achieved steady-state behavior. 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", + "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.ipynb), we can alter the flowsheet to use a plug flow reactor (PFR). 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 a CSTR which assumes well-mixed liquid behavior, the concentration profiles will vary spatially in one dimension. In actuality, following start-up flow reactor exhibit dynamic behavior as they approach a steady-state equilibrium; we will assume our system has already achieved steady-state behavior. 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", @@ -182,7 +182,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:" ] }, { @@ -232,9 +232,7 @@ { "cell_type": "code", "execution_count": 7, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "m.fs.R101 = PFR(\n", @@ -380,772 +378,999 @@ }, "outputs": [ { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "# Flowsheet Plug Flow Reactor (PFR) Simulation and Optimization of Ethylene Glycol Production\n", - "Author: Andrew Lee \n", - "Maintainer: Andrew Lee \n", - "\n", - "\n", - "## Learning Outcomes\n", - "\n", - "\n", - "- Call and implement the IDAES PFR unit model\n", - "- Construct a steady-state flowsheet using the IDAES unit model library\n", - "- Connecting unit models in a flowsheet using Arcs\n", - "- Fomulate 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 plug flow reactor (PFR). 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 a CSTR which assumes well-mixed liquid behavior, the concentration profiles will vary spatially in one dimension. In actuality, following start-up flow reactor exhibit dynamic behavior as they approach a steady-state equilibrium; we will assume our system has already achieved steady-state behavior. 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)" - ] - }, - { - "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", - "- PFR\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 import (\n", - " GenericParameterBlock,\n", - " GenericReactionParameterBlock,\n", - ")\n", - "from idaes.models.unit_models import Feed, Mixer, Heater, PFR, 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", - "The reaction package here assumes Arrhenius kinetic behavior for the PFR, for which $k_0$ and $E_a$ are known *a priori* (if unknown, they may be obtained using one of the parameter estimation tools within IDAES).\n", - "\n", - "$ r = -kVC_{EO} $, $ k = k_0 e^{(-E_a/RT)}$, with the variables as follows:\n", - "\n", - "$r$ - reaction rate extent in moles of ethylene oxide consumed per second; note that the traditional reaction rate would be given by $rate = r/V$ in moles per $m^3$ per second \n", - "$k$ - reaction rate constant per second \n", - "$V$ - volume of PFR in $m^3$, note that this is *liquid volume* and not the *total volume* of the reactor itself \n", - "$C_{EO}$ - bulk concentration of ethylene oxide in moles per $m^3$ (the limiting reagent, since we assume excess catalyst and water) \n", - "$k_0$ - pre-exponential Arrhenius factor per second \n", - "$E_a$ - reaction activation energy in kJ per mole of ethylene oxide consumed \n", - "$R$ - gas constant in J/mol-K \n", - "$T$ - reactor temperature in K\n", - "\n", - "These calculations are contained within the property, reaction and unit model packages, and do not need to be entered into the flowsheet. More information on property estimation may be found in the IDAES documentation on [Parameter Estimation](https://idaes-pse.readthedocs.io/en/stable/how_to_guides/workflow/data_rec_parmest.html).\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 `PFR`. 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 = PFR(\n", - " property_package=m.fs.thermo_params,\n", - " reaction_package=m.fs.reaction_params,\n", - " has_equilibrium_reactions=False,\n", - " has_heat_of_reaction=True,\n", - " has_heat_transfer=True,\n", - " has_pressure_change=False,\n", - " transformation_method=\"dae.finite_difference\",\n", - " transformation_scheme=\"BACKWARD\",\n", - " finite_elements=20,\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 specifed 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 `PFR`. 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 documentaiton]( 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. To calculate cooling cost, it is important to note that the heat duty is not constant throughout the reactor's length and is expressed in terms of heat per length (J/m/s). This is why we utilize the trapezoid rule to calculate the total heat duty of the reactor:$Q=\\Delta x\\big(\\sum_{k=1}^{N-1}(Q_k)+\\frac{Q_N+Q_0}{2}\\big)$ \n", - "where k is the subinterval in the length domain, N is the number of intervals, and $\\Delta x$ is the length of the interval.\n", - "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": { - "scrolled": true - }, - "outputs": [], - "source": [ - "m.fs.cooling_cost = Expression(\n", - " expr=2.12e-8\n", - " * m.fs.R101.length\n", - " / m.fs.R101.config.finite_elements\n", - " * (\n", - " -sum(\n", - " m.fs.R101.heat_duty[0, k]\n", - " for k in m.fs.R101.control_volume.length_domain\n", - " if 0.0 <= k < 1.0\n", - " )\n", - " )\n", - " - (\n", - " value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(1)])\n", - " - value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(-1)])\n", - " )\n", - " / 2\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 2 unit specifications and 1 specification for each finite element. Therefore, we have 35 degrees of freedom to specify: temperature, pressure and flow of all four components on both streams; outlet heater temperature; a reactor property such as conversion or heat duty at each finite element; reactor volume and reactor length." - ] - }, - { - "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": [ - "For the `PFR`, we have to define the conversion in terms of ethylene oxide. Note that the `PFR` reaction volume variable (m.fs.R101.volume) does not need to be defined here since it is internally defined by the `PFR` model. We'll estimate 50% conversion for our initial flowsheet:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.R101.conversion = Var(\n", - " bounds=(0, 1), initialize=0.80, 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", - "for x in m.fs.R101.control_volume.length_domain:\n", - " if x == 0:\n", - " continue\n", - " m.fs.R101.control_volume.properties[0, x].temperature.fix(\n", - " 328.15 * pyunits.K\n", - " ) # equal inlet reactor temperature\n", - "\n", - "m.fs.R101.conversion.fix(0.5)\n", - "\n", - "m.fs.R101.length.fix(1 * pyunits.m)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As we did not place a specification on reactor duty, the solver may try positive values to increase the reaction temperature and rate. To prevent the optimization from diverging, we need to set an upper bound restricting heat flow to cooling only:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m.fs.R101.heat_duty.setub(\n", - " 0 * pyunits.J / pyunits.m / pyunits.s\n", - ") # heat duty is only used for cooling" - ] - }, - { - "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.3f} 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%}\")\n", - "print()\n", - "print(\n", - " f\"Total heat duty required = \"\n", - " f\"\"\"{(value(m.fs.R101.length) / value(m.fs.R101.config.finite_elements) * \n", - " (value(sum(m.fs.R101.heat_duty[0, k] for k in m.fs.R101.control_volume.length_domain if 0.0 <= k < 1.0))\n", - " + (value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(1)])\n", - " + value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(-1)]))/2))/1e6:0.3f}\"\"\"\n", - " f\" MJ\"\n", - ")\n", - "print()\n", - "print(f\"Tube area required = {value(m.fs.R101.area):0.3f} m^2\")\n", - "print()\n", - "print(f\"Tube length required = {value(m.fs.R101.length):0.3f} m\")\n", - "print()\n", - "print(f\"Tube volume required = {value(m.fs.R101.volume):0.3f} m^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": 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:" - ] - }, + "name": "stdout", + "output_type": "stream", + "text": [ + "35\n" + ] + } + ], + "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) == 35" + ] + }, + { + "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": [ + "For the `PFR`, we have to define the conversion in terms of ethylene oxide. Note that the `PFR` reaction volume variable (m.fs.R101.volume) does not need to be defined here since it is internally defined by the `PFR` model. We'll estimate 50% conversion for our initial flowsheet:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.R101.conversion = Var(\n", + " bounds=(0, 1), initialize=0.80, 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", + "for x in m.fs.R101.control_volume.length_domain:\n", + " if x == 0:\n", + " continue\n", + " m.fs.R101.control_volume.properties[0, x].temperature.fix(\n", + " 328.15 * pyunits.K\n", + " ) # equal inlet reactor temperature\n", + "\n", + "m.fs.R101.conversion.fix(0.5)\n", + "\n", + "m.fs.R101.length.fix(1 * pyunits.m)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As we did not place a specification on reactor duty, the solver may try positive values to increase the reaction temperature and rate. To prevent the optimization from diverging, we need to set an upper bound restricting heat flow to cooling only:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "m.fs.R101.heat_duty.setub(\n", + " 0 * pyunits.J / pyunits.m / pyunits.s\n", + ") # heat duty is only used for cooling" + ] + }, + { + "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": 18, + "metadata": {}, + "outputs": [ { - "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.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.R101.length.unfix()\n", - "m.fs.R101.length.setlb(0 * pyunits.m)\n", - "m.fs.R101.length.setub(5 * pyunits.m)\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", - "for x in m.fs.R101.control_volume.length_domain:\n", - " if x == 0:\n", - " continue\n", - " m.fs.R101.control_volume.properties[\n", - " 0, x\n", - " ].temperature.unfix() # allow for temperature change in each finite element" - ] - }, + "name": "stdout", + "output_type": "stream", + "text": [ + "0\n" + ] + } + ], + "source": [ + "print(degrees_of_freedom(m))" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "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": 20, + "metadata": {}, + "outputs": [ { - "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" - ] - }, + "name": "stdout", + "output_type": "stream", + "text": [ + "2024-05-20 14:15:05 [INFO] idaes.init.fs.OXIDE.properties: Starting initialization\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.OXIDE.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.OXIDE.properties: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.OXIDE: Initialization Complete.\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.ACID.properties: Starting initialization\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.ACID.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.ACID.properties: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.ACID: Initialization Complete.\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.M101.reagent_feed_state: Starting initialization\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.M101.reagent_feed_state: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.M101.catalyst_feed_state: Starting initialization\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.M101.catalyst_feed_state: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:05 [INFO] idaes.init.fs.M101.mixed_state: Starting initialization\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.M101.mixed_state: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.M101.mixed_state: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.M101: Initialization Complete: optimal - Optimal Solution Found\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.H101.control_volume.properties_in: Starting initialization\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.H101.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.H101.control_volume.properties_out: Starting initialization\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.H101.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.H101.control_volume: Initialization Complete\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.H101: Initialization Complete: optimal - Optimal Solution Found\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.R101.control_volume.properties: Starting initialization\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.R101.control_volume.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.R101.control_volume.reactions: Initialization Complete.\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.R101.control_volume: Initialization Complete\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.R101: Initialization Complete: optimal - Optimal Solution Found\n", + "2024-05-20 14:15:06 [INFO] idaes.init.fs.PROD.properties: Starting initialization\n", + "2024-05-20 14:15:07 [INFO] idaes.init.fs.PROD.properties: Property initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:07 [INFO] idaes.init.fs.PROD.properties: Property package initialization: optimal - Optimal Solution Found.\n", + "2024-05-20 14:15:07 [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": 21, + "metadata": { + "scrolled": true + }, + "outputs": [ { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "results = solver.solve(m, tee=True)" - ] - }, + "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...: 1923\n", + "Number of nonzeros in inequality constraint Jacobian.: 0\n", + "Number of nonzeros in Lagrangian Hessian.............: 1323\n", + "\n", + "Total number of variables............................: 608\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 257\n", + " variables with only upper bounds: 20\n", + "Total number of equality constraints.................: 608\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 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 0.0000000e+00 1.55e+07 2.12e+05 -1.0 1.11e+08 - 2.85e-01 9.90e-01f 1\n", + " 2 0.0000000e+00 2.43e+05 5.35e+03 -1.0 1.11e+06 - 8.25e-01 9.90e-01h 1\n", + " 3 0.0000000e+00 2.41e+03 5.68e+01 -1.0 1.11e+04 - 9.90e-01 9.90e-01h 1\n", + " 4 0.0000000e+00 1.83e+01 3.21e+03 -1.0 1.10e+02 - 9.90e-01 9.92e-01h 1\n", + " 5 0.0000000e+00 4.61e-07 4.87e+03 -1.0 8.35e-01 - 9.91e-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.6686898166535613e+06 1.6686898166535613e+06\n", + "Constraint violation....: 4.6147033572196960e-07 4.6147033572196960e-07\n", + "Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Overall NLP error.......: 4.6147033572196960e-07 1.6686898166535613e+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.008\n", + "Total CPU secs in NLP function evaluations = 0.001\n", + "\n", + "EXIT: Optimal Solution Found.\n" + ] + } + ], + "source": [ + "# Solve the model\n", + "results = solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "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": 23, + "metadata": {}, + "outputs": [ { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "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(\"PFR reactor results\")\n", - "\n", - "m.fs.R101.report()" - ] - }, + "name": "stdout", + "output_type": "stream", + "text": [ + "operating cost = $6.589 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": 24, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "import pytest\n", + "\n", + "assert value(m.fs.operating_cost) / 1e6 == pytest.approx(6.589, 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": 25, + "metadata": {}, + "outputs": [ { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Display optimal values for the decision variables and design variables:" - ] - }, + "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", + " Area : 0.98560 : meter ** 2 : False : (None, None)\n", + "\n", + "------------------------------------------------------------------------------------\n", + " Stream Table\n", + " Units Inlet Outlet \n", + " Molar Flowrate ('Liq', 'ethylene_oxide') mole / second 58.000 29.000\n", + " Molar Flowrate ('Liq', 'water') mole / second 239.60 210.60\n", + " Molar Flowrate ('Liq', 'sulfuric_acid') mole / second 0.33401 0.33401\n", + " Molar Flowrate ('Liq', 'ethylene_glycol') mole / second 2.0000e-05 29.000\n", + " Temperature kelvin 328.15 328.15\n", + " Pressure pascal 1.0000e+05 1.0000e+05\n", + "====================================================================================\n", + "\n", + "Conversion achieved = 50.0%\n", + "\n", + "Total heat duty required = -3.617 MJ\n", + "\n", + "Tube area required = 0.986 m^2\n", + "\n", + "Tube length required = 1.000 m\n", + "\n", + "Tube volume required = 0.986 m^3\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\"Total heat duty required = \"\n", + " f\"\"\"{(value(m.fs.R101.length) / value(m.fs.R101.config.finite_elements) * \n", + " (value(sum(m.fs.R101.heat_duty[0, k] for k in m.fs.R101.control_volume.length_domain if 0.0 <= k < 1.0))\n", + " + (value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(1)])\n", + " + value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(-1)]))/2))/1e6:0.3f}\"\"\"\n", + " f\" MJ\"\n", + ")\n", + "print()\n", + "print(f\"Tube area required = {value(m.fs.R101.area):0.3f} m^2\")\n", + "print()\n", + "print(f\"Tube length required = {value(m.fs.R101.length):0.3f} m\")\n", + "print()\n", + "print(f\"Tube volume required = {value(m.fs.R101.volume):0.3f} m^3\")" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "assert value(m.fs.R101.conversion) == pytest.approx(0.5000, rel=1e-3)\n", + "assert value(m.fs.R101.area) == pytest.approx(0.9856, rel=1e-3)\n", + "assert (\n", + " value(m.fs.R101.length)\n", + " / value(m.fs.R101.config.finite_elements)\n", + " * value(\n", + " sum(\n", + " m.fs.R101.heat_duty[0, k]\n", + " for k in m.fs.R101.control_volume.length_domain\n", + " if 0.0 <= k < 1.0\n", + " )\n", + " )\n", + " + (\n", + " value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(1)])\n", + " + value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(-1)])\n", + " )\n", + " / 2\n", + ") / 1e6 == pytest.approx(-4.882, rel=1e-3)\n", + "assert value(m.fs.R101.outlet.temperature[0]) / 1e2 == pytest.approx(3.2815, 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": 27, + "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": 28, + "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.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.R101.length.unfix()\n", + "m.fs.R101.length.setlb(0 * pyunits.m)\n", + "m.fs.R101.length.setub(5 * pyunits.m)\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", + "for x in m.fs.R101.control_volume.length_domain:\n", + " if x == 0:\n", + " continue\n", + " m.fs.R101.control_volume.properties[\n", + " 0, x\n", + " ].temperature.unfix() # allow for temperature change in each finite element" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "assert degrees_of_freedom(m) == 22 # 2 unit variables and 20 finite element variables" + ] + }, + { + "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": 30, + "metadata": { + "scrolled": true + }, + "outputs": [ { - "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.3f} K\")\n", - "\n", - "print()\n", - "print(\n", - " \"Total heat duty required = \",\n", - " f\"\"\"{(value(m.fs.R101.length) / value(m.fs.R101.config.finite_elements) * (value(sum(m.fs.R101.heat_duty[0, k] for k in m.fs.R101.control_volume.length_domain if 0.0 <= k < 1.0))\n", - " + (value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(1)])\n", - " + value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(-1)]))/2))/1e6:0.3f}\"\"\"\n", - " f\" MJ\",\n", - ")\n", - "print()\n", - "print(f\"Tube area required = {value(m.fs.R101.area):0.3f} m^2\")\n", - "\n", - "print()\n", - "print(f\"Tube length required = {value(m.fs.R101.length):0.3f} m\")\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.3f}\"\n", - " f\" m^3 = {value(pyunits.convert(1.2*m.fs.R101.volume, 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%}\")" - ] - }, + "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...: 2067\n", + "Number of nonzeros in inequality constraint Jacobian.: 1\n", + "Number of nonzeros in Lagrangian Hessian.............: 1886\n", + "\n", + "Total number of variables............................: 631\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 280\n", + " variables with only upper bounds: 21\n", + "Total number of equality constraints.................: 608\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 6.5890113e+06 3.50e+06 1.00e+02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 6.7062650e+06 2.91e+06 8.32e+01 -1.0 2.08e+06 - 7.82e-02 1.70e-01h 1\n", + " 2 8.4753804e+06 8.06e+04 2.56e+02 -1.0 2.00e+06 - 3.84e-01 9.90e-01H 1\n", + " 3 8.5041953e+06 3.12e+03 7.64e+01 -1.0 3.91e+05 - 9.34e-01 9.91e-01h 1\n", + " 4 8.3381636e+06 4.86e+04 1.44e+04 -1.0 2.06e+07 - 1.90e-01 7.05e-02f 4\n", + " 5 8.3413021e+06 5.73e+00 1.76e+06 -1.0 1.24e+04 -4.0 4.88e-01 1.00e+00h 1\n", + " 6 8.3407522e+06 1.49e-02 1.60e+04 -1.0 1.37e+03 -4.5 9.90e-01 1.00e+00h 1\n", + " 7 8.3407580e+06 4.14e-07 1.58e+02 -1.0 2.07e+02 -5.0 9.90e-01 1.00e+00f 1\n", + " 8 8.1632997e+06 5.70e+04 1.51e+02 -1.0 5.98e+07 - 4.19e-02 2.31e-02f 3\n", + " 9 8.0620878e+06 8.17e+04 1.51e+02 -1.0 6.78e+08 - 4.18e-03 8.82e-04f 3\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 8.0790088e+06 1.82e+02 1.89e+05 -1.0 2.04e+04 -5.4 3.89e-01 1.00e+00h 1\n", + " 11 8.0777533e+06 7.66e-02 1.62e+03 -1.0 2.09e+03 -5.9 1.00e+00 1.00e+00h 1\n", + " 12 8.0775471e+06 5.88e-05 7.07e-01 -1.0 4.87e+02 -6.4 1.00e+00 1.00e+00f 1\n", + " 13 8.0767398e+06 9.74e-04 1.03e-02 -3.8 1.48e+02 -6.9 1.00e+00 1.00e+00f 1\n", + " 14 8.0742991e+06 8.88e-03 1.52e-03 -3.8 4.46e+02 -7.3 1.00e+00 1.00e+00f 1\n", + " 15 8.0669748e+06 8.00e-02 1.78e-04 -3.8 1.34e+03 -7.8 1.00e+00 1.00e+00f 1\n", + " 16 8.0633841e+06 8.60e-02 2.07e+03 -5.7 3.99e+03 -8.3 1.00e+00 1.65e-01f 1\n", + " 17 8.0632366e+06 5.83e-02 1.09e+03 -5.7 1.14e+03 -8.8 1.00e+00 4.74e-01h 1\n", + " 18 8.0624636e+06 5.24e-01 5.43e-02 -5.7 3.19e+03 -9.2 1.00e+00 1.00e+00f 1\n", + " 19 8.0601485e+06 4.71e+00 1.34e-03 -5.7 9.56e+03 -9.7 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.0532064e+06 4.24e+01 1.20e-02 -5.7 2.87e+04 -10.2 1.00e+00 1.00e+00f 1\n", + " 21 8.0324020e+06 3.82e+02 1.04e-01 -5.7 8.62e+04 -10.7 1.00e+00 1.00e+00f 1\n", + " 22 7.9701802e+06 3.45e+03 8.34e-01 -5.7 2.60e+05 -11.2 1.00e+00 1.00e+00f 1\n", + " 23 7.8922689e+06 7.47e+03 1.67e+00 -5.7 7.91e+05 -11.6 1.00e+00 4.19e-01f 1\n", + " 24 7.8922688e+06 7.47e+03 1.67e+00 -5.7 6.79e+05 -12.1 1.00e+00 1.46e-06h 2\n", + " 25 7.7985022e+06 9.96e+04 2.01e+01 -5.7 2.01e+06 -12.6 1.00e+00 1.00e+00f 1\n", + " 26 7.7060388e+06 1.87e+05 1.99e+01 -5.7 1.76e+07 -13.1 1.00e+00 7.81e-02f 1\n", + " 27 7.7060362e+06 1.87e+05 1.99e+01 -5.7 2.96e+06 -12.6 1.00e+00 2.16e-05h 1\n", + " 28 7.6629029e+06 1.88e+05 1.98e+01 -5.7 2.50e+07 -13.1 6.85e-01 2.99e-02f 1\n", + " 29 7.6320707e+06 1.29e+05 1.45e+01 -5.7 2.92e+06 -12.7 1.00e+00 3.07e-01h 1\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 30 7.6318689e+06 1.29e+05 1.44e+01 -5.7 1.56e+07 -13.2 1.00e+00 4.21e-04h 1\n", + " 31 7.5995354e+06 9.16e+04 9.56e+00 -5.7 2.81e+06 -12.7 1.00e+00 3.75e-01h 1\n", + " 32 7.5987331e+06 9.14e+04 9.54e+00 -5.7 1.49e+07 -13.2 1.00e+00 2.11e-03h 1\n", + " 33 7.5662181e+06 7.45e+04 5.70e+00 -5.7 2.73e+06 -12.8 1.00e+00 4.35e-01h 1\n", + " 34 7.5652516e+06 7.43e+04 5.68e+00 -5.7 1.39e+07 -13.3 1.00e+00 3.28e-03h 1\n", + " 35 7.5328904e+06 7.18e+04 2.96e+00 -5.7 2.65e+06 -12.8 1.00e+00 5.07e-01h 1\n", + " 36 7.5323971e+06 7.16e+04 2.96e+00 -5.7 1.25e+07 -13.3 1.00e+00 2.25e-03h 1\n", + " 37 7.4997231e+06 7.89e+04 1.21e+00 -5.7 2.54e+06 -12.9 1.00e+00 6.13e-01h 1\n", + " 38 7.4995708e+06 7.88e+04 1.21e+00 -5.7 1.09e+07 -13.4 1.00e+00 9.50e-04h 1\n", + " 39 7.4663333e+06 9.24e+04 1.24e+00 -5.7 2.41e+06 -12.9 1.00e+00 7.51e-01h 1\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 40 7.4662826e+06 9.24e+04 1.24e+00 -5.7 9.38e+06 -13.4 1.00e+00 4.25e-04h 1\n", + " 41 7.4314438e+06 1.12e+05 1.34e+00 -5.7 2.26e+06 -13.0 1.00e+00 9.30e-01h 1\n", + " 42 7.4314076e+06 1.12e+05 1.34e+00 -5.7 8.56e+06 -13.5 1.00e+00 3.73e-04h 1\n", + " 43 7.3989205e+06 1.15e+05 1.16e+00 -5.7 2.45e+06 -13.0 1.00e+00 1.00e+00h 1\n", + " 44 7.3930244e+06 1.17e+05 1.16e+00 -5.7 1.58e+07 -13.5 1.00e+00 2.51e-02h 1\n", + " 45 7.3589148e+06 1.48e+05 1.35e+00 -5.7 3.32e+06 -13.1 1.00e+00 1.00e+00h 1\n", + " 46 7.3243165e+06 2.28e+05 2.15e+00 -5.7 2.53e+07 -13.6 1.00e+00 1.03e-01h 1\n", + " 47 7.2636674e+06 3.65e+05 6.79e+00 -5.7 5.90e+06 -13.2 1.00e+00 1.00e+00h 1\n", + " 48 7.2517688e+06 7.53e+04 1.17e+00 -5.7 2.61e+06 -12.7 1.00e+00 9.04e-01h 1\n", + " 49 7.1085647e+06 5.44e+05 1.14e+01 -5.7 1.22e+07 -13.2 1.00e+00 6.21e-01f 1\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 50 7.0618984e+06 1.18e+05 1.10e+00 -5.7 4.67e+06 -12.8 1.00e+00 1.00e+00h 1\n", + " 51 6.4252245e+06 1.53e+06 9.14e+00 -5.7 2.47e+07 -13.3 1.00e+00 6.86e-01f 1\n", + " 52 6.4252097e+06 1.53e+06 9.14e+00 -5.7 3.23e+08 -13.7 4.17e-02 1.69e-06h 1\n", + " 53 5.7942297e+06 1.66e+05 8.10e+02 -5.7 3.24e+07 -13.3 8.57e-05 1.00e+00f 1\n", + " 54 5.7889255e+06 1.65e+05 8.09e+02 -5.7 8.47e+07 -13.8 1.21e-01 1.54e-03f 1\n", + " 55 5.0757844e+06 2.33e+05 7.26e+02 -5.7 2.54e+08 -14.3 2.12e-01 7.06e-02f 1\n", + " 56 4.9573062e+06 2.34e+05 7.23e+02 -5.7 7.32e+08 -14.7 1.60e-01 3.97e-03f 1\n", + " 57 4.9548257e+06 2.34e+05 7.23e+02 -5.7 3.95e+09 - 2.25e-02 1.62e-05f 1\n", + " 58 4.8544637e+06 2.35e+05 7.22e+02 -5.7 2.32e+09 - 4.85e-02 1.15e-03f 1\n", + " 59 4.7849161e+06 2.34e+05 7.21e+02 -5.7 1.53e+09 - 1.58e-02 1.26e-03f 1\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 60 4.5951938e+06 2.34e+05 7.18e+02 -5.7 1.17e+09 - 1.05e-02 4.72e-03f 1\n", + " 61 4.5927372e+06 2.34e+05 7.17e+02 -5.7 2.27e+08 -15.2 4.75e-03 3.54e-04f 1\n", + " 62 4.5520698e+06 2.33e+05 7.12e+02 -5.7 2.65e+08 - 8.03e-03 7.23e-03f 1\n", + " 63 4.5378438e+06 2.32e+05 7.10e+02 -5.7 2.29e+08 - 6.07e-03 2.98e-03f 1\n", + " 64 4.5086029e+06 2.30e+05 7.04e+02 -5.7 1.69e+08 - 7.78e-02 8.92e-03f 1\n", + " 65 4.4857030e+06 2.30e+05 7.03e+02 -5.7 1.31e+09 -14.8 6.05e-09 6.24e-04f 1\n", + " 66 4.4777329e+06 2.29e+05 6.99e+02 -5.7 4.81e+07 -13.5 2.28e-02 5.33e-03f 1\n", + " 67 4.4610632e+06 2.21e+05 6.75e+02 -5.7 1.30e+07 -13.0 2.06e-01 3.47e-02h 1\n", + " 68 4.4489784e+06 2.20e+05 6.72e+02 -5.7 9.26e+07 -13.5 1.17e-02 4.85e-03f 1\n", + " 69 4.4487175e+06 2.19e+05 6.70e+02 -5.7 1.42e+06 -12.2 7.76e-01 2.48e-03h 1\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 70 4.4443745e+06 2.14e+05 6.54e+02 -5.7 3.43e+06 -12.7 6.70e-01 2.33e-02h 1\n", + " 71 4.4231622e+06 2.01e+05 6.14e+02 -5.7 1.09e+07 -13.1 1.64e-01 6.19e-02h 1\n", + " 72 4.4231616e+06 2.01e+05 6.14e+02 -5.7 1.47e+06 -11.8 9.27e-02 1.29e-04h 1\n", + " 73 4.4225287e+06 1.55e+05 4.74e+02 -5.7 2.40e+06 - 9.20e-01 2.31e-01h 1\n", + " 74 4.4220119e+06 1.57e+05 2.39e+02 -5.7 2.14e+06 - 2.55e-05 5.03e-01h 1\n", + " 75 4.4215706e+06 5.20e+04 7.52e+01 -5.7 1.06e+06 -12.3 2.24e-01 6.64e-01h 1\n", + " 76 4.4215368e+06 3.93e+04 4.73e+00 -5.7 9.27e+05 -12.8 1.65e-01 1.00e+00h 1\n", + " 77 4.4215312e+06 1.77e+02 5.05e+00 -5.7 1.32e+05 -13.2 8.81e-01 1.00e+00h 1\n", + " 78 4.4215310e+06 7.99e-01 1.09e-02 -5.7 7.99e+03 -13.7 1.00e+00 1.00e+00h 1\n", + " 79 4.4215310e+06 1.36e+04 9.80e-03 -5.7 3.39e+06 - 1.00e+00 1.08e-01h 2\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 80 4.4215310e+06 9.90e+03 9.32e-04 -5.7 6.52e+05 - 1.00e+00 1.00e+00h 1\n", + " 81 4.4215310e+06 4.34e+00 2.21e-05 -5.7 3.57e+04 - 1.00e+00 1.00e+00h 1\n", + " 82 4.4215310e+06 2.75e+00 6.26e-08 -5.7 2.67e+04 - 1.00e+00 1.00e+00h 1\n", + " 83 4.4215310e+06 1.10e-03 1.84e-11 -5.7 4.99e+02 - 1.00e+00 1.00e+00h 1\n", + " 84 4.4215302e+06 3.05e-07 3.19e-10 -8.6 2.39e+01 - 1.00e+00 1.00e+00h 1\n", + "\n", + "Number of Iterations....: 84\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 2.0399399058705549e+02 4.4215302056112560e+06\n", + "Dual infeasibility......: 3.1883161994853613e-10 6.9106135629265704e-06\n", + "Constraint violation....: 3.0547380447387695e-07 3.0547380447387695e-07\n", + "Complementarity.........: 2.5107403576252549e-09 5.4419810592164126e-05\n", + "Overall NLP error.......: 3.0547380447387695e-07 5.4419810592164126e-05\n", + "\n", + "\n", + "Number of objective function evaluations = 103\n", + "Number of objective gradient evaluations = 85\n", + "Number of equality constraint evaluations = 103\n", + "Number of inequality constraint evaluations = 103\n", + "Number of equality constraint Jacobian evaluations = 85\n", + "Number of inequality constraint Jacobian evaluations = 85\n", + "Number of Lagrangian Hessian evaluations = 84\n", + "Total CPU secs in IPOPT (w/o function evaluations) = 0.234\n", + "Total CPU secs in NLP function evaluations = 0.020\n", + "\n", + "EXIT: Optimal Solution Found.\n" + ] + } + ], + "source": [ + "results = solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "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": 32, + "metadata": {}, + "outputs": [ { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] + "name": "stdout", + "output_type": "stream", + "text": [ + "operating cost = $4.422 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", + "PFR 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", + " Area : 2.9260 : meter ** 2 : False : (None, None)\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 286.11\n", + " Pressure pascal 1.0000e+05 1.0000e+05\n", + "====================================================================================\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(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(\"PFR reactor results\")\n", + "\n", + "m.fs.R101.report()" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "assert value(m.fs.operating_cost) / 1e6 == pytest.approx(4.422, rel=1e-3)\n", + "assert value(m.fs.R101.area) == pytest.approx(2.9260, rel=1e-3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Display optimal values for the decision variables and design variables:" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Optimal Values\n", + "\n", + "H101 outlet temperature = 328.150 K\n", + "\n", + "Total heat duty required = -3.790 MJ\n", + "\n", + "Tube area required = 2.926 m^2\n", + "\n", + "Tube length required = 4.982 m\n", + "\n", + "Assuming a 20% design factor for reactor volume,total CSTR volume required = 17.494 m^3 = 4621.519 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", + " \"Total heat duty required = \",\n", + " f\"\"\"{(value(m.fs.R101.length) / value(m.fs.R101.config.finite_elements) * (value(sum(m.fs.R101.heat_duty[0, k] for k in m.fs.R101.control_volume.length_domain if 0.0 <= k < 1.0))\n", + " + (value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(1)])\n", + " + value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(-1)]))/2))/1e6:0.3f}\"\"\"\n", + " f\" MJ\",\n", + ")\n", + "print()\n", + "print(f\"Tube area required = {value(m.fs.R101.area):0.3f} m^2\")\n", + "\n", + "print()\n", + "print(f\"Tube length required = {value(m.fs.R101.length):0.3f} m\")\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.3f}\"\n", + " f\" m^3 = {value(pyunits.convert(1.2*m.fs.R101.volume, 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": 35, + "metadata": { + "tags": [ + "testing" + ] + }, + "outputs": [], + "source": [ + "assert value(m.fs.H101.outlet.temperature[0]) / 100 == pytest.approx(3.2815, rel=1e-3)\n", + "assert (\n", + " value(m.fs.R101.length)\n", + " / value(m.fs.R101.config.finite_elements)\n", + " * (\n", + " value(\n", + " sum(\n", + " m.fs.R101.heat_duty[0, k]\n", + " for k in m.fs.R101.control_volume.length_domain\n", + " if 0.0 <= k < 1.0\n", + " )\n", + " )\n", + " + (\n", + " value(m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(1)])\n", + " + value(\n", + " m.fs.R101.heat_duty[0, m.fs.R101.control_volume.length_domain.at(-1)]\n", + " )\n", + " )\n", + " / 2\n", + " )\n", + ") / 1e6 == pytest.approx(-3.7892, rel=1e-3)\n", + "assert value(m.fs.R101.area) == pytest.approx(2.926, rel=1e-3)\n", + "assert value(m.fs.R101.control_volume.length) == pytest.approx(4.9788, rel=1e-3)\n", + "assert value(m.fs.R101.volume * 1.2) == pytest.approx(17.494, 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.000, rel=1e-3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "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": { @@ -1157,9 +1382,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 }