diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/__init__.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/example.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/example.py new file mode 100644 index 00000000..6aae5b6f --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/example.py @@ -0,0 +1,247 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the +# software owners: The Regents of the University of California, through +# Lawrence Berkeley National Laboratory, National Technology & Engineering +# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia +# University Research Corporation, et al. All rights reserved. +# +# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +"""Example of debugging a structural singularity using the IDAES diagnostics +toolbox. This script reproduces the functionality in the structural_singularity.ipynb +notebook, but breaks the example into helper functions that are easier to test. + +As a script, this module runs the example end-to-end. + +""" +from idaes_examples.mod.diagnostics.gas_solid_contactors.model import make_model +from idaes_examples.mod.diagnostics.util import get_subsystem_at_time +from idaes.core.util.model_statistics import degrees_of_freedom, large_residuals_set +from idaes.core.util.model_diagnostics import DiagnosticsToolbox +import pyomo.environ as pyo +from pyomo.core.expr import replace_expressions +from pyomo.util.subsystems import TemporarySubsystemManager +from pyomo.contrib.incidence_analysis import IncidenceGraphInterface +import logging + + +def check_dof_and_residuals(model): + dof = degrees_of_freedom(model) + has_large_residuals = bool(large_residuals_set(model, tol=1e-5)) + print(f"Degrees of freedom: {dof}") + print(f"Has large residuals: {has_large_residuals}") + return dof, has_large_residuals + + +def attempt_solve(model): + solver = pyo.SolverFactory("ipopt") + solver.options["max_iter"] = 20 + solver.options["print_user_options"] = "yes" + solver.options["OF_print_info_string"] = "yes" + res = solver.solve(model, tee=True) + return res + + +def fix_degrees_of_freedom(model): + model.fs.MB.gas_phase.properties[:, 0].flow_mol.fix() + model.fs.MB.solid_phase.properties[:, 1].flow_mass.fix() + model.piecewise_constant_constraints.deactivate() + + +def free_degrees_of_freedom(model): + model.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix() + model.fs.MB.gas_phase.properties[0, 0].flow_mol.fix() + model.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix() + model.fs.MB.solid_phase.properties[0, 1].flow_mass.fix() + model.piecewise_constant_constraints.activate() + + +def get_subsystem_at_t0(model): + t0 = model.fs.time.first() + t_block, inputs = get_subsystem_at_time(model, model.fs.time, t0) + return t_block, inputs + + +def add_particle_porosity_variable(model): + model.fs.MB.particle_porosity = pyo.Var( + model.fs.time, + model.fs.MB.length_domain, + initialize=model.fs.solid_properties.particle_porosity.value, + ) + + +def display_constraints_containing_variable(model, var): + igraph = IncidenceGraphInterface(model, include_fixed=True) + print(f"Constraints containing {var.name}:") + for con in igraph.get_adjacent_to(var): + print(f" {con.name}") + + +def replace_porosity_parameter_with_variable(model): + porosity_param = model.fs.solid_properties.particle_porosity + for t, x in model.fs.time * model.fs.MB.length_domain: + substitution_map = {id(porosity_param): model.fs.MB.particle_porosity[t, x]} + sp = model.fs.MB.solid_phase + cons = [ + sp.properties[t, x].density_particle_constraint, + sp.reactions[t, x].gen_rate_expression["R1"], + ] + for con in cons: + con.set_value( + replace_expressions( + con.expr, + substitution_map, + descend_into_named_expressions=True, + ) + ) + + +def add_density_flowrate_constraint(model): + @model.fs.MB.Constraint(model.fs.time, model.fs.MB.length_domain) + def density_flowrate_constraint(mb, t, x): + return ( + mb.velocity_superficial_solid[t] * mb.bed_area + * mb.solid_phase.properties[t, x].dens_mass_particle + == mb.solid_phase.properties[t, x].flow_mass + ) + + +def main(): + model = make_model() + # Before trying to solve the model, let's make sure it conforms to our + # expectations. I.e. it (a) has degrees of freedom and (b) is initialized to + # a feasible point. + check_dof_and_residuals(model) + # Looks good so far, let's try to solve! + attempt_solve(model) + + # Let's run the diagnostics toolbox on the model and see what it has to say + fix_degrees_of_freedom(model) + dt = DiagnosticsToolbox(model) + + # Before calling report_structural_issues, we'll effectively disable Pyomo + # logging messages. This is not recommended in general, but we do it here + # to suppress unit inconsistency errors that otherwise flood our screen. + # This model has unit inconsistency errors as it was created in IDAES 1.7, + # before we enforced that models use units. + logging.getLogger("pyomo").setLevel(logging.CRITICAL) + + # Now we can finally see what the diagnostics toolbox has to say + dt.report_structural_issues() + # We got the following warnings: + # - Inconsistent units + # - Structural singularity + # - Potential evaluation errors + # We'll ignore inconsistent units and potential evaluation errors, and focus on + # the structural singularity. + dt.display_underconstrained_set() + dt.display_overconstrained_set() + + # Suppose the above doesn't give us any leads. We'll try to break the problem + # down into subsystems at each point in time. These should individually be + # nonsingular. + t_block, inputs = get_subsystem_at_t0(model) + with TemporarySubsystemManager(to_fix=inputs): + dt = DiagnosticsToolbox(t_block) + dt.report_structural_issues() + dt.display_underconstrained_set() + dt.display_overconstrained_set() + # The overconstrained system decomposes into smaller independent blocks, which + # are easier to debug. + + # After some thought, we decide we need to make particle porosity a variable, + # and add an equation linking flow rate and density. We'll make these changes + # on a fresh copy of the model. + model2 = make_model() + fix_degrees_of_freedom(model2) + + # Add a new particle porosity variable + add_particle_porosity_variable(model2) + # Display the constraints containing our old porosity "parameter" + porosity_param = model2.fs.solid_properties.particle_porosity + display_constraints_containing_variable(model2, porosity_param) + # Replace the old porosity parameter with the new porosity variable + replace_porosity_parameter_with_variable(model2) + # Add density-flow rate constraint + add_density_flowrate_constraint(model2) + + # Re-check structural diagnostics + dt = DiagnosticsToolbox(model2) + dt.report_structural_issues() + + # The structural singularity appears to be gone. Let's try to solve. + free_degrees_of_freedom(model2) + attempt_solve(model2) + + # This doesn't look any better. Let's check for numerical issues. + fix_degrees_of_freedom(model2) + dt.report_numerical_issues() + + # We seem to have nearly parallel constraints. Let's see what they are. + dt.display_near_parallel_constraints() + + # What is this "solid_super_vel"? + model2.fs.MB.solid_super_vel[0].pprint() + + # This is the constraint we just added. Looks like it was already defined at + # the solid inlet. We'll just deactivate the new constraint here. + model2.fs.MB.density_flowrate_constraint[:, 1.0].deactivate() + + # But now we've added degrees of freedom. Let's re-check the structural + # diagnostics + dt = DiagnosticsToolbox(model2) + dt.report_structural_issues() + + # After some thought, we decide we need to fix particle porosity at the solid inlet + model2.fs.MB.particle_porosity[:, 1.0].fix() + + # Let's check the structural diagnostics again. + dt = DiagnosticsToolbox(model2) + dt.report_structural_issues() + # Looks good! + + # Now let's try to solve + free_degrees_of_freedom(model2) + attempt_solve(model2) + + +# Below are functions that can be used to construct the model at any point at +# which it may be interesting. These are used for testing. + + +def create_original_model(): + """Create the original model we attempt to solve""" + model = make_model() + return model + + +def create_original_square_model(): + """Create the model at the point at which the first diagnostic checks are run""" + model = make_model() + fix_degrees_of_freedom(model) + return model + + +def create_square_model_with_new_variable_and_constraint(): + """Create the model after the first attempt to fix the singularity""" + model = make_model() + fix_degrees_of_freedom(model) + add_particle_porosity_variable(model) + replace_porosity_parameter_with_variable(model) + add_density_flowrate_constraint(model) + return model + + +def create_corrected_square_model(): + """Create the model after correcting the singularity""" + model = create_square_model_with_new_variable_and_constraint() + model.fs.MB.density_flowrate_constraint[:, 1.0].deactivate() + model.fs.MB.particle_porosity[:, 1.0].fix() + return model + + +if __name__ == "__main__": + main() diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/model.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/model.py new file mode 100644 index 00000000..31832646 --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/model.py @@ -0,0 +1,338 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the +# software owners: The Regents of the University of California, through +# Lawrence Berkeley National Laboratory, National Technology & Engineering +# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia +# University Research Corporation, et al. All rights reserved. +# +# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +""" +import pyomo.environ as pyo +from idaes.core import FlowsheetBlock +import idaes.logger as idaeslog +from pyomo.contrib.mpc import DynamicModelInterface + +from idaes_examples.mod.diagnostics.gas_solid_contactors.unit_models.moving_bed import ( + MBR, + BiMBR, +) +from idaes_examples.mod.diagnostics.gas_solid_contactors.properties.methane_iron_OC_reduction.gas_phase_thermo import ( + GasPhaseThermoParameterBlock, +) +from idaes_examples.mod.diagnostics.gas_solid_contactors.properties.methane_iron_OC_reduction.solid_phase_thermo import ( + SolidPhaseThermoParameterBlock, +) +from idaes_examples.mod.diagnostics.gas_solid_contactors.properties.methane_iron_OC_reduction.hetero_reactions import ( + HeteroReactionParameterBlock, +) + + +def set_default_design_variables(m): + m.fs.MB.bed_diameter.fix(6.5) # m + m.fs.MB.bed_height.fix(5) # m + + +def set_default_inlet_conditions(m, fix_porosity=False): + m.fs.MB.gas_inlet.flow_mol[:].fix(128.20513) # mol/s + m.fs.MB.gas_inlet.temperature[:].fix(298.15) # K + m.fs.MB.gas_inlet.pressure[:].fix(2.00) # bar + m.fs.MB.gas_inlet.mole_frac_comp[:, "CO2"].fix(0.02499) + m.fs.MB.gas_inlet.mole_frac_comp[:, "H2O"].fix(0.00001) + m.fs.MB.gas_inlet.mole_frac_comp[:, "CH4"].fix(0.975) + + m.fs.MB.solid_inlet.flow_mass[:].fix(591.4) # kg/s + m.fs.MB.solid_inlet.temperature[:].fix(1183.15) # K + m.fs.MB.solid_inlet.mass_frac_comp[:, "Fe2O3"].fix(0.45) + m.fs.MB.solid_inlet.mass_frac_comp[:, "Fe3O4"].fix(1e-9) + m.fs.MB.solid_inlet.mass_frac_comp[:, "Al2O3"].fix(0.55) + # Only applicable in the patched version of the model + if fix_porosity: + m.fs.MB.solid_inlet.particle_porosity[:].fix(0.27) + + +def fix_initial_conditions(m): + t0 = m.fs.time.first() + x0 = m.fs.MB.gas_phase.length_domain.first() + xf = m.fs.MB.gas_phase.length_domain.last() + m.fs.MB.gas_phase.material_holdup[t0, ...].fix() + m.fs.MB.gas_phase.energy_holdup[t0, ...].fix() + m.fs.MB.gas_phase.material_holdup[t0, x0, ...].unfix() + m.fs.MB.gas_phase.energy_holdup[t0, x0, ...].unfix() + + m.fs.MB.solid_phase.material_holdup[t0, ...].fix() + m.fs.MB.solid_phase.energy_holdup[t0, ...].fix() + m.fs.MB.solid_phase.material_holdup[t0, xf, ...].unfix() + m.fs.MB.solid_phase.energy_holdup[t0, xf, ...].unfix() + + +def initialize_derivative_variables(m): + t0 = m.fs.time.first() + m.fs.MB.gas_phase.material_accumulation[...].set_value(0.0) + m.fs.MB.gas_phase.energy_accumulation[...].set_value(0.0) + m.fs.MB.solid_phase.material_accumulation[...].set_value(0.0) + m.fs.MB.solid_phase.energy_accumulation[...].set_value(0.0) + + +def add_piecewise_constant_constraints( + m, + sample_points=None, + n_samples=5, + sample_time=60.0, +): + if sample_points is None: + sample_points = [i*sample_time for i in range(n_samples + 1)] + inputs = [ + m.fs.MB.gas_inlet.flow_mol, + m.fs.MB.solid_inlet.flow_mass, + ] + dyninterface = DynamicModelInterface(m, m.fs.time) + input_set, pwc_con = dyninterface.get_piecewise_constant_constraints( + inputs, + sample_points, + tolerance=1e-8, + ) + m.input_set = input_set + m.piecewise_constant_constraints = pwc_con + return inputs + + +def get_state_variable_names(space): + setpoint_states = [] + setpoint_states.extend( + "fs.MB.gas_phase.properties[*,%s].flow_mol" % x + for x in space if x != space.first() + ) + setpoint_states.extend( + "fs.MB.gas_phase.properties[*,%s].temperature" % x + for x in space if x != space.first() + ) + setpoint_states.extend( + "fs.MB.gas_phase.properties[*,%s].pressure" % x + for x in space if x != space.first() + ) + setpoint_states.extend( + "fs.MB.gas_phase.properties[*,%s].mole_frac_comp[CH4]" % x + for x in space if x != space.first() + ) + setpoint_states.extend( + "fs.MB.gas_phase.properties[*,%s].mole_frac_comp[H2O]" % x + for x in space if x != space.first() + ) + setpoint_states.extend( + "fs.MB.gas_phase.properties[*,%s].mole_frac_comp[CO2]" % x + for x in space if x != space.first() + ) + setpoint_states.extend( + "fs.MB.solid_phase.properties[*,%s].flow_mass" % x + for x in space if x != space.last() + ) + setpoint_states.extend( + "fs.MB.solid_phase.properties[*,%s].temperature" % x + for x in space if x != space.last() + ) + setpoint_states.extend( + "fs.MB.solid_phase.properties[*,%s].mass_frac_comp[Fe2O3]" % x + for x in space if x != space.last() + ) + setpoint_states.extend( + "fs.MB.solid_phase.properties[*,%s].mass_frac_comp[Fe3O4]" % x + for x in space if x != space.last() + ) + setpoint_states.extend( + "fs.MB.solid_phase.properties[*,%s].mass_frac_comp[Al2O3]" % x + for x in space if x != space.last() + ) + return setpoint_states + + +def add_objective( + m, + inputs=None, + fix_porosity=False, +): + x0 = m.fs.MB.gas_phase.length_domain.first() + xf = m.fs.MB.gas_phase.length_domain.last() + if inputs is None: + slices = [ + m.fs.MB.gas_phase.properties[:, x0].flow_mol, + m.fs.MB.solid_phase.properties[:, xf].flow_mass, + ] + cuids = [pyo.ComponentUID(slice_) for slice_ in slices] + values = [128.20513, 591.4] + inputs = dict(zip(cuids, values)) + + nxfe = m.fs.MB.gas_phase.length_domain.get_discretization_info()["nfe"] + m_setpoint = make_model( + steady=True, + nxfe=nxfe, + square=True, + fix_porosity=fix_porosity, + ) + setpoint_helper = DynamicModelInterface(m_setpoint, m_setpoint.fs.time) + setpoint_helper.load_data(inputs) + + ipopt = pyo.SolverFactory("ipopt") + res = ipopt.solve(m_setpoint, tee=True) + pyo.assert_optimal_termination(res) + + setpoint_data = setpoint_helper.get_data_at_time() + + state_var_names = get_state_variable_names(m.fs.MB.gas_phase.length_domain) + differential_vars = [m.find_component(name) for name in state_var_names] + + dyninterface = DynamicModelInterface(m, m.fs.time) + m.tracking_var_set, m.tracking_cost = dyninterface.get_penalty_from_target( + setpoint_data, + variables=differential_vars, + ) + m.tracking_obj = pyo.Objective( + expr=sum( + m.tracking_cost[i, t] + for i in m.tracking_var_set + for t in m.fs.time if t != m.fs.time.first() + ) + ) + + +def make_model( + steady=False, + n_samples=5, + sample_time=60.0, + ntfe_per_sample=2, + nxfe=10, + t0=0.0, + initialize=True, + square=False, + fix_porosity=False, +): + dynamic = not steady + m = pyo.ConcreteModel() + fs_config = {"dynamic": dynamic} + if dynamic: + horizon = n_samples * sample_time + fs_config["time_set"] = [t0, horizon] + else: + fs_config["time_set"] = [t0] + fs_config["time_units"] = pyo.units.s + m.fs = FlowsheetBlock(**fs_config) + + # Set up thermo props and reaction props + m.fs.gas_properties = GasPhaseThermoParameterBlock() + m.fs.solid_properties = SolidPhaseThermoParameterBlock() + + m.fs.hetero_reactions = HeteroReactionParameterBlock( + **{ + "solid_property_package": m.fs.solid_properties, + "gas_property_package": m.fs.gas_properties, + } + ) + + m.fs.MB = BiMBR( + **{ + "transformation_method": "dae.finite_difference", + "finite_elements": nxfe, + "has_holdup": True, + "gas_phase_config": {"property_package": m.fs.gas_properties}, + "solid_phase_config": { + "property_package": m.fs.solid_properties, + "reaction_package": m.fs.hetero_reactions, + }, + } + ) + + set_default_design_variables(m) + set_default_inlet_conditions(m, fix_porosity=fix_porosity) + + solver = pyo.SolverFactory("ipopt") + if steady: + # If steady, initialize and solve to default steady state + + # State arguments for initializing property state blocks + # Gas phase temperature is initialized at solid + # temperature because thermal mass of solid >> thermal mass of gas + # Particularly useful for initialization if reaction takes place + blk = m.fs.MB + gas_phase_state_args = { + "flow_mol": blk.gas_inlet.flow_mol[t0].value, + "temperature": blk.solid_inlet.temperature[t0].value, + "pressure": blk.gas_inlet.pressure[t0].value, + "mole_frac": { + "CH4": blk.gas_inlet.mole_frac_comp[t0, "CH4"].value, + "CO2": blk.gas_inlet.mole_frac_comp[t0, "CO2"].value, + "H2O": blk.gas_inlet.mole_frac_comp[t0, "H2O"].value, + }, + } + solid_phase_state_args = { + "flow_mass": blk.solid_inlet.flow_mass[t0].value, + "temperature": blk.solid_inlet.temperature[t0].value, + "mass_frac": { + "Fe2O3": blk.solid_inlet.mass_frac_comp[t0, "Fe2O3"].value, + "Fe3O4": blk.solid_inlet.mass_frac_comp[t0, "Fe3O4"].value, + "Al2O3": blk.solid_inlet.mass_frac_comp[t0, "Al2O3"].value, + }, + } + + if initialize: + print("Initializing steady model") + m.fs.MB.initialize( + outlvl=idaeslog.INFO, + gas_phase_state_args=gas_phase_state_args, + solid_phase_state_args=solid_phase_state_args, + ) + # Create a solver + solver.solve(m.fs.MB, tee=True) + + if steady and not square: + raise ValueError("square=False, steady=True is not supported") + + if dynamic: + ntfe = ntfe_per_sample * n_samples + disc = pyo.TransformationFactory("dae.finite_difference") + disc.apply_to(m, wrt=m.fs.time, scheme="BACKWARD", nfe=ntfe) + + fix_initial_conditions(m) + # Set inlets again now that time has been discretized. + set_default_inlet_conditions(m, fix_porosity=fix_porosity) + initialize_derivative_variables(m) + + if initialize: + print("Constructing a steady model to initialize the dynamic model") + m_steady = make_model( + steady=True, + nxfe=nxfe, + square=True, + fix_porosity=fix_porosity, + ) + m_steady_helper = DynamicModelInterface(m_steady, m_steady.fs.time) + t0 = m_steady.fs.time.first() + steady_data = m_steady_helper.get_data_at_time() + steady_scalar_data = m_steady_helper.get_scalar_variable_data() + + m_helper = DynamicModelInterface(m, m.fs.time) + m_helper.load_data(steady_data) + m_helper.load_data(steady_scalar_data) + print("Solving square problem with dynamic model") + solver.solve(m, tee=True) + + if not square: + x0 = m.fs.MB.length_domain.first() + xf = m.fs.MB.length_domain.last() + inputs = { + m.fs.MB.gas_phase.properties[:, x0].flow_mol: 125.0, + m.fs.MB.solid_phase.properties[:, xf].flow_mass: 600.0, + } + + add_objective(m, inputs=inputs, fix_porosity=fix_porosity) + + m.fs.MB.gas_inlet.flow_mol[:].unfix() + m.fs.MB.gas_inlet.flow_mol[0].fix() + m.fs.MB.solid_inlet.flow_mass[:].unfix() + m.fs.MB.solid_inlet.flow_mass[0].fix() + add_piecewise_constant_constraints(m) + + return m diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/__init__.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/__init__.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/gas_phase_thermo.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/gas_phase_thermo.py new file mode 100644 index 00000000..f0eff91b --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/gas_phase_thermo.py @@ -0,0 +1,780 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the +# software owners: The Regents of the University of California, through +# Lawrence Berkeley National Laboratory, National Technology & Engineering +# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia +# University Research Corporation, et al. All rights reserved. +# +# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +This package provides the necessary constraints for gas phase properties for +the CLC of methane +Components - Methane (CH4), Carbon Dioxide (CO2), Water (H2O) + +Equations written in this model were derived from: +(1) B.E. Poling, J.M. Prausnitz, J.P. O'connell, The Properties of Gases and +Liquids, Mcgraw-Hill, New York, 2001. +(2) National Institute of Standards and Technology, NIST Chemistry WebBook, +https://webbook.nist.gov/chemistry/ (accessed March 10, 2018). + +""" + +# Import Pyomo libraries +from pyomo.environ import (Constraint, + Param, + PositiveReals, + Reals, + value, + Var, + units as pyunits, + ) +from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.opt import SolverFactory + +# Import IDAES cores +from idaes.core import (declare_process_block_class, + MaterialFlowBasis, + PhysicalParameterBlock, + StateBlockData, + StateBlock, + Component, + VaporPhase) +from idaes.core.util.initialization import (fix_state_vars, + revert_state_vars, + solve_indexed_blocks) +from idaes.core.util.misc import add_object_reference +from idaes.core.util.model_statistics import ( + degrees_of_freedom, + number_unfixed_variables_in_activated_equalities) +import idaes.logger as idaeslog + +# Some more information about this module +__author__ = "Chinedu Okoli" + + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +@declare_process_block_class("GasPhaseThermoParameterBlock") +class PhysicalParameterData(PhysicalParameterBlock): + """ + Property Parameter Block Class + + Contains parameters and indexing sets associated with properties for + methane CLC. + """ + + def build(self): + ''' + Callable method for Block construction. + ''' + super(PhysicalParameterData, self).build() + + self._state_block_class = GasPhaseThermoStateBlock + + # Create Phase object + self.Vap = VaporPhase() + + # Create Component objects + self.CH4 = Component() + self.CO2 = Component() + self.H2O = Component() + + # Thermodynamic reference state + self.pressure_ref = Param(within=PositiveReals, + mutable=True, + default=1.01325, + doc='Reference pressure [bar]') + self.temperature_ref = Param(within=PositiveReals, + mutable=True, + default=298.15, + doc='Thermodynamic Reference' + 'Temperature [K]') + + # Gas Constant + self.gas_const = Param(within=PositiveReals, + default=8.314459848e-3, + doc='Gas Constant [kJ/mol.K]') + + # ------------------------------------------------------------------------- + """ Pure gas component properties""" + + # Mol. weights of gas - units = kg/mol. ref: NIST webbook + mw_comp_dict = {'CH4': 0.016, 'CO2': 0.044, 'H2O': 0.018} + self.mw_comp = Param( + self.component_list, + mutable=False, + initialize=mw_comp_dict, + doc="Molecular weights of gas components [kg/mol]") + + # Std. heat of formation of comp. - units = kJ/(mol comp) - ref: NIST + enth_mol_form_comp_dict = {'CH4': -74.8731, 'CO2': -393.5224, + 'H2O': -241.8264} + self.enth_mol_form_comp = Param( + self.component_list, + mutable=False, + initialize=enth_mol_form_comp_dict, + doc="Component molar heats of formation [kJ/mol]") + + # Ideal gas spec. heat capacity parameters(Shomate) of + # components - ref: NIST webbook. Shomate equations from NIST. + # Parameters A-E are used for cp calcs while A-H are used for enthalpy + # calc. + # 1e3*cp_comp = A + B*T + C*T^2 + D*T^3 + E/(T^2) + # where T = Temperature (K)/1000, and cp_comp = (kJ/mol.K) + # H_comp = H - H(298.15) = A*T + B*T^2/2 + C*T^3/3 + + # D*T^4/4 - E/T + F - H where T = Temp (K)/1000 and H_comp = (kJ/mol) + cp_param_dict = { + ('CH4', 1): -0.7030290, + ('CH4', 2): 108.4773000, + ('CH4', 3): -42.5215700, + ('CH4', 4): 5.8627880, + ('CH4', 5): 0.6785650, + ('CH4', 6): -76.8437600, + ('CH4', 7): 158.7163000, + ('CH4', 8): -74.8731000, + ('CO2', 1): 24.9973500, + ('CO2', 2): 55.1869600, + ('CO2', 3): -33.6913700, + ('CO2', 4): 7.9483870, + ('CO2', 5): -0.1366380, + ('CO2', 6): -403.6075000, + ('CO2', 7): 228.2431000, + ('CO2', 8): -393.5224000, + ('H2O', 1): 30.0920000, + ('H2O', 2): 6.8325140, + ('H2O', 3): 6.7934350, + ('H2O', 4): -2.5344800, + ('H2O', 5): 0.0821390, + ('H2O', 6): -250.8810000, + ('H2O', 7): 223.3967000, + ('H2O', 8): -241.8264000 + } + self.cp_param = Param(self.component_list, + range(1, 10), + mutable=False, + initialize=cp_param_dict, + doc="Shomate equation heat capacity parameters") + + # Viscosity constants: + # Reference: Perry and Green Handbook; McGraw Hill, 2008 + visc_d_param_dict = {('CH4', 1): 5.2546e-7, ('CH4', 2): 0.59006, + ('CH4', 3): 105.67, ('CH4', 4): 0, + ('CO2', 1): 2.148e-6, ('CO2', 2): 0.46, + ('CO2', 3): 290, ('CO2', 4): 0, + ('H2O', 1): 1.7096e-8, ('H2O', 2): 1.1146, + ('H2O', 3): 0, ('H2O', 4): 0} + self.visc_d_param = Param(self.component_list, + range(1, 10), + mutable=True, + initialize=visc_d_param_dict, + doc="Dynamic viscosity constants") + + # Thermal conductivity constants: + # Reference: Perry and Green Handbook; McGraw Hill, 2008 + therm_cond_param_dict = {('CH4', 1): 8.3983e-6, ('CH4', 2): 1.4268, + ('CH4', 3): -49.654, ('CH4', 4): 0, + ('CO2', 1): 3.69, ('CO2', 2): -0.3838, + ('CO2', 3): 964, ('CO2', 4): 1.86e6, + ('H2O', 1): 6.204e-6, ('H2O', 2): 1.3973, + ('H2O', 3): 0, ('H2O', 4): 0} + self.therm_cond_param = Param(self.component_list, + range(1, 10), + mutable=True, + initialize=therm_cond_param_dict, + doc="Thermal conductivity constants") + + # Component diffusion volumes: + # Ref: (1) Prop gas & liquids (2) Fuller et al. IECR, 58(5), 19, 1966 + diff_vol_param_dict = {'CH4': 24.42, 'CO2': 26.9, 'H2O': 13.1} + self.diff_vol_param = Param(self.component_list, + mutable=True, + initialize=diff_vol_param_dict, + doc="Component diffusion volumes") + + @classmethod + def define_metadata(cls, obj): + obj.add_properties({ + 'flow_mol': {'method': None, 'units': 'mol/s'}, + 'pressure': {'method': None, 'units': 'bar'}, + 'temperature': {'method': None, 'units': 'K'}, + 'mole_frac_comp': {'method': None, 'units': None}, + 'mw': {'method': '_mw', 'units': 'kg/mol'}, + 'cp_mol': {'method': '_cp_mol', 'units': 'kJ/mol.K'}, + 'cp_mol_comp': {'method': '_cp_mol_comp', + 'units': 'kJ/mol.K'}, + 'cp_mass': {'method': '_cp_mass', 'units': 'kJ/kg.K'}, + 'dens_mol': {'method': '_dens_mol', + 'units': 'mol/m^3'}, + 'dens_mol_comp': {'method': '_dens_mol_comp', + 'units': 'mol/m^3'}, + 'dens_mass': {'method': '_dens_mass', + 'units': 'kg/m^3'}, + 'enth_mol': {'method': '_enth_mol', 'units': 'kJ/mol'}, + 'enth_mol_comp': {'method': '_enth_mol_comp', + 'units': 'kJ/mol'}, + 'visc_d': {'method': '_visc_d', 'units': 'kg/m.s'}, + 'therm_cond': {'method': '_therm_cond', 'units': 'kJ/m.K.s'}, + 'diffusion_comp': {'method': '_diffusion_comp', + 'units': 'cm2/s'}}) + + obj.add_default_units({'time': pyunits.s, + 'length': pyunits.m, + 'mass': pyunits.kg, + 'amount': pyunits.mol, + 'temperature': pyunits.K, + #'energy': pyunits.kJ, + #'holdup': pyunits.mol, + }) + + +class _GasPhaseThermoStateBlock(StateBlock): + """ + This Class contains methods which should be applied to Property Blocks as a + whole, rather than individual elements of indexed Property Blocks. + """ + def initialize(blk, state_args=None, hold_state=False, + state_vars_fixed=False, outlvl=idaeslog.NOTSET, + solver="ipopt", optarg={"tol": 1e-8}): + """ + Initialization routine for property package. + Keyword Arguments: + state_args : Dictionary with initial guesses for the state vars + chosen. Note that if this method is triggered + through the control volume, and if initial guesses + were not provided at the unit model level, the + control volume passes the inlet values as initial + guess. + Keys for the state_args dictionary are: + flow_mol, temperature, pressure and mole_frac_comp + outlvl : sets output level of initialization routine + optarg : solver options dictionary object (default=None) + solver : str indicating which solver to use during + initialization (default = "ipopt") + hold_state : flag indicating whether the initialization routine + should unfix any state variables fixed during + initialization (default=False). + - True - states variables are not unfixed, and + a dict of returned containing flags for + which states were fixed during + initialization. + - False - state variables are unfixed after + initialization by calling the + release_state method + Returns: + If hold_states is True, returns a dict containing flags for + which states were fixed during initialization. + """ + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties") + solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="properties") + + init_log.info_high('Starting initialization') + + # Deactivate the constraints specific for non-inlet blocks i.e. + # when defined state is False + for k in blk.keys(): + if blk[k].config.defined_state is False: + blk[k].sum_component_eqn.deactivate() + + # Fix state variables if not already fixed + if state_vars_fixed is False: + flags = fix_state_vars(blk, state_args) + else: + # Check when the state vars are fixed already result in dof 0 + for k in blk.keys(): + if degrees_of_freedom(blk[k]) != 0: + raise Exception("State vars fixed but degrees of freedom " + "for state block is not zero during " + "initialization.") + + # Set solver options + opt = SolverFactory(solver) + opt.options = optarg + + # --------------------------------------------------------------------- + # Initialise values + for k in blk.keys(): + + if hasattr(blk[k], "mw_eqn"): + calculate_variable_from_constraint( + blk[k].mw, + blk[k].mw_eqn) + + if hasattr(blk[k], "ideal_gas"): + calculate_variable_from_constraint( + blk[k].dens_mol, + blk[k].ideal_gas) + + if hasattr(blk[k], "dens_mass_basis"): + calculate_variable_from_constraint( + blk[k].dens_mass, + blk[k].dens_mass_basis) + + if hasattr(blk[k], "mixture_heat_capacity_eqn"): + calculate_variable_from_constraint( + blk[k].cp_mol, + blk[k].mixture_heat_capacity_eqn) + + if hasattr(blk[k], "cp_mass_basis"): + calculate_variable_from_constraint( + blk[k].cp_mass, + blk[k].cp_mass_basis) + + if hasattr(blk[k], "visc_d_constraint"): + calculate_variable_from_constraint( + blk[k].visc_d, + blk[k].visc_d_constraint) + + if hasattr(blk[k], "therm_cond_constraint"): + calculate_variable_from_constraint( + blk[k].therm_cond, + blk[k].therm_cond_constraint) + + if hasattr(blk[k], "mixture_enthalpy_eqn"): + calculate_variable_from_constraint( + blk[k].enth_mol, + blk[k].mixture_enthalpy_eqn) + + for j in blk[k]._params.component_list: + + if hasattr(blk[k], "comp_conc_eqn"): + calculate_variable_from_constraint( + blk[k].dens_mol_comp[j], + blk[k].comp_conc_eqn[j]) + + if hasattr(blk[k], "diffusion_comp_constraint"): + calculate_variable_from_constraint( + blk[k].diffusion_comp[j], + blk[k].diffusion_comp_constraint[j]) + + if hasattr(blk[k], "cp_shomate_eqn"): + calculate_variable_from_constraint( + blk[k].cp_mol_comp[j], + blk[k].cp_shomate_eqn[j]) + + if hasattr(blk[k], "enthalpy_shomate_eqn"): + calculate_variable_from_constraint( + blk[k].enth_mol_comp[j], + blk[k].enthalpy_shomate_eqn[j]) + + # Solve property block if non-empty + free_vars = 0 + for k in blk.keys(): + free_vars += number_unfixed_variables_in_activated_equalities( + blk[k]) + + if free_vars > 0: + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solve_indexed_blocks(opt, [blk], tee=slc.tee) + else: + res = "" + init_log.info_high("Initialization complete {}.".format( + idaeslog.condition(res)) + ) + + # --------------------------------------------------------------------- + if state_vars_fixed is False: + if hold_state is True: + return flags + else: + blk.release_state(flags) + + def release_state(blk, flags, outlvl=0): + """ + Method to release state variables fixed during initialization. + Keyword Arguments: + flags : dict containing information of which state variables + were fixed during initialization, and should now be + unfixed. This dict is returned by initialize if + hold_state=True. + outlvl : sets output level of logging + """ + if flags is None: + return + + # Unfix state variables + revert_state_vars(blk, flags) + + # Activate state variable related constraints + for k in blk.keys(): + if blk[k].config.defined_state is False: + blk[k].sum_component_eqn.activate() + + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties") + init_log.info_high('States released.') + + +@declare_process_block_class("GasPhaseThermoStateBlock", + block_class=_GasPhaseThermoStateBlock) +class GasPhaseThermoStateBlockData(StateBlockData): + """ + Property package for gas phase properties of methane combustion in CLC FR + """ + + def build(self): + """ + Callable method for Block construction + """ + super(GasPhaseThermoStateBlockData, self).build() + + # Object reference for molecular weight if needed by CV1D + # Molecular weights + add_object_reference(self, "mw_comp", + self.config.parameters.mw_comp) + + """List the necessary state variable objects.""" + self.flow_mol = Var(initialize=1.0, + domain=Reals, + doc='Component molar flowrate [mol/s]') + self.mole_frac_comp = Var( + self._params.component_list, + domain=Reals, + initialize=1 / len(self._params.component_list), + doc='State component mole fractions [-]') + self.pressure = Var(initialize=1.01325, + domain=Reals, + doc='State pressure [bar]') + self.temperature = Var(initialize=298.15, + domain=Reals, + doc='State temperature [K]') + + # Create standard constraints + # Sum mole fractions if not inlet block + if self.config.defined_state is False: + def sum_component_eqn(b): + return 1e2 == 1e2 * sum(b.mole_frac_comp[j] + for j in b._params.component_list) + self.sum_component_eqn = Constraint(rule=sum_component_eqn) + + def _mw(self): + # Molecular weight of gas mixture + self.mw = Var(domain=Reals, + initialize=1.0, + doc="Molecular weight of gas mixture [kg/mol]") + + def mw_eqn(b): + return (b.mw == + sum(b.mole_frac_comp[j]*b._params.mw_comp[j] + for j in b._params.component_list)) + try: + # Try to build constraint + self.mw_eqn = Constraint(rule=mw_eqn) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.mw) + self.del_component(self.mw_eqn) + raise + + def _dens_mol(self): + # Molar density + self.dens_mol = Var(domain=Reals, + initialize=1.0, + doc="Molar density/concentration [mol/m3]") + + def ideal_gas(b): + return (b.dens_mol*b._params.gas_const*b.temperature*1e-2 == + b.pressure) + try: + # Try to build constraint + self.ideal_gas = Constraint(rule=ideal_gas) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.dens_mol) + self.del_component(self.ideal_gas) + raise + + def _dens_mol_comp(self): + # Mixture heat capacities + self.dens_mol_comp = Var(self._params.component_list, + domain=Reals, + initialize=1.0, + doc='Component molar concentration' + '[mol/m3]') + + def comp_conc_eqn(b, j): + return (b.dens_mol_comp[j] == + b.dens_mol*b.mole_frac_comp[j]) + try: + # Try to build constraint + self.comp_conc_eqn = Constraint(self._params.component_list, + rule=comp_conc_eqn) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.dens_mol_comp) + self.del_component(self.comp_conc_eqn) + raise + + def _dens_mass(self): + # Mass density + self.dens_mass = Var(domain=Reals, + initialize=1.0, + doc="Mass density [kg/m3]") + + def dens_mass_basis(b): + return b.dens_mass == b.mw*b.dens_mol + try: + # Try to build constraint + self.dens_mass_basis = Constraint(rule=dens_mass_basis) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.dens_mass) + self.del_component(self.dens_mass_basis) + raise + + def _visc_d(self): + # Mixture dynamic viscosity + self.visc_d = Var(domain=Reals, + initialize=1e-5, + doc="Mixture dynamic viscosity [kg/m.s]") + + def visc_d_comp(i): + return self._params.visc_d_param[i, 1] * \ + (self.temperature**self._params.visc_d_param[i, 2]) \ + / ((1 + (self._params.visc_d_param[i, 3]/self.temperature)) + + (self._params.visc_d_param[i, 4] / + (self.temperature**2))) + + def visc_d_constraint(b): + return 1e6*b.visc_d == 1e6*sum(b.mole_frac_comp[i]*visc_d_comp(i) + / (sum(b.mole_frac_comp[j] + * (b._params.mw_comp[j] / + b._params.mw_comp[i])**0.5 + for j in + b._params.component_list)) + for i in + b._params.component_list) + try: + # Try to build constraint + self.visc_d_constraint = Constraint(rule=visc_d_constraint) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.visc_d) + self.del_component(self.visc_d_constraint) + raise + + def _diffusion_comp(self): + # Component diffusion in a gas mixture - units of cm2/s to help scaling + self.diffusion_comp = Var(self._params.component_list, + domain=Reals, + initialize=1e-5, + doc='Component diffusion in a gas mixture' + '[cm2/s]') + + def D_bin(i, j): + # 1e3 used to multiply MW to convert from kg/mol to kg/kmol + return ((1.43e-3*(self.temperature**1.75) * + ((1e3 * self._params.mw_comp[i] + + 1e3 * self._params.mw_comp[j]) + / (2 * (1e3 * self._params.mw_comp[i]) * + (1e3*self._params.mw_comp[j])))**0.5) + / ((self.pressure) + * ((self._params.diff_vol_param[i]**(1/3)) + + (self._params.diff_vol_param[j]**(1/3)))**2)) + + def diffusion_comp_constraint(b, i): + return (b.diffusion_comp[i] + * sum(b.mole_frac_comp[j]/D_bin(i, j) + for j in b._params.component_list if i != j) + == (1-b.mole_frac_comp[i])) + try: + # Try to build constraint + self.diffusion_comp_constraint \ + = Constraint(self._params.component_list, + rule=diffusion_comp_constraint) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.diffusion_comp) + self.del_component(self.diffusion_comp_constraint) + raise + + def _therm_cond(self): + # Thermal conductivity of gas + self.therm_cond = Var(domain=Reals, + initialize=1e-5, + doc="Thermal conductivity of gas [kJ/m.K.s]") + + def therm_cond_comp(i): + return self._params.therm_cond_param[i, 1] \ + * (self.temperature**self._params.therm_cond_param[i, 2]) \ + / ((1 + (self._params.therm_cond_param[i, 3] / + self.temperature)) + + (self._params.therm_cond_param[i, 4] / + (self.temperature**2))) + + def A_bin(i, j): + return (1 + ((therm_cond_comp(j)/therm_cond_comp(i))**0.5) + * ((self._params.mw_comp[j] / + self._params.mw_comp[i])**0.25))**2 \ + / (8*(1+(self._params.mw_comp[j] / + self._params.mw_comp[i])))**0.5 + + def therm_cond_constraint(b): + # The 1e-3 term is used as a conversion factor to a kJ basis + return 1e6*b.therm_cond == 1e6*(1e-3) * \ + sum(b.mole_frac_comp[i] + * therm_cond_comp(i) + / (sum(b.mole_frac_comp[j] * + A_bin(i, j)**0.5 + for j in + b._params.component_list)) + for i in + b._params.component_list) + try: + # Try to build constraint + self.therm_cond_constraint = Constraint(rule=therm_cond_constraint) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.therm_cond) + self.del_component(self.therm_cond_constraint) + raise + + def _cp_mol_comp(self): + # Pure component vapour heat capacities + self.cp_mol_comp = Var(self._params.component_list, + domain=Reals, + initialize=1.0, + doc="Pure component vapour heat capacities " + "[kJ/mol.K]") + + def pure_component_cp_mol(b, j): + return b.cp_mol_comp[j] == 1e-3*( + b._params.cp_param[j, 1] + + b._params.cp_param[j, 2]*(b.temperature*1e-3) + + b._params.cp_param[j, 3]*(b.temperature*1e-3)**2 + + b._params.cp_param[j, 4]*(b.temperature*1e-3)**3 + + b._params.cp_param[j, 5]/((b.temperature*1e-3)**2)) + try: + # Try to build constraint + self.cp_shomate_eqn = Constraint(self._params.component_list, + rule=pure_component_cp_mol) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.cp_mol_comp) + self.del_component(self.cp_shomate_eqn) + raise + + def _cp_mol(self): + # Mixture heat capacities + self.cp_mol = Var(domain=Reals, + initialize=1.0, + doc="Mixture heat capacity [kJ/mol.K]") + + def cp_mol(b): + return b.cp_mol == sum(b.cp_mol_comp[j]*b.mole_frac_comp[j] + for j in b._params.component_list) + try: + # Try to build constraint + self.mixture_heat_capacity_eqn = Constraint( + rule=cp_mol) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.cp_mol) + self.del_component(self.mixture_heat_capacity_eqn) + raise + + def _cp_mass(self): + # Mixture heat capacities + self.cp_mass = Var(domain=Reals, + initialize=1.0, + doc="Mixture heat capacity, mass-basis [kJ/kg.K]") + + def cp_mass(b): + return b.cp_mass*b.mw == b.cp_mol + try: + # Try to build constraint + self.cp_mass_basis = Constraint(rule=cp_mass) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.cp_mass) + self.del_component(self.cp_mass_basis) + raise + + def _enth_mol_comp(self): + # Pure component vapour enthalpies + self.enth_mol_comp = Var( + self._params.component_list, + domain=Reals, + initialize=1.0, + doc="Pure component enthalpies [kJ/mol]") + + def pure_comp_enthalpy(b, j): + return b.enth_mol_comp[j] == ( + b._params.cp_param[j, 1]*(b.temperature*1e-3) + + b._params.cp_param[j, 2]*((b.temperature*1e-3)**2)/2 + + b._params.cp_param[j, 3]*((b.temperature*1e-3)**3)/3 + + b._params.cp_param[j, 4]*((b.temperature*1e-3)**4)/4 - + b._params.cp_param[j, 5]/(b.temperature*1e-3) + + b._params.cp_param[j, 6] - + b._params.cp_param[j, 8]) + try: + # Try to build constraint + self.enthalpy_shomate_eqn = Constraint(self._params.component_list, + rule=pure_comp_enthalpy) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.enth_mol_comp) + self.del_component(self.enthalpy_shomate_eqn) + raise + + def _enth_mol(self): + # Mixture molar enthalpy + self.enth_mol = Var( + domain=Reals, + initialize=1.0, + doc='Mixture specific enthalpy [kJ/mol]') + try: + # Try to build constraint + self.mixture_enthalpy_eqn = Constraint(expr=( + self.enth_mol == sum(self.mole_frac_comp[j] * + self.enth_mol_comp[j] + for j in + self._params.component_list))) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.enth_mol) + self.del_component(self.mixture_enthalpy_eqn) + raise + + def get_material_flow_terms(b, p, j): + return b.flow_mol*b.mole_frac_comp[j] + + def get_enthalpy_flow_terms(b, p): + return b.flow_mol*b.enth_mol + + def get_material_density_terms(b, p, j): + return b.dens_mol_comp[j] + + def get_energy_density_terms(b, p): + return b.dens_mol*b.enth_mol + + def define_state_vars(b): + return {"flow_mol": b.flow_mol, + "temperature": b.temperature, + "pressure": b.pressure, + "mole_frac_comp": b.mole_frac_comp} + + def get_material_flow_basis(b): + return MaterialFlowBasis.molar + + def model_check(blk): + """ + Model checks for property block + """ + # Check temperature bounds + if value(blk.temperature) < blk.temperature.lb: + _log.error('{} Temperature set below lower bound.' + .format(blk.name)) + if value(blk.temperature) > blk.temperature.ub: + _log.error('{} Temperature set above upper bound.' + .format(blk.name)) + + # Check pressure bounds + if value(blk.pressure) < blk.pressure.lb: + _log.error('{} Pressure set below lower bound.'.format(blk.name)) + if value(blk.pressure) > blk.pressure.ub: + _log.error('{} Pressure set above upper bound.'.format(blk.name)) diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/hetero_reactions.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/hetero_reactions.py new file mode 100644 index 00000000..7a076cb6 --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/hetero_reactions.py @@ -0,0 +1,523 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the +# software owners: The Regents of the University of California, through +# Lawrence Berkeley National Laboratory, National Technology & Engineering +# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia +# University Research Corporation, et al. All rights reserved. +# +# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +Property package for the reaction of CH4 with an iron-based OC. +Overall reducer reactions for Methane combustion: +(1) CH4 + 12Fe2O3 => 8Fe3O4 + CO2 + 2H2O + +Parameters and equations written in this model were primarily derived from: +A. Abad, J. Adánez, F. García-Labiano, L.F. de Diego, P. Gayán, J. Celaya, +Mapping of the range of operational conditions for cu-, Fe-, and Ni-based +oxygen carriers in chemical-looping combustion, +Chem. Eng. Sci. 62 (2007) 533–549. + +""" + +# Import Pyomo libraries +from pyomo.environ import (Constraint, + exp, + Param, + PositiveReals, + Reals, + Set, + value, + Var) +import pyomo.environ as pyo +from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.opt import SolverFactory +from pyomo.common.config import ConfigBlock, ConfigValue, In + + +# Import IDAES cores +from idaes.core import (declare_process_block_class, + MaterialFlowBasis, + ReactionParameterBlock, + ReactionBlockDataBase, + ReactionBlockBase, + Component, + VaporPhase, + SolidPhase) +from idaes.core.util.misc import add_object_reference +from idaes.core.util.initialization import (fix_state_vars, + revert_state_vars, + solve_indexed_blocks) +from idaes.core.util.model_statistics import ( + number_unfixed_variables_in_activated_equalities) +from idaes.core.util.config import (is_state_block, + is_physical_parameter_block, + is_reaction_parameter_block) +import idaes.logger as idaeslog + +# Some more information about this module +__author__ = "Chinedu Okoli" + + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +@declare_process_block_class("HeteroReactionParameterBlock") +class ReactionParameterData(ReactionParameterBlock): + """ + Property Parameter Block Class + + Contains parameters and indexing sets associated with properties for + superheated steam. + + """ + # Create Class ConfigBlock + CONFIG = ConfigBlock() + CONFIG.declare("gas_property_package", ConfigValue( + description="Reference to associated PropertyPackageParameter " + "object for the gas phase.", + domain=is_physical_parameter_block)) + CONFIG.declare("solid_property_package", ConfigValue( + description="Reference to associated PropertyPackageParameter " + "object for the solid phase.", + domain=is_physical_parameter_block)) + CONFIG.declare("default_arguments", ConfigBlock( + description="Default arguments to use with Property Package", + implicit=True)) + + def build(self): + ''' + Callable method for Block construction. + ''' + super(ReactionParameterBlock, self).build() + + self._reaction_block_class = ReactionBlock + + # Create Phase objects + self.Vap = VaporPhase() + self.Sol = SolidPhase() + + # Create Component objects + self.CH4 = Component() + self.CO2 = Component() + self.H2O = Component() + self.Fe2O3 = Component() + self.Fe3O4 = Component() + self.Al2O3 = Component() + + # Component list subsets + self.gas_component_list = Set(initialize=['CO2', 'H2O', 'CH4']) + self.sol_component_list = Set(initialize=['Fe2O3', 'Fe3O4', 'Al2O3']) + + # Reaction Index + self.rate_reaction_idx = Set(initialize=["R1"]) + + # Gas Constant + self.gas_const = Param(within=PositiveReals, + mutable=False, + default=8.314459848e-3, + doc='Gas Constant [kJ/mol.K]') + + # Smoothing factor + self.eps = Param(mutable=True, + default=1e-8, + doc='Smoothing Factor') + # Reaction rate scale factor + self._scale_factor_rxn = Param(mutable=True, + default=1, + doc='Scale Factor for reaction eqn.' + 'Used to help initialization routine') + + # Reaction Stoichiometry + self.rate_reaction_stoichiometry = {("R1", "Vap", "CH4"): -1, + ("R1", "Vap", "CO2"): 1, + ("R1", "Vap", "H2O"): 2, + ("R1", "Sol", "Fe2O3"): -12, + ("R1", "Sol", "Fe3O4"): 8, + ("R1", "Sol", "Al2O3"): 0} + + # Reaction stoichiometric coefficient + self.rxn_stoich_coeff = Param(self.rate_reaction_idx, + default=12, + mutable=True, + doc='Reaction stoichiometric' + 'coefficient [-]') + + # Standard Heat of Reaction - kJ/mol_rxn + dh_rxn_dict = {"R1": 136.5843} + self.dh_rxn = Param(self.rate_reaction_idx, + initialize=dh_rxn_dict, + doc="Heat of reaction [kJ/mol]") + + # ------------------------------------------------------------------------- + """ Reaction properties that can be estimated""" + + # Particle grain radius within OC particle + self.grain_radius = Var(domain=Reals, + initialize=2.6e-7, + doc='Representative particle grain' + 'radius within OC particle [m]') + self.grain_radius.fix() + + # Molar density OC particle + self.dens_mol_sol = Var(domain=Reals, + initialize=32811, + doc='Molar density of OC particle [mol/m^3]') + self.dens_mol_sol.fix() + + # Available volume for reaction - from EPAT report (1-ep)' + self.a_vol = Var(domain=Reals, + initialize=0.28, + doc='Available reaction vol. per vol. of OC') + self.a_vol.fix() + + # Activation Energy + self.energy_activation = Var(self.rate_reaction_idx, + domain=Reals, + initialize=4.9e1, + doc='Activation energy [kJ/mol]') + self.energy_activation.fix() + + # Reaction order + self.rxn_order = Var(self.rate_reaction_idx, + domain=Reals, + initialize=1.3, + doc='Reaction order in gas species [-]') + self.rxn_order.fix() + + # Pre-exponential factor + self.k0_rxn = Var(self.rate_reaction_idx, + domain=Reals, + initialize=8e-4, + doc='Pre-exponential factor' + '[mol^(1-N_reaction)m^(3*N_reaction -2)/s]') + self.k0_rxn.fix() + + @classmethod + def define_metadata(cls, obj): + obj.define_custom_properties( + { + "OC_conv": {"method": "_OC_conv", "units": None}, + "OC_conv_temp": {"method": "_OC_conv_temp", "units": None}, + } + ) + obj.add_properties({ + 'k_rxn': {'method': '_k_rxn', + 'units': 'mol^(1-N_reaction)m^(3*N_reaction -2)/s]'}, + 'OC_conv': {'method': "_OC_conv", 'units': None}, + 'OC_conv_temp': {'method': "_OC_conv_temp", 'units': None}, + 'reaction_rate': {'method': "_reaction_rate", + 'units': 'mol_rxn/m3.s'} + }) + obj.add_default_units({'time': pyo.units.s, + 'length': pyo.units.m, + 'mass': pyo.units.kg, + 'amount': pyo.units.mol, + 'temperature': pyo.units.K, + #'energy': pyo.units.kJ + }) + + +class _ReactionBlock(ReactionBlockBase): + """ + This Class contains methods which should be applied to Reaction Blocks as a + whole, rather than individual elements of indexed Reaction Blocks. + """ + def initialize(blk, outlvl=idaeslog.NOTSET, + optarg={'tol': 1e-8}, solver='ipopt'): + ''' + Initialisation routine for reaction package. + + Keyword Arguments: + outlvl : sets output level of initialization routine + * 0 = Use default idaes.init logger setting + * 1 = Maximum output + * 2 = Include solver output + * 3 = Return solver state for each step in subroutines + * 4 = Return solver state for each step in routine + * 5 = Final initialization status and exceptions + * 6 = No output + optarg : solver options dictionary object (default=None) + solver : str indicating which solver to use during + initialization (default = "ipopt") + Returns: + None + ''' + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="reactions") + solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="reactions") + + init_log.info_high('Starting initialization') + + # TODO - Update in the future as needed + # Get a single representative block for getting config arguments + for k in blk.keys(): + break + + # Fix state variables if not already fixed + # Fix state variables of the primary (solid) state block + state_var_flags = fix_state_vars(blk[k].config.solid_state_block) + + # Fix values of secondary (gas) state block variables if not fixed, + # as well as the solid density variable. + # This is done to keep the initialization problem square + Cflag = {} # Gas concentration flag + Dflag = {} # Solid density flag + + for k in blk.keys(): + for j in blk[k]._params.gas_component_list: + if blk[k].gas_state_ref.dens_mol_comp[j].fixed is True: + Cflag[k, j] = True + else: + Cflag[k, j] = False + blk[k].gas_state_ref.dens_mol_comp[j].fix( + blk[k].gas_state_ref.dens_mol_comp[j].value) + if blk[k].solid_state_ref.dens_mass_skeletal.fixed is True: + Dflag[k] = True + else: + Dflag[k] = False + blk[k].solid_state_ref.dens_mass_skeletal.fix( + blk[k].solid_state_ref.dens_mass_skeletal.value) + + # Set solver options + opt = SolverFactory(solver) + opt.options = optarg + + # Initialise values + for k in blk.keys(): + if hasattr(blk[k], "OC_conv_eqn"): + calculate_variable_from_constraint( + blk[k].OC_conv, + blk[k].OC_conv_eqn) + + if hasattr(blk[k], "OC_conv_temp_eqn"): + calculate_variable_from_constraint( + blk[k].OC_conv_temp, + blk[k].OC_conv_temp_eqn) + + for j in blk[k]._params.rate_reaction_idx: + if hasattr(blk[k], "rate_constant_eqn"): + calculate_variable_from_constraint( + blk[k].k_rxn[j], + blk[k].rate_constant_eqn[j]) + + if hasattr(blk[k], "gen_rate_expression"): + calculate_variable_from_constraint( + blk[k].reaction_rate[j], + blk[k].gen_rate_expression[j]) + + # Solve property block if non-empty + free_vars = 0 + for k in blk.keys(): + free_vars += number_unfixed_variables_in_activated_equalities( + blk[k]) + + if free_vars > 0: + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solve_indexed_blocks(opt, [blk], tee=slc.tee) + else: + res = "" + init_log.info_high("reactions initialization complete {}.".format( + idaeslog.condition(res)) + ) + + # --------------------------------------------------------------------- + # Revert state vars and other variables to pre-initialization states + # Revert state variables of the primary (solid) state block + revert_state_vars(blk[k].config.solid_state_block, state_var_flags) + + for k in blk.keys(): + for j in blk[k]._params.gas_component_list: + if Cflag[k, j] is False: + blk[k].gas_state_ref.dens_mol_comp[j].unfix() + if Dflag[k] is False: + blk[k].solid_state_ref.dens_mass_skeletal.unfix() + + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="reactions") + init_log.info_high('States released.') + + +@declare_process_block_class("ReactionBlock", + block_class=_ReactionBlock) +class ReactionBlockData(ReactionBlockDataBase): + """ + Heterogeneous reaction package for methane reacting with Fe2O3 based OC + """ + # Create Class ConfigBlock + CONFIG = ConfigBlock() + CONFIG.declare("parameters", ConfigValue( + domain=is_reaction_parameter_block, + description=""" + A reference to an instance of the Reaction Parameter + Block associated with this property package. + """)) + CONFIG.declare("solid_state_block", ConfigValue( + domain=is_state_block, + description=""" + A reference to an instance of a StateBlock for the + solid phase with which this reaction block should be associated. + """)) + CONFIG.declare("gas_state_block", ConfigValue( + domain=is_state_block, + description=""" + A reference to an instance of a StateBlock for the + gas phase with which this reaction block should be associated. + """)) + CONFIG.declare("has_equilibrium", ConfigValue( + default=False, + domain=In([True, False]), + description="Equilibrium reaction construction flag", + doc=""" + Indicates whether terms for equilibrium controlled reactions + should be constructed, + **default** - True. + **Valid values:** { + **True** - include equilibrium reaction terms, + **False** - exclude equilibrium reaction terms.} + """)) + + def build(self): + """ + Callable method for Block construction + """ + super(ReactionBlockDataBase, self).build() + + # Object references to the corresponding state blocks and parameters + add_object_reference(self, "_params", self.config.parameters) + + add_object_reference(self, + "solid_state_ref", + self.config.solid_state_block[self.index()]) + add_object_reference(self, + "gas_state_ref", + self.config.gas_state_block[self.index()]) + + # Object reference for parameters if needed by CV1D + # Reaction stoichiometry + add_object_reference( + self, + "rate_reaction_stoichiometry", + self.config.parameters.rate_reaction_stoichiometry) + + # Heat of reaction + add_object_reference( + self, + "dh_rxn", + self.config.parameters.dh_rxn) + + # Rate constant method + def _k_rxn(self): + self.k_rxn = Var(self._params.rate_reaction_idx, + domain=Reals, + initialize=1, + doc='Rate constant ' + '[mol^(1-N_reaction)m^(3*N_reaction -2)/s]') + + def rate_constant_eqn(b, j): + if j == 'R1': + return 1e6 * self.k_rxn[j] == \ + 1e6 * (self._params.k0_rxn[j] * + exp(-self._params.energy_activation[j] / + (self._params.gas_const * + self.solid_state_ref.temperature))) + else: + return Constraint.Skip + try: + # Try to build constraint + self.rate_constant_eqn = Constraint(self._params.rate_reaction_idx, + rule=rate_constant_eqn) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.k_rxn) + self.del_component(self.rate_constant_eqn) + raise + + # Conversion of oxygen carrier + def _OC_conv(self): + self.OC_conv = Var(domain=Reals, initialize=0.0, + doc='Fraction of metal oxide converted') + + def OC_conv_eqn(b): + return 1e6 * b.OC_conv * \ + (b.solid_state_ref.mass_frac_comp['Fe3O4'] + + (b.solid_state_ref._params.mw_comp['Fe3O4'] / + b.solid_state_ref._params.mw_comp['Fe2O3']) * + (b._params.rate_reaction_stoichiometry + ['R1', 'Sol', 'Fe3O4'] + / -b._params.rate_reaction_stoichiometry + ['R1', 'Sol', 'Fe2O3']) * + b.solid_state_ref.mass_frac_comp['Fe2O3']) == \ + 1e6 * b.solid_state_ref.mass_frac_comp['Fe3O4'] + try: + # Try to build constraint + self.OC_conv_eqn = Constraint(rule=OC_conv_eqn) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.OC_conv) + self.del_component(self.OC_conv_eqn) + + # Conversion of oxygen carrier reformulated + def _OC_conv_temp(self): + self.OC_conv_temp = Var(domain=Reals, initialize=1.0, + doc='Reformulation term for' + 'X to help eqn scaling') + + def OC_conv_temp_eqn(b): + return 1e3*b.OC_conv_temp**3 == 1e3*(1-b.OC_conv)**2 + try: + # Try to build constraint + self.OC_conv_temp_eqn = Constraint(rule=OC_conv_temp_eqn) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.OC_conv_temp) + self.del_component(self.OC_conv_temp_eqn) + + # General rate of reaction method + def _reaction_rate(self): + self.reaction_rate = Var(self._params.rate_reaction_idx, + domain=Reals, + initialize=0, + doc="Gen. rate of reaction [mol_rxn/m3.s]") + + def rate_rule(b, r): + return b.reaction_rate[r]*1e4 == b._params._scale_factor_rxn*1e4*( + b.solid_state_ref.mass_frac_comp['Fe2O3'] * + (1 - b.solid_state_ref._params.particle_porosity) * + b.solid_state_ref.dens_mass_skeletal * + (b._params.a_vol / + (b.solid_state_ref._params.mw_comp['Fe2O3'])) * + 3*b._params.rxn_stoich_coeff[r]*b.k_rxn[r] * + (((b.gas_state_ref.dens_mol_comp['CH4']**2 + + b._params.eps**2)**0.5) ** + b._params.rxn_order[r]) * + b.OC_conv_temp/(b._params.dens_mol_sol * + b._params.grain_radius) / + (-b._params.rate_reaction_stoichiometry['R1', 'Sol', 'Fe2O3'])) + try: + # Try to build constraint + self.gen_rate_expression = Constraint( + self._params.rate_reaction_idx, + rule=rate_rule) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.reaction_rate) + self.del_component(self.gen_rate_expression) + raise + + def get_reaction_rate_basis(b): + return MaterialFlowBasis.molar + + def model_check(blk): + """ + Model checks for property block + """ + # Check temperature bounds + if value(blk.temperature) < blk.temperature.lb: + _log.error('{} Temperature set below lower bound.'.format(blk.name) + ) + if value(blk.temperature) > blk.temperature.ub: + _log.error('{} Temperature set above upper bound.'.format(blk.name) + ) diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/solid_phase_thermo.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/solid_phase_thermo.py new file mode 100644 index 00000000..2ee5f0db --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/solid_phase_thermo.py @@ -0,0 +1,580 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the +# software owners: The Regents of the University of California, through +# Lawrence Berkeley National Laboratory, National Technology & Engineering +# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia +# University Research Corporation, et al. All rights reserved. +# +# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +This package provides the necessary constraints for solid phase properties of +an iron-based oxygen carrier +Components - Fe2O3, Fe3O4, Al2O3 + +Equations written in this model were primarily derived from: +National Institute of Standards and Technology, NIST Chemistry WebBook, +https://webbook.nist.gov/chemistry/ (accessed March 10, 2018). + +""" + +# Import Pyomo libraries +from pyomo.environ import (Constraint, + Param, + Reals, + value, + Var) +import pyomo.environ as pyo +from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.opt import SolverFactory + +# Import IDAES cores +from idaes.core import (declare_process_block_class, + MaterialFlowBasis, + PhysicalParameterBlock, + StateBlockData, + StateBlock, + Component, + SolidPhase) +from idaes.core.util.initialization import (fix_state_vars, + revert_state_vars, + solve_indexed_blocks) +from idaes.core.util.misc import add_object_reference +from idaes.core.util.model_statistics import ( + degrees_of_freedom, + number_unfixed_variables_in_activated_equalities) +import idaes.logger as idaeslog + +# Some more information about this module +__author__ = "Chinedu Okoli" + + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +@declare_process_block_class("SolidPhaseThermoParameterBlock") +class PhysicalParameterData(PhysicalParameterBlock): + """ + Property Parameter Block Class + + Contains parameters and indexing sets associated with properties for + methane CLC. + """ + + def build(self): + ''' + Callable method for Block construction. + ''' + super(PhysicalParameterData, self).build() + + self._state_block_class = SolidPhaseThermoStateBlock + + # Create Phase object + self.Sol = SolidPhase() + + # Create Component objects + self.Fe2O3 = Component() + self.Fe3O4 = Component() + self.Al2O3 = Component() + + # ------------------------------------------------------------------------- + """ Pure solid component properties""" + + # Mol. weights of solid components - units = kg/mol. ref: NIST webbook + mw_comp_dict = {'Fe2O3': 0.15969, 'Fe3O4': 0.231533, 'Al2O3': 0.10196} + self.mw_comp = Param( + self.component_list, + mutable=False, + initialize=mw_comp_dict, + doc="Molecular weights of solid components [kg/mol]") + + # Skeletal density of solid components - units = kg/m3. ref: NIST + dens_mass_comp_skeletal_dict = { + 'Fe2O3': 5250, 'Fe3O4': 5000, 'Al2O3': 3987} + self.dens_mass_comp_skeletal = Param( + self.component_list, + mutable=False, + initialize=dens_mass_comp_skeletal_dict, + doc='Skeletal density of solid components' + '[kg/m3]') + + # Ideal gas spec. heat capacity parameters(Shomate) of + # components - ref: NIST webbook. Shomate equations from NIST. + # Parameters A-E are used for cp calcs while A-H are used for enthalpy + # calc. + # 1e3*cp_comp = A + B*T + C*T^2 + D*T^3 + E/(T^2) + # where T = Temperature (K)/1000, and cp_comp = (kJ/mol.K) + # H_comp = H - H(298.15) = A*T + B*T^2/2 + C*T^3/3 + + # D*T^4/4 - E/T + F - H where T = Temp (K)/1000 and H_comp = (kJ/mol) + cp_param_dict = { + ('Al2O3', 1): 102.4290, + ('Al2O3', 2): 38.74980, + ('Al2O3', 3): -15.91090, + ('Al2O3', 4): 2.628181, + ('Al2O3', 5): -3.007551, + ('Al2O3', 6): -1717.930, + ('Al2O3', 7): 146.9970, + ('Al2O3', 8): -1675.690, + ('Fe3O4', 1): 200.8320000, + ('Fe3O4', 2): 1.586435e-7, + ('Fe3O4', 3): -6.661682e-8, + ('Fe3O4', 4): 9.452452e-9, + ('Fe3O4', 5): 3.18602e-8, + ('Fe3O4', 6): -1174.1350000, + ('Fe3O4', 7): 388.0790000, + ('Fe3O4', 8): -1120.8940000, + ('Fe2O3', 1): 110.9362000, + ('Fe2O3', 2): 32.0471400, + ('Fe2O3', 3): -9.1923330, + ('Fe2O3', 4): 0.9015060, + ('Fe2O3', 5): 5.4336770, + ('Fe2O3', 6): -843.1471000, + ('Fe2O3', 7): 228.3548000, + ('Fe2O3', 8): -825.5032000} + self.cp_param = Param(self.component_list, + range(1, 10), + mutable=False, + initialize=cp_param_dict, + doc="Shomate equation heat capacity parameters") + + # Std. heat of formation of comp. - units = kJ/(mol comp) - ref: NIST + enth_mol_form_comp_dict = {'Fe2O3': -825.5032, 'Fe3O4': -1120.894, + 'Al2O3': -1675.690} + self.enth_mol_form_comp = Param( + self.component_list, + mutable=False, + initialize=enth_mol_form_comp_dict, + doc="Component molar heats of formation [kJ/mol]") + + # ------------------------------------------------------------------------- + """ Mixed solid properties""" + # These are setup as fixed vars to allow for parameter estimation + + # Particle size + self.particle_dia = Var(domain=Reals, + initialize=1.5e-3, + doc='Diameter of solid particles [m]') + self.particle_dia.fix() + + # Particle porosity: + # The porosity of the OC particle is assumed to be a known parameter, + # and it is calculated from the known bulk density of the fresh OC + # particle (3251.75 kg/m3), and the known skeletal density of the + # fresh OC particle (calculated from the known composition of the + # fresh particle, and the skeletal density of its components) + self.particle_porosity = Var(domain=Reals, + initialize=1.5e-3, + doc='Porosity of oxygen carrier [-]') + self.particle_porosity.fix() + + # TODO -provide reference + # Minimum fluidization velocity - EPAT value used for Davidson model + self.velocity_mf = Var(domain=Reals, + initialize=0.039624, + doc='Velocity at minimum fluidization [m/s]') + self.velocity_mf.fix() + + # Minimum fluidization voidage - educated guess as rough + # estimate from ergun equation results (0.4) are suspicious + self.voidage_mf = Var(domain=Reals, + initialize=0.45, + doc='Voidage at minimum fluidization [-]') + self.voidage_mf.fix() + + # Particle thermal conductivity + self.therm_cond_sol = Var(domain=Reals, + initialize=12.3e-3, + doc='Thermal conductivity of solid' + 'particles [kJ/m.K.s]') + self.therm_cond_sol.fix() + + @classmethod + def define_metadata(cls, obj): + obj.define_custom_properties( + { + "particle_porosity": {"method": None, "units": pyo.units.dimensionless}, + "dens_mass_skeletal": { + "method": "_dens_mass_skeletal", + "units": obj.derived_units.DENSITY_MASS, + }, + "dens_mass_particle": { + "method": "_dens_mass_particle", + "units": obj.derived_units.DENSITY_MASS, + }, + } + ) + obj.add_properties({ + 'flow_mass': {'method': None, 'units': 'kg/s'}, + 'temperature': {'method': None, 'units': 'K'}, + 'mass_frac_comp': {'method': None, 'units': None}, + 'dens_mass_skeletal': {'method': '_dens_mass_skeletal', + 'units': 'kg/m3'}, + 'dens_mass_particle': {'method': '_dens_mass_particle', + 'units': 'kg/m3'}, + 'cp_mol_comp': {'method': '_cp_mol_comp', + 'units': 'kJ/mol.K'}, + 'cp_mass': {'method': '_cp_mass', 'units': 'kJ/kg.K'}, + 'enth_mass': {'method': '_enth_mass', 'units': 'kJ/kg'}, + 'enth_mol_comp': {'method': '_enth_mol_comp', + 'units': 'kJ/mol'}}) + + obj.add_default_units({'time': pyo.units.s, + 'length': pyo.units.m, + 'mass': pyo.units.kg, + 'amount': pyo.units.mol, + 'temperature': pyo.units.K, + #'energy': pyo.units.kJ, + #'holdup': pyo.units.kg, + }) + + +class _SolidPhaseThermoStateBlock(StateBlock): + """ + This Class contains methods which should be applied to Property Blocks as a + whole, rather than individual elements of indexed Property Blocks. + """ + def initialize(blk, state_args=None, hold_state=False, + state_vars_fixed=False, outlvl=idaeslog.NOTSET, + solver="ipopt", optarg={"tol": 1e-8}): + """ + Initialization routine for property package. + Keyword Arguments: + state_args : Dictionary with initial guesses for the state vars + chosen. Note that if this method is triggered + through the control volume, and if initial guesses + were not provided at the unit model level, the + control volume passes the inlet values as initial + guess. + Keys for the state_args dictionary are: + flow_mass, temperature, and mass_frac_comp + outlvl : sets output level of initialization routine + optarg : solver options dictionary object (default=None) + solver : str indicating which solver to use during + initialization (default = "ipopt") + hold_state : flag indicating whether the initialization routine + should unfix any state variables fixed during + initialization (default=False). + - True - states variables are not unfixed, and + a dict of returned containing flags for + which states were fixed during + initialization. + - False - state variables are unfixed after + initialization by calling the + release_state method + Returns: + If hold_states is True, returns a dict containing flags for + which states were fixed during initialization. + """ + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties") + solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="properties") + + init_log.info_high('Starting initialization') + + # Deactivate the constraints specific for outlet block i.e. + # when defined state is False + for k in blk.keys(): + if blk[k].config.defined_state is False: + blk[k].sum_component_eqn.deactivate() + + # Fix state variables if not already fixed + if state_vars_fixed is False: + flags = fix_state_vars(blk, state_args) + else: + # Check when the state vars are fixed already result in dof 0 + for k in blk.keys(): + if degrees_of_freedom(blk[k]) != 0: + raise Exception("State vars fixed but degrees of freedom " + "for state block is not zero during " + "initialization.") + + # Set solver options + opt = SolverFactory(solver) + opt.options = optarg + + # --------------------------------------------------------------------- + # Initialise values + for k in blk.keys(): + if hasattr(blk[k], "density_skeletal_constraint"): + calculate_variable_from_constraint( + blk[k].dens_mass_skeletal, + blk[k].density_skeletal_constraint) + + if hasattr(blk[k], "mixture_heat_capacity_eqn"): + calculate_variable_from_constraint( + blk[k].cp_mass, + blk[k].mixture_heat_capacity_eqn) + + if hasattr(blk[k], "mixture_enthalpy_eqn"): + calculate_variable_from_constraint( + blk[k].enth_mass, + blk[k].mixture_enthalpy_eqn) + + for j in blk[k]._params.component_list: + + if hasattr(blk[k], "cp_shomate_eqn"): + calculate_variable_from_constraint(blk[k].cp_mol_comp[j], + blk[k].cp_shomate_eqn[j] + ) + + if hasattr(blk[k], "enthalpy_shomate_eqn"): + calculate_variable_from_constraint( + blk[k].enth_mol_comp[j], + blk[k].enthalpy_shomate_eqn[j]) + + # Solve property block if non-empty + free_vars = 0 + for k in blk.keys(): + free_vars += number_unfixed_variables_in_activated_equalities( + blk[k]) + + if free_vars > 0: + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solve_indexed_blocks(opt, [blk], tee=slc.tee) + else: + res = "" + init_log.info_high("Initialization complete {}.".format( + idaeslog.condition(res)) + ) + + # --------------------------------------------------------------------- + if state_vars_fixed is False: + if hold_state is True: + return flags + else: + blk.release_state(flags) + + def release_state(blk, flags, outlvl=0): + """ + Method to release state variables fixed during initialization. + Keyword Arguments: + flags : dict containing information of which state variables + were fixed during initialization, and should now be + unfixed. This dict is returned by initialize if + hold_state=True. + outlvl : sets output level of of logging + """ + if flags is None: + return + + # Unfix state variables + revert_state_vars(blk, flags) + + # Activate state variable related constraints + for k in blk.keys(): + if blk[k].config.defined_state is False: + blk[k].sum_component_eqn.activate() + + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties") + init_log.info_high('States released.') + + +@declare_process_block_class("SolidPhaseThermoStateBlock", + block_class=_SolidPhaseThermoStateBlock) +class SolidPhaseThermoStateBlockData(StateBlockData): + """ + Property package for gas phase properties of methane combustion in CLC FR + """ + + def build(self): + """ + Callable method for Block construction + """ + super(SolidPhaseThermoStateBlockData, self).build() + + # Object reference for molecular weight if needed by CV1D + # Molecular weights + add_object_reference(self, "mw_comp", + self.config.parameters.mw_comp) + + self._make_state_vars() + + def _make_state_vars(self): + """List the necessary state variable objects.""" + self.flow_mass = Var(initialize=1.0, + domain=Reals, + doc='Component mass flowrate [kg/s]') + self.mass_frac_comp = Var( + self._params.component_list, + initialize=1 / len(self._params.component_list), + doc='State component mass fractions [-]') + self.temperature = Var(initialize=298.15, + domain=Reals, + doc='State temperature [K]') + + # Create standard constraints + # Sum mass fractions if not inlet block + if self.config.defined_state is False: + def sum_component_eqn(b): + return 1e2 == 1e2 * sum(b.mass_frac_comp[j] + for j in b._params.component_list) + self.sum_component_eqn = Constraint(rule=sum_component_eqn) + + def _dens_mass_skeletal(self): + # Skeletal density of OC solid particles + self.dens_mass_skeletal = Var(domain=Reals, + initialize=3251.75, + doc='Skeletal density of OC' + '[kg/m3]') + + def density_skeletal_constraint(b): + return (b.dens_mass_skeletal * sum( + b.mass_frac_comp[j] / + b._params.dens_mass_comp_skeletal[j] + for j in b._params.component_list) == + 1) + try: + # Try to build constraint + self.density_skeletal_constraint = Constraint( + rule=density_skeletal_constraint) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.dens_mass_skeletal) + self.del_component(self.density_skeletal_constraint) + raise + + def _dens_mass_particle(self): + # Particle density of OC (includes the OC pores) + self.dens_mass_particle = Var(domain=Reals, + initialize=3251.75, + doc='Particle density of oxygen carrier ' + '[kg/m3]') + + def density_particle_constraint(b): + return (b.dens_mass_particle == (1 - b._params.particle_porosity) * + b.dens_mass_skeletal) + try: + # Try to build constraint + self.density_particle_constraint = Constraint( + rule=density_particle_constraint) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.dens_mass_particle) + self.del_component(self.density_particle_constraint) + raise + + def _cp_mol_comp(self): + # Pure component solid heat capacities + self.cp_mol_comp = Var(self._params.component_list, + domain=Reals, + initialize=1.0, + doc="Pure component solid heat capacities " + "[kJ/mol.K]") + + def pure_component_cp_mol(b, j): + return b.cp_mol_comp[j] == 1e-3*( + b._params.cp_param[j, 1] + + b._params.cp_param[j, 2]*(b.temperature*1e-3) + + b._params.cp_param[j, 3]*(b.temperature*1e-3)**2 + + b._params.cp_param[j, 4]*(b.temperature*1e-3)**3 + + b._params.cp_param[j, 5]/((b.temperature*1e-3)**2)) + try: + # Try to build constraint + self.cp_shomate_eqn = Constraint(self._params.component_list, + rule=pure_component_cp_mol) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.cp_mol_comp) + self.del_component(self.cp_shomate_eqn) + raise + + def _cp_mass(self): + # Mixture heat capacities + self.cp_mass = Var(domain=Reals, + initialize=1.0, + doc="Mixture heat capacity, mass-basis [kJ/kg.K]") + + def cp_mass(b): + return b.cp_mass == sum(b.cp_mol_comp[j]*b.mass_frac_comp[j] + * (1/b._params.mw_comp[j]) + for j in b._params.component_list) + try: + # Try to build constraint + self.mixture_heat_capacity_eqn = Constraint(rule=cp_mass) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.cp_mass) + self.del_component(self.mixture_heat_capacity_eqn) + raise + + def _enth_mol_comp(self): + # Pure component vapour enthalpies + self.enth_mol_comp = Var( + self._params.component_list, + domain=Reals, + initialize=1.0, + doc="Pure component enthalpies [kJ/mol]") + + def pure_comp_enthalpy(b, j): + return b.enth_mol_comp[j] == ( + b._params.cp_param[j, 1]*(b.temperature*1e-3) + + b._params.cp_param[j, 2]*((b.temperature*1e-3)**2)/2 + + b._params.cp_param[j, 3]*((b.temperature*1e-3)**3)/3 + + b._params.cp_param[j, 4]*((b.temperature*1e-3)**4)/4 - + b._params.cp_param[j, 5]/(b.temperature*1e-3) + + b._params.cp_param[j, 6] - + b._params.cp_param[j, 8]) + try: + # Try to build constraint + self.enthalpy_shomate_eqn = Constraint(self._params.component_list, + rule=pure_comp_enthalpy) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.enth_mol_comp) + self.del_component(self.enthalpy_shomate_eqn) + raise + + def _enth_mass(self): + # Mixture mass enthalpy + self.enth_mass = Var(domain=Reals, + initialize=0.0, + doc='Mixture specific enthalpy [kJ/kg]') + try: + # Try to build constraint + self.mixture_enthalpy_eqn = Constraint(expr=( + self.enth_mass == sum( + self.mass_frac_comp[j] * + self.enth_mol_comp[j] + * (1/self._params.mw_comp[j]) + for j in self._params.component_list + ))) + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.enth_mass) + self.del_component(self.mixture_enthalpy_eqn) + raise + + def get_material_flow_terms(b, p, j): + return b.flow_mass*b.mass_frac_comp[j] + + def get_enthalpy_flow_terms(b, p): + return b.flow_mass*b.enth_mass + + def get_material_density_terms(b, p, j): + return b.dens_mass_particle * b.mass_frac_comp[j] + + def get_energy_density_terms(b, p): + return b.dens_mass_particle * b.enth_mass + + def define_state_vars(b): + return {"flow_mass": b.flow_mass, + "temperature": b.temperature, + "mass_frac_comp": b.mass_frac_comp} + + def get_material_flow_basis(b): + return MaterialFlowBasis.mass + + def model_check(blk): + """ + Model checks for property block + """ + # Check temperature bounds + if value(blk.temperature) < blk.temperature.lb: + _log.error('{} Temperature set below lower bound.' + .format(blk.name)) + if value(blk.temperature) > blk.temperature.ub: + _log.error('{} Temperature set above upper bound.' + .format(blk.name)) diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/__init__.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/test_CLC_gas_prop.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/test_CLC_gas_prop.py new file mode 100644 index 00000000..635b8b12 --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/test_CLC_gas_prop.py @@ -0,0 +1,70 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the +# software owners: The Regents of the University of California, through +# Lawrence Berkeley National Laboratory, National Technology & Engineering +# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia +# University Research Corporation, et al. All rights reserved. +# +# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +Tests for CLC gas phase thermo state block; tests for construction and solve +Author: Chinedu Okoli +""" + +import pytest +from pyomo.environ import (ConcreteModel, Var) +from idaes.core import FlowsheetBlock +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.testing import initialization_tester +from idaes.core.solvers import get_solver +from idaes_examples.mod.diagnostics.gas_solid_contactors.properties.methane_iron_OC_reduction. \ + gas_phase_thermo import GasPhaseThermoParameterBlock + +# Get default solver for testing +solver = get_solver() + + +# ----------------------------------------------------------------------------- +@pytest.fixture(scope="class") +def gas_prop(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + # gas properties and state inlet block + m.fs.properties = GasPhaseThermoParameterBlock() + m.fs.unit = m.fs.properties.build_state_block( + parameters=m.fs.properties, defined_state=True + ) + + m.fs.unit.flow_mol.fix(1) + m.fs.unit.temperature.fix(450) + m.fs.unit.pressure.fix(1.60) + m.fs.unit.mole_frac_comp["CO2"].fix(0.4772) + m.fs.unit.mole_frac_comp["H2O"].fix(0.0646) + m.fs.unit.mole_frac_comp["CH4"].fix(0.4582) + + return m + + +def test_build_inlet_state_block(gas_prop): + assert isinstance(gas_prop.fs.unit.mw, Var) + assert isinstance(gas_prop.fs.unit.dens_mol, Var) + assert isinstance(gas_prop.fs.unit.dens_mol_comp, Var) + assert isinstance(gas_prop.fs.unit.dens_mass, Var) + assert isinstance(gas_prop.fs.unit.cp_mol_comp, Var) + assert isinstance(gas_prop.fs.unit.cp_mol, Var) + assert isinstance(gas_prop.fs.unit.cp_mass, Var) + assert isinstance(gas_prop.fs.unit.visc_d, Var) + + +def test_setInputs_state_block(gas_prop): + assert degrees_of_freedom(gas_prop.fs.unit) == 0 + + +def test_initialize(gas_prop): + initialization_tester( + gas_prop) diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/test_CLC_rxn_prop.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/test_CLC_rxn_prop.py new file mode 100644 index 00000000..da8ec67a --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/test_CLC_rxn_prop.py @@ -0,0 +1,91 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the +# software owners: The Regents of the University of California, through +# Lawrence Berkeley National Laboratory, National Technology & Engineering +# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia +# University Research Corporation, et al. All rights reserved. +# +# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +Tests for CLC heterogeneous reaction block; tests for construction and solve +Author: Chinedu Okoli +""" + +import pytest +from pyomo.environ import ConcreteModel, Var +from idaes.core import FlowsheetBlock +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.testing import initialization_tester +from idaes.core.solvers import get_solver + +from idaes_examples.mod.diagnostics.gas_solid_contactors.properties.methane_iron_OC_reduction. \ + gas_phase_thermo import GasPhaseThermoParameterBlock +from idaes_examples.mod.diagnostics.gas_solid_contactors.properties.methane_iron_OC_reduction. \ + solid_phase_thermo import SolidPhaseThermoParameterBlock +from idaes_examples.mod.diagnostics.gas_solid_contactors.properties.methane_iron_OC_reduction. \ + hetero_reactions import HeteroReactionParameterBlock + + +# Get default solver for testing +solver = get_solver() + + +# ----------------------------------------------------------------------------- +@pytest.fixture(scope="class") +def rxn_prop(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + # Set up thermo props and reaction props + m.fs.solid_properties = SolidPhaseThermoParameterBlock() + m.fs.solid_state_block = m.fs.solid_properties.build_state_block( + parameters=m.fs.solid_properties, defined_state=True + ) + + m.fs.gas_properties = GasPhaseThermoParameterBlock() + m.fs.gas_state_block = m.fs.gas_properties.build_state_block( + parameters=m.fs.gas_properties, defined_state=True + ) + + m.fs.reactions = HeteroReactionParameterBlock( + solid_property_package=m.fs.solid_properties, + gas_property_package=m.fs.gas_properties, + ) + m.fs.unit = m.fs.reactions.reaction_block_class( + parameters=m.fs.reactions, + solid_state_block=m.fs.solid_state_block, + gas_state_block=m.fs.gas_state_block, + has_equilibrium=False, + ) + + # Fix required variables to make reaction model square + # (gas mixture and component densities, + # solid density and component fractions) + m.fs.solid_state_block.mass_frac_comp["Fe2O3"].fix(0.45) + m.fs.solid_state_block.mass_frac_comp["Fe3O4"].fix(1e-9) + m.fs.solid_state_block.mass_frac_comp["Al2O3"].fix(0.55) + m.fs.gas_state_block.dens_mol.fix(10) + m.fs.gas_state_block.dens_mol_comp.fix(10) + m.fs.solid_state_block.dens_mass_skeletal.fix(1) + + return m + + +def test_build_reaction_block(rxn_prop): + assert isinstance(rxn_prop.fs.unit.k_rxn, Var) + assert isinstance(rxn_prop.fs.unit.OC_conv, Var) + assert isinstance(rxn_prop.fs.unit.OC_conv_temp, Var) + assert isinstance(rxn_prop.fs.unit.reaction_rate, Var) + + +def test_setInputs_reaction_block(rxn_prop): + assert degrees_of_freedom(rxn_prop.fs.unit) == 0 + + +def test_initialize(rxn_prop): + initialization_tester( + rxn_prop) diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/test_CLC_solid_prop.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/test_CLC_solid_prop.py new file mode 100644 index 00000000..46bb6ebc --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/properties/methane_iron_OC_reduction/tests/test_CLC_solid_prop.py @@ -0,0 +1,69 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the +# software owners: The Regents of the University of California, through +# Lawrence Berkeley National Laboratory, National Technology & Engineering +# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia +# University Research Corporation, et al. All rights reserved. +# +# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +Tests for CLC solid phase thermo state block; tests for construction and solve +Author: Chinedu Okoli +""" + +import pytest + +from pyomo.environ import ConcreteModel, Var +from idaes.core import FlowsheetBlock +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.testing import initialization_tester +from idaes.core.solvers import get_solver +from idaes_examples.mod.diagnostics.gas_solid_contactors.properties.methane_iron_OC_reduction. \ + solid_phase_thermo import SolidPhaseThermoParameterBlock + +# Get default solver for testing +solver = get_solver() + + +# ----------------------------------------------------------------------------- +@pytest.fixture(scope="class") +def solid_prop(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + # solid properties and state inlet block + m.fs.properties = SolidPhaseThermoParameterBlock() + + m.fs.unit = m.fs.properties.build_state_block( + parameters=m.fs.properties, + defined_state=True, + ) + + m.fs.unit.flow_mass.fix(1) + m.fs.unit.temperature.fix(1183.15) + m.fs.unit.mass_frac_comp["Fe2O3"].fix(0.45) + m.fs.unit.mass_frac_comp["Fe3O4"].fix(1e-9) + m.fs.unit.mass_frac_comp["Al2O3"].fix(0.55) + + return m + + +def test_build_inlet_state_block(solid_prop): + assert isinstance(solid_prop.fs.unit.dens_mass_skeletal, Var) + assert isinstance(solid_prop.fs.unit.enth_mol_comp, Var) + assert isinstance(solid_prop.fs.unit.enth_mass, Var) + assert isinstance(solid_prop.fs.unit.cp_mol_comp, Var) + assert isinstance(solid_prop.fs.unit.cp_mass, Var) + + +def test_setInputs_state_block(solid_prop): + assert degrees_of_freedom(solid_prop.fs.unit) == 0 + + +def test_initialize(solid_prop): + initialization_tester( + solid_prop) diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/tests/__init__.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/tests/test_moving_bed_example.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/tests/test_moving_bed_example.py new file mode 100644 index 00000000..a0ac34b2 --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/tests/test_moving_bed_example.py @@ -0,0 +1,83 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the +# software owners: The Regents of the University of California, through +# Lawrence Berkeley National Laboratory, National Technology & Engineering +# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia +# University Research Corporation, et al. All rights reserved. +# +# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +import pytest +import idaes_examples.mod.diagnostics.gas_solid_contactors.example as example +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.model_diagnostics import check_parallel_jacobian +import pyomo.environ as pyo +from pyomo.contrib.incidence_analysis import IncidenceGraphInterface + + +class TestMovingBedExample: + + def test_original_model(self): + model = example.create_original_model() + assert degrees_of_freedom(model) == 10 + # Make sure the model doesn't solve + solver = pyo.SolverFactory("ipopt") + solver.options["max_iter"] = 10 + res = solver.solve(model) + assert not pyo.check_optimal_termination(res) + + def test_original_square_model(self): + model = example.create_original_square_model() + assert degrees_of_freedom(model) == 0 + igraph = IncidenceGraphInterface(model) + vdm, cdm = igraph.dulmage_mendelsohn() + # Make sure the model is structurally singular + assert len(vdm.unmatched) == 10 + assert len(cdm.unmatched) == 10 + + def test_model_with_new_variable_and_constraint(self): + # First attempt to fix the model + model = example.create_square_model_with_new_variable_and_constraint() + assert degrees_of_freedom(model) == 0 + igraph = IncidenceGraphInterface(model) + vdm, cdm = igraph.dulmage_mendelsohn() + # Make sure the model is *not* structurally singular + assert not vdm.unmatched + assert not cdm.unmatched + + # Make sure the model still doesn't solve + example.free_degrees_of_freedom(model) + solver = pyo.SolverFactory("ipopt") + solver.options["max_iter"] = 10 + res = solver.solve(model) + assert not pyo.check_optimal_termination(res) + example.fix_degrees_of_freedom(model) + + # Make sure the model is numerically singular + parallel = check_parallel_jacobian(model, direction="row", tolerance=1e-8) + assert len(parallel) == 11 + + def test_corrected_model(self): + model = example.create_corrected_square_model() + + igraph = IncidenceGraphInterface(model) + vdm, cdm = igraph.dulmage_mendelsohn() + # Make sure the model is *not* structurally singular + assert not vdm.unmatched + assert not cdm.unmatched + + example.free_degrees_of_freedom(model) + assert degrees_of_freedom(model) == 10 + + # Make sure the model solves + solver = pyo.SolverFactory("ipopt") + solver.options["max_iter"] = 10 + res = solver.solve(model) + assert pyo.check_optimal_termination(res) + + +if __name__ == "__main__": + pytest.main([__file__]) diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/unit_models/__init__.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/unit_models/__init__.py new file mode 100644 index 00000000..0ce5f8e1 --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/unit_models/__init__.py @@ -0,0 +1 @@ +from .moving_bed import MBR diff --git a/idaes_examples/mod/diagnostics/gas_solid_contactors/unit_models/moving_bed.py b/idaes_examples/mod/diagnostics/gas_solid_contactors/unit_models/moving_bed.py new file mode 100644 index 00000000..a242e65e --- /dev/null +++ b/idaes_examples/mod/diagnostics/gas_solid_contactors/unit_models/moving_bed.py @@ -0,0 +1,1612 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the +# software owners: The Regents of the University of California, through +# Lawrence Berkeley National Laboratory, National Technology & Engineering +# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia +# University Research Corporation, et al. All rights reserved. +# +# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +IDAES Moving Bed Model. + +The moving bed model is a 1D axially discretized model with a gas and +solid phase and a counter-current flow direction. The model captures +the gas-solid interaction between both phases through reaction, mass +and heat transfer. + +Equations written in this model were derived from: +A. Ostace, A. Lee, C.O. Okoli, A.P. Burgard, D.C. Miller, D. Bhattacharyya, +Mathematical modeling of a moving-bed reactor for chemical looping combustion +of methane, in: M.R. Eden, M. Ierapetritou, G.P. Towler (Eds.),13th Int. Symp. +Process Syst. Eng. (PSE 2018), Computer-Aided Chemical Engineering 2018, +pp. 325–330 , San Diego, CA. + +Assumptions: +Property package contains temperature and pressure variables. +Property package contains minimum fluidization velocity. + +""" +from __future__ import division + +# Import Python libraries +import matplotlib.pyplot as plt + +# Import Pyomo libraries +from pyomo.environ import (SolverFactory, Var, Param, Reals, value, + TransformationFactory, Constraint, + TerminationCondition) +from pyomo.core.base.block import BlockData +from pyomo.common.config import ConfigBlock, ConfigValue, In +from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.dae import ContinuousSet + +# Import IDAES cores +from idaes.core import (ControlVolume1DBlock, + UnitModelBlockData, + declare_process_block_class, + MaterialBalanceType, + EnergyBalanceType, + MomentumBalanceType, + FlowDirection) +from idaes.core.util.config import (is_physical_parameter_block, + is_reaction_parameter_block) +from idaes.core.util.exceptions import (ConfigurationError, + BurntToast) +from idaes.core.util.tables import create_stream_table_dataframe +from idaes.core.base.control_volume1d import DistributedVars +from idaes.core.util.constants import Constants as constants +from idaes.core.util.math import smooth_abs +import idaes.logger as idaeslog + +__author__ = "Chinedu Okoli", "Anca Ostace" + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +@declare_process_block_class("MBR") +class MBRData(UnitModelBlockData): + """Standard Moving Bed Unit Model Class.""" + + # Create template for unit level config arguments + CONFIG = UnitModelBlockData.CONFIG() + + # Unit level config arguments + CONFIG.declare("finite_elements", ConfigValue( + default=10, + domain=int, + description="Number of finite elements length domain", + doc="""Number of finite elements to use when discretizing length +domain (default=20)""")) + CONFIG.declare("length_domain_set", ConfigValue( + default=[0.0, 1.0], + domain=list, + description="Number of finite elements length domain", + doc="""length_domain_set - (optional) list of point to use to +initialize a new ContinuousSet if length_domain is not +provided (default = [0.0, 1.0])""")) + CONFIG.declare("transformation_method", ConfigValue( + default="dae.finite_difference", + description="Method to use for DAE transformation", + doc="""Method to use to transform domain. Must be a method recognised +by the Pyomo TransformationFactory, +**default** - "dae.finite_difference". +**Valid values:** { +**"dae.finite_difference"** - Use a finite difference transformation method, +**"dae.collocation"** - use a collocation transformation method}""")) + CONFIG.declare("transformation_scheme", ConfigValue( + default=None, + domain=In([None, "BACKWARD", "FORWARD", "LAGRANGE-RADAU"]), + description="Scheme to use for DAE transformation", + doc="""Scheme to use when transforming domain. See Pyomo +documentation for supported schemes, +**default** - None. +**Valid values:** { +**None** - defaults to "BACKWARD" for finite difference transformation method, +and to "LAGRANGE-RADAU" for collocation transformation method, +**"BACKWARD"** - Use a finite difference transformation method, +**"FORWARD""** - use a finite difference transformation method, +**"LAGRANGE-RADAU""** - use a collocation transformation method}""")) + CONFIG.declare("collocation_points", ConfigValue( + default=3, + domain=int, + description="Number of collocation points per finite element", + doc="""Number of collocation points to use per finite element when +discretizing length domain (default=3)""")) + CONFIG.declare("flow_type", ConfigValue( + default="counter_current", + domain=In(['counter_current']), + description="Flow configuration of Moving Bed", + doc="""Flow configuration of Moving Bed +- counter_current: gas side flows from 0 to 1 +solid side flows from 1 to 0""")) + CONFIG.declare("material_balance_type", ConfigValue( + default=MaterialBalanceType.componentTotal, + domain=In(MaterialBalanceType), + description="Material balance construction flag", + doc="""Indicates what type of mass balance should be constructed, +**default** - MaterialBalanceType.componentTotal. +**Valid values:** { +**MaterialBalanceType.none** - exclude material balances, +**MaterialBalanceType.componentPhase** - use phase component balances, +**MaterialBalanceType.componentTotal** - use total component balances, +**MaterialBalanceType.elementTotal** - use total element balances, +**MaterialBalanceType.total** - use total material balance.}""")) + CONFIG.declare("energy_balance_type", ConfigValue( + default=EnergyBalanceType.enthalpyTotal, + domain=In(EnergyBalanceType), + description="Energy balance construction flag", + doc="""Indicates what type of energy balance should be constructed, +**default** - EnergyBalanceType.enthalpyTotal. +**Valid values:** { +**EnergyBalanceType.none** - exclude energy balances, +**EnergyBalanceType.enthalpyTotal** - single enthalpy balance for material, +**EnergyBalanceType.enthalpyPhase** - enthalpy balances for each phase, +**EnergyBalanceType.energyTotal** - single energy balance for material, +**EnergyBalanceType.energyPhase** - energy balances for each phase.}""")) + CONFIG.declare("momentum_balance_type", ConfigValue( + default=MomentumBalanceType.pressureTotal, + domain=In(MomentumBalanceType), + description="Momentum balance construction flag", + doc="""Indicates what type of momentum balance should be constructed, +**default** - MomentumBalanceType.pressureTotal. +**Valid values:** { +**MomentumBalanceType.none** - exclude momentum balances, +**MomentumBalanceType.pressureTotal** - single pressure balance for material, +**MomentumBalanceType.pressurePhase** - pressure balances for each phase, +**MomentumBalanceType.momentumTotal** - single momentum balance for material, +**MomentumBalanceType.momentumPhase** - momentum balances for each phase.}""")) + CONFIG.declare("has_pressure_change", ConfigValue( + default=True, + domain=In([True, False]), + description="Pressure change term construction flag", + doc="""Indicates whether terms for pressure change should be +constructed, +**default** - False. +**Valid values:** { +**True** - include pressure change terms, +**False** - exclude pressure change terms.}""")) + CONFIG.declare("pressure_drop_type", ConfigValue( + default="simple_correlation", + domain=In(["simple_correlation", "ergun_correlation"]), + description="Construction flag for type of pressure drop", + doc="""Indicates what type of pressure drop correlation should be used, +**default** - "simple_correlation". +**Valid values:** { +**"simple_correlation"** - Use a simplified pressure drop correlation, +**"ergun_correlation"** - Use the ergun equation.}""")) + + # Create template for phase specific config arguments + _PhaseTemplate = UnitModelBlockData.CONFIG() + _PhaseTemplate.declare("has_equilibrium_reactions", ConfigValue( + default=False, + domain=In([True, False]), + description="Equilibrium reaction construction flag", + doc="""Indicates whether terms for equilibrium controlled reactions +should be constructed, +**default** - True. +**Valid values:** { +**True** - include equilibrium reaction terms, +**False** - exclude equilibrium reaction terms.}""")) + _PhaseTemplate.declare("property_package", ConfigValue( + default=None, + domain=is_physical_parameter_block, + description="Property package to use for control volume", + doc="""Property parameter object used to define property calculations +(default = 'use_parent_value') +- 'use_parent_value' - get package from parent (default = None) +- a ParameterBlock object""")) + _PhaseTemplate.declare("property_package_args", ConfigValue( + default={}, + domain=dict, + description="Arguments for constructing gas property package", + doc="""A dict of arguments to be passed to the PropertyBlockData +and used when constructing these +(default = 'use_parent_value') +- 'use_parent_value' - get package from parent (default = None) +- a dict (see property package for documentation)""")) + _PhaseTemplate.declare("reaction_package", ConfigValue( + default=None, + domain=is_reaction_parameter_block, + description="Reaction package to use for control volume", + doc="""Reaction parameter object used to define reaction calculations, +**default** - None. +**Valid values:** { +**None** - no reaction package, +**ReactionParameterBlock** - a ReactionParameterBlock object.}""")) + _PhaseTemplate.declare("reaction_package_args", ConfigBlock( + implicit=True, + implicit_domain=ConfigBlock, + description="Arguments to use for constructing reaction packages", + doc="""A ConfigBlock with arguments to be passed to a reaction block(s) +and used when constructing these, +**default** - None. +**Valid values:** { +see reaction package for documentation.}""")) + + # Create individual config blocks for gas and solid sides + CONFIG.declare("gas_phase_config", + _PhaseTemplate(doc="gas phase config arguments")) + CONFIG.declare("solid_phase_config", + _PhaseTemplate(doc="solid phase config arguments")) + + # ========================================================================= + def build(self): + """ + Begin building model (pre-DAE transformation). + + Args: + None + + Returns: + None + """ + # Call UnitModel.build to build default attributes + super(MBRData, self).build() + + # Consistency check for transformation method and transformation scheme + if (self.config.transformation_method == "dae.finite_difference" and + self.config.transformation_scheme is None): + self.config.transformation_scheme = "BACKWARD" + elif (self.config.transformation_method == "dae.collocation" and + self.config.transformation_scheme is None): + self.config.transformation_scheme = "LAGRANGE-RADAU" + elif (self.config.transformation_method == "dae.finite_difference" and + self.config.transformation_scheme != "BACKWARD" and + self.config.transformation_scheme != "FORWARD"): + raise ConfigurationError("{} invalid value for " + "transformation_scheme argument. " + "Must be ""BACKWARD"" or ""FORWARD"" " + "if transformation_method is" + " ""dae.finite_difference""." + .format(self.name)) + elif (self.config.transformation_method == "dae.collocation" and + self.config.transformation_scheme != "LAGRANGE-RADAU"): + raise ConfigurationError("{} invalid value for " + "transformation_scheme argument." + "Must be ""LAGRANGE-RADAU"" if " + "transformation_method is" + " ""dae.collocation""." + .format(self.name)) + + # Set flow directions for the control volume blocks + # Gas flows from 0 to 1, solid flows from 1 to 0 + # An if statement is used here despite only one option to allow for + # future extensions to other flow configurations + if self.config.flow_type == "counter_current": + set_direction_gas = FlowDirection.forward + set_direction_solid = FlowDirection.backward + else: + raise BurntToast( + "{} encountered unrecognized argument " + "for flow type. Please contact the IDAES" + " developers with this bug.".format(self.name)) + + # Set arguments for gas sides if homoogeneous reaction block + if self.config.gas_phase_config.reaction_package is not None: + has_rate_reaction_gas_phase = True + else: + has_rate_reaction_gas_phase = False + + # Set arguments for gas and solid sides if heterogeneous reaction block + if self.config.solid_phase_config.reaction_package is not None: + has_rate_reaction_solid_phase = True + has_mass_transfer_gas_phase = True + else: + has_rate_reaction_solid_phase = False + has_mass_transfer_gas_phase = False + + # Set heat transfer terms + if self.config.energy_balance_type != EnergyBalanceType.none: + has_heat_transfer = True + else: + has_heat_transfer = False + + # Set heat of reaction terms + if (self.config.energy_balance_type != EnergyBalanceType.none + and self.config.gas_phase_config.reaction_package is not None): + has_heat_of_reaction_gas_phase = True + else: + has_heat_of_reaction_gas_phase = False + + if (self.config.energy_balance_type != EnergyBalanceType.none + and self.config.solid_phase_config. + reaction_package is not None): + has_heat_of_reaction_solid_phase = True + else: + has_heat_of_reaction_solid_phase = False + + # Create a unit model length domain + self.length_domain = ContinuousSet( + bounds=(0.0, 1.0), + initialize=self.config.length_domain_set, + doc="Normalized length domain") + + # ========================================================================= + """ Build Control volume 1D for gas phase and + populate gas control volume""" + + self.gas_phase = ControlVolume1DBlock(**{ + "transformation_method": self.config.transformation_method, + "transformation_scheme": self.config.transformation_scheme, + "finite_elements": self.config.finite_elements, + "collocation_points": self.config.collocation_points, + "dynamic": self.config.dynamic, + "has_holdup": self.config.has_holdup, + "area_definition": DistributedVars.variant, + "property_package": self.config.gas_phase_config.property_package, + "property_package_args": + self.config.gas_phase_config.property_package_args, + "reaction_package": self.config.gas_phase_config.reaction_package, + "reaction_package_args": + self.config.gas_phase_config.reaction_package_args}) + + self.gas_phase.add_geometry( + length_domain=self.length_domain, + length_domain_set=self.config.length_domain_set, + flow_direction=set_direction_gas) + + self.gas_phase.add_state_blocks( + information_flow=set_direction_gas, + has_phase_equilibrium=False) + + if self.config.gas_phase_config.reaction_package is not None: + self.gas_phase.add_reaction_blocks( + has_equilibrium=self.config.gas_phase_config. + has_equilibrium_reactions) + + self.gas_phase.add_material_balances( + balance_type=self.config.material_balance_type, + has_phase_equilibrium=False, + has_mass_transfer=has_mass_transfer_gas_phase, + has_rate_reactions=has_rate_reaction_gas_phase) + + self.gas_phase.add_energy_balances( + balance_type=self.config.energy_balance_type, + has_heat_transfer=has_heat_transfer, + has_heat_of_reaction=has_heat_of_reaction_gas_phase) + + self.gas_phase.add_momentum_balances( + balance_type=self.config.momentum_balance_type, + has_pressure_change=self.config.has_pressure_change) + + # ========================================================================= + """ Build Control volume 1D for solid phase and + populate solid control volume""" + + # Set argument for heterogeneous reaction block + self.solid_phase = ControlVolume1DBlock(**{ + "transformation_method": self.config.transformation_method, + "transformation_scheme": self.config.transformation_scheme, + "finite_elements": self.config.finite_elements, + "collocation_points": self.config.collocation_points, + "dynamic": self.config.dynamic, + "has_holdup": self.config.has_holdup, + "area_definition": DistributedVars.variant, + "property_package": + self.config.solid_phase_config.property_package, + "property_package_args": + self.config.solid_phase_config.property_package_args, + "reaction_package": + self.config.solid_phase_config.reaction_package, + "reaction_package_args": + self.config.solid_phase_config.reaction_package_args}) + + self.solid_phase.add_geometry( + length_domain=self.length_domain, + length_domain_set=self.config.length_domain_set, + flow_direction=set_direction_solid) + + self.solid_phase.add_state_blocks( + information_flow=set_direction_solid, + has_phase_equilibrium=False) + + if self.config.solid_phase_config.reaction_package is not None: + # TODO - a generalization of the heterogeneous reaction block + # The heterogeneous reaction block does not use the + # add_reaction_blocks in control volumes as control volumes are + # currently setup to handle only homogeneous reaction properties. + # Thus appending the heterogeneous reaction block to the + # solid state block is currently hard coded here. + + tmp_dict = dict(**self.config.solid_phase_config. + reaction_package_args) + tmp_dict["gas_state_block"] = self.gas_phase.properties + tmp_dict["solid_state_block"] = self.solid_phase.properties + tmp_dict["has_equilibrium"] = (self.config.solid_phase_config. + has_equilibrium_reactions) + tmp_dict["parameters"] = (self.config.solid_phase_config. + reaction_package) + self.solid_phase.reactions = ( + self.config.solid_phase_config.reaction_package. + reaction_block_class( + self.flowsheet().config.time, + self.length_domain, + doc="Reaction properties in control volume", + **tmp_dict)) + + self.solid_phase.add_material_balances( + balance_type=self.config.material_balance_type, + has_phase_equilibrium=False, + has_mass_transfer=False, + has_rate_reactions=has_rate_reaction_solid_phase) + + self.solid_phase.add_energy_balances( + balance_type=self.config.energy_balance_type, + has_heat_transfer=has_heat_transfer, + has_heat_of_reaction=has_heat_of_reaction_solid_phase) + + self.solid_phase.add_momentum_balances( + balance_type=MomentumBalanceType.none, + has_pressure_change=False) + + # ========================================================================= + """ Add ports""" + # Add Ports for gas side + self.add_inlet_port(name="gas_inlet", block=self.gas_phase) + self.add_outlet_port(name="gas_outlet", block=self.gas_phase) + + # Add Ports for solid side + self.add_inlet_port(name="solid_inlet", block=self.solid_phase) + self.add_outlet_port(name="solid_outlet", block=self.solid_phase) + + # ========================================================================= + """ Add performance equation method""" + self._apply_transformation() + self._make_performance() + + # ========================================================================= + def _apply_transformation(self): + """ + Method to apply DAE transformation to the Control Volume length domain. + Transformation applied will be based on the Control Volume + configuration arguments. + """ + if self.config.finite_elements is None: + raise ConfigurationError( + "{} was not provided a value for the finite_elements" + " configuration argument. Please provide a valid value." + .format(self.name)) + + if self.config.transformation_method == "dae.finite_difference": + self.discretizer = TransformationFactory( + self.config.transformation_method) + self.discretizer.apply_to(self, + wrt=self.length_domain, + nfe=self.config.finite_elements, + scheme=self.config.transformation_scheme) + elif self.config.transformation_method == "dae.collocation": + self.discretizer = TransformationFactory( + self.config.transformation_method) + self.discretizer.apply_to( + self, + wrt=self.length_domain, + nfe=self.config.finite_elements, + ncp=self.config.collocation_points, + scheme=self.config.transformation_scheme) + + def _make_performance(self): + """ + Constraints for unit model. + + Args: + None + + Returns: + None + """ + # local aliases used to shorten object names + gas_phase = self.config.gas_phase_config + solid_phase = self.config.solid_phase_config + + # Declare Mutable Parameters + self.eps = Param(mutable=True, + default=1e-8, + doc='Smoothing Factor for Smooth IF Statements') + + # Unit Model variables + self.bed_diameter = Var(domain=Reals, + initialize=1, + doc='Reactor diameter [m]') + self.bed_area = Var(domain=Reals, + initialize=1, + doc='Reactor cross-sectional area [m2]') + self.bed_height = Var(domain=Reals, initialize=1, + doc='Bed length [m]') + + # Phase specific variables + self.velocity_superficial_gas = Var( + self.flowsheet().config.time, + self.length_domain, + domain=Reals, initialize=0.05, + doc='Gas superficial velocity [m/s]') + self.velocity_superficial_solid = Var( + self.flowsheet().config.time, + domain=Reals, initialize=0.005, + doc='Solid superficial velocity [m/s]') + + # Dimensionless numbers, mass and heat transfer coefficients + self.Re_particle = Var(self.flowsheet().config.time, + self.length_domain, + domain=Reals, initialize=1.0, + doc='Particle Reynolds number [-]') + + self.Pr = Var(self.flowsheet().config.time, + self.length_domain, + domain=Reals, initialize=1.0, + doc='Prandtl number of gas in bed [-]') + + self.Nu_particle = Var(self.flowsheet().config.time, + self.length_domain, + domain=Reals, initialize=1.0, + doc='Particle Nusselt number [-]') + self.gas_solid_htc = Var(self.flowsheet().config.time, + self.length_domain, + domain=Reals, initialize=1.0, + doc='Gas-solid heat transfer coefficient' + '[kJ/(m2Ks)]') + + # Fixed variables (these are parameters that can be estimated) + self.bed_voidage = Var(domain=Reals, + initialize=0.4, + doc="Bed voidage [-]") + self.bed_voidage.fix() + + # ========================================================================= + # Add performance equations + + # --------------------------------------------------------------------- + # Geometry constraints + + # Bed area + @self.Constraint(doc="Bed area") + def bed_area_eqn(b): + return b.bed_area == ( + constants.pi*(0.5*b.bed_diameter)**2) + + # Area of gas side, and solid side + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Gas side area") + def gas_phase_area(b, t, x): + return (b.gas_phase.area[t, x] == + b.bed_area*b.bed_voidage) + + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Solid side area") + def solid_phase_area(b, t, x): + return (b.solid_phase.area[t, x] == + b.bed_area*(1-b.bed_voidage)) + + # Length of gas side, and solid side + @self.Constraint(doc="Gas side length") + def gas_phase_length(b): + return (b.gas_phase.length == b.bed_height) + + @self.Constraint(doc="Solid side length") + def solid_phase_length(b): + return (b.solid_phase.length == b.bed_height) + + # --------------------------------------------------------------------- + # Hydrodynamic constraints + + # Gas superficial velocity + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Gas superficial velocity") + def gas_super_vel(b, t, x): + return (b.velocity_superficial_gas[t, x] * b.bed_area * + b.gas_phase.properties[t, x].dens_mol == + b.gas_phase.properties[t, x].flow_mol) + + # Solid superficial velocity + @self.Constraint(self.flowsheet().config.time, + doc="Solid superficial velocity") + def solid_super_vel(b, t): + return (b.velocity_superficial_solid[t] * b.bed_area * + b.solid_phase.properties[t, 1].dens_mass_particle == + b.solid_phase.properties[t, 1].flow_mass) + + # Gas side pressure drop calculation + if (self.config.has_pressure_change and + self.config.pressure_drop_type == + "simple_correlation"): + # Simplified pressure drop + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Gas side pressure drop calculation -" + "simplified pressure drop") + def gas_phase_config_pressure_drop(b, t, x): + return b.gas_phase.deltaP[t, x]*1e5 == -0.2*( + b.velocity_superficial_gas[t, x] * + (b.solid_phase.properties[t, x].dens_mass_particle - + b.gas_phase.properties[t, x].dens_mass)) + elif (self.config.has_pressure_change and + self.config.pressure_drop_type == "ergun_correlation"): + # Ergun equation + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Gas side pressure drop calculation -" + "ergun equation") + def gas_phase_config_pressure_drop(b, t, x): + return (1e2*-b.gas_phase.deltaP[t, x]*1e5 == + 1e2*( + 150*(1 - b.bed_voidage) ** 2 * + b.gas_phase.properties[t, x].visc_d * + (b.velocity_superficial_gas[t, x] + + b.velocity_superficial_solid[t]) / + (b.solid_phase.properties[t, x]. + _params.particle_dia ** 2 * b.bed_voidage ** 3)) + + 1e2*( + 1.75*b.gas_phase.properties[t, x].dens_mass * + (1 - b.bed_voidage) * + (b.velocity_superficial_gas[t, x] + + b.velocity_superficial_solid[t]) ** 2 / + (b.solid_phase.properties[t, x]._params.particle_dia * + b.bed_voidage**3))) + # The above expression has no absolute values - assumes: + # (velocity_superficial_gas + velocity_superficial_solid) > 0 + else: + raise BurntToast( + "{} encountered unrecognized argument for " + "the pressure drop correlation. Please contact the IDAES" + " developers with this bug.".format(self.name)) + # --------------------------------------------------------------------- + # Reaction constraints + + # Build homogeneous reaction constraints + if gas_phase.reaction_package is not None: + # Gas side rate reaction extent + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + gas_phase.reaction_package.rate_reaction_idx, + doc="Gas side rate reaction extent") + def gas_phase_config_rxn_ext(b, t, x, r): + return 1e3*b.gas_phase.rate_reaction_extent[t, x, r] == 1e3*( + b.gas_phase.reactions[t, x].reaction_rate[r] * + b.gas_phase.area[t, x]) + + # Build hetereogeneous reaction constraints + if solid_phase.reaction_package is not None: + # Solid side rate reaction extent + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + solid_phase.reaction_package.rate_reaction_idx, + doc="Solid side rate reaction extent") + def solid_phase_config_rxn_ext(b, t, x, r): + return 1e3*b.solid_phase.rate_reaction_extent[t, x, r] == 1e3*( + b.solid_phase.reactions[t, x].reaction_rate[r] * + b.solid_phase.area[t, x]) + + # Gas side heterogeneous rate reaction generation + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + gas_phase.property_package.phase_list, + gas_phase.property_package.component_list, + doc='Gas side heterogeneous' + 'rate reaction generation') + def gas_comp_hetero_rxn(b, t, x, p, j): + return 1e3*b.gas_phase.mass_transfer_term[t, x, p, j] == 1e3*( + sum(b.solid_phase.reactions[t, x]. + rate_reaction_stoichiometry[r, p, j] * + b.solid_phase.reactions[t, x].reaction_rate[r] + for r in ( + solid_phase.reaction_package.rate_reaction_idx) + ) * b.solid_phase.area[t, x]) + + # --------------------------------------------------------------------- + if self.config.energy_balance_type != EnergyBalanceType.none: + # Solid phase - gas to solid heat transfer + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Solid phase - gas to solid heat transfer") + def solid_phase_heat_transfer(b, t, x): + return (b.solid_phase.heat[t, x] * + b.solid_phase.properties[t, x]._params.particle_dia == + 6 * b.gas_solid_htc[t, x] * + (b.gas_phase.properties[t, x].temperature - + b.solid_phase.properties[t, x].temperature) * + b.solid_phase.area[t, x]) + + # Dimensionless numbers, mass and heat transfer coefficients + # Particle Reynolds number + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Particle Reynolds number") + def reynolds_number_particle(b, t, x): + return (b.Re_particle[t, x] * + b.gas_phase.properties[t, x].visc_d == + b.velocity_superficial_gas[t, x] * + b.solid_phase.properties[t, x]._params.particle_dia * + b.gas_phase.properties[t, x].dens_mass) + + # Prandtl number + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Prandtl number of gas in bed") + def prandtl_number(b, t, x): + return (b.Pr[t, x] * + b.gas_phase.properties[t, x].therm_cond == + b.solid_phase.properties[t, x].cp_mass * + b.gas_phase.properties[t, x].visc_d) + + # Particle Nusselt number + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Particle Nusselt number") + def nusselt_number_particle(b, t, x): + return (b.Nu_particle[t, x] ** 3 == + (2.0 + 1.1 * (smooth_abs(b.Re_particle[t, x], b.eps) ** + 0.6) ** 3) * + b.Pr[t, x]) + + # Gas-solid heat transfer coefficient + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Gas-solid heat transfer coefficient") + def gas_solid_htc_eqn(b, t, x): + return (1e-3*b.gas_solid_htc[t, x] * + b.solid_phase.properties[t, x]._params.particle_dia == + 1e-3 * b.Nu_particle[t, x] * + b.gas_phase.properties[t, x].therm_cond) + + # Gas phase - gas to solid heat transfer + @self.Constraint(self.flowsheet().config.time, + self.length_domain, + doc="Gas phase - gas to solid heat transfer") + def gas_phase_heat_transfer(b, t, x): + return (b.gas_phase.heat[t, x] * + b.solid_phase.properties[t, x]._params.particle_dia == + -6 * b.gas_solid_htc[t, x] * + (b.gas_phase.properties[t, x].temperature - + b.solid_phase.properties[t, x].temperature) * + b.solid_phase.area[t, x]) + + elif self.config.energy_balance_type == EnergyBalanceType.none: + # If energy balance is none fix gas and solid temperatures to inlet + @self.Constraint( + self.flowsheet().config.time, + self.length_domain, + doc="Isothermal gas phase constraint") + def isothermal_gas_phase(b, t, x): + if x == self.length_domain.first(): + return Constraint.Skip + else: + return ( + b.gas_phase.properties[t, x].temperature == + b.gas_inlet.temperature[t]) + + @self.Constraint( + self.flowsheet().config.time, + self.length_domain, + doc="Isothermal solid phase constraint") + def isothermal_solid_phase(b, t, x): + if x == self.length_domain.last(): + return Constraint.Skip + else: + return ( + b.solid_phase.properties[t, x].temperature == + b.solid_inlet.temperature[t]) + + # ========================================================================= + # Model initialization routine + + def initialize(blk, gas_phase_state_args={}, solid_phase_state_args={}, + outlvl=idaeslog.NOTSET, + solver='ipopt', optarg={'tol': 1e-6}): + """ + Initialisation routine for MB unit (default solver ipopt). + + Keyword Arguments: + state_args : a dict of arguments to be passed to the property + package(s) to provide an initial state for + initialization (see documentation of the specific + property package) (default = {}). + outlvl : sets output level of initialisation routine + optarg : solver options dictionary object (default={'tol': 1e-6}) + solver : str indicating which solver to use during + initialization (default = 'ipopt') + + Returns: + None + """ + + # Set up logger for initialization and solve + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit") + solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit") + + # Set solver options + opt = SolverFactory(solver) + opt.options = optarg + + # --------------------------------------------------------------------- + # local aliases used to shorten object names + gas_phase = blk.config.gas_phase_config + solid_phase = blk.config.solid_phase_config + + # Keep all unit model geometry constraints, derivative_var constraints, + # and property block constraints active. Additionally, in control + # volumes - keep conservation linking constraints and + # holdup calculation (for dynamic flowsheets) constraints active + + geometry_constraints_terms = ["bed_area_eqn", + "solid_phase_area", + "gas_phase_area", + "gas_phase_length", + "solid_phase_length"] + endswith_terms = ("_disc_eq", "linking_constraint", + "linking_constraints", "_holdup_calculation") + startswith_terms = ("properties") + + for c in blk.component_objects(Constraint, descend_into=True): + if not c.parent_block().local_name.startswith(startswith_terms) \ + and not c.local_name.endswith(endswith_terms) \ + and c.local_name not in geometry_constraints_terms: + c.deactivate() + + # --------------------------------------------------------------------- + # Initialize thermophysical property constraints + init_log.info('Initialize Thermophysical Properties') + # Initialize gas_phase block + gas_phase_flags = blk.gas_phase.properties.initialize( + state_args=gas_phase_state_args, + outlvl=outlvl, + optarg=optarg, + solver=solver, + hold_state=True) + + # Initialize solid_phase properties block + solid_phase_flags = blk.solid_phase.properties.initialize( + state_args=solid_phase_state_args, + outlvl=outlvl, + optarg=optarg, + solver=solver, + hold_state=True) + + init_log.info_high("Initialization Step 1 Complete.") + + # --------------------------------------------------------------------- + # Initialize hydrodynamics (velocities) + for t in blk.flowsheet().config.time: + calculate_variable_from_constraint( + blk.velocity_superficial_solid[t], + blk.solid_super_vel[t]) + for x in blk.length_domain: + calculate_variable_from_constraint( + blk.velocity_superficial_gas[t, x], + blk.gas_super_vel[t, x]) + + blk.gas_super_vel.activate() + blk.solid_super_vel.activate() + + init_log.info('Initialize Hydrodynamics') + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + results = opt.solve(blk, tee=slc.tee) + if results.solver.termination_condition \ + == TerminationCondition.optimal: + init_log.info_high( + "Initialization Step 2 {}.".format( + idaeslog.condition(results)) + ) + else: + _log.warning('{} Initialisation Step 2 Failed.' + .format(blk.name)) + + # --------------------------------------------------------------------- + # Initialize mass balance - no reaction and no pressure drop + + # Unfix material balance state variables but keep other states fixed + blk.gas_phase.properties.release_state( + gas_phase_flags) + blk.solid_phase.properties.release_state( + solid_phase_flags) + for t in blk.flowsheet().config.time: + for x in blk.length_domain: + blk.gas_phase.properties[t, x].pressure.fix() + blk.gas_phase.properties[t, x].temperature.fix() + blk.solid_phase.properties[t, x].temperature.fix() + + blk.gas_phase.material_balances.activate() + + if gas_phase.reaction_package is not None: + for t in blk.flowsheet().config.time: + gas_rxn_gen = blk.gas_phase.rate_reaction_generation + for x in blk.length_domain: + for p in gas_phase.property_package.phase_list: + for j in gas_phase.property_package.component_list: + (gas_rxn_gen[t, x, p, j].fix(0.0)) + + blk.solid_phase.material_balances.activate() + + if solid_phase.reaction_package is not None: + for t in blk.flowsheet().config.time: + solid_rxn_gen = blk.solid_phase.rate_reaction_generation + for x in blk.length_domain: + for p in solid_phase.property_package.phase_list: + for j in solid_phase.property_package.component_list: + (solid_rxn_gen[t, x, p, j].fix(0.0)) + + # Gas side heterogeneous rate reaction generation + for x in blk.length_domain: + for p in gas_phase.property_package.phase_list: + for j in gas_phase.property_package.component_list: + (blk.gas_phase.mass_transfer_term[t, x, p, j].fix( + 0.0)) + + init_log.info('Initialize Mass Balances') + init_log.info_high('initialize mass balances - no reactions ' + 'and no pressure drop') + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + results = opt.solve(blk, tee=slc.tee) + if results.solver.termination_condition \ + == TerminationCondition.optimal: + init_log.info_high( + "Initialization Step 3a {}.".format( + idaeslog.condition(results)) + ) + else: + _log.warning('{} Initialisation Step 3a Failed.' + .format(blk.name)) + + # Initialize mass balance - with reaction and no pressure drop + if gas_phase.reaction_package is not None: + # local aliases used to shorten object names + gas_rxn_gen = blk.gas_phase.rate_reaction_generation + gas_phase_stoichiometry_eqn = ( + blk.gas_phase.rate_reaction_stoichiometry_constraint) + + # Initialize reaction property package + blk.gas_phase.reactions.activate() + for t in blk.flowsheet().config.time: + for x in blk.length_domain: + obj = blk.gas_phase.reactions[t, x] + for c in obj.component_objects( + Constraint, descend_into=False): + c.activate() + + blk.gas_phase.reactions.initialize(outlvl=outlvl, + optarg=optarg, + solver=solver) + + for t in blk.flowsheet().config.time: + for x in blk.length_domain: + for r in gas_phase.reaction_package.rate_reaction_idx: + calculate_variable_from_constraint( + blk.gas_phase.rate_reaction_extent[t, x, r], + blk.gas_phase_config_rxn_ext[t, x, r]) + for p in gas_phase.property_package.phase_list: + for j in gas_phase.property_package.component_list: + (gas_rxn_gen[t, x, p, j].unfix()) + calculate_variable_from_constraint( + gas_rxn_gen[t, x, p, j], + gas_phase_stoichiometry_eqn[t, x, p, j]) + + gas_phase_stoichiometry_eqn.activate() + blk.gas_phase_config_rxn_ext.activate() + + if solid_phase.reaction_package is not None: + # local aliases used to shorten object names + solid_rxn_gen = blk.solid_phase.rate_reaction_generation + solid_phase_stoichiometry_eqn = ( + blk.solid_phase.rate_reaction_stoichiometry_constraint) + gas_mass_transfer_term = blk.gas_phase.mass_transfer_term + + # Initialize reaction property package + blk.solid_phase.reactions.activate() + for t in blk.flowsheet().config.time: + for x in blk.length_domain: + obj = blk.solid_phase.reactions[t, x] + for c in obj.component_objects( + Constraint, descend_into=False): + c.activate() + + blk.solid_phase.reactions.initialize(outlvl=outlvl, + optarg=optarg, + solver=solver) + + for t in blk.flowsheet().config.time: + for x in blk.length_domain: + for p in gas_phase.property_package.phase_list: + for j in gas_phase.property_package.component_list: + (gas_mass_transfer_term[t, x, p, j].unfix()) + calculate_variable_from_constraint( + gas_mass_transfer_term[t, x, p, j], + blk.gas_comp_hetero_rxn[t, x, p, j]) + for x in blk.length_domain: + for r in solid_phase.reaction_package.rate_reaction_idx: + calculate_variable_from_constraint( + blk.solid_phase.rate_reaction_extent[t, x, r], + blk.solid_phase_config_rxn_ext[t, x, r]) + for p in solid_phase.property_package.phase_list: + for j in solid_phase.property_package.component_list: + (solid_rxn_gen[t, x, p, j].unfix()) + if (t, x, p, j) in solid_phase_stoichiometry_eqn: + calculate_variable_from_constraint( + solid_rxn_gen[t, x, p, j], + solid_phase_stoichiometry_eqn[t, x, p, j], + ) + + blk.gas_comp_hetero_rxn.activate() + blk.solid_phase.rate_reaction_stoichiometry_constraint.activate() + blk.solid_phase_config_rxn_ext.activate() + + if (gas_phase.reaction_package is not None or + solid_phase.reaction_package is not None): + init_log.info_high('initialize mass balances - with reactions ' + 'and no pressure drop') + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + results = opt.solve(blk, tee=slc.tee) + if results.solver.termination_condition \ + == TerminationCondition.optimal: + init_log.info_high( + "Initialization Step 3b {}.".format( + idaeslog.condition(results)) + ) + else: + _log.warning('{} Initialisation Step 3b Failed.' + .format(blk.name)) + + # Initialize mass balance - with pressure drop + for t in blk.flowsheet().config.time: + for x in blk.length_domain: + # Unfix all pressure variables except at the inlet + if (blk.gas_phase.properties[t, x].config.defined_state + is False): + blk.gas_phase.properties[t, x].pressure.unfix() + + blk.gas_phase.pressure_balance.activate() + + # Set scaling factors for pressure balance equation + blk.gas_phase.scaling_factor_pressure = 1e2 + + if blk.config.has_pressure_change: + blk.gas_phase_config_pressure_drop.activate() + + for t in blk.flowsheet().config.time: + for x in blk.length_domain: + calculate_variable_from_constraint( + blk.gas_phase.deltaP[t, x], + blk.gas_phase_config_pressure_drop[t, x]) + + init_log.info_high('initialize mass balances - with reactions ' + 'and pressure drop') + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + results = opt.solve(blk, tee=slc.tee) + if results.solver.termination_condition \ + == TerminationCondition.optimal: + init_log.info_high( + "Initialization Step 3c {}.".format( + idaeslog.condition(results)) + ) + else: + _log.warning('{} Initialisation Step 3c Failed.' + .format(blk.name)) + # --------------------------------------------------------------------- + # Initialize energy balance + if blk.config.energy_balance_type != EnergyBalanceType.none: + # Initialize dimensionless numbers, + # mass and heat transfer coefficients + for t in blk.flowsheet().config.time: + for x in blk.length_domain: + calculate_variable_from_constraint( + blk.Re_particle[t, x], + blk.reynolds_number_particle[t, x]) + calculate_variable_from_constraint( + blk.Pr[t, x], + blk.prandtl_number[t, x]) + calculate_variable_from_constraint( + blk.Nu_particle[t, x], + blk.nusselt_number_particle[t, x]) + calculate_variable_from_constraint( + blk.gas_solid_htc[t, x], + blk.gas_solid_htc_eqn[t, x]) + calculate_variable_from_constraint( + blk.gas_phase.heat[t, x], + blk.gas_phase_heat_transfer[t, x]) + calculate_variable_from_constraint( + blk.solid_phase.heat[t, x], + blk.solid_phase_heat_transfer[t, x]) + + # Unfix temperatures + for t in blk.flowsheet().config.time: + for x in blk.length_domain: + # Unfix all gas temperature variables except at the inlet + if (blk.gas_phase.properties[t, x].config.defined_state + is False): + blk.gas_phase.properties[t, x].temperature.unfix() + for x in blk.length_domain: + # Unfix all solid temperature variables except at the inlet + if (blk.solid_phase.properties[t, x].config.defined_state + is False): + blk.solid_phase.properties[t, x].temperature.unfix() + + blk.reynolds_number_particle.activate() + blk.prandtl_number.activate() + blk.nusselt_number_particle.activate() + blk.gas_solid_htc_eqn.activate() + + # Activate gas phase energy balance equations + blk.gas_phase_heat_transfer.activate() + blk.gas_phase.enthalpy_balances.activate() + + # Activate solid phase energy balance equations + blk.solid_phase_heat_transfer.activate() + blk.solid_phase.enthalpy_balances.activate() + + init_log.info('Initialize Energy Balances') + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + results = opt.solve(blk, tee=slc.tee) + if results.solver.termination_condition \ + == TerminationCondition.optimal: + init_log.info_high( + "Initialization Step 4 {}.".format( + idaeslog.condition(results)) + ) + else: + _log.warning('{} Initialisation Step 4 Failed.' + .format(blk.name)) + + # Initialize energy balance + if blk.config.energy_balance_type == EnergyBalanceType.none: + for t in blk.flowsheet().config.time: + for x in blk.length_domain: + # Unfix all gas temperature variables except at the inlet + if (blk.gas_phase.properties[t, x].config.defined_state + is False): + blk.gas_phase.properties[t, x].temperature.unfix() + for x in blk.length_domain: + # Unfix all solid temperature variables except at the inlet + if (blk.solid_phase.properties[t, x].config.defined_state + is False): + blk.solid_phase.properties[t, x].temperature.unfix() + + blk.isothermal_gas_phase.activate() + blk.isothermal_solid_phase.activate() + + init_log.info('Initialize Energy Balances') + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + results = opt.solve(blk, tee=slc.tee) + if results.solver.termination_condition \ + == TerminationCondition.optimal: + init_log.info_high( + "Initialization Step 4 {}.".format( + idaeslog.condition(results)) + ) + else: + _log.warning('{} Initialisation Step 4 Failed.' + .format(blk.name)) + + def results_plot(blk): + ''' + Plot method for common moving bed variables + + Variables plotted: + Tg : Temperature in gas phase + Ts : Temperature in solid phase + vg : Superficial gas velocity + P : Pressure in gas phase + Ftotal : Total molar flowrate of gas + Mtotal : Total mass flowrate of solid + Cg : Concentration of gas components in the gas phase + y_frac : Mole fraction of gas components in the gas phase + x_frac : Mass fraction of solid components in the solid phase + ''' + print() + print('================================= Reactor plots ===============' + '==================') + # local aliases used to shorten object names + gas_phase = blk.config.gas_phase_config + solid_phase = blk.config.solid_phase_config + + Tg = [] + Ts = [] + P = [] + Ftotal = [] + Mtotal = [] + vg = [] + + for t in blk.flowsheet().config.time: + for x in blk.gas_phase.length_domain: + vg.append(value(blk.velocity_superficial_gas[t, x])) + + fig_vg = plt.figure(1) + plt.plot(blk.gas_phase.length_domain, vg, label='vg') + plt.legend(loc=0, ncol=2) + plt.grid() + plt.xlabel("Bed height [-]") + plt.ylabel("Superficial gas velocity [m/s]") + fig_vg.savefig('superficial_vel.png') + + # Pressure profile + for t in blk.flowsheet().config.time: + for x in blk.gas_phase.length_domain: + P.append(blk.gas_phase.properties[t, x].pressure.value) + + fig_P = plt.figure(2) + plt.plot(blk.gas_phase.length_domain, P, label='P') + plt.legend(loc=0, ncol=2) + plt.grid() + plt.xlabel("Bed height [-]") + plt.ylabel("Total Pressure [bar]") + fig_P.savefig('Pressure.png') + + # Temperature profile + for t in blk.flowsheet().config.time: + for x in blk.gas_phase.length_domain: + Tg.append(blk.gas_phase.properties[t, x].temperature.value) + for x in blk.solid_phase.length_domain: + Ts.append(blk.solid_phase.properties[t, x].temperature.value) + fig_T = plt.figure(3) + plt.plot(blk.gas_phase.length_domain, Tg, label='Tg') + plt.plot(blk.solid_phase.length_domain, Ts, label='Ts') + plt.legend(loc=0, ncol=2) + plt.grid() + plt.xlabel("Bed height [-]") + plt.ylabel("Temperature [K]") + fig_T.savefig('Temperature.png') + + # Bulk gas phase total molar flow rate + for t in blk.flowsheet().config.time: + for x in blk.gas_phase.length_domain: + Ftotal.append(blk.gas_phase.properties[t, x].flow_mol.value) + fig_Ftotal = plt.figure(4) + plt.plot(blk.gas_phase.length_domain, Ftotal) + plt.grid() + plt.xlabel("Bed height [-]") + plt.ylabel("Total molar gas flow rate [mol/s]") + fig_Ftotal.savefig('Total_gas_flow.png') + + # Bulk solid phase total mass flow rate + for t in blk.flowsheet().config.time: + for x in blk.solid_phase.length_domain: + Mtotal.append(blk.solid_phase.properties[t, x].flow_mass.value) + fig_Mtotal = plt.figure(5) + plt.plot(blk.solid_phase.length_domain, Mtotal) + plt.grid() + plt.xlabel("Bed height [-]") + plt.ylabel("Solid total mass flow rate [kg/s]") + fig_Mtotal.savefig('Total_solid_flow.png') + + # Gas phase mole fractions + for t in blk.flowsheet().config.time: + for j in gas_phase.property_package.component_list: + y_frac = [] + for x in blk.gas_phase.length_domain: + y_frac.append(value( + blk.gas_phase.properties[t, x].mole_frac_comp[j])) + fig_y = plt.figure(6) + plt.plot(blk.gas_phase.length_domain, y_frac, label=j) + plt.legend(loc=0, ncol=len(gas_phase.property_package.component_list)) + plt.grid() + plt.xlabel("Bed height [-]") + plt.ylabel("y_frac [-]") + fig_y.savefig('Gas_mole_fractions.png') + + # Solid phase mass fractions + for t in blk.flowsheet().config.time: + for j in solid_phase.property_package.component_list: + x_frac = [] + for x in blk.solid_phase.length_domain: + x_frac.append(value( + blk.solid_phase.properties[t, x].mass_frac_comp[j])) + fig_x = plt.figure(7) + plt.plot(blk.solid_phase.length_domain, x_frac, label=j) + plt.legend(loc=0, ncol=len( + solid_phase.property_package.component_list)) + plt.grid() + plt.xlabel("Bed height [-]") + plt.ylabel("x_frac [-]") + fig_x.savefig('Solid_mass_fractions.png') + + # Gas phase concentrations + for t in blk.flowsheet().config.time: + for j in gas_phase.property_package.component_list: + Cg = [] + for x in blk.gas_phase.length_domain: + Cg.append( + blk.gas_phase.properties[t, x]. + dens_mol_comp[j].value) + fig_Cg = plt.figure(8) + plt.plot(blk.gas_phase.length_domain, Cg, label=j) + plt.legend(loc=0, ncol=len(gas_phase.property_package.component_list)) + plt.grid() + plt.xlabel("Bed height [-]") + plt.ylabel("Concentration [mol/m3]") + fig_Cg.savefig('Gas_concentration.png') + + def _get_stream_table_contents(self, time_point=0): + return create_stream_table_dataframe( + { + "Gas Inlet": self.gas_inlet, + "Gas Outlet": self.gas_outlet, + "Solid Inlet": self.solid_inlet, + "Solid Outlet": self.solid_outlet, + }, + time_point=time_point, + ) + + +@declare_process_block_class("BiMBR") +class BidirectionalMBData(MBRData): + """ + All this subclass does is override the build method + to let me construct gas and solid control volumes + with different length domains. + """ + + def build(self): + """ + Begin building model (pre-DAE transformation). + + Args: + None + + Returns: + None + """ + # Call UnitModel.build to build default attributes + super(MBRData, self).build() + + # Set flow directions for the control volume blocks + # Gas flows from 0 to 1, solid flows from 1 to 0 + # An if statement is used here despite only one option to allow for + # future extensions to other flow configurations + if self.config.flow_type == "counter_current": + set_direction_gas = FlowDirection.forward + set_direction_solid = FlowDirection.backward + + # Set transformation scheme to be in the "opposite + # direction" as flow. + self.GAS_TRANSFORM_SCHEME = "BACKWARD" + self.SOLID_TRANSFORM_SCHEME = "FORWARD" + else: + raise BurntToast( + "{} encountered unrecognized argument " + "for flow type. Please contact the IDAES" + " developers with this bug.".format(self.name)) + # Set arguments for gas sides if homoogeneous reaction block + if self.config.gas_phase_config.reaction_package is not None: + has_rate_reaction_gas_phase = True + else: + has_rate_reaction_gas_phase = False + + # Set arguments for gas and solid sides if heterogeneous reaction block + if self.config.solid_phase_config.reaction_package is not None: + has_rate_reaction_solid_phase = True + has_mass_transfer_gas_phase = True + else: + has_rate_reaction_solid_phase = False + has_mass_transfer_gas_phase = False + + # Set heat transfer terms + if self.config.energy_balance_type != EnergyBalanceType.none: + has_heat_transfer = True + else: + has_heat_transfer = False + + # Set heat of reaction terms + if (self.config.energy_balance_type != EnergyBalanceType.none + and self.config.gas_phase_config.reaction_package is not None): + has_heat_of_reaction_gas_phase = True + else: + has_heat_of_reaction_gas_phase = False + + if (self.config.energy_balance_type != EnergyBalanceType.none + and self.config.solid_phase_config. + reaction_package is not None): + has_heat_of_reaction_solid_phase = True + else: + has_heat_of_reaction_solid_phase = False + + # Create two different length domains; one for each phase. + # Add them to this block so I can use one of them to index + # all the variables on this block. + # I can then call discretization on this block, which will + # discretize the variables on the control volume blocks. + self.solid_length_domain = ContinuousSet( + bounds=(0.0, 1.0), + initialize=self.config.length_domain_set, + doc="Normalized length domain", + ) + self.gas_length_domain = ContinuousSet( + bounds=(0.0, 1.0), + initialize=self.config.length_domain_set, + doc="Normalized length domain", + ) + self.bed_height = Var(domain=Reals, initialize=1, doc="Bed length [m]") + + # This looks like a hack... Why am I getting around the ProcessBlock + # __setattr__? + super(BlockData, self).__setattr__( + 'length_domain', + self.solid_length_domain, + ) + + # ========================================================================= + """ Build Control volume 1D for gas phase and + populate gas control volume""" + + self.gas_phase = ControlVolume1DBlock(**{ + "transformation_method": self.config.transformation_method, + #"transformation_scheme": self.config.transformation_scheme, + "transformation_scheme": self.GAS_TRANSFORM_SCHEME, + "finite_elements": self.config.finite_elements, + "collocation_points": self.config.collocation_points, + "dynamic": self.config.dynamic, + "has_holdup": self.config.has_holdup, + "area_definition": DistributedVars.variant, + "property_package": self.config.gas_phase_config.property_package, + "property_package_args": + self.config.gas_phase_config.property_package_args, + "reaction_package": self.config.gas_phase_config.reaction_package, + "reaction_package_args": + self.config.gas_phase_config.reaction_package_args}) + + # Pass gas_length_domain to the gas phase control volume + # Note that length_domain_set is redundant as the set is + # already initialized. + self.gas_phase.add_geometry( + length_domain=self.gas_length_domain, + length_domain_set=self.config.length_domain_set, + flow_direction=set_direction_gas, + ) + + self.gas_phase.add_state_blocks( + information_flow=set_direction_gas, + has_phase_equilibrium=False) + + if self.config.gas_phase_config.reaction_package is not None: + self.gas_phase.add_reaction_blocks( + has_equilibrium=self.config.gas_phase_config. + has_equilibrium_reactions) + + self.gas_phase.add_material_balances( + balance_type=self.config.material_balance_type, + has_phase_equilibrium=False, + has_mass_transfer=has_mass_transfer_gas_phase, + has_rate_reactions=has_rate_reaction_gas_phase) + + self.gas_phase.add_energy_balances( + balance_type=self.config.energy_balance_type, + has_heat_transfer=has_heat_transfer, + has_heat_of_reaction=has_heat_of_reaction_gas_phase) + + self.gas_phase.add_momentum_balances( + balance_type=self.config.momentum_balance_type, + has_pressure_change=self.config.has_pressure_change) + + # ========================================================================= + """ Build Control volume 1D for solid phase and + populate solid control volume""" + + # Set argument for heterogeneous reaction block + self.solid_phase = ControlVolume1DBlock(**{ + "transformation_method": self.config.transformation_method, + "transformation_scheme": self.SOLID_TRANSFORM_SCHEME, + "finite_elements": self.config.finite_elements, + "collocation_points": self.config.collocation_points, + # ^ These arguments have no effect as the transformation + # is applied in this class. + "dynamic": self.config.dynamic, + "has_holdup": self.config.has_holdup, + "area_definition": DistributedVars.variant, + "property_package": + self.config.solid_phase_config.property_package, + "property_package_args": + self.config.solid_phase_config.property_package_args, + "reaction_package": + self.config.solid_phase_config.reaction_package, + "reaction_package_args": + self.config.solid_phase_config.reaction_package_args}) + + # Same comment as made for the gas phase. + # Pass in the set we've constructed for this purpose. + self.solid_phase.add_geometry( + length_domain=self.solid_length_domain, + length_domain_set=self.config.length_domain_set, + flow_direction=set_direction_solid) + + # These constraints no longer are created by make_performance, + # So I make them here... + # Length of gas side, and solid side + @self.Constraint(doc="Gas side length") + def gas_phase_length(b): + return (b.gas_phase.length == b.bed_height) + + @self.Constraint(doc="Solid side length") + def solid_phase_length(b): + return (b.solid_phase.length == b.bed_height) + + # Many other methods of the MBR base class actually rely on this + # attribute, so here I slap on a reference. The particular set + # I use shouldn't matter, as long as the derivatives are constructed + # wrt the correct set. + + self.solid_phase.add_state_blocks( + information_flow=set_direction_solid, + has_phase_equilibrium=False) + + if self.config.solid_phase_config.reaction_package is not None: + # TODO - a generalization of the heterogeneous reaction block + # The heterogeneous reaction block does not use the + # add_reaction_blocks in control volumes as control volumes are + # currently setup to handle only homogeneous reaction properties. + # Thus appending the heterogeneous reaction block to the + # solid state block is currently hard coded here. + + tmp_dict = dict(**self.config.solid_phase_config. + reaction_package_args) + tmp_dict["gas_state_block"] = self.gas_phase.properties + tmp_dict["solid_state_block"] = self.solid_phase.properties + tmp_dict["has_equilibrium"] = (self.config.solid_phase_config. + has_equilibrium_reactions) + tmp_dict["parameters"] = (self.config.solid_phase_config. + reaction_package) + self.solid_phase.reactions = ( + self.config.solid_phase_config.reaction_package. + reaction_block_class( + self.flowsheet().config.time, + self.length_domain, + doc="Reaction properties in control volume", + **tmp_dict)) + + self.solid_phase.add_material_balances( + balance_type=self.config.material_balance_type, + has_phase_equilibrium=False, + has_mass_transfer=False, + has_rate_reactions=has_rate_reaction_solid_phase) + + self.solid_phase.add_energy_balances( + balance_type=self.config.energy_balance_type, + has_heat_transfer=has_heat_transfer, + has_heat_of_reaction=has_heat_of_reaction_solid_phase) + + self.solid_phase.add_momentum_balances( + balance_type=MomentumBalanceType.none, + has_pressure_change=False) + + # ========================================================================= + """ Add ports""" + # Add Ports for gas side + self.add_inlet_port(name="gas_inlet", block=self.gas_phase) + self.add_outlet_port(name="gas_outlet", block=self.gas_phase) + + # Add Ports for solid side + self.add_inlet_port(name="solid_inlet", block=self.solid_phase) + self.add_outlet_port(name="solid_outlet", block=self.solid_phase) + + # ========================================================================= + """ Add performance equation method""" + self._apply_transformation() + self._make_performance() + + def _apply_transformation(self): + if self.config.finite_elements is None: + raise ConfigurationError( + 'PROVIDE FINITE_ELEMENTS!!!' + ) + if self.config.transformation_method == 'dae.finite_difference': + self.discretizer = TransformationFactory( + self.config.transformation_method) + # Apply discretization to gas and solid phases separately + self.discretizer.apply_to( + self, + wrt=self.gas_length_domain, + nfe=self.config.finite_elements, + scheme=self.GAS_TRANSFORM_SCHEME, + ) + self.discretizer.apply_to( + self, + wrt=self.solid_length_domain, + nfe=self.config.finite_elements, + scheme=self.SOLID_TRANSFORM_SCHEME, + ) + else: + raise ConfigurationError( + 'THIS SUBCLASS ONLY SUPPORTS FINITE DIFFERENCE!!!' + ) diff --git a/idaes_examples/mod/diagnostics/util.py b/idaes_examples/mod/diagnostics/util.py new file mode 100644 index 00000000..872980c3 --- /dev/null +++ b/idaes_examples/mod/diagnostics/util.py @@ -0,0 +1,44 @@ +import pyomo.environ as pyo +from pyomo.core.expr import EqualityExpression, identify_variables +from pyomo.util.subsystems import create_subsystem_block +from pyomo.common.collections import ComponentMap, ComponentSet +from pyomo.dae.flatten import flatten_dae_components + + +def _get_variables(cons, include_fixed=False): + seen = set() + variables = [] + for con in cons: + for var in identify_variables(con.expr, include_fixed=include_fixed): + if id(var) not in seen: + seen.add(id(var)) + variables.append(var) + return variables + + +def _remove_duplicates(items): + seen = set() + filtered = [] + for item in items: + if id(item) not in seen: + seen.add(id(item)) + filtered.append(item) + return filtered + + +def get_subsystem_at_time(model, time, t): + t0 = time.first() + t1 = time.next(t0) + indices = ComponentMap([(time, t1)]) + scalar_vars, dae_vars = flatten_dae_components(model, time, pyo.Var, indices=indices) + scalar_cons, dae_cons = flatten_dae_components(model, time, pyo.Constraint, indices=indices) + + constraints = _remove_duplicates([ + con[t] for con in dae_cons + if t in con and con[t].active and isinstance(con[t].expr, EqualityExpression) + ]) + var_set = ComponentSet(_get_variables(constraints)) + variables = _remove_duplicates([var[t] for var in dae_vars if var[t] in var_set]) + + subsystem = create_subsystem_block(constraints, variables) + return subsystem, list(subsystem.input_vars.values()) diff --git a/idaes_examples/notebooks/_toc.yml b/idaes_examples/notebooks/_toc.yml index 42a02238..54b632e0 100644 --- a/idaes_examples/notebooks/_toc.yml +++ b/idaes_examples/notebooks/_toc.yml @@ -19,6 +19,7 @@ parts: sections: - file: docs/diagnostics/diagnostics_toolbox_doc - file: docs/diagnostics/degeneracy_hunter_doc + - file: docs/diagnostics/structural_singularity_doc - file: docs/param_est/index sections: - file: docs/param_est/parameter_estimation_nrtl_using_state_block_doc @@ -117,4 +118,4 @@ parts: # Moved this one to 'held' # - file: active/power_gen/ngcc/ngcc_soec_doc - file: active/power_gen/ngcc/ngcc_doc -root: index \ No newline at end of file +root: index diff --git a/idaes_examples/notebooks/docs/diagnostics/index.md b/idaes_examples/notebooks/docs/diagnostics/index.md index b38fb2fb..4defbd31 100644 --- a/idaes_examples/notebooks/docs/diagnostics/index.md +++ b/idaes_examples/notebooks/docs/diagnostics/index.md @@ -3,3 +3,4 @@ IDAES tools for diagnostics: * Introduction to the Diagnostics Toolbox * Solver degeneracy debugging (Degeneracy Hunter) +* Debugging a structural singularity diff --git a/idaes_examples/notebooks/docs/diagnostics/structural_singularity.ipynb b/idaes_examples/notebooks/docs/diagnostics/structural_singularity.ipynb new file mode 100644 index 00000000..565758df --- /dev/null +++ b/idaes_examples/notebooks/docs/diagnostics/structural_singularity.ipynb @@ -0,0 +1,746 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "a261df3e-3957-4a85-ac6d-1bbe9eb053f8", + "metadata": {}, + "outputs": [], + "source": [ + "###############################################################################\n", + "# The Institute for the Design of Advanced Energy Systems Integrated Platform\n", + "# Framework (IDAES IP) was produced under the DOE Institute for the\n", + "# Design of Advanced Energy Systems (IDAES).\n", + "#\n", + "# Copyright (c) 2018-2023 by the software owners: The Regents of the\n", + "# University of California, through Lawrence Berkeley National Laboratory,\n", + "# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon\n", + "# University, West Virginia University Research Corporation, et al.\n", + "# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md\n", + "# for full copyright and license information.\n", + "###############################################################################" + ] + }, + { + "cell_type": "markdown", + "id": "86cdef66", + "metadata": {}, + "source": [ + "Debugging a Structural Singularity\n", + "===========================\n", + "Author: Robert Parker\\\n", + "Maintainer: Robert Parker\\\n", + "Updated: 2024-06-10\n", + "\n", + "In this tutorial, we will use the [IDAES Diagnostics Toolbox](https://idaes-pse.readthedocs.io/en/2.4.0/explanations/model_diagnostics/index.html#diagnostics-toolbox)\n", + "to diagnose and fix a structural singularity that is preventing us from solving an optimization problem." + ] + }, + { + "cell_type": "markdown", + "id": "f456b3f7", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "# Constructing the model\n", + "\n", + "Suppose a collaborator has given us a model to work with. They give us a square model and tell us what the degrees of freedom are. We construct an optimization problem and try to solve it. In this tutorial, we don't want to worry too much about the details that go into constructing the model. This has been provided in the `idaes_examples.mod.diagnostics.gas_solid_contactors.model` module.\n", + "\n", + "## Model details (OKAY TO SKIP)\n", + "\n", + "The model we are trying to optimize is a dynamic model of a moving bed chemical looping combustion reactor. The model has been described by [Okoli et al.][1] and [Parker and Biegler][2]. This is a gas-solid reactor with counter-current flow. The degrees of freedom are gas and solid inlet flow rates, and we are trying to minimize the deviation from a desired operating point via a least-squares objective function.\n", + "\n", + "[1]: https://www.sciencedirect.com/science/article/pii/S0032591019302803\n", + "[2]: https://www.sciencedirect.com/science/article/pii/S2405896322008825\n", + "\n", + "Again, we don't want to worry too much about the model. The `make_model` function will construct the optimization problem that we want to solve, and whenever we do something model-specific, we will explicitly make note of it." + ] + }, + { + "cell_type": "markdown", + "id": "c276aa4a", + "metadata": {}, + "source": [ + "# Trying to solve the original model\n", + "\n", + "With that out of the way, let's construct the model and try to solve it!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "201f63f2", + "metadata": {}, + "outputs": [], + "source": [ + "from idaes_examples.mod.diagnostics.gas_solid_contactors.model import make_model\n", + "import logging\n", + "\n", + "# We'll turn off IDAES logging. This is not recommended in general, but this is an old model\n", + "# (from IDAES 1.7) that has been ported to work with the current version of IDAES. It generates\n", + "# a lot of warnings.\n", + "logging.getLogger(\"idaes\").setLevel(logging.CRITICAL)\n", + "# We'll also turn off Pyomo logging. This will suppress unit inconsistency warnings later,\n", + "# which otherwise flood our console and slow down this notebook. We have unit inconsistencies\n", + "# as, in IDAES 1.7, we didn't rigorously enforce that models use units.\n", + "logging.getLogger(\"pyomo\").setLevel(logging.CRITICAL)\n", + "\n", + "# This constructs a dynamic model with degrees of freedom and an objective function.\n", + "model = make_model()" + ] + }, + { + "cell_type": "markdown", + "id": "2cac2868", + "metadata": {}, + "source": [ + "Before trying to solve the model, let's make sure it conforms to our expectations, i.e. it (a) has degrees of freedom and (b) is well-initialized to a feasible point." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8133d68f", + "metadata": {}, + "outputs": [], + "source": [ + "# Import some useful utilities from the model_statistics module.\n", + "# Degrees of freedom and constraint residuals are always good things to check before\n", + "# trying to solve a simulation or optimization problem.\n", + "from idaes.core.util.model_statistics import degrees_of_freedom, large_residuals_set\n", + "\n", + "dof = degrees_of_freedom(model)\n", + "print(f\"Degrees of freedom: {dof}\")\n", + "has_large_residuals = bool(large_residuals_set(model, tol=1e-5))\n", + "print(f\"Has large residuals: {has_large_residuals}\")" + ] + }, + { + "cell_type": "markdown", + "id": "86e1a118", + "metadata": {}, + "source": [ + "In the above `make_model` function, the model has been \"solved\" to arrive at a feasible point, then degrees of freedom have been unfixed and an objective function has been added to give us an optimization problem. This looks good so far, so let's try to solve the optimization problem." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7b8ddf55", + "metadata": {}, + "outputs": [], + "source": [ + "# Import pyomo.environ for access to solvers\n", + "import pyomo.environ as pyo\n", + "\n", + "solver = pyo.SolverFactory(\"ipopt\")\n", + "solver.options[\"max_iter\"] = 20\n", + "solver.options[\"print_user_options\"] = \"yes\"\n", + "solver.options[\"OF_print_info_string\"] = \"yes\"\n", + "res = solver.solve(model, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "d3513d46", + "metadata": {}, + "source": [ + "IPOPT fails to solve the optimization problem... You can try increasing the iteration limit, but it is very unlikely that this model will ever solve. A telltale sign that something is wrong with our model is the persistence of regularization coefficients, that is, numbers in the `lg(rg)` column of the IPOPT log. These coefficients can have multiple causes. One is that the constraint Jacobian (partial derivative matrix) is singular, which indicates a problem with our model. We have set the `print_info_string` option in IPOPT to display \"diagnostic tags\" to help interpret these regularization coefficients. The \"L\" and \"l\" diagnostic tags, which appear repeatedly, indicate that the Jacobian is singular. For more information on IPOPT diagnostic tags, see the IPOPT [documentation](https://coin-or.github.io/Ipopt/OUTPUT.html)." + ] + }, + { + "cell_type": "markdown", + "id": "b01c2de5", + "metadata": {}, + "source": [ + "# Debugging the original model\n", + "\n", + "Let's run the diagnostics toolbox on the model and see what it has to say.\n", + "\n", + "For good practice, we'll first make sure the model we're debugging is square. Remember that we're assuming we already know how to toggle degrees of freedom in our model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "107471c0", + "metadata": {}, + "outputs": [], + "source": [ + "# Fix gas and solid flow rates at their respective inlets\n", + "model.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "# Part of our optimization problem was a set of constraints to enforce piecewise\n", + "# constant control inputs. We need to deactivate these as well.\n", + "model.piecewise_constant_constraints.deactivate()" + ] + }, + { + "cell_type": "markdown", + "id": "36f82a0a", + "metadata": {}, + "source": [ + "Now we can run the diagnostics toolbox." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7667a9f", + "metadata": {}, + "outputs": [], + "source": [ + "from idaes.core.util.model_diagnostics import DiagnosticsToolbox\n", + "\n", + "dt = DiagnosticsToolbox(model)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "id": "3b4530cf", + "metadata": {}, + "source": [ + "Let's look at the warnings we got:\n", + "- Inconsistent units\n", + "- Structural singularity\n", + "- Potential evaluation errors\n", + "\n", + "We'll ignore the inconsistent units. The property package and unit model here were extracted from IDAES 1.7, before we rigorously enforced that all models use units. The potential evaluation errors we see here may be worth looking into, but looking at the failing IPOPT log above, we don't notice any evaluation errors. (If evaluation errors occurred in IPOPT, we would see a message like \"Error in AMPL evaluation\" in the IPOPT iteration log, which we don't see here.) The structural singularity looks like the most promising avenue to debug, especially as the IPOPT log displays persistent regularization coefficients that appear to be caused by a singular Jacobian.\n", + "\n", + "Let's follow the toolbox's advice and display the under and over-constrained sets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9976bc61", + "metadata": {}, + "outputs": [], + "source": [ + "dt.display_underconstrained_set()\n", + "dt.display_overconstrained_set()" + ] + }, + { + "cell_type": "markdown", + "id": "eb67e99c", + "metadata": {}, + "source": [ + "## Over and under-constrained subsystems\n", + "\n", + "Structural singularities are characterized by the [Dulmage-Mendelson decomposition][3], which partitions a system into minimal over and under-constrained subsystems. These subsystems contain the potentially unmatched constraints and variables, respectively. Here, \"unmatched\" effectively means \"causing a singularity\". [Pothen and Fan][4] give a good overview of the Dulmage-Mendelsohn decomposition and [Parker et al.][5] give several examples.\n", + "\n", + "[3]: https://www.cambridge.org/core/journals/canadian-journal-of-mathematics/article/coverings-of-bipartite-graphs/413735C5888AB542B92D0C4F402800B1\n", + "[4]: https://dl.acm.org/doi/10.1145/98267.98287\n", + "[5]: https://www.sciencedirect.com/science/article/pii/S0098135423002533\n", + "\n", + "The most straightforward way to fix a structural singularity is to fix variables that are in the under-constrained system and deactivate constraints in the over-constrained subsystem. However, this may not be applicable for every model. For example, we may need to add variables and constraints instead. What over and under-constrained subsystems are telling us is that something is wrong with our modeling assumptions. The particular fix that is appropriate will depend heavily on the model.\n", + "\n", + "If the above output gives us any clues, we can go ahead and start trying to fix things. However, suppose it doesn't. A good strategy is to try to break down the model into smaller, square subsystems that we think should be nonsingular. For a dynamic model like this one, a good candidate is the subsystem of variables and equations at each point in time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea05b00e", + "metadata": {}, + "outputs": [], + "source": [ + "# We've included a utility function to extract the subsystem of variables and equations\n", + "# at a specified point in time. If you are dealing with a process flowsheet, here you\n", + "# may want to extract each unit model individually.\n", + "from idaes_examples.mod.diagnostics.util import get_subsystem_at_time\n", + "# TemporarySubsystemManager is used to temporarily fix some variables to make sure\n", + "# we're debugging a square subsystem.\n", + "from pyomo.util.subsystems import TemporarySubsystemManager\n", + "\n", + "# Let's start with t=0. Really, we'd probably want to do this in a loop and try all time points.\n", + "t0 = model.fs.time.first()\n", + "t_block, inputs = get_subsystem_at_time(model, model.fs.time, t0)\n", + "# We'll temporarily fix the \"inputs\" to make sure we have a square system while debugging\n", + "with TemporarySubsystemManager(to_fix=inputs):\n", + " dt = DiagnosticsToolbox(t_block)\n", + " dt.report_structural_issues()\n", + " dt.display_underconstrained_set()\n", + " dt.display_overconstrained_set()" + ] + }, + { + "cell_type": "markdown", + "id": "986b5113", + "metadata": {}, + "source": [ + "These over and under-constrained subsystems aren't much smaller, but now the over-constrained system decomposes into 10 small, independent blocks. These should be easier to debug.\n", + "\n", + "## Debugging the over-constrained subsystem\n", + "\n", + "To debug the over-constrained subsystem, we look for a constraint that is not calculating any of the variables in the subsystem. The \"odd constraint out\" here seems to be the mass fraction sum, `sum_component_eqn`. This must \"solve for\" one of the mass fractions, which means one of the `material_holdup_calculation` equations must \"solve for\" particle density rather than mass fraction. If we want to see what variables are contained in one of these constraints, we can always `pprint` it:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4029972c", + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].sum_component_eqn.pprint()\n", + "model.fs.MB.solid_phase.material_holdup_calculation[0, 0.9, \"Sol\", \"Fe3O4\"].pprint()" + ] + }, + { + "cell_type": "markdown", + "id": "a11d6609", + "metadata": {}, + "source": [ + "If one of these `material_holdup_calculation` equations is solving for particle density, then that means that `density_particle_constraint` is not actually solving for density. Maybe `density_particle_constraint` is over-determining our system?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d1836161", + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].density_particle_constraint.pprint()" + ] + }, + { + "cell_type": "markdown", + "id": "c8e9abf7", + "metadata": {}, + "source": [ + "But this looks like a very reasonable constraint. After some thought, which admittedly requires some knowledge of the process we are modeling, we decide that the right approach is to make particle porosity a variable. We have assumed that porosity is constant, but this overconstrained subsystem is telling us that this assumption is not valid.\n", + "\n", + "### How did we figure this out? (OKAY TO SKIP)\n", + "Adding a variable (including by unfixing a parameter) to an over-constraining constraint will often remove that constraint from the over-constrained subsystem. But how did we know that this was the right thing to do? If you just care about using the diagnostics toolbox to extract as much information about a singularity as possible, you can skip this section. But if you are curious how we determined that particle porosity should not be constant, read on.\n", + "\n", + "`dens_mass_skeletal` is determined purely by the composition of solid, which is made up of Fe2O3, Fe3O4, and inert Ti2O3. We can view the `density_skeletal_constraint` as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c13a3a8-8e77-498d-b6fe-ad1f88413cbf", + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].density_skeletal_constraint.pprint()" + ] + }, + { + "cell_type": "markdown", + "id": "c42000a6-e0f1-4118-b67f-d878bcf777a2", + "metadata": {}, + "source": [ + "If we assume a constant particle porosity, this gives us a particle porosity that is also uniquely determined by the solid composition by the above `density_particle_constraint`:\n", + "```\n", + "dens_mass_particle = (1 - porosity) * dens_mass_skeletal\n", + "```\n", + "But the composition of the solid is determined by the (somewhat misnamed) `material_holdup_calculation` constraints. While the name of these constraints implies they \"calculate holdups,\" material holdups at $t=0$ are fixed as initial conditions (because holdups are the differential variables with respect to time in this model). At other time points, we assume that holdups are specified by differential and discretization equations of the model. This means that the `material_holdup_calculation` constraints actually calculate the solid phase mass fractions *from* the holdups. But as we hinted at above, the 4-by-4 system of holdup calculation constraints, `sum_component_eqn` (which simply constrains the sum of mass fractions to be one), mass fractions, and `dens_mass_particle`, uniquely solve for `dens_mass_particle` *as well as* the mass fractions. But if the holdup variables can be used to solve for the mass fractions, they *also* solve for `dens_mass_skeletal`. So both sides of `density_particle_constraint` are already uniquely determined! This implies that we don't need this constraint at all, but we also know that this constraint has to hold. Something has to give. With this in mind, we actually have several options for how to resolve this overspecification:\n", + "1. Remove `density_particle_constraint`. Then we would have `dens_mass_particle` and `dens_mass_skeletal`, with no relationship between them. This would leave us with a mathematically sound model, but with densities that contradict constant particle porosity that we have assumed (which is used elsewhere in the reaction rate calculation equations).\n", + "2. Remove the constraints that calculate skeletal density from composition.\n", + "3. Relax particle porosity from a parameter to a variable.\n", + "\n", + "Options 2 and 3 are equally valid. We've chosen option 3, meaning we assume that the particle \"evolves\" with a density that is well determined from its constituent species, rather than changing density to accommodate whatever mass it accumulates via reaction without altering its volume. This exercise should remind us that all mathematical modeling is somewhat of an art. In the process of choosing the \"least bad\" model, it is fairly easy to over or under-specify something by making the wrong combination of assumptions, and the Dulmage-Mendelsohn decomposition is a great tool for detecting when this has happened.\n", + "\n", + "## Debugging the under-constrained subsystem\n", + "\n", + "The under-constrained system does not decompose into independent subsystems, making it more difficult to debug. However, by inspection, we notice that the same constraints and variables seem to be repeated at each point in the length domain. For each point in space, the \"odd variable out\" seems to be the total flow rate `flow_mass`. Using some intuition about this particular process model, we may conclude that this variable should be calculated from the solid phase velocity, which is constant. We expect an equation that looks like\n", + "```\n", + "flow_mass == velocity * area * density\n", + "```\n", + "\n", + "But this equation isn't here... so we need to add it.\n", + "\n", + "# Fixing the model\n", + "\n", + "We'll start by creating a fresh copy of the model, so we don't accidentally rely on IPOPT's point of termination." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01cd1929", + "metadata": {}, + "outputs": [], + "source": [ + "model2 = make_model()\n", + "# Make the model square while we try to fix the structural singularity\n", + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.deactivate()" + ] + }, + { + "cell_type": "markdown", + "id": "565599a6", + "metadata": {}, + "source": [ + "## Adding a new particle porosity variable" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "479f3067", + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.particle_porosity = pyo.Var(\n", + " model2.fs.time, model2.fs.MB.length_domain, initialize=model2.fs.solid_properties.particle_porosity.value\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "b614d7a0", + "metadata": {}, + "source": [ + "Now we need to replace the old particle porosity parameter with this new variable. Luckily, the old parameter is actually implemented as a fixed variable, so we can easily identify all the constraints it participates in with `IncidenceGraphInterface`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d45b4a9", + "metadata": {}, + "outputs": [], + "source": [ + "from pyomo.contrib.incidence_analysis import IncidenceGraphInterface\n", + "\n", + "igraph = IncidenceGraphInterface(model2, include_fixed=True)\n", + "porosity_param = model2.fs.solid_properties.particle_porosity\n", + "print(f\"Constraints containing {porosity_param.name}:\")\n", + "for con in igraph.get_adjacent_to(porosity_param):\n", + " print(f\" {con.name}\")" + ] + }, + { + "cell_type": "markdown", + "id": "0e6fee8f", + "metadata": {}, + "source": [ + "Particle porosity only appears in two constraints: the density constraint we saw above, and the reaction rate equation. We can replace particle porosity in these constraints using Pyomo's `replace_expressions` function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d62bd7a", + "metadata": {}, + "outputs": [], + "source": [ + "from pyomo.core.expr import replace_expressions\n", + "\n", + "for t, x in model2.fs.time * model2.fs.MB.length_domain:\n", + " substitution_map = {id(porosity_param): model2.fs.MB.particle_porosity[t, x]}\n", + " sp = model2.fs.MB.solid_phase\n", + " cons = [sp.properties[t, x].density_particle_constraint, sp.reactions[t, x].gen_rate_expression[\"R1\"]]\n", + " for con in cons:\n", + " con.set_value(replace_expressions(con.expr, substitution_map, descend_into_named_expressions=True))" + ] + }, + { + "cell_type": "markdown", + "id": "cff3dc24", + "metadata": {}, + "source": [ + "We have added a new `particle_porosity` variable, and are using it in the relevant locations. Now we can move on to adding the missing constraint.\n", + "\n", + "## Adding a new density-flow rate constraint" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2cd11e5d", + "metadata": {}, + "outputs": [], + "source": [ + "@model2.fs.MB.Constraint(model2.fs.time, model2.fs.MB.length_domain)\n", + "def density_flowrate_constraint(mb, t, x):\n", + " return (\n", + " mb.velocity_superficial_solid[t] * mb.bed_area\n", + " * mb.solid_phase.properties[t, x].dens_mass_particle\n", + " == mb.solid_phase.properties[t, x].flow_mass\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "eca5108e", + "metadata": {}, + "source": [ + "## Testing the new model\n", + "\n", + "Let's see if these changes have fixed our model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc391ba8", + "metadata": {}, + "outputs": [], + "source": [ + "# Construct a new diagnostics toolbox\n", + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "id": "de070064", + "metadata": {}, + "source": [ + "The structural singularity seems to be gone! Let's unfix our degrees of freedom and see if we can solve." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fc14d9e2", + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix()\n", + "model2.fs.MB.gas_phase.properties[0, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix()\n", + "model2.fs.MB.solid_phase.properties[0, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.activate()\n", + "\n", + "res = solver.solve(model2, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "fd8748b9", + "metadata": {}, + "source": [ + "This doesn't look much better. What's going on? I thought we just fixed the issue?\n", + "\n", + "# Debugging the model, take two\n", + "\n", + "Let's check the diagnostics toolbox for numerical issues." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "199e993e", + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.deactivate()\n", + "dt.report_numerical_issues()" + ] + }, + { + "cell_type": "markdown", + "id": "a2040b5b", + "metadata": {}, + "source": [ + "Looks like we have \"parallel constraints\", which are another form of singularity. Let's follow the toolbox's advice to see what they are." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b8c21fd9-510c-455d-9360-9005710986b6", + "metadata": {}, + "outputs": [], + "source": [ + "dt.display_near_parallel_constraints()" + ] + }, + { + "cell_type": "markdown", + "id": "e2b8d8d3-344f-496e-9128-0aac83a445d1", + "metadata": {}, + "source": [ + "`density_flowrate_constraint` is the constraint that we added. What is `solid_super_vel`?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f845733e-a315-4833-b0a9-012cc75a3e48", + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.solid_super_vel[0].pprint()" + ] + }, + { + "cell_type": "markdown", + "id": "0559c8c4", + "metadata": {}, + "source": [ + "This is the same as the constraint we just added! Looks like that constraint already existed at the solid inlet. We can easily deactivate the new constraints at this point in the length domain:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8f086a50", + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.density_flowrate_constraint[:, 1.0].deactivate();" + ] + }, + { + "cell_type": "markdown", + "id": "9d257074", + "metadata": {}, + "source": [ + "But now we have removed constraints from a square model, and expect to have degrees of freedom. Let's see what the diagnostics toolbox has to say." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c1af8b6c", + "metadata": {}, + "outputs": [], + "source": [ + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "id": "0b27aee5", + "metadata": {}, + "source": [ + "But this doesn't help us very much. We have some extraneous degrees of freedom, but with 8881 variables in the under-constrained subsystem, it will be difficult to tell what they are. After some thought (and model-specific intuition), we land on the conclusion that maybe we need to fix particle porosity at the solid inlet. Here, total flow rate is specified, and the `solid_super_vel` equation is using it to compute velocity. So we need `dens_mass_particle` to be known, which means we need `particle_porosity` to be fixed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "656ae921", + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.particle_porosity[:, 1.0].fix();" + ] + }, + { + "cell_type": "markdown", + "id": "cbb89807", + "metadata": {}, + "source": [ + "Let's run the diagnostics toolbox as a sanity check." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6cf86e28", + "metadata": {}, + "outputs": [], + "source": [ + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()\n", + "dt.report_numerical_issues()" + ] + }, + { + "cell_type": "markdown", + "id": "016f8799", + "metadata": {}, + "source": [ + "Looks good! Now we can release our degrees of freedom and try to solve again." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2f83e2e", + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix()\n", + "model2.fs.MB.gas_phase.properties[0, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix()\n", + "model2.fs.MB.solid_phase.properties[0, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.activate()\n", + "\n", + "res = solver.solve(model2, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "b8587acd", + "metadata": {}, + "source": [ + "It worked! For the simple optimization problem we have set up, this solve looks a lot more like what we expect.\n", + "\n", + "# Takeaways from this tutorial\n", + "What have we learned?\n", + "1. IPOPT using non-zero regularization coefficients hints at a singular Jacobian (especially when \"L\"/\"l\" diagnostic tags are present).\n", + "2. When this happens, start by calling `report_structural_issues` to check for a structural singularity. If this looks good, call `report_numerical_issues` to check for a numerical singularity.\n", + "3. When debugging a structural singularity, decomposing a problem into subsystems that each should be nonsingular (e.g. unit models or points in time) is very useful.\n", + "4. The solution to a structural singularity is often to relax a fixed parameter, add a constraint that was forgotten, remove a constraint that was redundant, or fix an extraneous degree of freedom.\n", + "5. Model-specific intuition is usually necessary to diagnose and fix modeling issues. (If you're an algorithm developer, learn about the models you're using! If you don't understand your models, you don't understand your algorithms!)\n", + "6. A modeling issue doesn't necessarily have a unique solution. This is especially true when the issue involves invalid assumptions.\n", + "7. Debugging is an iterative process — fixing one issue can introduce another." + ] + }, + { + "cell_type": "markdown", + "id": "33983d41", + "metadata": {}, + "source": [ + "# References\n", + "\n", + "[[1]] Okoli et al., \"A framework for the optimization of chemical looping combustion processes\". *Powder Tech*, 2020.\n", + "\n", + "[[2]] Parker and Biegler, \"Dynamic modeling and nonlinear model predictive control of a moving bed chemical looping combustion reactor\". *IFAC PapersOnline*, 2022.\n", + "\n", + "[[3]] Dulmage and Mendelsohn, \"Coverings of bipartite graphs\". *Can J. Math.*, 1958.\n", + "\n", + "[[4]] Pothen and Fan, \"Computing the block triangular form of a sparse matrix\". *ACM Trans. Math. Softw.*, 1990.\n", + "\n", + "[[5]] Parker et al., \"Applications of the Dulmage-Mendelsohn decomposition for debugging nonlinear optimization problems\". *Comp. Chem. Eng.*, 2023.\n", + "\n", + "[1]: https://www.sciencedirect.com/science/article/pii/S0032591019302803\n", + "[2]: https://www.sciencedirect.com/science/article/pii/S2405896322008825\n", + "[3]: https://www.cambridge.org/core/journals/canadian-journal-of-mathematics/article/coverings-of-bipartite-graphs/413735C5888AB542B92D0C4F402800B1\n", + "[4]: https://dl.acm.org/doi/10.1145/98267.98287\n", + "[5]: https://www.sciencedirect.com/science/article/pii/S0098135423002533\n" + ] + } + ], + "metadata": { + "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.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/idaes_examples/notebooks/docs/diagnostics/structural_singularity_doc.ipynb b/idaes_examples/notebooks/docs/diagnostics/structural_singularity_doc.ipynb new file mode 100644 index 00000000..46b7fb0b --- /dev/null +++ b/idaes_examples/notebooks/docs/diagnostics/structural_singularity_doc.ipynb @@ -0,0 +1,690 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "###############################################################################\n", + "# The Institute for the Design of Advanced Energy Systems Integrated Platform\n", + "# Framework (IDAES IP) was produced under the DOE Institute for the\n", + "# Design of Advanced Energy Systems (IDAES).\n", + "#\n", + "# Copyright (c) 2018-2023 by the software owners: The Regents of the\n", + "# University of California, through Lawrence Berkeley National Laboratory,\n", + "# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon\n", + "# University, West Virginia University Research Corporation, et al.\n", + "# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md\n", + "# for full copyright and license information.\n", + "###############################################################################" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Debugging a Structural Singularity\n", + "===========================\n", + "Author: Robert Parker\\\n", + "Maintainer: Robert Parker\\\n", + "Updated: 2024-06-10\n", + "\n", + "In this tutorial, we will use the [IDAES Diagnostics Toolbox](https://idaes-pse.readthedocs.io/en/2.4.0/explanations/model_diagnostics/index.html#diagnostics-toolbox)\n", + "to diagnose and fix a structural singularity that is preventing us from solving an optimization problem." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "# Constructing the model\n", + "\n", + "Suppose a collaborator has given us a model to work with. They give us a square model and tell us what the degrees of freedom are. We construct an optimization problem and try to solve it. In this tutorial, we don't want to worry too much about the details that go into constructing the model. This has been provided in the `idaes_examples.mod.diagnostics.gas_solid_contactors.model` module.\n", + "\n", + "## Model details (OKAY TO SKIP)\n", + "\n", + "The model we are trying to optimize is a dynamic model of a moving bed chemical looping combustion reactor. The model has been described by [Okoli et al.][1] and [Parker and Biegler][2]. This is a gas-solid reactor with counter-current flow. The degrees of freedom are gas and solid inlet flow rates, and we are trying to minimize the deviation from a desired operating point via a least-squares objective function.\n", + "\n", + "[1]: https://www.sciencedirect.com/science/article/pii/S0032591019302803\n", + "[2]: https://www.sciencedirect.com/science/article/pii/S2405896322008825\n", + "\n", + "Again, we don't want to worry too much about the model. The `make_model` function will construct the optimization problem that we want to solve, and whenever we do something model-specific, we will explicitly make note of it." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Trying to solve the original model\n", + "\n", + "With that out of the way, let's construct the model and try to solve it!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from idaes_examples.mod.diagnostics.gas_solid_contactors.model import make_model\n", + "import logging\n", + "\n", + "# We'll turn off IDAES logging. This is not recommended in general, but this is an old model\n", + "# (from IDAES 1.7) that has been ported to work with the current version of IDAES. It generates\n", + "# a lot of warnings.\n", + "logging.getLogger(\"idaes\").setLevel(logging.CRITICAL)\n", + "# We'll also turn off Pyomo logging. This will suppress unit inconsistency warnings later,\n", + "# which otherwise flood our console and slow down this notebook. We have unit inconsistencies\n", + "# as, in IDAES 1.7, we didn't rigorously enforce that models use units.\n", + "logging.getLogger(\"pyomo\").setLevel(logging.CRITICAL)\n", + "\n", + "# This constructs a dynamic model with degrees of freedom and an objective function.\n", + "model = make_model()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Before trying to solve the model, let's make sure it conforms to our expectations, i.e. it (a) has degrees of freedom and (b) is well-initialized to a feasible point." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import some useful utilities from the model_statistics module.\n", + "# Degrees of freedom and constraint residuals are always good things to check before\n", + "# trying to solve a simulation or optimization problem.\n", + "from idaes.core.util.model_statistics import degrees_of_freedom, large_residuals_set\n", + "\n", + "dof = degrees_of_freedom(model)\n", + "print(f\"Degrees of freedom: {dof}\")\n", + "has_large_residuals = bool(large_residuals_set(model, tol=1e-5))\n", + "print(f\"Has large residuals: {has_large_residuals}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the above `make_model` function, the model has been \"solved\" to arrive at a feasible point, then degrees of freedom have been unfixed and an objective function has been added to give us an optimization problem. This looks good so far, so let's try to solve the optimization problem." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import pyomo.environ for access to solvers\n", + "import pyomo.environ as pyo\n", + "\n", + "solver = pyo.SolverFactory(\"ipopt\")\n", + "solver.options[\"max_iter\"] = 20\n", + "solver.options[\"print_user_options\"] = \"yes\"\n", + "solver.options[\"OF_print_info_string\"] = \"yes\"\n", + "res = solver.solve(model, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "IPOPT fails to solve the optimization problem... You can try increasing the iteration limit, but it is very unlikely that this model will ever solve. A telltale sign that something is wrong with our model is the persistence of regularization coefficients, that is, numbers in the `lg(rg)` column of the IPOPT log. These coefficients can have multiple causes. One is that the constraint Jacobian (partial derivative matrix) is singular, which indicates a problem with our model. We have set the `print_info_string` option in IPOPT to display \"diagnostic tags\" to help interpret these regularization coefficients. The \"L\" and \"l\" diagnostic tags, which appear repeatedly, indicate that the Jacobian is singular. For more information on IPOPT diagnostic tags, see the IPOPT [documentation](https://coin-or.github.io/Ipopt/OUTPUT.html)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Debugging the original model\n", + "\n", + "Let's run the diagnostics toolbox on the model and see what it has to say.\n", + "\n", + "For good practice, we'll first make sure the model we're debugging is square. Remember that we're assuming we already know how to toggle degrees of freedom in our model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Fix gas and solid flow rates at their respective inlets\n", + "model.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "# Part of our optimization problem was a set of constraints to enforce piecewise\n", + "# constant control inputs. We need to deactivate these as well.\n", + "model.piecewise_constant_constraints.deactivate()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can run the diagnostics toolbox." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from idaes.core.util.model_diagnostics import DiagnosticsToolbox\n", + "\n", + "dt = DiagnosticsToolbox(model)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's look at the warnings we got:\n", + "- Inconsistent units\n", + "- Structural singularity\n", + "- Potential evaluation errors\n", + "\n", + "We'll ignore the inconsistent units. The property package and unit model here were extracted from IDAES 1.7, before we rigorously enforced that all models use units. The potential evaluation errors we see here may be worth looking into, but looking at the failing IPOPT log above, we don't notice any evaluation errors. (If evaluation errors occurred in IPOPT, we would see a message like \"Error in AMPL evaluation\" in the IPOPT iteration log, which we don't see here.) The structural singularity looks like the most promising avenue to debug, especially as the IPOPT log displays persistent regularization coefficients that appear to be caused by a singular Jacobian.\n", + "\n", + "Let's follow the toolbox's advice and display the under and over-constrained sets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt.display_underconstrained_set()\n", + "dt.display_overconstrained_set()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Over and under-constrained subsystems\n", + "\n", + "Structural singularities are characterized by the [Dulmage-Mendelson decomposition][3], which partitions a system into minimal over and under-constrained subsystems. These subsystems contain the potentially unmatched constraints and variables, respectively. Here, \"unmatched\" effectively means \"causing a singularity\". [Pothen and Fan][4] give a good overview of the Dulmage-Mendelsohn decomposition and [Parker et al.][5] give several examples.\n", + "\n", + "[3]: https://www.cambridge.org/core/journals/canadian-journal-of-mathematics/article/coverings-of-bipartite-graphs/413735C5888AB542B92D0C4F402800B1\n", + "[4]: https://dl.acm.org/doi/10.1145/98267.98287\n", + "[5]: https://www.sciencedirect.com/science/article/pii/S0098135423002533\n", + "\n", + "The most straightforward way to fix a structural singularity is to fix variables that are in the under-constrained system and deactivate constraints in the over-constrained subsystem. However, this may not be applicable for every model. For example, we may need to add variables and constraints instead. What over and under-constrained subsystems are telling us is that something is wrong with our modeling assumptions. The particular fix that is appropriate will depend heavily on the model.\n", + "\n", + "If the above output gives us any clues, we can go ahead and start trying to fix things. However, suppose it doesn't. A good strategy is to try to break down the model into smaller, square subsystems that we think should be nonsingular. For a dynamic model like this one, a good candidate is the subsystem of variables and equations at each point in time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# We've included a utility function to extract the subsystem of variables and equations\n", + "# at a specified point in time. If you are dealing with a process flowsheet, here you\n", + "# may want to extract each unit model individually.\n", + "from idaes_examples.mod.diagnostics.util import get_subsystem_at_time\n", + "# TemporarySubsystemManager is used to temporarily fix some variables to make sure\n", + "# we're debugging a square subsystem.\n", + "from pyomo.util.subsystems import TemporarySubsystemManager\n", + "\n", + "# Let's start with t=0. Really, we'd probably want to do this in a loop and try all time points.\n", + "t0 = model.fs.time.first()\n", + "t_block, inputs = get_subsystem_at_time(model, model.fs.time, t0)\n", + "# We'll temporarily fix the \"inputs\" to make sure we have a square system while debugging\n", + "with TemporarySubsystemManager(to_fix=inputs):\n", + " dt = DiagnosticsToolbox(t_block)\n", + " dt.report_structural_issues()\n", + " dt.display_underconstrained_set()\n", + " dt.display_overconstrained_set()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These over and under-constrained subsystems aren't much smaller, but now the over-constrained system decomposes into 10 small, independent blocks. These should be easier to debug.\n", + "\n", + "## Debugging the over-constrained subsystem\n", + "\n", + "To debug the over-constrained subsystem, we look for a constraint that is not calculating any of the variables in the subsystem. The \"odd constraint out\" here seems to be the mass fraction sum, `sum_component_eqn`. This must \"solve for\" one of the mass fractions, which means one of the `material_holdup_calculation` equations must \"solve for\" particle density rather than mass fraction. If we want to see what variables are contained in one of these constraints, we can always `pprint` it:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].sum_component_eqn.pprint()\n", + "model.fs.MB.solid_phase.material_holdup_calculation[0, 0.9, \"Sol\", \"Fe3O4\"].pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If one of these `material_holdup_calculation` equations is solving for particle density, then that means that `density_particle_constraint` is not actually solving for density. Maybe `density_particle_constraint` is over-determining our system?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].density_particle_constraint.pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "But this looks like a very reasonable constraint. After some thought, which admittedly requires some knowledge of the process we are modeling, we decide that the right approach is to make particle porosity a variable. We have assumed that porosity is constant, but this overconstrained subsystem is telling us that this assumption is not valid.\n", + "\n", + "### How did we figure this out? (OKAY TO SKIP)\n", + "Adding a variable (including by unfixing a parameter) to an over-constraining constraint will often remove that constraint from the over-constrained subsystem. But how did we know that this was the right thing to do? If you just care about using the diagnostics toolbox to extract as much information about a singularity as possible, you can skip this section. But if you are curious how we determined that particle porosity should not be constant, read on.\n", + "\n", + "`dens_mass_skeletal` is determined purely by the composition of solid, which is made up of Fe2O3, Fe3O4, and inert Ti2O3. We can view the `density_skeletal_constraint` as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].density_skeletal_constraint.pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we assume a constant particle porosity, this gives us a particle porosity that is also uniquely determined by the solid composition by the above `density_particle_constraint`:\n", + "```\n", + "dens_mass_particle = (1 - porosity) * dens_mass_skeletal\n", + "```\n", + "But the composition of the solid is determined by the (somewhat misnamed) `material_holdup_calculation` constraints. While the name of these constraints implies they \"calculate holdups,\" material holdups at $t=0$ are fixed as initial conditions (because holdups are the differential variables with respect to time in this model). At other time points, we assume that holdups are specified by differential and discretization equations of the model. This means that the `material_holdup_calculation` constraints actually calculate the solid phase mass fractions *from* the holdups. But as we hinted at above, the 4-by-4 system of holdup calculation constraints, `sum_component_eqn` (which simply constrains the sum of mass fractions to be one), mass fractions, and `dens_mass_particle`, uniquely solve for `dens_mass_particle` *as well as* the mass fractions. But if the holdup variables can be used to solve for the mass fractions, they *also* solve for `dens_mass_skeletal`. So both sides of `density_particle_constraint` are already uniquely determined! This implies that we don't need this constraint at all, but we also know that this constraint has to hold. Something has to give. With this in mind, we actually have several options for how to resolve this overspecification:\n", + "1. Remove `density_particle_constraint`. Then we would have `dens_mass_particle` and `dens_mass_skeletal`, with no relationship between them. This would leave us with a mathematically sound model, but with densities that contradict constant particle porosity that we have assumed (which is used elsewhere in the reaction rate calculation equations).\n", + "2. Remove the constraints that calculate skeletal density from composition.\n", + "3. Relax particle porosity from a parameter to a variable.\n", + "\n", + "Options 2 and 3 are equally valid. We've chosen option 3, meaning we assume that the particle \"evolves\" with a density that is well determined from its constituent species, rather than changing density to accommodate whatever mass it accumulates via reaction without altering its volume. This exercise should remind us that all mathematical modeling is somewhat of an art. In the process of choosing the \"least bad\" model, it is fairly easy to over or under-specify something by making the wrong combination of assumptions, and the Dulmage-Mendelsohn decomposition is a great tool for detecting when this has happened.\n", + "\n", + "## Debugging the under-constrained subsystem\n", + "\n", + "The under-constrained system does not decompose into independent subsystems, making it more difficult to debug. However, by inspection, we notice that the same constraints and variables seem to be repeated at each point in the length domain. For each point in space, the \"odd variable out\" seems to be the total flow rate `flow_mass`. Using some intuition about this particular process model, we may conclude that this variable should be calculated from the solid phase velocity, which is constant. We expect an equation that looks like\n", + "```\n", + "flow_mass == velocity * area * density\n", + "```\n", + "\n", + "But this equation isn't here... so we need to add it.\n", + "\n", + "# Fixing the model\n", + "\n", + "We'll start by creating a fresh copy of the model, so we don't accidentally rely on IPOPT's point of termination." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2 = make_model()\n", + "# Make the model square while we try to fix the structural singularity\n", + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.deactivate()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Adding a new particle porosity variable" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.particle_porosity = pyo.Var(\n", + " model2.fs.time, model2.fs.MB.length_domain, initialize=model2.fs.solid_properties.particle_porosity.value\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we need to replace the old particle porosity parameter with this new variable. Luckily, the old parameter is actually implemented as a fixed variable, so we can easily identify all the constraints it participates in with `IncidenceGraphInterface`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyomo.contrib.incidence_analysis import IncidenceGraphInterface\n", + "\n", + "igraph = IncidenceGraphInterface(model2, include_fixed=True)\n", + "porosity_param = model2.fs.solid_properties.particle_porosity\n", + "print(f\"Constraints containing {porosity_param.name}:\")\n", + "for con in igraph.get_adjacent_to(porosity_param):\n", + " print(f\" {con.name}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Particle porosity only appears in two constraints: the density constraint we saw above, and the reaction rate equation. We can replace particle porosity in these constraints using Pyomo's `replace_expressions` function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyomo.core.expr import replace_expressions\n", + "\n", + "for t, x in model2.fs.time * model2.fs.MB.length_domain:\n", + " substitution_map = {id(porosity_param): model2.fs.MB.particle_porosity[t, x]}\n", + " sp = model2.fs.MB.solid_phase\n", + " cons = [sp.properties[t, x].density_particle_constraint, sp.reactions[t, x].gen_rate_expression[\"R1\"]]\n", + " for con in cons:\n", + " con.set_value(replace_expressions(con.expr, substitution_map, descend_into_named_expressions=True))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have added a new `particle_porosity` variable, and are using it in the relevant locations. Now we can move on to adding the missing constraint.\n", + "\n", + "## Adding a new density-flow rate constraint" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@model2.fs.MB.Constraint(model2.fs.time, model2.fs.MB.length_domain)\n", + "def density_flowrate_constraint(mb, t, x):\n", + " return (\n", + " mb.velocity_superficial_solid[t] * mb.bed_area\n", + " * mb.solid_phase.properties[t, x].dens_mass_particle\n", + " == mb.solid_phase.properties[t, x].flow_mass\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Testing the new model\n", + "\n", + "Let's see if these changes have fixed our model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Construct a new diagnostics toolbox\n", + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The structural singularity seems to be gone! Let's unfix our degrees of freedom and see if we can solve." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix()\n", + "model2.fs.MB.gas_phase.properties[0, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix()\n", + "model2.fs.MB.solid_phase.properties[0, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.activate()\n", + "\n", + "res = solver.solve(model2, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This doesn't look much better. What's going on? I thought we just fixed the issue?\n", + "\n", + "# Debugging the model, take two\n", + "\n", + "Let's check the diagnostics toolbox for numerical issues." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.deactivate()\n", + "dt.report_numerical_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looks like we have \"parallel constraints\", which are another form of singularity. Let's follow the toolbox's advice to see what they are." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt.display_near_parallel_constraints()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`density_flowrate_constraint` is the constraint that we added. What is `solid_super_vel`?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.solid_super_vel[0].pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is the same as the constraint we just added! Looks like that constraint already existed at the solid inlet. We can easily deactivate the new constraints at this point in the length domain:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.density_flowrate_constraint[:, 1.0].deactivate();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "But now we have removed constraints from a square model, and expect to have degrees of freedom. Let's see what the diagnostics toolbox has to say." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "But this doesn't help us very much. We have some extraneous degrees of freedom, but with 8881 variables in the under-constrained subsystem, it will be difficult to tell what they are. After some thought (and model-specific intuition), we land on the conclusion that maybe we need to fix particle porosity at the solid inlet. Here, total flow rate is specified, and the `solid_super_vel` equation is using it to compute velocity. So we need `dens_mass_particle` to be known, which means we need `particle_porosity` to be fixed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.particle_porosity[:, 1.0].fix();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's run the diagnostics toolbox as a sanity check." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()\n", + "dt.report_numerical_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looks good! Now we can release our degrees of freedom and try to solve again." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix()\n", + "model2.fs.MB.gas_phase.properties[0, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix()\n", + "model2.fs.MB.solid_phase.properties[0, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.activate()\n", + "\n", + "res = solver.solve(model2, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It worked! For the simple optimization problem we have set up, this solve looks a lot more like what we expect.\n", + "\n", + "# Takeaways from this tutorial\n", + "What have we learned?\n", + "1. IPOPT using non-zero regularization coefficients hints at a singular Jacobian (especially when \"L\"/\"l\" diagnostic tags are present).\n", + "2. When this happens, start by calling `report_structural_issues` to check for a structural singularity. If this looks good, call `report_numerical_issues` to check for a numerical singularity.\n", + "3. When debugging a structural singularity, decomposing a problem into subsystems that each should be nonsingular (e.g. unit models or points in time) is very useful.\n", + "4. The solution to a structural singularity is often to relax a fixed parameter, add a constraint that was forgotten, remove a constraint that was redundant, or fix an extraneous degree of freedom.\n", + "5. Model-specific intuition is usually necessary to diagnose and fix modeling issues. (If you're an algorithm developer, learn about the models you're using! If you don't understand your models, you don't understand your algorithms!)\n", + "6. A modeling issue doesn't necessarily have a unique solution. This is especially true when the issue involves invalid assumptions.\n", + "7. Debugging is an iterative process — fixing one issue can introduce another." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# References\n", + "\n", + "[[1]] Okoli et al., \"A framework for the optimization of chemical looping combustion processes\". *Powder Tech*, 2020.\n", + "\n", + "[[2]] Parker and Biegler, \"Dynamic modeling and nonlinear model predictive control of a moving bed chemical looping combustion reactor\". *IFAC PapersOnline*, 2022.\n", + "\n", + "[[3]] Dulmage and Mendelsohn, \"Coverings of bipartite graphs\". *Can J. Math.*, 1958.\n", + "\n", + "[[4]] Pothen and Fan, \"Computing the block triangular form of a sparse matrix\". *ACM Trans. Math. Softw.*, 1990.\n", + "\n", + "[[5]] Parker et al., \"Applications of the Dulmage-Mendelsohn decomposition for debugging nonlinear optimization problems\". *Comp. Chem. Eng.*, 2023.\n", + "\n", + "[1]: https://www.sciencedirect.com/science/article/pii/S0032591019302803\n", + "[2]: https://www.sciencedirect.com/science/article/pii/S2405896322008825\n", + "[3]: https://www.cambridge.org/core/journals/canadian-journal-of-mathematics/article/coverings-of-bipartite-graphs/413735C5888AB542B92D0C4F402800B1\n", + "[4]: https://dl.acm.org/doi/10.1145/98267.98287\n", + "[5]: https://www.sciencedirect.com/science/article/pii/S0098135423002533\n" + ] + } + ], + "metadata": { + "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.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 3 +} \ No newline at end of file diff --git a/idaes_examples/notebooks/docs/diagnostics/structural_singularity_test.ipynb b/idaes_examples/notebooks/docs/diagnostics/structural_singularity_test.ipynb new file mode 100644 index 00000000..c49f13f8 --- /dev/null +++ b/idaes_examples/notebooks/docs/diagnostics/structural_singularity_test.ipynb @@ -0,0 +1,690 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "###############################################################################\n", + "# The Institute for the Design of Advanced Energy Systems Integrated Platform\n", + "# Framework (IDAES IP) was produced under the DOE Institute for the\n", + "# Design of Advanced Energy Systems (IDAES).\n", + "#\n", + "# Copyright (c) 2018-2023 by the software owners: The Regents of the\n", + "# University of California, through Lawrence Berkeley National Laboratory,\n", + "# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon\n", + "# University, West Virginia University Research Corporation, et al.\n", + "# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md\n", + "# for full copyright and license information.\n", + "###############################################################################" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Debugging a Structural Singularity\n", + "===========================\n", + "Author: Robert Parker\\\n", + "Maintainer: Robert Parker\\\n", + "Updated: 2024-06-10\n", + "\n", + "In this tutorial, we will use the [IDAES Diagnostics Toolbox](https://idaes-pse.readthedocs.io/en/2.4.0/explanations/model_diagnostics/index.html#diagnostics-toolbox)\n", + "to diagnose and fix a structural singularity that is preventing us from solving an optimization problem." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "# Constructing the model\n", + "\n", + "Suppose a collaborator has given us a model to work with. They give us a square model and tell us what the degrees of freedom are. We construct an optimization problem and try to solve it. In this tutorial, we don't want to worry too much about the details that go into constructing the model. This has been provided in the `idaes_examples.mod.diagnostics.gas_solid_contactors.model` module.\n", + "\n", + "## Model details (OKAY TO SKIP)\n", + "\n", + "The model we are trying to optimize is a dynamic model of a moving bed chemical looping combustion reactor. The model has been described by [Okoli et al.][1] and [Parker and Biegler][2]. This is a gas-solid reactor with counter-current flow. The degrees of freedom are gas and solid inlet flow rates, and we are trying to minimize the deviation from a desired operating point via a least-squares objective function.\n", + "\n", + "[1]: https://www.sciencedirect.com/science/article/pii/S0032591019302803\n", + "[2]: https://www.sciencedirect.com/science/article/pii/S2405896322008825\n", + "\n", + "Again, we don't want to worry too much about the model. The `make_model` function will construct the optimization problem that we want to solve, and whenever we do something model-specific, we will explicitly make note of it." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Trying to solve the original model\n", + "\n", + "With that out of the way, let's construct the model and try to solve it!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from idaes_examples.mod.diagnostics.gas_solid_contactors.model import make_model\n", + "import logging\n", + "\n", + "# We'll turn off IDAES logging. This is not recommended in general, but this is an old model\n", + "# (from IDAES 1.7) that has been ported to work with the current version of IDAES. It generates\n", + "# a lot of warnings.\n", + "logging.getLogger(\"idaes\").setLevel(logging.CRITICAL)\n", + "# We'll also turn off Pyomo logging. This will suppress unit inconsistency warnings later,\n", + "# which otherwise flood our console and slow down this notebook. We have unit inconsistencies\n", + "# as, in IDAES 1.7, we didn't rigorously enforce that models use units.\n", + "logging.getLogger(\"pyomo\").setLevel(logging.CRITICAL)\n", + "\n", + "# This constructs a dynamic model with degrees of freedom and an objective function.\n", + "model = make_model()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Before trying to solve the model, let's make sure it conforms to our expectations, i.e. it (a) has degrees of freedom and (b) is well-initialized to a feasible point." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import some useful utilities from the model_statistics module.\n", + "# Degrees of freedom and constraint residuals are always good things to check before\n", + "# trying to solve a simulation or optimization problem.\n", + "from idaes.core.util.model_statistics import degrees_of_freedom, large_residuals_set\n", + "\n", + "dof = degrees_of_freedom(model)\n", + "print(f\"Degrees of freedom: {dof}\")\n", + "has_large_residuals = bool(large_residuals_set(model, tol=1e-5))\n", + "print(f\"Has large residuals: {has_large_residuals}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the above `make_model` function, the model has been \"solved\" to arrive at a feasible point, then degrees of freedom have been unfixed and an objective function has been added to give us an optimization problem. This looks good so far, so let's try to solve the optimization problem." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import pyomo.environ for access to solvers\n", + "import pyomo.environ as pyo\n", + "\n", + "solver = pyo.SolverFactory(\"ipopt\")\n", + "solver.options[\"max_iter\"] = 20\n", + "solver.options[\"print_user_options\"] = \"yes\"\n", + "solver.options[\"OF_print_info_string\"] = \"yes\"\n", + "res = solver.solve(model, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "IPOPT fails to solve the optimization problem... You can try increasing the iteration limit, but it is very unlikely that this model will ever solve. A telltale sign that something is wrong with our model is the persistence of regularization coefficients, that is, numbers in the `lg(rg)` column of the IPOPT log. These coefficients can have multiple causes. One is that the constraint Jacobian (partial derivative matrix) is singular, which indicates a problem with our model. We have set the `print_info_string` option in IPOPT to display \"diagnostic tags\" to help interpret these regularization coefficients. The \"L\" and \"l\" diagnostic tags, which appear repeatedly, indicate that the Jacobian is singular. For more information on IPOPT diagnostic tags, see the IPOPT [documentation](https://coin-or.github.io/Ipopt/OUTPUT.html)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Debugging the original model\n", + "\n", + "Let's run the diagnostics toolbox on the model and see what it has to say.\n", + "\n", + "For good practice, we'll first make sure the model we're debugging is square. Remember that we're assuming we already know how to toggle degrees of freedom in our model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Fix gas and solid flow rates at their respective inlets\n", + "model.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "# Part of our optimization problem was a set of constraints to enforce piecewise\n", + "# constant control inputs. We need to deactivate these as well.\n", + "model.piecewise_constant_constraints.deactivate()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can run the diagnostics toolbox." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from idaes.core.util.model_diagnostics import DiagnosticsToolbox\n", + "\n", + "dt = DiagnosticsToolbox(model)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's look at the warnings we got:\n", + "- Inconsistent units\n", + "- Structural singularity\n", + "- Potential evaluation errors\n", + "\n", + "We'll ignore the inconsistent units. The property package and unit model here were extracted from IDAES 1.7, before we rigorously enforced that all models use units. The potential evaluation errors we see here may be worth looking into, but looking at the failing IPOPT log above, we don't notice any evaluation errors. (If evaluation errors occurred in IPOPT, we would see a message like \"Error in AMPL evaluation\" in the IPOPT iteration log, which we don't see here.) The structural singularity looks like the most promising avenue to debug, especially as the IPOPT log displays persistent regularization coefficients that appear to be caused by a singular Jacobian.\n", + "\n", + "Let's follow the toolbox's advice and display the under and over-constrained sets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt.display_underconstrained_set()\n", + "dt.display_overconstrained_set()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Over and under-constrained subsystems\n", + "\n", + "Structural singularities are characterized by the [Dulmage-Mendelson decomposition][3], which partitions a system into minimal over and under-constrained subsystems. These subsystems contain the potentially unmatched constraints and variables, respectively. Here, \"unmatched\" effectively means \"causing a singularity\". [Pothen and Fan][4] give a good overview of the Dulmage-Mendelsohn decomposition and [Parker et al.][5] give several examples.\n", + "\n", + "[3]: https://www.cambridge.org/core/journals/canadian-journal-of-mathematics/article/coverings-of-bipartite-graphs/413735C5888AB542B92D0C4F402800B1\n", + "[4]: https://dl.acm.org/doi/10.1145/98267.98287\n", + "[5]: https://www.sciencedirect.com/science/article/pii/S0098135423002533\n", + "\n", + "The most straightforward way to fix a structural singularity is to fix variables that are in the under-constrained system and deactivate constraints in the over-constrained subsystem. However, this may not be applicable for every model. For example, we may need to add variables and constraints instead. What over and under-constrained subsystems are telling us is that something is wrong with our modeling assumptions. The particular fix that is appropriate will depend heavily on the model.\n", + "\n", + "If the above output gives us any clues, we can go ahead and start trying to fix things. However, suppose it doesn't. A good strategy is to try to break down the model into smaller, square subsystems that we think should be nonsingular. For a dynamic model like this one, a good candidate is the subsystem of variables and equations at each point in time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# We've included a utility function to extract the subsystem of variables and equations\n", + "# at a specified point in time. If you are dealing with a process flowsheet, here you\n", + "# may want to extract each unit model individually.\n", + "from idaes_examples.mod.diagnostics.util import get_subsystem_at_time\n", + "# TemporarySubsystemManager is used to temporarily fix some variables to make sure\n", + "# we're debugging a square subsystem.\n", + "from pyomo.util.subsystems import TemporarySubsystemManager\n", + "\n", + "# Let's start with t=0. Really, we'd probably want to do this in a loop and try all time points.\n", + "t0 = model.fs.time.first()\n", + "t_block, inputs = get_subsystem_at_time(model, model.fs.time, t0)\n", + "# We'll temporarily fix the \"inputs\" to make sure we have a square system while debugging\n", + "with TemporarySubsystemManager(to_fix=inputs):\n", + " dt = DiagnosticsToolbox(t_block)\n", + " dt.report_structural_issues()\n", + " dt.display_underconstrained_set()\n", + " dt.display_overconstrained_set()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These over and under-constrained subsystems aren't much smaller, but now the over-constrained system decomposes into 10 small, independent blocks. These should be easier to debug.\n", + "\n", + "## Debugging the over-constrained subsystem\n", + "\n", + "To debug the over-constrained subsystem, we look for a constraint that is not calculating any of the variables in the subsystem. The \"odd constraint out\" here seems to be the mass fraction sum, `sum_component_eqn`. This must \"solve for\" one of the mass fractions, which means one of the `material_holdup_calculation` equations must \"solve for\" particle density rather than mass fraction. If we want to see what variables are contained in one of these constraints, we can always `pprint` it:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].sum_component_eqn.pprint()\n", + "model.fs.MB.solid_phase.material_holdup_calculation[0, 0.9, \"Sol\", \"Fe3O4\"].pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If one of these `material_holdup_calculation` equations is solving for particle density, then that means that `density_particle_constraint` is not actually solving for density. Maybe `density_particle_constraint` is over-determining our system?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].density_particle_constraint.pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "But this looks like a very reasonable constraint. After some thought, which admittedly requires some knowledge of the process we are modeling, we decide that the right approach is to make particle porosity a variable. We have assumed that porosity is constant, but this overconstrained subsystem is telling us that this assumption is not valid.\n", + "\n", + "### How did we figure this out? (OKAY TO SKIP)\n", + "Adding a variable (including by unfixing a parameter) to an over-constraining constraint will often remove that constraint from the over-constrained subsystem. But how did we know that this was the right thing to do? If you just care about using the diagnostics toolbox to extract as much information about a singularity as possible, you can skip this section. But if you are curious how we determined that particle porosity should not be constant, read on.\n", + "\n", + "`dens_mass_skeletal` is determined purely by the composition of solid, which is made up of Fe2O3, Fe3O4, and inert Ti2O3. We can view the `density_skeletal_constraint` as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].density_skeletal_constraint.pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we assume a constant particle porosity, this gives us a particle porosity that is also uniquely determined by the solid composition by the above `density_particle_constraint`:\n", + "```\n", + "dens_mass_particle = (1 - porosity) * dens_mass_skeletal\n", + "```\n", + "But the composition of the solid is determined by the (somewhat misnamed) `material_holdup_calculation` constraints. While the name of these constraints implies they \"calculate holdups,\" material holdups at $t=0$ are fixed as initial conditions (because holdups are the differential variables with respect to time in this model). At other time points, we assume that holdups are specified by differential and discretization equations of the model. This means that the `material_holdup_calculation` constraints actually calculate the solid phase mass fractions *from* the holdups. But as we hinted at above, the 4-by-4 system of holdup calculation constraints, `sum_component_eqn` (which simply constrains the sum of mass fractions to be one), mass fractions, and `dens_mass_particle`, uniquely solve for `dens_mass_particle` *as well as* the mass fractions. But if the holdup variables can be used to solve for the mass fractions, they *also* solve for `dens_mass_skeletal`. So both sides of `density_particle_constraint` are already uniquely determined! This implies that we don't need this constraint at all, but we also know that this constraint has to hold. Something has to give. With this in mind, we actually have several options for how to resolve this overspecification:\n", + "1. Remove `density_particle_constraint`. Then we would have `dens_mass_particle` and `dens_mass_skeletal`, with no relationship between them. This would leave us with a mathematically sound model, but with densities that contradict constant particle porosity that we have assumed (which is used elsewhere in the reaction rate calculation equations).\n", + "2. Remove the constraints that calculate skeletal density from composition.\n", + "3. Relax particle porosity from a parameter to a variable.\n", + "\n", + "Options 2 and 3 are equally valid. We've chosen option 3, meaning we assume that the particle \"evolves\" with a density that is well determined from its constituent species, rather than changing density to accommodate whatever mass it accumulates via reaction without altering its volume. This exercise should remind us that all mathematical modeling is somewhat of an art. In the process of choosing the \"least bad\" model, it is fairly easy to over or under-specify something by making the wrong combination of assumptions, and the Dulmage-Mendelsohn decomposition is a great tool for detecting when this has happened.\n", + "\n", + "## Debugging the under-constrained subsystem\n", + "\n", + "The under-constrained system does not decompose into independent subsystems, making it more difficult to debug. However, by inspection, we notice that the same constraints and variables seem to be repeated at each point in the length domain. For each point in space, the \"odd variable out\" seems to be the total flow rate `flow_mass`. Using some intuition about this particular process model, we may conclude that this variable should be calculated from the solid phase velocity, which is constant. We expect an equation that looks like\n", + "```\n", + "flow_mass == velocity * area * density\n", + "```\n", + "\n", + "But this equation isn't here... so we need to add it.\n", + "\n", + "# Fixing the model\n", + "\n", + "We'll start by creating a fresh copy of the model, so we don't accidentally rely on IPOPT's point of termination." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2 = make_model()\n", + "# Make the model square while we try to fix the structural singularity\n", + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.deactivate()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Adding a new particle porosity variable" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.particle_porosity = pyo.Var(\n", + " model2.fs.time, model2.fs.MB.length_domain, initialize=model2.fs.solid_properties.particle_porosity.value\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we need to replace the old particle porosity parameter with this new variable. Luckily, the old parameter is actually implemented as a fixed variable, so we can easily identify all the constraints it participates in with `IncidenceGraphInterface`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyomo.contrib.incidence_analysis import IncidenceGraphInterface\n", + "\n", + "igraph = IncidenceGraphInterface(model2, include_fixed=True)\n", + "porosity_param = model2.fs.solid_properties.particle_porosity\n", + "print(f\"Constraints containing {porosity_param.name}:\")\n", + "for con in igraph.get_adjacent_to(porosity_param):\n", + " print(f\" {con.name}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Particle porosity only appears in two constraints: the density constraint we saw above, and the reaction rate equation. We can replace particle porosity in these constraints using Pyomo's `replace_expressions` function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyomo.core.expr import replace_expressions\n", + "\n", + "for t, x in model2.fs.time * model2.fs.MB.length_domain:\n", + " substitution_map = {id(porosity_param): model2.fs.MB.particle_porosity[t, x]}\n", + " sp = model2.fs.MB.solid_phase\n", + " cons = [sp.properties[t, x].density_particle_constraint, sp.reactions[t, x].gen_rate_expression[\"R1\"]]\n", + " for con in cons:\n", + " con.set_value(replace_expressions(con.expr, substitution_map, descend_into_named_expressions=True))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have added a new `particle_porosity` variable, and are using it in the relevant locations. Now we can move on to adding the missing constraint.\n", + "\n", + "## Adding a new density-flow rate constraint" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@model2.fs.MB.Constraint(model2.fs.time, model2.fs.MB.length_domain)\n", + "def density_flowrate_constraint(mb, t, x):\n", + " return (\n", + " mb.velocity_superficial_solid[t] * mb.bed_area\n", + " * mb.solid_phase.properties[t, x].dens_mass_particle\n", + " == mb.solid_phase.properties[t, x].flow_mass\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Testing the new model\n", + "\n", + "Let's see if these changes have fixed our model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Construct a new diagnostics toolbox\n", + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The structural singularity seems to be gone! Let's unfix our degrees of freedom and see if we can solve." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix()\n", + "model2.fs.MB.gas_phase.properties[0, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix()\n", + "model2.fs.MB.solid_phase.properties[0, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.activate()\n", + "\n", + "res = solver.solve(model2, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This doesn't look much better. What's going on? I thought we just fixed the issue?\n", + "\n", + "# Debugging the model, take two\n", + "\n", + "Let's check the diagnostics toolbox for numerical issues." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.deactivate()\n", + "dt.report_numerical_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looks like we have \"parallel constraints\", which are another form of singularity. Let's follow the toolbox's advice to see what they are." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt.display_near_parallel_constraints()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`density_flowrate_constraint` is the constraint that we added. What is `solid_super_vel`?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.solid_super_vel[0].pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is the same as the constraint we just added! Looks like that constraint already existed at the solid inlet. We can easily deactivate the new constraints at this point in the length domain:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.density_flowrate_constraint[:, 1.0].deactivate();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "But now we have removed constraints from a square model, and expect to have degrees of freedom. Let's see what the diagnostics toolbox has to say." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "But this doesn't help us very much. We have some extraneous degrees of freedom, but with 8881 variables in the under-constrained subsystem, it will be difficult to tell what they are. After some thought (and model-specific intuition), we land on the conclusion that maybe we need to fix particle porosity at the solid inlet. Here, total flow rate is specified, and the `solid_super_vel` equation is using it to compute velocity. So we need `dens_mass_particle` to be known, which means we need `particle_porosity` to be fixed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.particle_porosity[:, 1.0].fix();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's run the diagnostics toolbox as a sanity check." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()\n", + "dt.report_numerical_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looks good! Now we can release our degrees of freedom and try to solve again." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix()\n", + "model2.fs.MB.gas_phase.properties[0, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix()\n", + "model2.fs.MB.solid_phase.properties[0, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.activate()\n", + "\n", + "res = solver.solve(model2, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It worked! For the simple optimization problem we have set up, this solve looks a lot more like what we expect.\n", + "\n", + "# Takeaways from this tutorial\n", + "What have we learned?\n", + "1. IPOPT using non-zero regularization coefficients hints at a singular Jacobian (especially when \"L\"/\"l\" diagnostic tags are present).\n", + "2. When this happens, start by calling `report_structural_issues` to check for a structural singularity. If this looks good, call `report_numerical_issues` to check for a numerical singularity.\n", + "3. When debugging a structural singularity, decomposing a problem into subsystems that each should be nonsingular (e.g. unit models or points in time) is very useful.\n", + "4. The solution to a structural singularity is often to relax a fixed parameter, add a constraint that was forgotten, remove a constraint that was redundant, or fix an extraneous degree of freedom.\n", + "5. Model-specific intuition is usually necessary to diagnose and fix modeling issues. (If you're an algorithm developer, learn about the models you're using! If you don't understand your models, you don't understand your algorithms!)\n", + "6. A modeling issue doesn't necessarily have a unique solution. This is especially true when the issue involves invalid assumptions.\n", + "7. Debugging is an iterative process — fixing one issue can introduce another." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# References\n", + "\n", + "[[1]] Okoli et al., \"A framework for the optimization of chemical looping combustion processes\". *Powder Tech*, 2020.\n", + "\n", + "[[2]] Parker and Biegler, \"Dynamic modeling and nonlinear model predictive control of a moving bed chemical looping combustion reactor\". *IFAC PapersOnline*, 2022.\n", + "\n", + "[[3]] Dulmage and Mendelsohn, \"Coverings of bipartite graphs\". *Can J. Math.*, 1958.\n", + "\n", + "[[4]] Pothen and Fan, \"Computing the block triangular form of a sparse matrix\". *ACM Trans. Math. Softw.*, 1990.\n", + "\n", + "[[5]] Parker et al., \"Applications of the Dulmage-Mendelsohn decomposition for debugging nonlinear optimization problems\". *Comp. Chem. Eng.*, 2023.\n", + "\n", + "[1]: https://www.sciencedirect.com/science/article/pii/S0032591019302803\n", + "[2]: https://www.sciencedirect.com/science/article/pii/S2405896322008825\n", + "[3]: https://www.cambridge.org/core/journals/canadian-journal-of-mathematics/article/coverings-of-bipartite-graphs/413735C5888AB542B92D0C4F402800B1\n", + "[4]: https://dl.acm.org/doi/10.1145/98267.98287\n", + "[5]: https://www.sciencedirect.com/science/article/pii/S0098135423002533\n" + ] + } + ], + "metadata": { + "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.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 3 +} \ No newline at end of file diff --git a/idaes_examples/notebooks/docs/diagnostics/structural_singularity_usr.ipynb b/idaes_examples/notebooks/docs/diagnostics/structural_singularity_usr.ipynb new file mode 100644 index 00000000..c49f13f8 --- /dev/null +++ b/idaes_examples/notebooks/docs/diagnostics/structural_singularity_usr.ipynb @@ -0,0 +1,690 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "###############################################################################\n", + "# The Institute for the Design of Advanced Energy Systems Integrated Platform\n", + "# Framework (IDAES IP) was produced under the DOE Institute for the\n", + "# Design of Advanced Energy Systems (IDAES).\n", + "#\n", + "# Copyright (c) 2018-2023 by the software owners: The Regents of the\n", + "# University of California, through Lawrence Berkeley National Laboratory,\n", + "# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon\n", + "# University, West Virginia University Research Corporation, et al.\n", + "# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md\n", + "# for full copyright and license information.\n", + "###############################################################################" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Debugging a Structural Singularity\n", + "===========================\n", + "Author: Robert Parker\\\n", + "Maintainer: Robert Parker\\\n", + "Updated: 2024-06-10\n", + "\n", + "In this tutorial, we will use the [IDAES Diagnostics Toolbox](https://idaes-pse.readthedocs.io/en/2.4.0/explanations/model_diagnostics/index.html#diagnostics-toolbox)\n", + "to diagnose and fix a structural singularity that is preventing us from solving an optimization problem." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "# Constructing the model\n", + "\n", + "Suppose a collaborator has given us a model to work with. They give us a square model and tell us what the degrees of freedom are. We construct an optimization problem and try to solve it. In this tutorial, we don't want to worry too much about the details that go into constructing the model. This has been provided in the `idaes_examples.mod.diagnostics.gas_solid_contactors.model` module.\n", + "\n", + "## Model details (OKAY TO SKIP)\n", + "\n", + "The model we are trying to optimize is a dynamic model of a moving bed chemical looping combustion reactor. The model has been described by [Okoli et al.][1] and [Parker and Biegler][2]. This is a gas-solid reactor with counter-current flow. The degrees of freedom are gas and solid inlet flow rates, and we are trying to minimize the deviation from a desired operating point via a least-squares objective function.\n", + "\n", + "[1]: https://www.sciencedirect.com/science/article/pii/S0032591019302803\n", + "[2]: https://www.sciencedirect.com/science/article/pii/S2405896322008825\n", + "\n", + "Again, we don't want to worry too much about the model. The `make_model` function will construct the optimization problem that we want to solve, and whenever we do something model-specific, we will explicitly make note of it." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Trying to solve the original model\n", + "\n", + "With that out of the way, let's construct the model and try to solve it!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from idaes_examples.mod.diagnostics.gas_solid_contactors.model import make_model\n", + "import logging\n", + "\n", + "# We'll turn off IDAES logging. This is not recommended in general, but this is an old model\n", + "# (from IDAES 1.7) that has been ported to work with the current version of IDAES. It generates\n", + "# a lot of warnings.\n", + "logging.getLogger(\"idaes\").setLevel(logging.CRITICAL)\n", + "# We'll also turn off Pyomo logging. This will suppress unit inconsistency warnings later,\n", + "# which otherwise flood our console and slow down this notebook. We have unit inconsistencies\n", + "# as, in IDAES 1.7, we didn't rigorously enforce that models use units.\n", + "logging.getLogger(\"pyomo\").setLevel(logging.CRITICAL)\n", + "\n", + "# This constructs a dynamic model with degrees of freedom and an objective function.\n", + "model = make_model()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Before trying to solve the model, let's make sure it conforms to our expectations, i.e. it (a) has degrees of freedom and (b) is well-initialized to a feasible point." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import some useful utilities from the model_statistics module.\n", + "# Degrees of freedom and constraint residuals are always good things to check before\n", + "# trying to solve a simulation or optimization problem.\n", + "from idaes.core.util.model_statistics import degrees_of_freedom, large_residuals_set\n", + "\n", + "dof = degrees_of_freedom(model)\n", + "print(f\"Degrees of freedom: {dof}\")\n", + "has_large_residuals = bool(large_residuals_set(model, tol=1e-5))\n", + "print(f\"Has large residuals: {has_large_residuals}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the above `make_model` function, the model has been \"solved\" to arrive at a feasible point, then degrees of freedom have been unfixed and an objective function has been added to give us an optimization problem. This looks good so far, so let's try to solve the optimization problem." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import pyomo.environ for access to solvers\n", + "import pyomo.environ as pyo\n", + "\n", + "solver = pyo.SolverFactory(\"ipopt\")\n", + "solver.options[\"max_iter\"] = 20\n", + "solver.options[\"print_user_options\"] = \"yes\"\n", + "solver.options[\"OF_print_info_string\"] = \"yes\"\n", + "res = solver.solve(model, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "IPOPT fails to solve the optimization problem... You can try increasing the iteration limit, but it is very unlikely that this model will ever solve. A telltale sign that something is wrong with our model is the persistence of regularization coefficients, that is, numbers in the `lg(rg)` column of the IPOPT log. These coefficients can have multiple causes. One is that the constraint Jacobian (partial derivative matrix) is singular, which indicates a problem with our model. We have set the `print_info_string` option in IPOPT to display \"diagnostic tags\" to help interpret these regularization coefficients. The \"L\" and \"l\" diagnostic tags, which appear repeatedly, indicate that the Jacobian is singular. For more information on IPOPT diagnostic tags, see the IPOPT [documentation](https://coin-or.github.io/Ipopt/OUTPUT.html)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Debugging the original model\n", + "\n", + "Let's run the diagnostics toolbox on the model and see what it has to say.\n", + "\n", + "For good practice, we'll first make sure the model we're debugging is square. Remember that we're assuming we already know how to toggle degrees of freedom in our model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Fix gas and solid flow rates at their respective inlets\n", + "model.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "# Part of our optimization problem was a set of constraints to enforce piecewise\n", + "# constant control inputs. We need to deactivate these as well.\n", + "model.piecewise_constant_constraints.deactivate()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can run the diagnostics toolbox." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from idaes.core.util.model_diagnostics import DiagnosticsToolbox\n", + "\n", + "dt = DiagnosticsToolbox(model)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's look at the warnings we got:\n", + "- Inconsistent units\n", + "- Structural singularity\n", + "- Potential evaluation errors\n", + "\n", + "We'll ignore the inconsistent units. The property package and unit model here were extracted from IDAES 1.7, before we rigorously enforced that all models use units. The potential evaluation errors we see here may be worth looking into, but looking at the failing IPOPT log above, we don't notice any evaluation errors. (If evaluation errors occurred in IPOPT, we would see a message like \"Error in AMPL evaluation\" in the IPOPT iteration log, which we don't see here.) The structural singularity looks like the most promising avenue to debug, especially as the IPOPT log displays persistent regularization coefficients that appear to be caused by a singular Jacobian.\n", + "\n", + "Let's follow the toolbox's advice and display the under and over-constrained sets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt.display_underconstrained_set()\n", + "dt.display_overconstrained_set()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Over and under-constrained subsystems\n", + "\n", + "Structural singularities are characterized by the [Dulmage-Mendelson decomposition][3], which partitions a system into minimal over and under-constrained subsystems. These subsystems contain the potentially unmatched constraints and variables, respectively. Here, \"unmatched\" effectively means \"causing a singularity\". [Pothen and Fan][4] give a good overview of the Dulmage-Mendelsohn decomposition and [Parker et al.][5] give several examples.\n", + "\n", + "[3]: https://www.cambridge.org/core/journals/canadian-journal-of-mathematics/article/coverings-of-bipartite-graphs/413735C5888AB542B92D0C4F402800B1\n", + "[4]: https://dl.acm.org/doi/10.1145/98267.98287\n", + "[5]: https://www.sciencedirect.com/science/article/pii/S0098135423002533\n", + "\n", + "The most straightforward way to fix a structural singularity is to fix variables that are in the under-constrained system and deactivate constraints in the over-constrained subsystem. However, this may not be applicable for every model. For example, we may need to add variables and constraints instead. What over and under-constrained subsystems are telling us is that something is wrong with our modeling assumptions. The particular fix that is appropriate will depend heavily on the model.\n", + "\n", + "If the above output gives us any clues, we can go ahead and start trying to fix things. However, suppose it doesn't. A good strategy is to try to break down the model into smaller, square subsystems that we think should be nonsingular. For a dynamic model like this one, a good candidate is the subsystem of variables and equations at each point in time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# We've included a utility function to extract the subsystem of variables and equations\n", + "# at a specified point in time. If you are dealing with a process flowsheet, here you\n", + "# may want to extract each unit model individually.\n", + "from idaes_examples.mod.diagnostics.util import get_subsystem_at_time\n", + "# TemporarySubsystemManager is used to temporarily fix some variables to make sure\n", + "# we're debugging a square subsystem.\n", + "from pyomo.util.subsystems import TemporarySubsystemManager\n", + "\n", + "# Let's start with t=0. Really, we'd probably want to do this in a loop and try all time points.\n", + "t0 = model.fs.time.first()\n", + "t_block, inputs = get_subsystem_at_time(model, model.fs.time, t0)\n", + "# We'll temporarily fix the \"inputs\" to make sure we have a square system while debugging\n", + "with TemporarySubsystemManager(to_fix=inputs):\n", + " dt = DiagnosticsToolbox(t_block)\n", + " dt.report_structural_issues()\n", + " dt.display_underconstrained_set()\n", + " dt.display_overconstrained_set()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These over and under-constrained subsystems aren't much smaller, but now the over-constrained system decomposes into 10 small, independent blocks. These should be easier to debug.\n", + "\n", + "## Debugging the over-constrained subsystem\n", + "\n", + "To debug the over-constrained subsystem, we look for a constraint that is not calculating any of the variables in the subsystem. The \"odd constraint out\" here seems to be the mass fraction sum, `sum_component_eqn`. This must \"solve for\" one of the mass fractions, which means one of the `material_holdup_calculation` equations must \"solve for\" particle density rather than mass fraction. If we want to see what variables are contained in one of these constraints, we can always `pprint` it:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].sum_component_eqn.pprint()\n", + "model.fs.MB.solid_phase.material_holdup_calculation[0, 0.9, \"Sol\", \"Fe3O4\"].pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If one of these `material_holdup_calculation` equations is solving for particle density, then that means that `density_particle_constraint` is not actually solving for density. Maybe `density_particle_constraint` is over-determining our system?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].density_particle_constraint.pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "But this looks like a very reasonable constraint. After some thought, which admittedly requires some knowledge of the process we are modeling, we decide that the right approach is to make particle porosity a variable. We have assumed that porosity is constant, but this overconstrained subsystem is telling us that this assumption is not valid.\n", + "\n", + "### How did we figure this out? (OKAY TO SKIP)\n", + "Adding a variable (including by unfixing a parameter) to an over-constraining constraint will often remove that constraint from the over-constrained subsystem. But how did we know that this was the right thing to do? If you just care about using the diagnostics toolbox to extract as much information about a singularity as possible, you can skip this section. But if you are curious how we determined that particle porosity should not be constant, read on.\n", + "\n", + "`dens_mass_skeletal` is determined purely by the composition of solid, which is made up of Fe2O3, Fe3O4, and inert Ti2O3. We can view the `density_skeletal_constraint` as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.fs.MB.solid_phase.properties[0, 0.9].density_skeletal_constraint.pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we assume a constant particle porosity, this gives us a particle porosity that is also uniquely determined by the solid composition by the above `density_particle_constraint`:\n", + "```\n", + "dens_mass_particle = (1 - porosity) * dens_mass_skeletal\n", + "```\n", + "But the composition of the solid is determined by the (somewhat misnamed) `material_holdup_calculation` constraints. While the name of these constraints implies they \"calculate holdups,\" material holdups at $t=0$ are fixed as initial conditions (because holdups are the differential variables with respect to time in this model). At other time points, we assume that holdups are specified by differential and discretization equations of the model. This means that the `material_holdup_calculation` constraints actually calculate the solid phase mass fractions *from* the holdups. But as we hinted at above, the 4-by-4 system of holdup calculation constraints, `sum_component_eqn` (which simply constrains the sum of mass fractions to be one), mass fractions, and `dens_mass_particle`, uniquely solve for `dens_mass_particle` *as well as* the mass fractions. But if the holdup variables can be used to solve for the mass fractions, they *also* solve for `dens_mass_skeletal`. So both sides of `density_particle_constraint` are already uniquely determined! This implies that we don't need this constraint at all, but we also know that this constraint has to hold. Something has to give. With this in mind, we actually have several options for how to resolve this overspecification:\n", + "1. Remove `density_particle_constraint`. Then we would have `dens_mass_particle` and `dens_mass_skeletal`, with no relationship between them. This would leave us with a mathematically sound model, but with densities that contradict constant particle porosity that we have assumed (which is used elsewhere in the reaction rate calculation equations).\n", + "2. Remove the constraints that calculate skeletal density from composition.\n", + "3. Relax particle porosity from a parameter to a variable.\n", + "\n", + "Options 2 and 3 are equally valid. We've chosen option 3, meaning we assume that the particle \"evolves\" with a density that is well determined from its constituent species, rather than changing density to accommodate whatever mass it accumulates via reaction without altering its volume. This exercise should remind us that all mathematical modeling is somewhat of an art. In the process of choosing the \"least bad\" model, it is fairly easy to over or under-specify something by making the wrong combination of assumptions, and the Dulmage-Mendelsohn decomposition is a great tool for detecting when this has happened.\n", + "\n", + "## Debugging the under-constrained subsystem\n", + "\n", + "The under-constrained system does not decompose into independent subsystems, making it more difficult to debug. However, by inspection, we notice that the same constraints and variables seem to be repeated at each point in the length domain. For each point in space, the \"odd variable out\" seems to be the total flow rate `flow_mass`. Using some intuition about this particular process model, we may conclude that this variable should be calculated from the solid phase velocity, which is constant. We expect an equation that looks like\n", + "```\n", + "flow_mass == velocity * area * density\n", + "```\n", + "\n", + "But this equation isn't here... so we need to add it.\n", + "\n", + "# Fixing the model\n", + "\n", + "We'll start by creating a fresh copy of the model, so we don't accidentally rely on IPOPT's point of termination." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2 = make_model()\n", + "# Make the model square while we try to fix the structural singularity\n", + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.deactivate()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Adding a new particle porosity variable" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.particle_porosity = pyo.Var(\n", + " model2.fs.time, model2.fs.MB.length_domain, initialize=model2.fs.solid_properties.particle_porosity.value\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we need to replace the old particle porosity parameter with this new variable. Luckily, the old parameter is actually implemented as a fixed variable, so we can easily identify all the constraints it participates in with `IncidenceGraphInterface`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyomo.contrib.incidence_analysis import IncidenceGraphInterface\n", + "\n", + "igraph = IncidenceGraphInterface(model2, include_fixed=True)\n", + "porosity_param = model2.fs.solid_properties.particle_porosity\n", + "print(f\"Constraints containing {porosity_param.name}:\")\n", + "for con in igraph.get_adjacent_to(porosity_param):\n", + " print(f\" {con.name}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Particle porosity only appears in two constraints: the density constraint we saw above, and the reaction rate equation. We can replace particle porosity in these constraints using Pyomo's `replace_expressions` function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyomo.core.expr import replace_expressions\n", + "\n", + "for t, x in model2.fs.time * model2.fs.MB.length_domain:\n", + " substitution_map = {id(porosity_param): model2.fs.MB.particle_porosity[t, x]}\n", + " sp = model2.fs.MB.solid_phase\n", + " cons = [sp.properties[t, x].density_particle_constraint, sp.reactions[t, x].gen_rate_expression[\"R1\"]]\n", + " for con in cons:\n", + " con.set_value(replace_expressions(con.expr, substitution_map, descend_into_named_expressions=True))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have added a new `particle_porosity` variable, and are using it in the relevant locations. Now we can move on to adding the missing constraint.\n", + "\n", + "## Adding a new density-flow rate constraint" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@model2.fs.MB.Constraint(model2.fs.time, model2.fs.MB.length_domain)\n", + "def density_flowrate_constraint(mb, t, x):\n", + " return (\n", + " mb.velocity_superficial_solid[t] * mb.bed_area\n", + " * mb.solid_phase.properties[t, x].dens_mass_particle\n", + " == mb.solid_phase.properties[t, x].flow_mass\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Testing the new model\n", + "\n", + "Let's see if these changes have fixed our model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Construct a new diagnostics toolbox\n", + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The structural singularity seems to be gone! Let's unfix our degrees of freedom and see if we can solve." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix()\n", + "model2.fs.MB.gas_phase.properties[0, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix()\n", + "model2.fs.MB.solid_phase.properties[0, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.activate()\n", + "\n", + "res = solver.solve(model2, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This doesn't look much better. What's going on? I thought we just fixed the issue?\n", + "\n", + "# Debugging the model, take two\n", + "\n", + "Let's check the diagnostics toolbox for numerical issues." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.deactivate()\n", + "dt.report_numerical_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looks like we have \"parallel constraints\", which are another form of singularity. Let's follow the toolbox's advice to see what they are." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt.display_near_parallel_constraints()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`density_flowrate_constraint` is the constraint that we added. What is `solid_super_vel`?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.solid_super_vel[0].pprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is the same as the constraint we just added! Looks like that constraint already existed at the solid inlet. We can easily deactivate the new constraints at this point in the length domain:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.density_flowrate_constraint[:, 1.0].deactivate();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "But now we have removed constraints from a square model, and expect to have degrees of freedom. Let's see what the diagnostics toolbox has to say." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "But this doesn't help us very much. We have some extraneous degrees of freedom, but with 8881 variables in the under-constrained subsystem, it will be difficult to tell what they are. After some thought (and model-specific intuition), we land on the conclusion that maybe we need to fix particle porosity at the solid inlet. Here, total flow rate is specified, and the `solid_super_vel` equation is using it to compute velocity. So we need `dens_mass_particle` to be known, which means we need `particle_porosity` to be fixed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.particle_porosity[:, 1.0].fix();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's run the diagnostics toolbox as a sanity check." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = DiagnosticsToolbox(model2)\n", + "dt.report_structural_issues()\n", + "dt.report_numerical_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looks good! Now we can release our degrees of freedom and try to solve again." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model2.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix()\n", + "model2.fs.MB.gas_phase.properties[0, 0].flow_mol.fix()\n", + "model2.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix()\n", + "model2.fs.MB.solid_phase.properties[0, 1].flow_mass.fix()\n", + "model2.piecewise_constant_constraints.activate()\n", + "\n", + "res = solver.solve(model2, tee=True)\n", + "print(f\"Converged successfully: {pyo.check_optimal_termination(res)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It worked! For the simple optimization problem we have set up, this solve looks a lot more like what we expect.\n", + "\n", + "# Takeaways from this tutorial\n", + "What have we learned?\n", + "1. IPOPT using non-zero regularization coefficients hints at a singular Jacobian (especially when \"L\"/\"l\" diagnostic tags are present).\n", + "2. When this happens, start by calling `report_structural_issues` to check for a structural singularity. If this looks good, call `report_numerical_issues` to check for a numerical singularity.\n", + "3. When debugging a structural singularity, decomposing a problem into subsystems that each should be nonsingular (e.g. unit models or points in time) is very useful.\n", + "4. The solution to a structural singularity is often to relax a fixed parameter, add a constraint that was forgotten, remove a constraint that was redundant, or fix an extraneous degree of freedom.\n", + "5. Model-specific intuition is usually necessary to diagnose and fix modeling issues. (If you're an algorithm developer, learn about the models you're using! If you don't understand your models, you don't understand your algorithms!)\n", + "6. A modeling issue doesn't necessarily have a unique solution. This is especially true when the issue involves invalid assumptions.\n", + "7. Debugging is an iterative process — fixing one issue can introduce another." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# References\n", + "\n", + "[[1]] Okoli et al., \"A framework for the optimization of chemical looping combustion processes\". *Powder Tech*, 2020.\n", + "\n", + "[[2]] Parker and Biegler, \"Dynamic modeling and nonlinear model predictive control of a moving bed chemical looping combustion reactor\". *IFAC PapersOnline*, 2022.\n", + "\n", + "[[3]] Dulmage and Mendelsohn, \"Coverings of bipartite graphs\". *Can J. Math.*, 1958.\n", + "\n", + "[[4]] Pothen and Fan, \"Computing the block triangular form of a sparse matrix\". *ACM Trans. Math. Softw.*, 1990.\n", + "\n", + "[[5]] Parker et al., \"Applications of the Dulmage-Mendelsohn decomposition for debugging nonlinear optimization problems\". *Comp. Chem. Eng.*, 2023.\n", + "\n", + "[1]: https://www.sciencedirect.com/science/article/pii/S0032591019302803\n", + "[2]: https://www.sciencedirect.com/science/article/pii/S2405896322008825\n", + "[3]: https://www.cambridge.org/core/journals/canadian-journal-of-mathematics/article/coverings-of-bipartite-graphs/413735C5888AB542B92D0C4F402800B1\n", + "[4]: https://dl.acm.org/doi/10.1145/98267.98287\n", + "[5]: https://www.sciencedirect.com/science/article/pii/S0098135423002533\n" + ] + } + ], + "metadata": { + "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.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 3 +} \ No newline at end of file