###############################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES).
#
# Copyright (c) 2018-2023 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.md and LICENSE.md
# for full copyright and license information.
###############################################################################
Creating Custom Unit Model#
Author: Javal Vyas
Maintainer: Javal Vyas
Updated: 2023-02-20
This tutorial is a comprehensive step-wise procedure to build a custom unit model from scratch. This tutorial will include creating a property package, a custom unit model and testing them. For this tutorial we shall create a custom unit model for Liquid - Liquid Extraction.
The Liquid - Liquid Extractor model contains two immiscible fluids forming the two phases. One of the phases, say phase_1 has a high concentration of solutes which is to be separated. A mass transfer happens between the two phases and the solute is transferred from phase_1 to phase_2. This mass transfer is governed by a parameter called the distribution coefficient.
After reviewing the working principles of the Liquid - Liquid Extractor, we shall proceed to create a custom unit model. We will require a property package for each phase, a custom unit model class and tests for the model and property packages.
Before commencing the development of the model, we need to state some assumptions which the following unit model will be using.
Steady-state only
Organic phase property package has a single phase named Org
Aqueous phase property package has a single phase named Aq
Organic and Aqueous phase properties need not have the same component list.
Thus as per the assumptions, we will be creating one property package for the aqueous phase (Aq), and the other for the Organic phase (Org).
1. Creating Organic Property Package#
Creating a property package is a 4 step process
Import necessary libraries
Creating Physical Parameter Data Block
Define State Block
Define State Block Data
1.1 Importing necessary packages#
Let us begin with the importing the necessary libraries where we will be using functionalities from IDAES and Pyomo.
# Import Python libraries
import logging
import idaes.logger as idaeslog
from idaes.core.util.initialization import fix_state_vars
# Import Pyomo libraries
from pyomo.environ import (
Param,
Set,
Var,
NonNegativeReals,
units,
Expression,
PositiveReals,
)
# Import IDAES cores
from idaes.core import (
declare_process_block_class,
MaterialFlowBasis,
PhysicalParameterBlock,
StateBlockData,
StateBlock,
MaterialBalanceType,
EnergyBalanceType,
Solute,
Solvent,
LiquidPhase,
)
from idaes.core.util.model_statistics import degrees_of_freedom
1.2 Physical Parameter Data Block#
A PhysicalParameterBlock serves as the central point of reference for all aspects of the property package and needs to define several things about the package. These are summarized below:
Units of measurement
What properties are supported and how they are implemented
What components and phases are included in the packages
All the global parameters necessary for calculating properties
A reference to the associated State Block class, so that construction of the State Block components can be automated from the Physical Parameter Block
To construct this block, we begin by declaring a process block class using a Python decorator. One can learn more about declare_process_block_class here. After constructing the process block, we define a build function which contains all the components that the property package would have. super function here is used to give access to methods and properties of a parent or sibling class and since this is used on the class PhysicalParameterData class, build has access to all the parent and sibling class methods.
The PhysicalParameterBlock then refers to the state block, in this case OrgPhaseStateBlock (which will be declared later), so that we can build a state block instance by only knowing the PhysicalParameterBlock we wish to use. Then we move on to list the number of phases in this property package. Then we assign the variable to the phase which follows a naming convention. Like here since the solvent is in the Organic phase, we will assign the Phase as OrganicPhase and the variable will be named Org as per the naming convention. The details of naming conventions can be found here. We will be following the same convention throughout the example.
After defining the list of the phases, we move on to list the components and their type in the phase. It can be a solute or a solvent in the Organic phase. Thus, we define the component and assign it to either being a solute or a solvent. In this case, the salts are the solutes and Ethylene dibromide is the solvent. Next, we define the physical properties involved in the package, like the heat capacity and density of the solvent, the reference temperature, and the distribution factor that would govern the mass transfer from one phase into another. Additionally, a parameter, the diffusion_factor, is introduced. This factor plays a crucial role in governing mass transfer between phases, necessitating its definition within the state block.
The final step in creating the Physical Parameter Block is to declare a classmethod named define_metadata, which takes two arguments: a class (cls) and an instance of that class (obj). In this method, we will call the predefined method add_default_units().
obj.add_default_units()sets the default units metadata for the property package, and here we define units to be used with this property package as default.
@declare_process_block_class("OrgPhase")
class PhysicalParameterData(PhysicalParameterBlock):
"""
Property Parameter Block Class
Contains parameters and indexing sets associated with properties for
organic Phase
"""
def build(self):
"""
Callable method for Block construction.
"""
super().build()
self._state_block_class = OrgPhaseStateBlock
# List of valid phases in property package
self.Org = LiquidPhase()
# Component list - a list of component identifiers
self.NaCl = Solute()
self.KNO3 = Solute()
self.CaSO4 = Solute()
self.solvent = (
Solvent()
) # Solvent used here is ethylene dibromide (Organic Polar)
# Heat capacity of solvent
self.cp_mass = Param(
mutable=True,
initialize=717.01,
doc="Specific heat capacity of solvent",
units=units.J / units.kg / units.K,
)
self.dens_mass = Param(
mutable=True,
initialize=2170,
doc="Density of ethylene dibromide",
units=units.kg / units.m**3,
)
self.temperature_ref = Param(
within=PositiveReals,
mutable=True,
default=298.15,
doc="Reference temperature",
units=units.K,
)
self.diffusion_factor = Param(
self.solute_set,
initialize={"NaCl": 2.15, "KNO3": 3, "CaSO4": 1.5},
within=PositiveReals,
mutable=True,
)
@classmethod
def define_metadata(cls, obj):
obj.add_default_units(
{
"time": units.hour,
"length": units.m,
"mass": units.g,
"amount": units.mol,
"temperature": units.K,
}
)
1.3 State Block#
After the PhysicalParameterBlock class has been created, the next step is to write the code necessary to create the State Blocks that will be used throughout the flowsheet. StateBlock contains all the information necessary to define the state of the system. This includes the state variables and constraints on those variables which are used to describe a state property like the enthalpy, material balance, etc.
Creating a State Block requires us to write two classes. The reason we write two classes is because of the inherent nature of how declare_process_block_data works. declare_process_block_data facilitates creating an IndexedComponent object which can handle multiple ComponentData objects which represent the component at each point in the indexing set. This makes it easier to build an instance of the model at each indexed point. However, State Blocks are slightly different, as they are always indexed (at least by time). Due to this, we often want to perform actions on all the elements of the indexed StateBlock all at once (rather than element by element).
The class _OrganicStateBlock is defined without the declare_process_block_data decorator and thus works as a traditional class and this facilitates performing a method on the class as a whole rather than individual elements of the indexed property blocks. In this class we define the fix_initialization_states function. fix_initialization_states function is to used to fix the state variable within the state block with the provided initial values (usually inlet conditions). It takes a block as the argument in which the state variables are to be fixed. It also takes state_args as an optional argument. state_args is a dictionary with the value for the state variables to be fixed. This function returns a dictionary indexed by the block, state variables and variable index indicating the fixed status of each variable before applying the function.
The above function comprise of the _OrganicStateBlock, next we shall see the construction of the OrgPhaseStateBlockData class.
class _OrganicStateBlock(StateBlock):
"""
This Class contains methods which should be applied to Property Blocks as a
whole, rather than individual elements of indexed Property Blocks.
"""
def fix_initialization_states(self):
fix_state_vars(self)
The class OrgPhaseStateBlockData is designated with the declare_process_block_class decorator, named OrgPhaseStateBlock, and inherits the block class from _OrganicStateBlock. This inheritance allows OrgPhaseStateBlockData to leverage functions from _OrganicStateBlock. Following the class definition, a build function similar to the one used in the PhysicalParameterData block is employed. The super function is utilized to enable the utilization of functions from the parent or sibling class.
The subsequent objective is to delineate the state variables, accomplished through the _make_state_vars method. This method encompasses all the essential state variables and associated data. For this particular property package, the required state variables are:
flow_vol- volumetric flow rateconc_mass_comp- mass fractionspressure- state pressuretemperature- state temperature
After establishing the state variables, the subsequent step involves setting up state properties as constraints. This includes specifying the relationships and limitations that dictate the system’s behavior. The following properties need to be articulated:
-get_material_flow_terms: quantifies the amount of material flow.
get_enthalpy_flow_terms: quantifies the amount of enthalpy flow.get_flow_rate: details volumetric flow rates.default_material_balance_type: defines the kind of material balance to be used.default_energy_balance_type: defines the kind of energy balance to be used.define_state_vars: involves defining state variables with units, akin to the define_metadata function in the PhysicalParameterData block.get_material_flow_basis: establishes the basis on which state variables are measured, whether in mass or molar terms.
These definitions mark the conclusion of the state block construction and thus the property package. For additional details on creating a property package, please refer to this resource.
@declare_process_block_class("OrgPhaseStateBlock", block_class=_OrganicStateBlock)
class OrgPhaseStateBlockData(StateBlockData):
"""
An example property package for Organic phzase for liquid liquid extraction
"""
def build(self):
"""
Callable method for Block construction
"""
super().build()
self._make_state_vars()
def _make_state_vars(self):
self.flow_vol = Var(
initialize=1,
domain=NonNegativeReals,
doc="Total volumetric flowrate",
units=units.L / units.hour,
)
self.conc_mass_comp = Var(
self.params.solute_set,
domain=NonNegativeReals,
initialize=1,
doc="Component mass concentrations",
units=units.g / units.L,
)
self.pressure = Var(
domain=NonNegativeReals,
initialize=1,
bounds=(1, 5),
units=units.atm,
doc="State pressure [atm]",
)
self.temperature = Var(
domain=NonNegativeReals,
initialize=300,
bounds=(273, 373),
units=units.K,
doc="State temperature [K]",
)
def material_flow_expression(self, j):
if j == "solvent":
return self.flow_vol * self.params.dens_mass
else:
return self.flow_vol * self.conc_mass_comp[j]
self.material_flow_expression = Expression(
self.component_list,
rule=material_flow_expression,
doc="Material flow terms",
)
def enthalpy_flow_expression(self):
return (
self.flow_vol
* self.params.dens_mass
* self.params.cp_mass
* (self.temperature - self.params.temperature_ref)
)
self.enthalpy_flow_expression = Expression(
rule=enthalpy_flow_expression, doc="Enthalpy flow term"
)
def get_flow_rate(self):
return self.flow_vol
def get_material_flow_terms(self, p, j):
return self.material_flow_expression[j]
def get_enthalpy_flow_terms(self, p):
return self.enthalpy_flow_expression
def default_material_balance_type(self):
return MaterialBalanceType.componentTotal
def default_energy_balance_type(self):
return EnergyBalanceType.enthalpyTotal
def define_state_vars(self):
return {
"flow_vol": self.flow_vol,
"conc_mass_comp": self.conc_mass_comp,
"temperature": self.temperature,
"pressure": self.pressure,
}
def get_material_flow_basis(self):
return MaterialFlowBasis.mass
2. Creating Aqueous Property Package#
The structure of Aqueous Property Package mirrors that of the Organic Property Package we previously developed. We’ll commence with an overview, importing the required libraries, followed by the creation of the physical property block and two state blocks. The distinctions in this package lie in the physical parameter values, and notably, the absence of the diffusion factor term, differentiating it from the prior package. The following code snippet should provide clarity on these distinctions.
# Import Python libraries
import logging
from idaes.core.util.initialization import fix_state_vars
# Import Pyomo libraries
from pyomo.environ import (
Param,
Var,
NonNegativeReals,
units,
Expression,
PositiveReals,
)
# Import IDAES cores
from idaes.core import (
declare_process_block_class,
MaterialFlowBasis,
PhysicalParameterBlock,
StateBlockData,
StateBlock,
MaterialBalanceType,
EnergyBalanceType,
Solute,
Solvent,
LiquidPhase,
)
# Some more information about this module
__author__ = "Javal Vyas"
# Set up logger
_log = logging.getLogger(__name__)
@declare_process_block_class("AqPhase")
class AqPhaseData(PhysicalParameterBlock):
"""
Property Parameter Block Class
Contains parameters and indexing sets associated with properties for
aqueous Phase
"""
def build(self):
"""
Callable method for Block construction.
"""
super().build()
self._state_block_class = AqPhaseStateBlock
# List of valid phases in property package
self.Aq = LiquidPhase()
# Component list - a list of component identifiers
self.NaCl = Solute()
self.KNO3 = Solute()
self.CaSO4 = Solute()
self.H2O = Solvent()
# Heat capacity of solvent
self.cp_mass = Param(
mutable=True,
initialize=4182,
doc="Specific heat capacity of solvent",
units=units.J / units.kg / units.K,
)
self.dens_mass = Param(
mutable=True,
initialize=997,
doc="Density of ethylene dibromide",
units=units.kg / units.m**3,
)
self.temperature_ref = Param(
within=PositiveReals,
mutable=True,
default=298.15,
doc="Reference temperature",
units=units.K,
)
@classmethod
def define_metadata(cls, obj):
obj.add_default_units(
{
"time": units.hour,
"length": units.m,
"mass": units.g,
"amount": units.mol,
"temperature": units.K,
}
)
class _AqueousStateBlock(StateBlock):
"""
This Class contains methods which should be applied to Property Blocks as a
whole, rather than individual elements of indexed Property Blocks.
"""
def fix_initialization_states(self):
fix_state_vars(self)
@declare_process_block_class("AqPhaseStateBlock", block_class=_AqueousStateBlock)
class AqPhaseStateBlockData(StateBlockData):
"""
An example property package for ideal gas properties with Gibbs energy
"""
def build(self):
"""
Callable method for Block construction
"""
super().build()
self._make_state_vars()
def _make_state_vars(self):
self.flow_vol = Var(
initialize=1,
domain=NonNegativeReals,
doc="Total volumetric flowrate",
units=units.L / units.hour,
)
self.conc_mass_comp = Var(
self.params.solute_set,
domain=NonNegativeReals,
initialize={"NaCl": 0.15, "KNO3": 0.2, "CaSO4": 0.1},
doc="Component mass concentrations",
units=units.g / units.L,
)
self.pressure = Var(
domain=NonNegativeReals,
initialize=1,
bounds=(1, 5),
units=units.atm,
doc="State pressure [atm]",
)
self.temperature = Var(
domain=NonNegativeReals,
initialize=300,
bounds=(273, 373),
units=units.K,
doc="State temperature [K]",
)
def material_flow_expression(self, j):
if j == "H2O":
return self.flow_vol * self.params.dens_mass
else:
return self.conc_mass_comp[j] * self.flow_vol
self.material_flow_expression = Expression(
self.component_list,
rule=material_flow_expression,
doc="Material flow terms",
)
def enthalpy_flow_expression(self):
return (
self.flow_vol
* self.params.dens_mass
* self.params.cp_mass
* (self.temperature - self.params.temperature_ref)
)
self.enthalpy_flow_expression = Expression(
rule=enthalpy_flow_expression, doc="Enthalpy flow term"
)
def get_flow_rate(self):
return self.flow_vol
def get_material_flow_terms(self, p, j):
return self.material_flow_expression[j]
def get_enthalpy_flow_terms(self, p):
return self.enthalpy_flow_expression
def default_material_balance_type(self):
return MaterialBalanceType.componentTotal
def default_energy_balance_type(self):
return EnergyBalanceType.enthalpyTotal
def define_state_vars(self):
return {
"flow_vol": self.flow_vol,
"conc_mass_comp": self.conc_mass_comp,
"temperature": self.temperature,
"pressure": self.pressure,
}
def get_material_flow_basis(self):
return MaterialFlowBasis.mass
3. Liquid Liquid Extractor Unit Model#
Following the creation of property packages, our next step is to develop a unit model that facilitates the mass transfer of solutes between phases. This involves importing necessary libraries, building the unit model, defining auxiliary functions, and establishing the initialization routine for the unit model.
3.1 Importing necessary libraries#
Let’s commence by importing the essential libraries from Pyomo and IDAES.
# Import Pyomo libraries
from pyomo.common.config import ConfigBlock, ConfigValue, In, Bool
from pyomo.environ import (
value,
Constraint,
check_optimal_termination,
)
# Import IDAES cores
from idaes.core import (
ControlVolume0DBlock,
declare_process_block_class,
MaterialBalanceType,
EnergyBalanceType,
MaterialFlowBasis,
MomentumBalanceType,
UnitModelBlockData,
useDefault,
)
from idaes.core.util.config import (
is_physical_parameter_block,
is_reaction_parameter_block,
)
import idaes.logger as idaeslog
from idaes.core.solvers import get_solver
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.core.util.exceptions import ConfigurationError, InitializationError
3.2 Creating the unit model#
Creating a unit model starts by creating a class called LiqExtractionData and use the declare_process_block_class decorator. The LiqExtractionData inherits the properties of UnitModelBlockData class, which allows us to create a control volume which is necessary for the unit model. After declaration of the class we proceed to define the relevant config arguments for the control volume. The config arguments includes the following properties:
material_balance_type- Indicates what type of mass balance should be constructedhas_pressure_change- Indicates whether terms for pressure change should be constructedhas_phase_equilibrium- Indicates whether terms for phase equilibrium should be constructedOrganic Property- Property parameter object used to define property calculations for the Organic phaseOrganic Property Arguments- Arguments to use for constructing Organic phase propertiesAqueous Property- Property parameter object used to define property calculations for the aqueous phaseAqueous Property Arguments- Arguments to use for constructing aqueous phase properties
As there are no pressure changes or reactions in this scenario, configuration arguments for these aspects are not included. However, additional details on configuration arguments can be found here.
@declare_process_block_class("LiqExtraction")
class LiqExtractionData(UnitModelBlockData):
"""
LiqExtraction Unit Model Class
"""
CONFIG = UnitModelBlockData.CONFIG()
CONFIG.declare(
"material_balance_type",
ConfigValue(
default=MaterialBalanceType.useDefault,
domain=In(MaterialBalanceType),
description="Material balance construction flag",
doc="""Indicates what type of mass balance should be constructed,
**default** - MaterialBalanceType.useDefault.
**Valid values:** {
**MaterialBalanceType.useDefault - refer to property package for default
balance type
**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(
"has_pressure_change",
ConfigValue(
default=False,
domain=Bool,
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(
"has_phase_equilibrium",
ConfigValue(
default=False,
domain=Bool,
description="Phase equilibrium construction flag",
doc="""Indicates whether terms for phase equilibrium should be
constructed,
**default** = False.
**Valid values:** {
**True** - include phase equilibrium terms
**False** - exclude phase equilibrium terms.}""",
),
)
CONFIG.declare(
"organic_property_package",
ConfigValue(
default=useDefault,
domain=is_physical_parameter_block,
description="Property package to use for organic phase",
doc="""Property parameter object used to define property calculations
for the organic phase,
**default** - useDefault.
**Valid values:** {
**useDefault** - use default package from parent model or flowsheet,
**PropertyParameterObject** - a PropertyParameterBlock object.}""",
),
)
CONFIG.declare(
"organic_property_package_args",
ConfigBlock(
implicit=True,
description="Arguments to use for constructing organic phase properties",
doc="""A ConfigBlock with arguments to be passed to organic phase
property block(s) and used when constructing these,
**default** - None.
**Valid values:** {
see property package for documentation.}""",
),
)
CONFIG.declare(
"aqueous_property_package",
ConfigValue(
default=useDefault,
domain=is_physical_parameter_block,
description="Property package to use for aqueous phase",
doc="""Property parameter object used to define property calculations
for the aqueous phase,
**default** - useDefault.
**Valid values:** {
**useDefault** - use default package from parent model or flowsheet,
**PropertyParameterObject** - a PropertyParameterBlock object.}""",
),
)
CONFIG.declare(
"aqueous_property_package_args",
ConfigBlock(
implicit=True,
description="Arguments to use for constructing aqueous phase properties",
doc="""A ConfigBlock with arguments to be passed to aqueous phase
property block(s) and used when constructing these,
**default** - None.
**Valid values:** {
see property package for documentation.}""",
),
)
Building the model#
After constructing the LiqExtractionData block and defining the config arguments for the control block, the next step is to write a build function that incorporates control volume and establishes constraints on the control volume to achieve the desired mass transfer. The control volume serves as a pivotal component in the unit model construction, representing the volume in which the process unfolds.
IDAES provides flexibility in choosing control volumes based on geometry, with options including 0D or 1D. In this instance, we opt for a 0D control volume, the most commonly used control volume. This choice is suitable for systems where there is a well-mixed volume of fluid or where spatial variations are deemed negligible.
The control volume encompasses parameters from (1-8), and its equations are configured to satisfy the specified config arguments. For a more in-depth understanding, users are encouraged to refer to this resource.
The build function is initiated using the super function to gain access to methods and properties of a parent or sibling class, in this case, the LiqExtractionData class. Following the super function, checks are performed on the property packages to ensure the appropriate names for the solvents, such as ‘Aq’ for the aqueous phase and ‘Org’ for the Organic phase. An error is raised if these conditions are not met. Subsequently, a check is performed to ensure there is at least one common component between the two property packages that can be transferred from one phase to another.
After these checks are completed without any exceptions raised, it is ensured that the property packages have the desired components with appropriate names. The next step is to create a control volume and assign it to a property package. Here, we initiate with the Organic phase and attach a 0D control volume to it. The control volume takes arguments about the dynamics of the block, and the property package, along with property package arguments.
The subsequent steps involve adding inlet and outlet state blocks to the control volume using the add_state_blocks function. This function takes arguments about the flow direction (defaulted to forward) and a flag for has_phase_equilibrium, which is read from the config. The control volume is now equipped with the inlet and outlet state blocks and has access to the Organic property package
Next, material balance equations are added to the control volume using the add_material_balance function, taking into account the type of material balance, has_phase_equilibrium, and the presence of has_mass_transfer. To understand this arguments further let us have a look at the material balance equation and how it is implemented in control volume.
\(\frac{\partial M_{t, p, j}}{\partial t} = F_{in, t, p, j} - F_{out, t, p, j} + N_{kinetic, t, p, j} + N_{equilibrium, t, p, j} + N_{pe, t, p, j} + N_{transfer, t, p, j} + N_{custom, t, p, j}\)
\(\frac{\partial M_{t, p, j}}{\partial t}\) - Material accumulation
\(F_{in, t, p, j}\) - Flow into the control volume
\(F_{out, t, p, j}\) - Flow out of the control volume
\(N_{kinetic, t, p, j}\) - Rate of reaction generation
\(N_{equilibrium, t, p, j}\) - Equilibrium reaction generation
\(N_{pe, t, p, j}\) - Equilibrium reaction extent
\(N_{transfer, t, p, j}\) - Mass transfer
\(N_{custom, t, p, j}\) - User defined terms in material balance
t indicates time index
p indicates phase index
j indicates component index
e indicates element index
r indicates reaction name index
Here we shall see that \(N_{transfer, t, p, j}\) is the term in the equation which is responsible for the mass transfer and the mass_transfer_term should only be equal to the amount being transferred and not include a material balance on our own. For a detailed description of the terms one should refer to the following resource
This concludes the creation of organic phase control volume. Similar procedure is done for the aqueous phase control volume with aqueous property package.
Now, the unit model has two control volumes with appropriate configurations and material, momentum and energy balances. The next step is to check the basis of the two property packages. They should both have the same flow basis, and an error is raised if this is not the case.
Following this, the add_inlet_ports and add_outlet_ports functions are used to create inlet and outlet ports. These ports are named and assigned to each control volume, resulting in labeled inlet and outlet ports for each control volume.
The subsequent steps involve writing unit-level constraints. A check if the basis is either molar or mass, and unit-level constraints are written accordingly. The first constraint pertains to the mass transfer term for the aqueous phase. The mass transfer term is equal to \(mass\_transfer\_term_{aq} = (D_{i})\frac{mass_{i}~in~aq~phase}{flowrate~of~aq~phase}\). The second constraint relates to the mass transfer term in the organic phase, which is the negative of the mass transfer term in the aqueous phase: \(mass\_transfer\_term_{org} = - mass\_transfer\_term_{aq} \)
Here \(mass\_transfer\_term_{p}\) is the term indicating the amount of material being transferred from/to the phase and \(D_{i}\) is the Distribution co-efficient for component i.
This marks the completion of the build function, and the unit model is now equipped with the necessary process constraints. The subsequent steps involve writing the initialization routine.
def build(self):
"""
Begin building model (pre-DAE transformation).
Args:
None
Returns:
None
"""
# Call UnitModel.build to setup dynamics
super().build()
# Check phase lists match assumptions
if self.config.aqueous_property_package.phase_list != ["Aq"]:
raise ConfigurationError(
f"{self.name} Liquid-Liquid Extractor model requires that the aquoues "
f"phase property package have a single phase named 'Aq'"
)
if self.config.organic_property_package.phase_list != ["Org"]:
raise ConfigurationError(
f"{self.name} Liquid-Liquid Extractor model requires that the organic "
f"phase property package have a single phase named 'Org'"
)
# Check for at least one common component in component lists
if not any(
j in self.config.aqueous_property_package.component_list
for j in self.config.organic_property_package.component_list
):
raise ConfigurationError(
f"{self.name} Liquid-Liquid Extractor model requires that the organic "
f"and aqueous phase property packages have at least one "
f"common component."
)
self.organic_phase = ControlVolume0DBlock(
dynamic=self.config.dynamic,
property_package=self.config.organic_property_package,
property_package_args=self.config.organic_property_package_args,
)
self.organic_phase.add_state_blocks(
has_phase_equilibrium=self.config.has_phase_equilibrium
)
# Separate organic and aqueous phases means that phase equilibrium will
# be handled at the unit model level, thus has_phase_equilibrium is
# False, but has_mass_transfer is True.
self.organic_phase.add_material_balances(
balance_type=self.config.material_balance_type,
has_phase_equilibrium=self.config.has_phase_equilibrium,
has_mass_transfer=True,
)
# ---------------------------------------------------------------------
self.aqueous_phase = ControlVolume0DBlock(
dynamic=self.config.dynamic,
property_package=self.config.aqueous_property_package,
property_package_args=self.config.aqueous_property_package_args,
)
self.aqueous_phase.add_state_blocks(
has_phase_equilibrium=self.config.has_phase_equilibrium
)
# Separate liquid and aqueous phases means that phase equilibrium will
# be handled at the unit model level, thus has_phase_equilibrium is
# False, but has_mass_transfer is True.
self.aqueous_phase.add_material_balances(
balance_type=self.config.material_balance_type,
# has_rate_reactions=False,
has_phase_equilibrium=self.config.has_phase_equilibrium,
has_mass_transfer=True,
)
self.aqueous_phase.add_geometry()
# ---------------------------------------------------------------------
# Check flow basis is compatible
t_init = self.flowsheet().time.first()
if (
self.aqueous_phase.properties_out[t_init].get_material_flow_basis()
!= self.organic_phase.properties_out[t_init].get_material_flow_basis()
):
raise ConfigurationError(
f"{self.name} aqueous and organic property packages must use the "
f"same material flow basis."
)
self.organic_phase.add_geometry()
# Add Ports
self.add_inlet_port(
name="organic_inlet", block=self.organic_phase, doc="Organic feed"
)
self.add_inlet_port(
name="aqueous_inlet", block=self.aqueous_phase, doc="Aqueous feed"
)
self.add_outlet_port(
name="organic_outlet", block=self.organic_phase, doc="Organic outlet"
)
self.add_outlet_port(
name="aqueous_outlet",
block=self.aqueous_phase,
doc="Aqueous outlet",
)
# ---------------------------------------------------------------------
# Add unit level constraints
# First, need the union and intersection of component lists
all_comps = (
self.aqueous_phase.properties_out.component_list
| self.organic_phase.properties_out.component_list
)
common_comps = (
self.aqueous_phase.properties_out.component_list
& self.organic_phase.properties_out.component_list
)
# Get units for unit conversion
aunits = self.config.aqueous_property_package.get_metadata().get_derived_units
lunits = self.config.organic_property_package.get_metadata().get_derived_units
flow_basis = self.aqueous_phase.properties_out[t_init].get_material_flow_basis()
if flow_basis == MaterialFlowBasis.mass:
fb = "flow_mass"
elif flow_basis == MaterialFlowBasis.molar:
fb = "flow_mole"
else:
raise ConfigurationError(
f"{self.name} Liquid-Liquid Extractor only supports mass "
f"basis for MaterialFlowBasis."
)
# Material balances
def rule_material_aq_balance(self, t, j):
if j in common_comps:
return self.aqueous_phase.mass_transfer_term[
t, "Aq", j
] == -self.organic_phase.config.property_package.diffusion_factor[j] * (
self.aqueous_phase.properties_in[t].get_material_flow_terms("Aq", j)
)
elif j in self.organic_phase.properties_out.component_list:
# No mass transfer term
# Set organic flowrate to an arbitrary small value
return self.organic_phase.mass_transfer_term[t, "Org", j] == 0 * lunits(fb)
elif j in self.aqueous_phase.properties_out.component_list:
# No mass transfer term
# Set aqueous flowrate to an arbitrary small value
return self.aqueous_phase.mass_transfer_term[t, "Aq", j] == 0 * aunits(fb)
self.material_aq_balance = Constraint(
self.flowsheet().time,
self.aqueous_phase.properties_out.component_list,
rule=rule_material_aq_balance,
doc="Unit level material balances for Aq",
)
def rule_material_liq_balance(self, t, j):
if j in common_comps:
return (
self.organic_phase.mass_transfer_term[t, "Org", j]
== -self.aqueous_phase.mass_transfer_term[t, "Aq", j]
)
else:
# No mass transfer term
# Set organic flowrate to an arbitrary small value
return self.organic_phase.mass_transfer_term[t, "Org", j] == 0 * aunits(fb)
self.material_org_balance = Constraint(
self.flowsheet().time,
self.organic_phase.properties_out.component_list,
rule=rule_material_liq_balance,
doc="Unit level material balances Org",
)
Initialization Routine#
After writing the unit model it is crucial to initialize the model properly, as non-linear models may encounter local minima or infeasibility if not initialized properly. IDAES provides us with a few initialization routines which may not work for all the models, and in such cases the developer will have to define their own initialization routines.
To create a custom initialization routine, model developers must create an initialize method as part of their model, and provide a sequence of steps intended to build up a feasible solution. Initialization routines generally make use of Pyomo’s tools for activating and deactivating constraints and often involve solving multiple sub-problems whilst building up an initial state.
For this tutorial we would use the pre-defined initialization routine of BlockTriangularizationInitializer when initializing the model in the flowsheet. This Initializer should be suitable for most models, but may struggle to initialize
tightly coupled systems of equations. This method of initialization will follow the following workflow.
Have precheck for structural singularity
Run incidence analysis on given block data and check matching.
Call Block Triangularization solver on model.
Call solve_strongly_connected_components on a given BlockData.
For more details about this initialization routine can be found here.
This marks the conclusion of creating a custom unit model, for a more detailed explanation on creating a unit model refer this resource. The next sections will deal with the diagonistics and testing of the property package and unit model.
3.3 Building a Flowsheet#
Once we have set up the unit model and its property packages, we can start building a flowsheet using them. In this tutorial, we’re focusing on a simple flowsheet with just a liquid-liquid extractor. To create the flowsheet we follow the following steps:
Import necessary libraries
Create a Pyomo model.
Inside the model, create a flowsheet block.
Assign property packages to the flowsheet block.
Add the liquid-liquid extractor to the flowsheet block.
Fix variable to make it a square problem
Run an initialization process.
Solve the flowsheet.
Following these steps, we’ve built a basic flowsheet using Pyomo. For more details, refer to the documentation.
import pyomo.environ as pyo
import idaes.core
import idaes.models.unit_models
from idaes.core.solvers import get_solver
import idaes.logger as idaeslog
from pyomo.network import Arc
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.core.initialization import InitializationStatus
from idaes.core.initialization.block_triangularization import (
BlockTriangularizationInitializer,
)
from liquid_extraction.organic_property import OrgPhase
from liquid_extraction.aqueous_property import AqPhase
from liquid_extraction.liquid_liquid_extractor import LiqExtraction
def build_model():
m = pyo.ConcreteModel()
m.fs = idaes.core.FlowsheetBlock(dynamic=False)
m.fs.org_properties = OrgPhase()
m.fs.aq_properties = AqPhase()
m.fs.lex = LiqExtraction(
dynamic=False,
has_pressure_change=False,
organic_property_package=m.fs.org_properties,
aqueous_property_package=m.fs.aq_properties,
)
return m
def fix_state_variables(m):
m.fs.lex.organic_inlet.flow_vol.fix(80 * pyo.units.L / pyo.units.hour)
m.fs.lex.organic_inlet.temperature.fix(300 * pyo.units.K)
m.fs.lex.organic_inlet.pressure.fix(1 * pyo.units.atm)
m.fs.lex.organic_inlet.conc_mass_comp[0, "NaCl"].fix(
1e-5 * pyo.units.g / pyo.units.L
)
m.fs.lex.organic_inlet.conc_mass_comp[0, "KNO3"].fix(
1e-5 * pyo.units.g / pyo.units.L
)
m.fs.lex.organic_inlet.conc_mass_comp[0, "CaSO4"].fix(
1e-5 * pyo.units.g / pyo.units.L
)
m.fs.lex.aqueous_inlet.flow_vol.fix(100 * pyo.units.L / pyo.units.hour)
m.fs.lex.aqueous_inlet.temperature.fix(300 * pyo.units.K)
m.fs.lex.aqueous_inlet.pressure.fix(1 * pyo.units.atm)
m.fs.lex.aqueous_inlet.conc_mass_comp[0, "NaCl"].fix(
0.15 * pyo.units.g / pyo.units.L
)
m.fs.lex.aqueous_inlet.conc_mass_comp[0, "KNO3"].fix(
0.2 * pyo.units.g / pyo.units.L
)
m.fs.lex.aqueous_inlet.conc_mass_comp[0, "CaSO4"].fix(
0.1 * pyo.units.g / pyo.units.L
)
return m
def initialize_model(m):
initializer = BlockTriangularizationInitializer()
initializer.initialize(m.fs.lex)
return m
def main():
m = build_model()
m = fix_state_variables(m)
m = initialize_model(m)
return m
if __name__ == main:
main()
4. Model Diagnostics using DiagnosticsToolbox#
Here, during initialization, we encounter warnings indicating that variables are being set to negative values, which is not expected behavior. These warnings suggest that there may be flaws in the model that require further investigation using the DiagnosticsToolbox from IDAES. A detailed notebook on using DiagnosticsToolbox can be found here.
To proceed with investigating these issues, we need to import the DiagnosticsToolbox. We can gain a better understanding of its functionality by running the help function on it.
from idaes.core.util import DiagnosticsToolbox
The help() function provides comprehensive information on the DiagnosticsToolbox and all its supported methods. However, it’s essential to focus on the initial steps outlined at the beginning of the docstring to get started effectively.
Here’s a breakdown of the steps to start with:
Instantiate Model:Ensure you have an instance of the model with a degrees of freedom equal to 0.Create DiagnosticsToolbox Instance:Next, instantiate a DiagnosticsToolbox object.Provide Model to DiagnosticsToolbox:Pass the model instance to the DiagnosticsToolbox.Call report_structural_issues() Function:Finally, call the report_structural_issues() function. This function will highlight any warnings in the model’s structure, such as unit inconsistencies or other issues related to variables in the caution section.
By following these steps, you can efficiently utilize the DiagnosticsToolbox to identify and address any structural issues or warnings in your model.
m = main()
dt = DiagnosticsToolbox(m)
dt.report_structural_issues()
WARNING (W1001): Setting Var
'fs.lex.aqueous_phase.properties_out[0.0].conc_mass_comp[NaCl]' to a value
`-0.1725` (float) not in domain NonNegativeReals.
See also https://pyomo.readthedocs.io/en/stable/errors.html#w1001
WARNING (W1001): Setting Var
'fs.lex.aqueous_phase.properties_out[0.0].conc_mass_comp[KNO3]' to a value
`-0.4` (float) not in domain NonNegativeReals.
See also https://pyomo.readthedocs.io/en/stable/errors.html#w1001
WARNING (W1001): Setting Var
'fs.lex.aqueous_phase.properties_out[0.0].conc_mass_comp[CaSO4]' to a value
`-0.05` (float) not in domain NonNegativeReals.
See also https://pyomo.readthedocs.io/en/stable/errors.html#w1001
====================================================================================
Model Statistics
Activated Blocks: 21 (Deactivated: 0)
Free Variables in Activated Constraints: 16 (External: 0)
Free Variables with only lower bounds: 8
Free Variables with only upper bounds: 0
Free Variables with upper and lower bounds: 0
Fixed Variables in Activated Constraints: 8 (External: 0)
Activated Equality Constraints: 16 (Deactivated: 0)
Activated Inequality Constraints: 0 (Deactivated: 0)
Activated Objectives: 0 (Deactivated: 0)
------------------------------------------------------------------------------------
0 WARNINGS
No warnings found!
------------------------------------------------------------------------------------
1 Cautions
Caution: 10 unused variables (4 fixed)
------------------------------------------------------------------------------------
Suggested next steps:
Try to initialize/solve your model and then call report_numerical_issues()
====================================================================================
Although no warnings were reported, it’s important to note that there are 3 variables fixed to 0 and 10 unused variables, out of which 4 are fixed. As indicated in the output, the next step is to solve the model. After solving, you should call the report_numerical_issues() function. This function will help identify any numerical issues that may arise during the solution process.
solver = pyo.SolverFactory("ipopt")
solver.solve(m, tee=True)
Ipopt 3.13.2:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 33
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 14
Total number of variables............................: 16
variables with only lower bounds: 8
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 16
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 4.10e+01 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 4.00e+01 4.93e+01 -1.0 4.10e-01 - 9.91e-01 2.41e-02h 1
2 0.0000000e+00 4.00e+01 2.03e+05 -1.0 4.00e-01 - 1.00e+00 2.47e-04h 1
3r 0.0000000e+00 4.00e+01 1.00e+03 1.6 0.00e+00 - 0.00e+00 3.09e-07R 4
4r 0.0000000e+00 4.00e+01 9.88e+04 1.6 3.68e+02 - 9.92e-01 2.29e-03f 1
5r 0.0000000e+00 3.60e+01 3.03e+00 1.6 4.01e+00 - 1.00e+00 1.00e+00f 1
6r 0.0000000e+00 3.69e+01 1.21e+01 -1.2 9.24e-01 - 9.69e-01 9.78e-01f 1
7r 0.0000000e+00 3.70e+01 2.11e-01 -1.9 1.00e-01 - 9.97e-01 1.00e+00f 1
8r 0.0000000e+00 3.78e+01 2.03e-02 -4.3 8.71e-01 - 9.71e-01 1.00e+00f 1
9r 0.0000000e+00 3.80e+01 2.62e-04 -6.4 1.24e-01 - 9.99e-01 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r 0.0000000e+00 3.81e+01 5.87e-09 -6.4 1.58e-01 - 1.00e+00 1.00e+00f 1
11r 0.0000000e+00 3.91e+01 1.09e-05 -9.0 9.35e-01 - 9.68e-01 1.00e+00f 1
Number of Iterations....: 11
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 5.1393961893966849e-07 5.1393961893966849e-07
Constraint violation....: 3.9105165554489545e+01 3.9105165554489545e+01
Complementarity.........: 9.0909090910996620e-10 9.0909090910996620e-10
Overall NLP error.......: 3.9105165554489545e+01 3.9105165554489545e+01
Number of objective function evaluations = 17
Number of objective gradient evaluations = 5
Number of equality constraint evaluations = 17
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 14
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 12
Total CPU secs in IPOPT (w/o function evaluations) = 0.004
Total CPU secs in NLP function evaluations = 0.000
EXIT: Converged to a point of local infeasibility. Problem may be infeasible.
WARNING: Loading a SolverResults object with a warning status into
model.name="unknown";
- termination condition: infeasible
- message from solver: Ipopt 3.13.2\x3a Converged to a locally infeasible
point. Problem may be infeasible.
{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 16, 'Number of variables': 16, 'Sense': 'unknown'}], 'Solver': [{'Status': 'warning', 'Message': 'Ipopt 3.13.2\\x3a Converged to a locally infeasible point. Problem may be infeasible.', 'Termination condition': 'infeasible', 'Id': 200, 'Error rc': 0, 'Time': 0.06552338600158691}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}
The model is probably infeasible thus indicating numerical issues with the model. We should call the report_numerical_issues() function and check what the constraints/variables causing this issue.
dt.report_numerical_issues()
====================================================================================
Model Statistics
Jacobian Condition Number: 7.955E+03
------------------------------------------------------------------------------------
2 WARNINGS
WARNING: 6 Constraints with large residuals (>1.0E-05)
WARNING: 5 Variables at or outside bounds (tol=0.0E+00)
------------------------------------------------------------------------------------
3 Cautions
Caution: 8 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04)
Caution: 5 Variables with value close to zero (tol=1.0E-08)
Caution: 3 Variables with extreme value (<1.0E-04 or >1.0E+04)
------------------------------------------------------------------------------------
Suggested next steps:
display_constraints_with_large_residuals()
display_variables_at_or_outside_bounds()
====================================================================================
In this scenario, it’s observed that the condition number of the Jacobian is high, indicating that the Jacobian is ill-conditioned. Additionally, there are 2 warnings related to constraints with large residuals and variables at or outside the bounds. The cautions mentioned in the output are also related to these warnings.
As suggested, the next steps would be to:
Call the
display_variables_at_or_outside_bounds()function to investigate variables at or outside the bounds.Call the
display_constraints_with_large_residuals()function to examine constraints with large residuals.
These steps will help identify the underlying causes of the numerical issues and constraints violations, allowing for further analysis and potential resolution.
dt.display_variables_at_or_outside_bounds()
====================================================================================
The following variable(s) have values at or outside their bounds (tol=0.0E+00):
fs.lex.organic_phase.properties_in[0.0].pressure (fixed): value=1.0 bounds=(1, 5)
fs.lex.organic_phase.properties_out[0.0].pressure (free): value=1 bounds=(1, 5)
fs.lex.aqueous_phase.properties_out[0.0].conc_mass_comp[NaCl] (free): value=0.0 bounds=(0, None)
fs.lex.aqueous_phase.properties_out[0.0].conc_mass_comp[KNO3] (free): value=0.0 bounds=(0, None)
fs.lex.aqueous_phase.properties_out[0.0].conc_mass_comp[CaSO4] (free): value=0.0 bounds=(0, None)
====================================================================================
In this scenario, there are a couple of issues to address:
The pressure variable is fixed to 1, which is its lower bound. This could potentially lead to numerical issues, although it may not affect the model significantly since there is no pressure change in the model. To mitigate this, consider adjusting the lower bound of the pressure variable to avoid having its value at or outside the bounds.
The more concerning issue is with the
conc_mass_compvariable attempting to go below 0 in the output. This suggests that there may be constraints involvingconc_mass_compin the aqueous phase causing this behavior. To investigate further, it’s recommended to call thedisplay_constraints_with_large_residuals()function. This will provide insights into whether constraints involvingconc_mass_compare contributing to the convergence issue.
dt.display_constraints_with_large_residuals()
====================================================================================
The following constraint(s) have large residuals (>1.0E-05):
fs.lex.material_aq_balance[0.0,NaCl]: 5.49716E-01
fs.lex.material_aq_balance[0.0,KNO3]: 8.94833E-01
fs.lex.material_aq_balance[0.0,CaSO4]: 5.48843E-02
fs.lex.aqueous_phase.material_balances[0.0,NaCl]: 1.67003E+01
fs.lex.aqueous_phase.material_balances[0.0,KNO3]: 3.91052E+01
fs.lex.aqueous_phase.material_balances[0.0,CaSO4]: 4.94512E+00
====================================================================================
As expected there are convergence issues with the constraints which have conc_mass_comp variable in them specifically in the aqeous phase. Now, let us investigate further by printing this constraints and checking the value of each term. Since this is an persistent issue across the components, we can focus on just one of the component to identify the issue.
m.fs.lex.aqueous_phase.material_balances[0.0, "NaCl"].pprint()
{Member of material_balances} : Material balances
Size=4, Index=fs._time*fs.aq_properties.component_list, Active=True
Key : Lower : Body : Upper : Active
(0.0, 'NaCl') : 0.0 : (fs.lex.aqueous_phase.properties_in[0.0].conc_mass_comp[NaCl]*fs.lex.aqueous_phase.properties_in[0.0].flow_vol) - (fs.lex.aqueous_phase.properties_out[0.0].conc_mass_comp[NaCl]*fs.lex.aqueous_phase.properties_out[0.0].flow_vol) + fs.lex.aqueous_phase.mass_transfer_term[0.0,Aq,NaCl] : 0.0 : True
m.fs.lex.aqueous_phase.properties_in[0.0].conc_mass_comp["NaCl"].pprint()
m.fs.lex.aqueous_phase.properties_in[0.0].flow_vol.pprint()
m.fs.lex.aqueous_phase.properties_out[0.0].conc_mass_comp["NaCl"].pprint()
m.fs.lex.aqueous_phase.properties_out[0.0].flow_vol.pprint()
m.fs.lex.aqueous_phase.mass_transfer_term[0.0, "Aq", "NaCl"].pprint()
{Member of conc_mass_comp} : Component mass concentrations
Size=3, Index=fs.aq_properties.solutes, Units=g/l
Key : Lower : Value : Upper : Fixed : Stale : Domain
NaCl : 0 : 0.15 : None : True : True : NonNegativeReals
flow_vol : Total volumetric flowrate
Size=1, Index=None, Units=l/h
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 100.0 : None : True : True : NonNegativeReals
{Member of conc_mass_comp} : Component mass concentrations
Size=3, Index=fs.aq_properties.solutes, Units=g/l
Key : Lower : Value : Upper : Fixed : Stale : Domain
NaCl : 0 : 0.0 : None : False : False : NonNegativeReals
flow_vol : Total volumetric flowrate
Size=1, Index=None, Units=l/h
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 100.0 : None : False : False : NonNegativeReals
{Member of mass_transfer_term} : Component material transfer into unit
Size=4, Index=fs._time*fs.aq_properties._phase_component_set, Units=g/h
Key : Lower : Value : Upper : Fixed : Stale : Domain
(0.0, 'Aq', 'NaCl') : None : -31.700284300098897 : None : False : False : Reals
It seems there is a discrepancy between the mass transfer term and the amount of input of NaCl. This can be inferred from the values where the input equals 15g/h and the mass_transfer_term equals -31.706g/h.
To further investigate this issue, it’s advisable to examine the material_aq_balance constraint within the unit model where the mass_transfer_term is defined. By printing out this constraint and analyzing its components, you can gain a better understanding of the discrepancy and take appropriate corrective actions.
m.fs.lex.material_aq_balance[0.0, "NaCl"].pprint()
{Member of material_aq_balance} : Unit level material balances for Aq
Size=4, Index=fs._time*fs.aq_properties.component_list, Active=True
Key : Lower : Body : Upper : Active
(0.0, 'NaCl') : 0.0 : fs.lex.aqueous_phase.mass_transfer_term[0.0,Aq,NaCl] + fs.org_properties.diffusion_factor[NaCl]*(fs.lex.aqueous_phase.properties_in[0.0].conc_mass_comp[NaCl]*fs.lex.aqueous_phase.properties_in[0.0].flow_vol) : 0.0 : True
Here the problem can be tracked down easily as there being a typing error while recording the distribution factor. The distribution factor here was wrongly written ignoring its magnitude which should have been 1e-2, but that was missed, thus adjusting the distribution factor parameter we should have this issue resolved.
m.fs.org_properties.diffusion_factor["NaCl"] = (
m.fs.org_properties.diffusion_factor["NaCl"] / 100
)
m.fs.org_properties.diffusion_factor["KNO3"] = (
m.fs.org_properties.diffusion_factor["KNO3"] / 100
)
m.fs.org_properties.diffusion_factor["CaSO4"] = (
m.fs.org_properties.diffusion_factor["CaSO4"] / 100
)
m.fs.lex.organic_phase.properties_in[0.0].pressure.setlb(0.5)
m.fs.lex.organic_phase.properties_out[0.0].pressure.setlb(0.5)
After the corrective actions, we should check if this have made any structural issues, for this we would call report_structural_issues()
dt.report_structural_issues()
====================================================================================
Model Statistics
Activated Blocks: 21 (Deactivated: 0)
Free Variables in Activated Constraints: 16 (External: 0)
Free Variables with only lower bounds: 8
Free Variables with only upper bounds: 0
Free Variables with upper and lower bounds: 0
Fixed Variables in Activated Constraints: 8 (External: 0)
Activated Equality Constraints: 16 (Deactivated: 0)
Activated Inequality Constraints: 0 (Deactivated: 0)
Activated Objectives: 0 (Deactivated: 0)
------------------------------------------------------------------------------------
0 WARNINGS
No warnings found!
------------------------------------------------------------------------------------
1 Cautions
Caution: 10 unused variables (4 fixed)
------------------------------------------------------------------------------------
Suggested next steps:
Try to initialize/solve your model and then call report_numerical_issues()
====================================================================================
Now since there are no warnings we can go ahead and solve the model and see if the results are optimal.
solver.solve(m, tee=True)
Ipopt 3.13.2:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 33
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 14
Total number of variables............................: 16
variables with only lower bounds: 8
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 16
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 5.85e+01 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 3.55e-15 8.41e+00 -1.0 5.85e+01 - 1.05e-01 1.00e+00h 1
Number of Iterations....: 1
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 3.5527136788005009e-15 3.5527136788005009e-15
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 3.5527136788005009e-15 3.5527136788005009e-15
Number of objective function evaluations = 2
Number of objective gradient evaluations = 2
Number of equality constraint evaluations = 2
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 1
Total CPU secs in IPOPT (w/o function evaluations) = 0.001
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 16, 'Number of variables': 16, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.07779264450073242}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}
This is a good sign that the model solved optimally and a solution was found.
NOTE: It is a good practice to run the model through DiagnosticsToolbox regardless of the solver termination status.
The next section we shall focus on testing the unit model.
5. Testing#
Testing is a crucial part of model development to ensure that the model works as expected, and remains reliable. Here’s an overview of why we conduct testing:
Verify Correctness: Testing ensure that the model works as expected and meets the specified requirements.Detect Bugs and Issues: Testing helps in identifying bugs, errors, or unexpected behaviors in the code or model, allowing for timely fixes.Ensure Reliability: Testing improves the reliability and robustness of the software, reducing the risk of failures when the user uses it.Support Changes: Tests provide confidence when making changes or adding new features, ensuring that existing functionalities are not affected and work as they should.
There are typically 3 types of tests:
Unit tests: Test runs quickly (under 2 seconds) and has no network/system dependencies. Uses only libraries installed by default with the softwareComponent test: Test may run more slowly (under 10 seconds, or so), e.g. it may run a solver or create a bunch of files. Like unit tests, it still shouldn’t depend on special libraries or dependencies.Integration test: Test may take a long time to run, and may have complex dependencies.
The expectation is that unit tests should be run by developers rather frequently, component tests should be run by the continuous integration system before running code, and integration tests are run across the codebase regularly, but infrequently (e.g. daily).
As a developer, testing is a crucial aspect of ensuring the reliability and correctness of the unit model. The testing process involves both Unit tests and Component tests, and pytest is used as the testing framework. A typical test is marked with @pytest.mark.level, where the level indicates the depth or specificity of the testing. This is written in a file usually named as test_*.py or *_test.py. The test files have functions written in them with the appropriate level of test being conducted.
For more detailed information on testing methodologies and procedures, developers are encouraged to refer to this resource. The resource provides comprehensive guidance on the testing process and ensures that the unit model meets the required standards and functionality.
5.1 Property package#
Unit Tests#
When writing tests for the Aqueous property phase package, it’s essential to focus on key aspects to ensure the correctness and robustness of the implementation. Here are the areas to cover in the unit tests:
Number of Config Dictionaries: Verify that the property phase package has the expected number of configuration dictionaries.
State Block Class Name: Confirm that the correct state block class is associated with the Aqueous property phase package.
Number of Phases: Check that the Aqueous property phase package defines the expected number of phases.
Components in the Phase and Physical Parameter Values: Test that the components present in the Aqueous phase match the anticipated list. Additionally, validate that the physical parameter values (such as density, viscosity, etc.) are correctly defined.
import pytest
from pyomo.environ import ConcreteModel, Param, value, Var
from pyomo.util.check_units import assert_units_consistent
from idaes.core import MaterialBalanceType, EnergyBalanceType
from liquid_extraction.organic_property import OrgPhase
from liquid_extraction.aqueous_property import AqPhase
from liquid_extraction.liquid_liquid_extractor import LiqExtraction
from idaes.core.solvers import get_solver
solver = get_solver()
class TestParamBlock(object):
@pytest.fixture(scope="class")
def model(self):
model = ConcreteModel()
model.params = AqPhase()
return model
@pytest.mark.unit
def test_config(self, model):
assert len(model.params.config) == 1
@pytest.mark.unit
def test_build(self, model):
assert len(model.params.phase_list) == 1
for i in model.params.phase_list:
assert i == "Aq"
assert len(model.params.component_list) == 4
for i in model.params.component_list:
assert i in ["H2O", "NaCl", "KNO3", "CaSO4"]
assert isinstance(model.params.cp_mass, Param)
assert value(model.params.cp_mass) == 4182
assert isinstance(model.params.dens_mass, Param)
assert value(model.params.dens_mass) == 997
assert isinstance(model.params.temperature_ref, Param)
assert value(model.params.temperature_ref) == 298.15
The next set of unit tests focuses on testing the build function in the state block. Here are the key aspects to cover in these tests:
Existence and Initialized Values of State Variables: Verify that the state variables are correctly defined and initialized within the state block. This ensures that the state block is properly constructed and ready for initialization.
Initialization Function Test: Check that state variables are not fixed before initialization and are released after initialization. This test ensures that the initialization process occurs as expected and that the state variables are appropriately managed throughout.
These unit tests provide comprehensive coverage for validating the functionality and behavior of the state block in the Aqueous property phase package. Similar tests can be written for the organic property package to ensure consistency and reliability across both packages.
class TestStateBlock(object):
@pytest.fixture(scope="class")
def model(self):
model = ConcreteModel()
model.params = AqPhase()
model.props = model.params.build_state_block([1])
return model
@pytest.mark.unit
def test_build(self, model):
assert isinstance(model.props[1].flow_vol, Var)
assert value(model.props[1].flow_vol) == 1
assert isinstance(model.props[1].temperature, Var)
assert value(model.props[1].temperature) == 300
assert isinstance(model.props[1].conc_mass_comp, Var)
assert len(model.props[1].conc_mass_comp) == 3
@pytest.mark.unit
def test_initialize(self, model):
assert not model.props[1].flow_vol.fixed
assert not model.props[1].temperature.fixed
assert not model.props[1].pressure.fixed
for i in model.props[1].conc_mass_comp:
assert not model.props[1].conc_mass_comp[i].fixed
model.props.initialize(hold_state=False, outlvl=1)
assert not model.props[1].flow_vol.fixed
assert not model.props[1].temperature.fixed
assert not model.props[1].pressure.fixed
for i in model.props[1].conc_mass_comp:
assert not model.props[1].conc_mass_comp[i].fixed
Component Tests#
In the component test, we aim to ensure unit consistency across the entire property package. Unlike unit tests that focus on individual functions, component tests assess the coherence and consistency of the entire package. Here’s what the component test will entail:
Unit Consistency Check: Verify that all units used within the property package are consistent throughout. This involves checking that all parameters, variables, and equations within the package adhere to the same unit system, ensuring compatibility.
By conducting a comprehensive component test, we can ensure that the property package functions as a cohesive unit, maintaining consistency and reliability across its entirety. This concludes our tests on the property package. Next we shall test the unit model.
@pytest.mark.component
def check_units(model):
model = ConcreteModel()
model.params = AqPhase()
assert_units_consistent(model)
5.2 Unit Model#
Unit tests#
Unit tests for the unit model encompass verifying the configuration arguments and the build function, similar to the approach taken for the property package. When testing the config arguments, we ensure that the correct number of arguments is provided and then match each argument with the expected one. This ensures that the unit model is properly configured and ready to operate as intended.
import pytest
import idaes.core
import idaes.models.unit_models
from idaes.core.solvers import get_solver
import idaes.logger as idaeslog
from pyomo.environ import value, check_optimal_termination, units
from pyomo.util.check_units import assert_units_consistent
from idaes.core.util.model_statistics import (
number_variables,
number_total_constraints,
)
from idaes.core.solvers import get_solver
from idaes.core.initialization import (
SingleControlVolumeUnitInitializer,
)
solver = get_solver()
@pytest.mark.unit
def test_config():
m = ConcreteModel()
m.fs = idaes.core.FlowsheetBlock(dynamic=False)
m.fs.org_properties = OrgPhase()
m.fs.aq_properties = AqPhase()
m.fs.unit = LiqExtraction(
dynamic=False,
has_pressure_change=False,
organic_property_package=m.fs.org_properties,
aqueous_property_package=m.fs.aq_properties,
)
# Check unit config arguments
assert len(m.fs.unit.config) == 9
# Check for config arguments
assert m.fs.unit.config.material_balance_type == MaterialBalanceType.useDefault
assert not m.fs.unit.config.has_pressure_change
assert not m.fs.unit.config.has_phase_equilibrium
assert m.fs.unit.config.organic_property_package is m.fs.org_properties
assert m.fs.unit.config.aqueous_property_package is m.fs.aq_properties
# Check for unit initializer
assert m.fs.unit.default_initializer is SingleControlVolumeUnitInitializer
In testing the build function, we verify whether the number of variables aligns with the intended values and also check for the existence of desired constraints within the unit model. This ensures that the unit model is constructed accurately and includes all the necessary variables and constraints required for its proper functioning.
class TestBuild(object):
@pytest.fixture(scope="class")
def model(self):
m = ConcreteModel()
m.fs = idaes.core.FlowsheetBlock(dynamic=False)
m.fs.org_properties = OrgPhase()
m.fs.aq_properties = AqPhase()
m.fs.unit = LiqExtraction(
dynamic=False,
has_pressure_change=False,
organic_property_package=m.fs.org_properties,
aqueous_property_package=m.fs.aq_properties,
)
m.fs.unit.organic_inlet.flow_vol.fix(80 * units.l / units.h)
m.fs.unit.organic_inlet.temperature.fix(300 * units.K)
m.fs.unit.organic_inlet.pressure.fix(1 * units.atm)
m.fs.unit.organic_inlet.conc_mass_comp[0, "NaCl"].fix(0 * units.g / units.l)
m.fs.unit.organic_inlet.conc_mass_comp[0, "KNO3"].fix(0 * units.g / units.l)
m.fs.unit.organic_inlet.conc_mass_comp[0, "CaSO4"].fix(0 * units.g / units.l)
m.fs.unit.aqueous_inlet.flow_vol.fix(10 * units.l / units.h)
m.fs.unit.aqueous_inlet.temperature.fix(300 * units.K)
m.fs.unit.aqueous_inlet.pressure.fix(1 * units.atm)
m.fs.unit.aqueous_inlet.conc_mass_comp[0, "NaCl"].fix(0.15 * units.g / units.l)
m.fs.unit.aqueous_inlet.conc_mass_comp[0, "KNO3"].fix(0.2 * units.g / units.l)
m.fs.unit.aqueous_inlet.conc_mass_comp[0, "CaSO4"].fix(0.1 * units.g / units.l)
return m
@pytest.mark.build
@pytest.mark.unit
def test_build(self, model):
assert hasattr(model.fs.unit, "aqueous_inlet")
assert len(model.fs.unit.aqueous_inlet.vars) == 4
assert hasattr(model.fs.unit.aqueous_inlet, "flow_vol")
assert hasattr(model.fs.unit.aqueous_inlet, "conc_mass_comp")
assert hasattr(model.fs.unit.aqueous_inlet, "temperature")
assert hasattr(model.fs.unit.aqueous_inlet, "pressure")
assert hasattr(model.fs.unit, "organic_inlet")
assert len(model.fs.unit.organic_inlet.vars) == 4
assert hasattr(model.fs.unit.organic_inlet, "flow_vol")
assert hasattr(model.fs.unit.organic_inlet, "conc_mass_comp")
assert hasattr(model.fs.unit.organic_inlet, "temperature")
assert hasattr(model.fs.unit.organic_inlet, "pressure")
assert hasattr(model.fs.unit, "aqueous_outlet")
assert len(model.fs.unit.aqueous_outlet.vars) == 4
assert hasattr(model.fs.unit.aqueous_outlet, "flow_vol")
assert hasattr(model.fs.unit.aqueous_outlet, "conc_mass_comp")
assert hasattr(model.fs.unit.aqueous_outlet, "temperature")
assert hasattr(model.fs.unit.aqueous_outlet, "pressure")
assert hasattr(model.fs.unit, "organic_outlet")
assert len(model.fs.unit.organic_outlet.vars) == 4
assert hasattr(model.fs.unit.organic_outlet, "flow_vol")
assert hasattr(model.fs.unit.organic_outlet, "conc_mass_comp")
assert hasattr(model.fs.unit.organic_outlet, "temperature")
assert hasattr(model.fs.unit.organic_outlet, "pressure")
assert hasattr(model.fs.unit, "material_aq_balance")
assert hasattr(model.fs.unit, "material_org_balance")
assert number_variables(model) == 34
assert number_total_constraints(model) == 16
Component tests#
During the component tests, we evaluate the performance of the unit model when integrated with the property package. This evaluation process typically involves several steps:
Unit Consistency Check: Verify that the unit model maintains consistency in its units throughout the model. This ensures that all variables and constraints within the model adhere to the same unit system, guaranteeing compatibility.
Termination Condition Verification: This involves checking whether the model terminates optimally with the given inlet conditions.
Variable Value Assessment: Check the values of outlet variables against the expected values. To account for the numerical tolerance of the solvers, the values are compared using the approx function with a relative tolerance.
Input Variable Stability Test: Verify that input variables, which should remain fixed during model operation, are not inadvertently unfixed or altered.
Structural Issues: Verify that there are no structural issues with the model.
By performing these checks, we conclude the testing for the unit model.
class TestFlowsheet:
@pytest.fixture
def model(self):
m = ConcreteModel()
m.fs = idaes.core.FlowsheetBlock(dynamic=False)
m.fs.org_properties = OrgPhase()
m.fs.aq_properties = AqPhase()
m.fs.unit = LiqExtraction(
dynamic=False,
has_pressure_change=False,
organic_property_package=m.fs.org_properties,
aqueous_property_package=m.fs.aq_properties,
)
m.fs.org_properties.diffusion_factor["NaCl"] = (
m.fs.org_properties.diffusion_factor["NaCl"] / 100
)
m.fs.org_properties.diffusion_factor["KNO3"] = (
m.fs.org_properties.diffusion_factor["KNO3"] / 100
)
m.fs.org_properties.diffusion_factor["CaSO4"] = (
m.fs.org_properties.diffusion_factor["CaSO4"] / 100
)
m.fs.unit.organic_inlet.flow_vol.fix(80 * units.ml / units.min)
m.fs.unit.organic_inlet.temperature.fix(300 * units.K)
m.fs.unit.organic_inlet.pressure.fix(1 * units.atm)
m.fs.unit.organic_inlet.conc_mass_comp[0, "NaCl"].fix(0 * units.g / units.kg)
m.fs.unit.organic_inlet.conc_mass_comp[0, "KNO3"].fix(0 * units.g / units.kg)
m.fs.unit.organic_inlet.conc_mass_comp[0, "CaSO4"].fix(0 * units.g / units.kg)
m.fs.unit.aqueous_inlet.flow_vol.fix(10 * units.ml / units.min)
m.fs.unit.aqueous_inlet.temperature.fix(300 * units.K)
m.fs.unit.aqueous_inlet.pressure.fix(1 * units.atm)
m.fs.unit.aqueous_inlet.conc_mass_comp[0, "NaCl"].fix(0.15 * units.g / units.kg)
m.fs.unit.aqueous_inlet.conc_mass_comp[0, "KNO3"].fix(0.2 * units.g / units.kg)
m.fs.unit.aqueous_inlet.conc_mass_comp[0, "CaSO4"].fix(0.1 * units.g / units.kg)
return m
@pytest.mark.component
def test_unit_model(self, model):
assert_units_consistent(model)
solver = get_solver()
results = solver.solve(model, tee=False)
# Check for optimal termination
assert check_optimal_termination(results)
# Checking for outlet flows
assert value(model.fs.unit.organic_outlet.flow_vol[0]) == pytest.approx(
80.0, rel=1e-5
)
assert value(model.fs.unit.aqueous_outlet.flow_vol[0]) == pytest.approx(
10.0, rel=1e-5
)
# Checking for outlet mass_comp
assert value(
model.fs.unit.organic_outlet.conc_mass_comp[0, "CaSO4"]
) == pytest.approx(0.000187499, rel=1e-5)
assert value(
model.fs.unit.organic_outlet.conc_mass_comp[0, "KNO3"]
) == pytest.approx(0.000749999, rel=1e-5)
assert value(
model.fs.unit.organic_outlet.conc_mass_comp[0, "NaCl"]
) == pytest.approx(0.000403124, rel=1e-5)
assert value(
model.fs.unit.aqueous_outlet.conc_mass_comp[0, "CaSO4"]
) == pytest.approx(0.0985, rel=1e-5)
assert value(
model.fs.unit.aqueous_outlet.conc_mass_comp[0, "KNO3"]
) == pytest.approx(0.194, rel=1e-5)
assert value(
model.fs.unit.aqueous_outlet.conc_mass_comp[0, "NaCl"]
) == pytest.approx(0.146775, rel=1e-5)
# Checking for outlet temperature
assert value(model.fs.unit.organic_outlet.temperature[0]) == pytest.approx(
300, rel=1e-5
)
assert value(model.fs.unit.aqueous_outlet.temperature[0]) == pytest.approx(
300, rel=1e-5
)
# Checking for outlet pressure
assert value(model.fs.unit.organic_outlet.pressure[0]) == pytest.approx(
1, rel=1e-5
)
assert value(model.fs.unit.aqueous_outlet.pressure[0]) == pytest.approx(
1, rel=1e-5
)
# Fixed state variables
assert model.fs.unit.organic_inlet.flow_vol[0].fixed
assert model.fs.unit.organic_inlet.conc_mass_comp[0, "NaCl"].fixed
assert model.fs.unit.organic_inlet.conc_mass_comp[0, "KNO3"].fixed
assert model.fs.unit.organic_inlet.conc_mass_comp[0, "CaSO4"].fixed
assert model.fs.unit.organic_inlet.temperature[0].fixed
assert model.fs.unit.organic_inlet.pressure[0].fixed
assert model.fs.unit.aqueous_inlet.flow_vol[0].fixed
assert model.fs.unit.aqueous_inlet.conc_mass_comp[0, "NaCl"].fixed
assert model.fs.unit.aqueous_inlet.conc_mass_comp[0, "KNO3"].fixed
assert model.fs.unit.aqueous_inlet.conc_mass_comp[0, "CaSO4"].fixed
assert model.fs.unit.aqueous_inlet.temperature[0].fixed
assert model.fs.unit.aqueous_inlet.pressure[0].fixed
@pytest.mark.component
def test_structural_issues(self, model):
dt = DiagnosticsToolbox(model)
dt.assert_no_structural_warnings()
In this tutorial, we have covered the comprehensive process of creating a custom unit model from scratch. Let’s recap the key steps we have undertaken:
Developing property package
Constructing the unit model
Creating a Flowsheet
Debugging the model using DiagnosticsToolbox
Writing tests for the unit model
By following the aforementioned procedure, one can create their own custom unit model. This would conclude the tutorial on creating custom unit model.