Custom Physical Property Packages

Jupyter Notebook source file

Introduction to Property Packages in IDAES¶

Calculation of thermophysical, transport and reaction properties form a key part of any process model, and it is important that these calculations are both accurate and tractable in order for the overall problem to be solved correctly. One of the features of the IDAES Integrated Platform is the ability for modelers to create their own property “packages” to calculate these properties, allowing them to customize the level of complexity and rigor to suit each application. This tutorial will introduce you to the basics of creating property packages for calculating thermophysical and transport properties within the IDAES Core Modeling Framework.

What is a Property?¶

Within the context of the IDAES Core Modeling Framework, a property is considered to be any variable describing the state of a material at a given point in space and time (both intensive and extensive), and anything which can be calculated directly from the state variables. Some common examples of properties include:

• material flowrates
• material compositions (mass and/or mole fractions, concentrations, etc.)
• temperature
• pressure
• density and specific volume
• specific heat capacity, enthalpy and entropy
• phase equilibria
• transport properties such as viscosity and thermal conductivity
• rates of reaction and chemical equilibria

The definition and calculation of all of these is defined via “property packages”, which contain all the variables and constraints associated with calculating these properties.

Note: It is important to emphasize that extensive flow information is considered part of the material state within the IDAES Core Modeling Framework. This is done to allow the framework to adapt to how the material flow (e.g. by volumetric flow and concentrations, or molar flow and mole fractions) is defined and automatically incorporate this into the material, energy and momentum balances.

How are Properties Used?¶

Before we begin discussing how to build a property package, it is important to understand where and how properties are used within the greater modeling framework. Hopefully by now you are familiar with the IDAES modeling hierarchy or flowsheets, unit models and state blocks.

A process flowsheet consists of a network of unit models connected together by Arcs (i.e. streams). Each unit model contains a set of material, energy and momentum balance equations as well as some performance equations that describe how the unit behaves. Each of these equations in turn depends on the properties of the material at some point in the unit operations (generally an inlet or an outlet), and it is here that the unit model calls upon a property package to obtain these properties.

For each point in space and time included in the model, a StateBlock is created which contains a full description of the material state (i.e. the full set of necessary state variables), along with calculations for all the properties required at that point. The unit model (and flowsheet) can then refer to the property in this StateBlock wherever it is required in a calculation.

What Properties do I Need?¶

An important aspect of the IDAES Core Modeling Framework is that a modeler only needs to provide calculations for those properties that they will use within their process. Put another way, modelers do not need to include calculations for a property that they are not going to use in their model – this allows modelers to avoid introducing unnecessary complexity into their models to calculate a property they do not actually need. When combined with flexibility elsewhere in the modeling framework to control which equations are written in the unit models, this can even allow users to avoid calculating properties that would normally be considered mandatory – for example, a property package for a conceptual design flowsheet which does not include energy or momentum balances would not need to define specific enthalpy or even temperature and pressure as these will not be required by the unit models.

This then raises the question of how do you know what properties you will need, especially if you are using models from a library you did not write yourself. To answer this, you should refer to the model documentation, and you can also use the IDAES Properties Interrogator tool to analyze your flowsheet and determine what properties are required.

Thermophysical Properties and Reaction Properties¶

Within the IDAES Core Modeling Framework, properties are divided into two classifications; thermophysical properties and reaction properties. Reaction properties are those properties related to chemical reactions (both equilibrium and rate-based, but not phase equilibrium) that occur within the system , whilst thermophysical properties include those properties related to thermodynamic relationships (including phase equilibrium) and transport properties. The reason for this separation is that thermophysical properties are required by all unit operations in a process (and need to be consistent with each other), whilst reaction properties are generally only required in specific unit operations identified as “reactors” (and each reactor may have a different set of chemical reactions occurring in it). Thus, reaction properties are separated from the thermophysical property calculations to allow for modular implementation in only specific reactor units. This tutorial only deals with thermophysical properties, and reaction properties will be dealt with in a later tutorial.

What is a Property “Package”?¶

Generally, properties (both thermophysical and reaction) are calculated using correlations that depend on some set of parameters (be they physical constants or empirical parameters). These parameters are constant across all instances of a property calculation in a flowsheet (i.e. each StateBlock uses the same parameters), it makes sense to store these parameters in a single, central location that all StateBlocks can refer to as necessary. Thus, the IDAES modeling framework has “Parameter Blocks” which are attached to the flowsheet to contain all the global parameters associated with a given set of property calculations.

Thus, the calculations of thermophysical properties within the IDAES modeling framework is achieved using a “package” of three related modeling components (or classes); the Physical Parameter Block, the State Block and the State Block Data classes. Each of these will be discussed further in the next section as we develop an example property package.

At a deeper level, the calculation for many thermophysical properties is a self-contained correlation that is more or less independent of the other properties around it. Thus, each set of thermophysical property calculations is a package of user-chosen sub-models for each property of interest to the user.

Steps in Creating a Property Package¶

Creating a new property package can be broken down into the following steps, which will be demonstrated in the next part of this tutorial.

1. Defining the units of measurement for the property package.
2. Defining the properties supported by the property package and the associated metadata.
3. Defining the phases and components of interest.
4. Defining the necessary parameters required to calculate the properties of interest.
5. Declaring the state variables to be used for the property package.
6. Creating variables and constraints to describe the properties of interest.
7. Creating an initialization routine for the property package.
8. Defining interface methods used to couple the property package with unit models.

Tutorial Example¶

For this tutorial, we will be building a simple property package that could be used to model a process for the hydrodealkylation of toluene to form benzene. This process involves five key chemical species:

• toluene
• benzene
• hydrogen
• methane
• diphenyl

The process also involves two reactions; a primary reaction of toluene and hydrogen reacting to form benzene and methane and a secondary reaction of benzene reacting to form diphenyl and hydrogen. These reactions occur within the vapor phase at high temperatures, however vapor-liquid equilibrium is used to separate the benzene product from the toluene and diphenyl. For this tutorial however, we will focus only on the vapor phase as implementing the equations for vapor-liquid equilibrium requires a separate tutorial of its own. Additionally, as this tutorial is only looking at writing the thermophysical property package, we do not need to concern our selves with the reaction properties (this will be the subject of another tutorial).

For this tutorial we will use total molar flowrate, component mole fractions, temperature and pressure as the state variables, and implement methods to calculate the following properties using the assumption of ideal gas behavior;

• molecular weight of each component
• density of the mixture using the ideal gas equation
• specific enthalpy of the mixture using the equations from "The Properties of Gases and Liquids 4th edition", by Reid, Prausnitz and Polling (1987).

This is a rather limited set of properties, but it is sufficient to demonstrate the different ways that properties can be implemented and to solve a simple example. Users can implement additional properties as required using these as examples.

A Note on this Tutorial¶

The build methods in the property package classes are generally written as a single, long method. However, to break the code into manageable pieces for discussion, in this tutorial we will create a number of smaller sub-methods that will then be called as part of the build method. This is done entirely for presentation purposes, and model developers should not feel compelled to write their models this way.

An example of how the example in this tutorial would be written without sub-methods can be found in the same folder with the name thermophysical_property_example.py.

Components of a Property Package¶

As mentioned previously, thermophysical property packages consist of three parts, which are written as Python classes. These components are:

• The Physical Parameter Block class, which contains all the global parameters associated with the property package,
• The State Block Data class, which contains the instructions on how to calculate all the properties at a given state, and,
• The State Block class, which is used to construct indexed sets of State Block Data objects and contains methods for acting on all multiple State Block Data objects at once (such as initialization).

It is not necessary to understand the reason for the distinction between the State Block and State Block Data classes. Suffice to say that this is due to the need to replicate the underlying component structure of Pyomo, and that the State Block represents the indexed Block representing a set of states across a given indexing set (most often time), and the State Block Data represents the individual elements of the indexed Block.

Importing Libraries¶

Before we begin writing the actual classes however, we need to import all the necessary components from the Pyomo and IDAES modeling libraries. To begin with, we are going to need a number of components from the Pyomo modeling environment to construct the variables, constraints and parameters that will make up the property package, and we will also make use of the Pyomo units of measurement tools to define the units of our properties. We will also make use of a number of components and supporting methods from the IDAES modeling framework and libraries.

Rather than describe the purpose of all of these here, we shall just import all of them here and discuss their use as they arise in the example.

In [1]:

# Import components from Pyomo
from pyomo.environ import Constraint, Expression, Reference, Param, units as pyunits, Var

# Import IDAES cores
from idaes.core import (declare_process_block_class,
                        MaterialFlowBasis,
                        PhysicalParameterBlock,
                        StateBlockData,
                        StateBlock,
                        MaterialBalanceType,
                        EnergyBalanceType,
                        Component,
                        VaporPhase)
from idaes.core.solvers import get_solver
from idaes.core.util.initialization import (fix_state_vars,
                                            revert_state_vars,
                                            solve_indexed_blocks)
from idaes.core.util.model_statistics import degrees_of_freedom, \
                                             number_unfixed_variables
from idaes.core.util.constants import Constants as const
import idaes.logger as idaeslog

The Physical Parameter Block¶

We will begin by constructing the Physical Parameter Block for our property example. This serves as the central point of reference for all aspects of the property package, and needs to define a number of 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

Step 1: Define Units of Measurement¶

The first step is to define the units of measurement for the property package, which will in turn be inherited by any unit model using this property package. Units of measurement for the property package are defined by setting units for the 7 base measurement quantities; time, length, mass, amount, temperature, current and luminous intensity (as current and luminous intensity are generally of lesser importance in process systems engineering, specifying units for these base quantities is optional). Within IDAES, units are specified using Pyomo's units of measurement features, which can be imported from pyomo.environ. For this example, the units of measurement features were given the name pyunits for clarity.

The units of measurement for all other quantities in the model can then be derived from these base quantities; for example the units of energy are mass*length^2/time^2. The framework expects all quantities in the property package to use these base units – the Pyomo units of measurement conversion tools can be used if conversion between different sets of units are required.

In order to set the base units, we need to create a dictionary which has each of the base quantities as a key, and provide a Pyomo recognized unit as the value as shown below.

In [2]:

units_metadata = {'time': pyunits.s,
                  'length': pyunits.m,
                  'mass': pyunits.kg,
                  'amount': pyunits.mol,
                  'temperature': pyunits.K}

Step 2: Define Supported Properties¶

The next step is to provide some metadata defining what properties are supported by the property package (including state variables). The first purpose of this metadata is to record a summary of what properties are supported to help user identify whether a given property package is suitable for their needs. The second purpose of the metadata is to allow us to simplify our property calculations by only construction those properties that are actually required by a given unit operation – a property package needs to support all the properties required by a process flowsheet, but not all of those properties are required in every unit operation. Thus, the IDAES modeling framework supports a “build-on-demand” approach for properties, such that only those properties that are required are constructed at any given point.

This is achieved through the use of the properties metadata, where for each property supported by the property package the user needs to define a method argument. This argument can take one of two forms:

• None indicates that the specified property is always constructed by the property package (as part of the State Block Data build method), and that if it can not be found then something has gone wrong. These types of properties would typically be the state variables (which must always be constructed), and any property that is expected to be used in most unit operations (such as specific enthalpy that will likely be used in the energy balances).
• for those properties that should only be constructed when required, the method argument should be the name of a method in the State Block Data class (as a string) which contains the instructions on building the necessary variables, expressions and constraints for calculating this property. If a unit model calls for this property and it has not been constructed yet, this method will be called in order to construct it.

In this tutorial, we will demonstrate both approaches (users can freely mix both approaches to suit their needs). For this example, the state variables, molecular weight and density calculations will be constructed for all states, whilst the molar enthalpy calculations will be built-on-demand. The cell, below shows the property metadata dictionary for this example - note that the method for enth_mol is a string which refers to a method we will construct later).

In [3]:

properties_metadata = {'flow_mol': {'method': None},
                       'mole_frac_comp': {'method': None},
                       'temperature': {'method': None},
                       'pressure': {'method': None},
                       'mw_comp': {'method': None},
                       'dens_mol': {'method': None},
                       'enth_mol': {'method': '_enth_mol'}}

Step 3: Define Component and Phase Lists¶

The next step in writing the Physical Parameter Block class is to define the phases and components present in the mixture. These are defined using Phase and Component objects which are imported from idaes.core. As Phase and Component objects are added to the proeprty package, the phase_list and component_list Sets required by the modeling framework are automatically populated. Even for systems where there is only a single phase or component, it is necessary to define phase and component objects as these are used to construct the necessary indexing sets used when building unit models.

For this example, we have 5 components of interest; benzene, toluene, hydrogen, methane and diphenyl. We define these in the property package by adding a generic Component object to the Physical Parameter Block; for example self.benzene = Component(). For more complex systems, IDAES supports a number of different component types and components can be assigned a number of arguments at constructions but these will not be discussed here.

Similarly, we also need to define the phases of interest in our system; in this case we have a single vapor phase. IDAES supports a number of different phase types, and the choice of phase type is used to determine how a phase behaves when a phase separation occurs. For example, in a flash separator, all vapor phases exit via the vapor outlet whilst all other phases exit via the liquid outlet. For this example, we add a single VaporPhase object to the Physical Parameter Block wit the name Vap as shown below.

In [4]:

def define_components_and_phases(self):
    # Define Component objects for all species
    self.benzene = Component()
    self.toluene = Component()
    self.methane = Component()
    self.hydrogen = Component()
    self.diphenyl = Component()

    # Define Phase objects for all phases
    self.Vap = VaporPhase()

Step 4: Define Parameters¶

Once the phases and components of interest have been defined, the next step is to begin defining the parameters necessary for the property calculations.

Param or Var:

The most obvious way to declare a "parameter" in a model would appear to be to use the Pyomo Param object. However, modelers should be aware the Param objects are never seen by the solver (they are converted to fixed floating point numbers by the solver writer). This means that it is not possible to use a solver to find the value for a parameter – i.e., it is not possible to use Param objects in a parameter estimation problem.

Instead, modelers should use fixed Var objects for any parameter that may need to be estimated at some point. Within IDAES, this means that most “parameters” are in fact declared as Var objects, with Param objects used only for parameters with well-known values (for example critical pressures and temperatures or molecular weights).

For this example, the first parameters we need to define are the reference state for our property calculations along with the molecular weights of each of the components of interest. These are fixed parameters that should not be estimated by parameter estimation, so we create Pyomo Param objects to represent each of these, as shown below. When we declare a Param, we also need to define a default value and the units of measurement for each parameter. Note that the units of measurement for these parameters does not necessarily need to match those defined in the properties metadata, but if they are not consistent then a unit conversion will be required at some point when calculating property values.

Within the IDAES Modeling Libraries, we generally define parameters and property calculations in the units of the source material, and then convert the calculated value into the desired units of measurement for the property package. This is done to minimize the chance of mistakes when entering parameter values as they can be taken directly from the source material. For example, Antoine coefficients are generally defined in units of bar and Kelvins and the calculation for saturation pressure is performed in these units with only the final value for saturation pressure converted to the desired units for pressure.

In [5]:

def define_basic_parameters(self):
    # Thermodynamic reference state
    self.pressure_ref = Param(mutable=True,
                              default=101325,
                              units=pyunits.Pa,
                              doc='Reference pressure')
    self.temperature_ref = Param(mutable=True,
                                default=298.15,
                                units=pyunits.K,
                                doc='Reference temperature')

    self.mw_comp = Param(self.component_list,
                        mutable=False,
                        initialize={'benzene': 78.1136E-3,
                                    'toluene': 92.1405E-3,
                                    'hydrogen': 2.016e-3,
                                    'methane': 16.043e-3,
                                    'diphenyl': 154.212e-4},
                        units=pyunits.kg/pyunits.mol,
                        doc="Molecular weight")

For this example, we also need to define the parameter associated with calculating the specific enthalpy of each component. As mentioned before, we will use the correlation proposed in “The Properties of Gases and Liquids, 4th Edition” by Reid, Prausnitz and Polling (1987), which has the form:

\begin{equation*} h_j – h_{j, ref}= A_j \times (T-T_{ref}) + \frac{B_j}{2}\times (T^2-T_{ref}^2) + \frac{C_j}{3}\times (T^3-T_{ref}^3) + \frac{D_j}{4}\times (T^4-T_{ref}^4) \end{equation*}

where $h_{j, ref}$ is the standard heat of formation of component $j$ in the vapor phase, and $A_j$, $B_j$, $C_j$, and $D_j$ are component-specific parameters in the correlation. At first glance, one might ask if we could declare a single object indexed by the list [“A”, “B”, “C”, “D”] and component to represent all the parameters as a single object; however it must be noted that the parameters $A$, $B$, $C$, and $D$ all have different units. Thus, we need to declare separate objects for each of $A$, $B$, $C$, and $D$ (along with $h_{ref}$) which are indexed by component so that we can assign the correct units to each.

However, these parameters are mostly empirical and are values that we may wish to estimate at some point, thus we will declare these as Pyomo Var objects rather than Param objects, which also means that we must fix the value of these parameters when we construct the property package. This is shown in the code below – note that each parameters (Var object is fixed immediately after it is declared).

In [6]:

def define_specific_heat_parameters(self):
    # Constants for specific heat capacity, enthalpy
    self.cp_mol_ig_comp_coeff_A = Var(
            self.component_list,
            initialize={"benzene": -3.392E1,
                        "toluene": -2.435E1,
                        "hydrogen": 2.714e1,
                        "methane": 1.925e1,
                        "diphenyl": -9.707e1},
            units=pyunits.J/pyunits.mol/pyunits.K,
            doc="Parameter A for ideal gas molar heat capacity")
    self.cp_mol_ig_comp_coeff_A.fix()

    self.cp_mol_ig_comp_coeff_B = Var(
            self.component_list,
            initialize={"benzene": 4.739E-1,
                        "toluene": 5.125E-1,
                        "hydrogen": 9.274e-3,
                        "methane": 5.213e-2,
                        "diphenyl": 1.106e0},
            units=pyunits.J/pyunits.mol/pyunits.K**2,
            doc="Parameter B for ideal gas molar heat capacity")
    self.cp_mol_ig_comp_coeff_B.fix()

    self.cp_mol_ig_comp_coeff_C = Var(
            self.component_list,
            initialize={"benzene": -3.017E-4,
                        "toluene": -2.765E-4,
                        "hydrogen": -1.381e-5,
                        "methane": -8.855e-4,
                        "diphenyl": -8.855e-4},
            units=pyunits.J/pyunits.mol/pyunits.K**3,
            doc="Parameter C for ideal gas molar heat capacity")
    self.cp_mol_ig_comp_coeff_C.fix()

    self.cp_mol_ig_comp_coeff_D = Var(
            self.component_list,
            initialize={"benzene": 7.130E-8,
                        "toluene": 4.911E-8,
                        "hydrogen": 7.645e-9,
                        "methane": -1.132e-8,
                        "diphenyl": 2.790e-7},
            units=pyunits.J/pyunits.mol/pyunits.K**4,
            doc="Parameter D for ideal gas molar heat capacity")
    self.cp_mol_ig_comp_coeff_D.fix()

    self.enth_mol_form_vap_comp_ref = Var(
            self.component_list,
            initialize={"benzene": -82.9e3,
                        "toluene": -50.1e3,
                        "hydrogen": 0,
                        "methane": -75e3,
                        "diphenyl": -180e3},
            units=pyunits.J/pyunits.mol,
            doc="Standard heat of formation at reference state")
    self.enth_mol_form_vap_comp_ref.fix()

Declaring the Physical Parameter Block¶

Now that the various parts of the Physical Parameter Block have been declared, we can assemble the actual class that will assemble these components in a flowsheet. There are four parts to declaring our new PhysicalParameterBlock class, which are:

1. Declaring the new class and inheriting from the PhysicalParameterBlock base class
2. Declaring any necessary configuration arguments
3. Writing the build method for our class
4. Creating a define_metadata method for the class.

Each of these steps are shown in the code example below.

First, we need to declare our new class and give it a unique name. In this example, we will call our new class HDAParameterBlock. The first two lines of the example below show how we declare our new class using the declare_process_block_decorator and inheriting from the PhysicalParameterBlock base class from the IDAES Core model libraries. Inheriting from the PhysicalParameterBlock brings us access to all the necessary features required by the IDAES modeling framework, whilst the declare_process_block_class decorator performs some boilerplate operations to replicate the expected object structure of Pyomo. Further details on these components can be found in the IDAES documentation.

Next, we need to set up any configuration arguments we need for the property package. This is done using Pyomo “Config Blocks” which provide a convenient way of declaring, organizing and documenting configuration arguments. To begin with, we can inherit from the CONFIG block declared in the PhysicalParameterBlock base class, which provides all the arguments that the IDAES modeling framework expects to be present. Modelers can then add additional configuration arguments to provide users with options when constructing their property packages, however we will not cover that in this tutorial.

The most significant part of any IDAES model class is the build method, which contains the instructions on how to construct an instance of the desired model and all IDAES models are expected to have a build method. The first step in any build method is to call super().build(), which will trigger the build method of the base class that the current class inherits from – this is important since this is how we automate construction of any underlying components required by the modeling framework and ensure that everything integrates smoothly. Next, a PhysicalParameterBlock needs to contain a pointer to the related StateBlock (which we will look at next) – this is used to allow us to build instances of the StateBlock by only knowing the PhysicalParameterBlock we wish to use. To do this, we create an attribute named _state_block_class attached to our class with a pointer to the StateBlock class; in this case self._state_block_class = HDAStateBlock, where HDAStateBlock is the name of the yet to be declared StateBlock. Finally, the build method needs to construct the actual parameters required for the property package, which we do here by calling the sub-methods written previously.

The final step in creating the PhysicalParameterBlock class is to declare a classmethod named define_metadata which takes two arguments; a class (cls) and an instance of that class (obj). This method in turn needs to call two pre-defined methods (inherited from the underlying base classes):

• obj.add_properties() is used to set the metadata regarding the supported properties, and here we pass the properties_metadata dict we created earlier as an argument.
• obj.add_default_units() sets the default units metadata for the property package, and here we pass the units_metadata dict we created earlier as an argument.

In [7]:

@declare_process_block_class("HDAParameterBlock")
class HDAParameterData(PhysicalParameterBlock):
    CONFIG = PhysicalParameterBlock.CONFIG()

    def build(self):
        '''
        Callable method for Block construction.
        '''
        super(HDAParameterData, self).build()

        self._state_block_class = HDAStateBlock
        
        define_components_and_phases(self)
        define_basic_parameters(self)
        define_specific_heat_parameters(self)

    @classmethod
    def define_metadata(cls, obj):
        """Define properties supported and units."""
        obj.add_properties(properties_metadata)

        obj.add_default_units(units_metadata)

State Block¶

After the Physical Parameter Block class has been created, the next step is to write the code necessary to create the State Blocks that will be used through out the flowsheet. Unlike other models however, creating a State Block actually required us to write two classes. In short, indexed Pyomo object components (e.g. Vars and Blocks) actually consist of two objects: an IndexedComponent object which serves as a container for multiple ComponentData objects which represent the component at each point in the indexing set. For example, a Var indexed by the Set [1, 2, 3, 4] actually consists of a single IndexedVar object which contains 4 VarData objects. Normally this behavior is hidden behind the declare_process_block_data decorator which handles the details of this structure (as a side note, unindexed components similarly involve two classes but this is hidden by the use of multiple inheritance.)

Normally, when we write models in IDAES, we are concerned only with the ComponentData object – i.e., the instructions on how to build an instance of the model at each indexed point (hence the naming convention used when declaring classes). However, State Blocks are slightly different in that we always expect State Blocks to be indexed (they will always be indexed by time, at a minimum). Due to this, we often want to perform actions on all the elements of the indexed State Block at once (rather than element by element), such as during initialization. Thus, we have a need to write methods that are attached to the IndexedStateBlock in addition to the normal methods for the StateBlockData object. Fortunately, the declare_process_block_data decorator facilitates this for us, but it does mean we need to declare two classes when creating State Blocks.

For this example, we will begin by describing the content of the StateBlockData objects, as this is where we create the variables and constraints that describe how to calculate the thermophysical properties of the material. After that, we will discuss how to create the class that contains methods to be applied to the IndexedStateBlock as a whole.

Step 5: Declare State Variables¶

The first step in defining a State Block is to create the “state variables” which will be used to define the “state” of the material at any given point. The concept of a “state variable” in IDAES is much the same as the concept in thermodynamics, with the exception that we include extensive flow information in the state definition in IDAES. In short, the “state variables” should be sufficient to fully define the state of the material (both extensive and intensive), and should result in a State Block with zero degrees of freedom if all the state variables are fixed.

For this example, our state variables will be:

• total molar flow rate (flow_mol),
• overall mole fractions of each component (mole_frac_comp),
• temperature (temperature), and
• pressure (pressure).

We declare these as Pyomo Var components which we will attach to our State Block Data object as shown below. Note that we need to define an initial (default value) for each variable along with its units of measurement (which should be consistent with those declared in the property package metadata). Additionally, modelers should provide bounds for the state variables (based on expected ranges of operation, known limitations of the system or fitted region of property correlations), and it is good practice to provide a documentation string for each.

When declaring variables (or Expressions) for physical properties, it is important to follow the IDAES naming conventions. The standard naming conventions allow unit model developers to know what name will be used for a given property in any property package, meaning that any unit model can work with any property package that uses the naming conventions.

In [8]:

def add_state_variables(self):
    self.flow_mol = Var(
            initialize=1,
            bounds=(1e-8, 1000),
            units=pyunits.mol/pyunits.s,
            doc='Molar flow rate')
    self.mole_frac_comp = Var(self.component_list,
                              initialize=0.2,
                              bounds=(0, None),
                              units=pyunits.dimensionless,
                              doc="Component mole fractions")
    self.pressure = Var(initialize=101325,
                        bounds=(101325, 400000),
                        units=pyunits.Pa,
                        doc='State pressure')
    self.temperature = Var(initialize=298.15,
                           bounds=(298.15, 1500),
                           units=pyunits.K,
                           doc='State temperature')

However, simply adding the state variables to our model is not sufficient, as the framework has no way of distinguishing these variables from any other variables in the property package. In order to provide a way for the framework to identify which variables should be considered to be the state variables, we create a dict which contains a pointer to all the state variables, along with their names (note that the name (key) should match the local name of the variable object).

In [9]:

def return_state_var_dict(self):
    return {"flow_mol": self.flow_mol,
            "mole_frac_comp": self.mole_frac_comp,
            "temperature": self.temperature,
            "pressure": self.pressure}

Step 6: Write Constraints and/or Expressions for Properties of Interest¶

Next, we need to write the equations which we will use to calculate the properties of interest as functions of the state variables. Within Pyomo, there are three ways in which this can be done:

1. by using a Var and a Constraint,
2. by using an Expression, or,
3. by using a Reference.

The different between the first two options is that an Expression does not appear in the problem passed to the solver – the Expression can be evaluated by the user and included in constraints in the same way as a variable, but when the problem is passed to the solver the Expression object is substituted for the expression it represents wherever it appears in the model. This means that there are fewer variables and constraints in the problem the solver sees, but that the constraints that do appear are more complex. There is no simple answer to which approach is best, and different applications may see better results with one form or the other. The third option, using a Reference is for cases where a property already exists elsewhere in the model, and we just want to create a local copy of the same object. In terms of properties, this most often occurs with fixed quantities which are declared in the Physical Parameter Block such as molecular weights. For the purposes of this example, we will demonstrate all of these approaches.

You may recall from the initial problem statement that we have three properties of interest in this example:

• component molecular weights,
• mixture density, and
• mixture specific enthalpy.

Let us start with the component molecular weights. These are of course global parameters for each component which are not going to change, and thus we declared these in the Physical Parameter Block. However, it is possible that a unit model may need these values in some correlation, and it is likely that the modeler would ask the State Block for these values. Thus, to make sure these properties are available in the State Block if needed, we will create a Pyomo Reference to the parameters in the Physical Parameter Block using the following code: self.mw_comp = Reference(self.params.mw_comp).

Next, we will create a Var for the mixture density and a Constraint to calculate it from the state variable. For this example, we are using the ideal gas equation, so $P = \rho \times R \times T$. First, we need to declare a variable for the density of the mixture, which is named dens_mol according to the naming standards. We then need to create the constraint with the ideal gas equation, which requires temperature and pressure (which we have already declared on the State Block), and the gas constant. The IDAES modeling framework contains a library of common constants, which we imported earlier using from idaes.core.util.constants import Constants as const and we can get the gas constant from there using const.gas_constant.

In [10]:

def add_molecular_weight_and_density(self):
    self.mw_comp = Reference(self.params.mw_comp)

    self.dens_mol = Var(initialize=1,
                        units=pyunits.mol/pyunits.m**3,
                        doc="Mixture density")

    self.ideal_gas_eq = Constraint(
        expr=self.pressure ==
        const.gas_constant*self.temperature*self.dens_mol)

The last property of interest we need to declare is the mixture specific enthalpy. For an ideal gas, the mixture specific enthalpy can be calculated from the component specific enthalpies using the following equation:

\begin{equation*} h_{mix}= \sum{x_j \times h_j} \end{equation*}

where $x_j$ is the mole fraction of component $j$. Recall that for this example we are using the following correlation for the component specific enthalpies.

\begin{equation*} h_j – h_{j, ref}= A_j \times (T-T_{ref}) + \frac{B_j}{2}\times (T^2-T_{ref}^2) + \frac{C_j}{3}\times (T^3-T_{ref}^3) + \frac{D_j}{4}\times (T^4-T_{ref}^4) \end{equation*}

For the specific enthalpy, we will create a Pyomo Expression rather than a Var and Constraint. In practice, this is much like creating a Constraint. However, rather than returning an equality between two expressions, an Expression requires a single numerical expression that can be used to compute the quantity of interest.

When writing this Expression, we need a few things from different parts of the property package. First, we need to get the parameters ($A_j$, $B_j$, $C_j$, $D_j$, $h_{j, ref}$ and $T_{ref}$) from the Physical Parameter Block. The State Block base class contains a pointer to the associated Physical Parameter Block named self.params which we can use to directly access these, and we will create a shorthand alias to keep the length of our equation down. We will also create a specific shorthand for the reference temperature, as the full name can get rather long. Then, we need to write out the full expression for the specific molar enthalpy, which also requires the temperature and mole fraction state variables from the State Block. Creating the Expression for the specific molar enthalpy is shown in the code below.

In [11]:

def add_enth_mol(self):
    def enth_rule(b):
        params = self.params
        T = self.temperature
        Tr = params.temperature_ref
        return sum(self.mole_frac_comp[j] * (
                        (params.cp_mol_ig_comp_coeff_D[j]/4)*(T**4-Tr**4) +
                        (params.cp_mol_ig_comp_coeff_C[j]/3)*(T**3-Tr**3) +
                        (params.cp_mol_ig_comp_coeff_B[j]/2)*(T**2-Tr**2) +
                        params.cp_mol_ig_comp_coeff_A[j]*(T-Tr) +
                        params.enth_mol_form_vap_comp_ref[j])
                    for j in self.component_list)

    self.enth_mol = Expression(rule=enth_rule)

We have now created the necessary calculations for the properties of interest in our property package, but we are not quite finished with equations yet. We also need to provide an equation that limits the sum of the component mole fractions to be equal to 1, i.e. $1 = \sum{x_j}$. However, there is a catch here; the mole fractions are part of the state variables, and thus they end up in the streams connecting unit operations together. Thus, at the inlet to a unit model, all of the mole fractions will be defined by the outlet of the previous unit operation. That is, if we add a constraint for the sum of mole fractions at the inlet of a unit operation, we will find we have over-specified the problem and it cannot be solved.

To deal with this issue, the each State Block includes a configuration argument (inherited from the base class) named defined_state. This argument is set by the parent unit model when it creates each instance of s State Block, and is set to True wherever the material state is fully defined (most commonly at the inlets). Thus, in our State Block we can add a check for the defined_state configuration argument (self.config.defined_state), and only write the sum of mole fractions constraint if this is False as shown below. Note also that we have added a scaling factor of 1,000 to the constraint to try to improve the robustness of the model. The IDAES modeling framework includes tools for automatically scaling variables and constraints, but these are beyond the scope of this example.

In [12]:

def add_mole_fraction_constraint(self):
    if self.config.defined_state is False:
        self.mole_fraction_constraint = Constraint(
            expr=1e3 == sum(1e3*self.mole_frac_comp[j]
                            for j in self.component_list))

Step 7: Define an Initialization Routine¶

Now that we have written the variables and constraints required to describe the properties of interest for our material, it would be tempting to think that most of the work is now done. Unfortunately, this is only half the work required as we still need to develop an initialization routine for our property package. Whilst we have provided default values for all the variables in our model as we created them, these are most likely not the actual conditions that the end user will wish to simulate. Given the often highly non-linear nature of physical property calculations, it is more often than not impossible to solve a State Block without providing a set of good initial values for all the variables we have declared; this is the purpose of the initialization routine.

The art of writing good, robust initialization routines is enough for a series of tutorials in itself, so this example will focus on the mechanics of how to write an initialization routine. Fortunately, our current example is simple enough that we do not need to do much to get a set of good initial guesses. Generally speaking, an initialization routine consists of three stages (these stages are general and apply to any model)):

1. Fix the state of the model such that there are no degrees of freedom. For State Blocks, it should only be necessary to fix the state variables to a set of initial guesses provided by the user or unit model, as well as deactivating any constraints like the sum of mole fractions.
2. Iteratively build up a solution for the full model. This often involves multiple steps and can involve deactivating constraints and fixing some variables to reduce complexity, as well as analytically calculating values for variables based on the known state (and any previously calculated variables). Solvers can be called as part of any step to efficiently initialize large numbers of variables simultaneously.
3. Return the state of the model to where it originally started (with the exception of variable values). Any variable that was fixed or constraint that was deactivated during initialization should be unfixed or reactivated, so that the degrees of freedom are restored to what they were before the initialization began.

State Blocks and Indexing¶

You should be aware by now that State Blocks are always indexed, at least by time and possibly by space as well. The methods we wrote above for creating variables and constraints all looked at a single element within the IndexedStateBlock, i.e. a single StateBlockData. This is because we need to write these variables and constraints for each state within the model. However, when it comes to initialization, we generally want to perform the same operations on each StateBlockData - i.e. the initialization for each element is the same. Thus, rather than writing an initialization method for a single element, and then calling it repeatedly as we iterate over all the elements, it would be convenient to write a single method that could initialize all the StateBlockData within an IndexedStateBlock simultaneously. This would greatly reduce the number of solver calls we would need to make, thus reducing the overhead associated with those solver calls.

Fortunately, the IDAES modeling framework supports this and the following code will show how to do this. The important thing to note is that the initialization methods are written to act on an arbitrary IndexedBlock which we refer to as blk (as opposed to the self we used previously). This blk represents the IndexedStateBlock and contains a set of elements which we can iterate over as needed (the StateBlockData objects). By having access to the IndexedStateBlock we can perform actions both on individual StateBlockData and on the overall IndexedStateBlock as needed.

Writing the Initialization Routine¶

For initializing State Blocks, the first step is to get our model to a state where it has no degrees of freedom. As mentioned earlier, fixing all of the state variables should be sufficient to fully define the state of the material – i.e. no degrees of freedom. Additionally, we want to initialize our State Block at a set of conditions that are good initial guesses for the final state of the model once it is finally solved. Of course, the State Block has no way of knowing what these initial values should be, so we depend on the unit model (or the end-user) to provide us with a set of initial values to use – this is done through a dict which we generally call state_args where the keys are the names of the state variables and the values are the initial guesses.

Before we start fixing the state variables, there is the possibility that all the state variables have already been fixed (e.g. by a the unit model during its own initialization routine). To allow us to save some time, we include a state_vars_fixed argument in our State Block initialization methods that lets the unit model tell us if the state variables are already fixed – if this is True then we know we can skip the step of checking the state variables ourselves. If state_vars_fixed is False however, then we need to go and fix all the state variables as part of our initialization routine. To save us the effort of having to code all of this ourselves, the IDAES toolkit contains a utility method named fix_state_vars (which we imported earlier), which takes the state_args dict and then iterates through all the state variables as defined by the State Block (using the dict we declared earlier in the return_state_var_dict sub-method). This method iterates through all the defined state variables and does the following:

1. If the variable is already fixed, it records this and does nothing. If a variable is already fixed, we assume it was fixed for a reason and that we should not change its value.
2. If the variable is not fixed, this is recorded and the method then checks the state_args dict for an initial guess for the variable. If a value is found, the variable is fixed to this value; otherwise, the variable is fixed to its current value.

Finally, the fix_state_vars method returns a dict that records which variables were fixed by the method, so that we can later reverse these changes. In the example below, we refer to this dict as flags.

At this point, all the state variables should now be fixed, but once again we have a small catch – if we fix all the state variables then we have a situation similar to the inlet of a unit where we cannot write a constraint on the sum of mole fractions and still solve the model. Thus we need to deactivate this constraint if it exists (remembering that this constraint will not exist in all State Blocks). We know that this constraint will only exist if defined_state is False, so we start by writing an IF statement to check for this and then use the Pyomo deactivate() method to deactivate the constraint (remembering that we will need to reactivate it later).

Before we move on however, it is probably a good idea to check the degrees of freedom to be sure that they are zero (as expected). There are a number of ways things could go wrong (e.g., the unit model said the state variables were fixed when not all of them were, or we missed a constraint we need to deactivate), so a quick check now might save someone a lot of pain in the future. We can use the IDAES degrees_of_freedom method to check the degrees of freedom of our State Block, and if this is not zero raise an Exception to let the user know something went wrong.

In [13]:

def prepare_state(blk, state_args, state_vars_fixed):
    # Fix state variables if not already fixed
    if state_vars_fixed is False:
        flags = fix_state_vars(blk, state_args)
    else:
        flags = None

    # Deactivate sum of mole fractions constraint
    for k in blk.keys():
        if blk[k].config.defined_state is False:
            blk[k].mole_fraction_constraint.deactivate()

    # Check that degrees of freedom are zero after fixing state vars
    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.")

    return flags

Now that we know we have no degrees of freedom in our State Blocks, we can begin the process of building up initial values for all the other variables in our model. The current example is simple enough that it can be solved directly once the state variables are fixed, but more rigorous models (such as those with phase equilibria) will likely require multiple steps to initialize. The code below shows how we call a solver to solve the entire IndexedStateBlock in one go.

Before we call the solver, we first need to make sure there is actually something to solve; depending on how the State Block is written (e.g. build-on-demand properties and use of Expressions), it is sometimes possible that there are actually no Constraints to be solved in the model. If we try to send a problem like that to a solver, we will likely get back an error message which is not what we want to see. Whilst we know that our State Block will always contain at least one constraint (for mixture density), we will add a check here anyway to show how it is done. First ,we create a counter to keep track of the number of unfixed variables in the system, free_vars. Then we iterate over all the elements of the blk (the IndexedStateBlock) and check how many free variables are in each. We use the number_unfixed_variables() method from the idaes.core.util.model_statistics module to do this, and add the result free_vars for each element. If the final value of free_vars is not zero, then we know there is something to solve for and we can proceed to call a solver; otherwise we know that we can skip this step.

In order to solve the entire IndexedStateBlock, we need to do things slightly differently than normal. The standard Pyomo SolverFactory cannot be applied to indexed blocks, so instead we use the IDAES solve_indexed_block method (imported from idaes.core.initialization) which puts a wrapper around the indexed block so that we can use Pyomo’s solver interface. In order to use this method, we need to provide a Pyomo SolverFactory object (called solver here, which also includes any attached solver options) along with the blk we wish to solve and where to send the solver results (the tee argument). Additionally, we want the user to have the ability to control the output from the solver through the IDAES logger interface, which we do by wrapping the solver call with the following line of code:

with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:

where idaeslog is an instance of the IDAES logger. Note that we send the solver output to this logger by setting tee=slc. In this way, all the output from the solver is passed into the logger allowing users to easily control the output level without needing to send additional arguments to the initialization methods.

If all goes well, the solver will successfully initialize our model and we can move on. However, sometimes the solver will fail catastrophically, in which case we need to make sure that our initialization routine can attempt to recover. In order to do this, we wrap the solver call within a Python try/except statement. This way, if the solver fails badly and returns an Exception, we can capture this and decide how to process – otherwise the execution of our model would terminate with the exception from the solver. In the case we encounter an Exception here, we will record results=None and try to continue with initializing our model in the hope that we can recover.

Finally, it is useful to provide the user with some feedback on how the initialization is proceeding. In the last lines below, we send a message to the IDAES logger with a message saying that the initialization step has been completed and append the final solver status.

In [14]:

def initialize_state(blk, solver, init_log, solve_log):
    # Check that there is something to solve for
    free_vars = 0
    for k in blk.keys():
        free_vars += number_unfixed_variables(blk[k])
    if free_vars > 0:
        # If there are free variables, call the solver to initialize
        try:
            with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
                res = solve_indexed_blocks(solver, [blk], tee=True)#slc.tee)
        except:
            res = None
    else:
        res = None

    init_log.info("Properties Initialized {}.".format(
        idaeslog.condition(res)))

One we have finished initializing all the variables in our State Block, the last step is to return the state of all variables and constraints to what they were before we started; i.e. unfixing any variables that were fixed as part of the initialization, and reactivating any constraints that were deactivated. However, at this point we need to be aware of the fact that the State Block initialization is generally occurring as part of the initialization of a unit model, and that there are often cases where the unit model might want us to delay restoring the state of the State Block. The most common occurrence of this is for State Blocks which represent inlets to unit models; State Block initialization is generally the first step in initializing a unit model, after which the unit model generally needs to do some additional steps to fully initialize itself. In order for the unit model to do this, the state variables for the inlet generally need to be fixed. Thus, rather than having to fix the state variables itself, the unit model would often like the inlet state to remain fixed until the unit model finishes initializing itself.

In order to achieve this, all State Block initialization methods get passed a hold_state argument by the unit model, and if this is True it indicates that the unit model wants the state to remain fixed until it is ready. Within the State Block then, we need to separate the code for restoring the model state into a separate release_state method which takes the flags dict from earlier as an argument. Then, in the State Block initialization routine, we need to check the hold_state argument and do the following:

• if hold_state is True, return the flags dict to the unit model so that it can call the restore_state method later when it is done, or
• if hold_state is False, call the restore_state method immediately, passing the flag dict as an argument.

This is shown in the code below.

In [15]:

def restore_state(blk, flags, hold_state):
    # Return state to initial conditions
    if hold_state is True:
        return flags
    else:
        blk.release_state(flags)

Now, we need to write the restore_state method for our State Block, which needs to perform two tasks. First, it needs to reactivate the sum of mole fractions constraint if it exists, and then it needs to unfix the state variables base on the information stored in the flags dict. Looking back to when we fixed the state variables initially, we will note that it is possible for flags = None (indicating that the state_Vars_fixed argument was True). In this case, we must assume the unit model knows what it is doing and we can skip unfixing the state variables here. Otherwise, we can call the revert_state_vars method (imported from idaes.core.util.initialization) and pass blk and flags as arguments, which will automatically restore all the state variables to their original state. Finally, it is helpful to log a message saying that the state has been restored to assist with debugging in the future.

In [16]:

def unfix_state(blk, flags, outlvl):
    init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties")

    # Reactivate sum of mole fractions constraint
    for k in blk.keys():
        if blk[k].config.defined_state is False:
            blk[k].mole_fraction_constraint.activate()

    if flags is not None:
        # Unfix state variables
        revert_state_vars(blk, flags)

    init_log.info_high("State Released.")

Creating the StateBlock and StateBlockData classes¶

We have now defined almost all the of methods we will need to construct the actual State Block classes, and those that remain will be built directly as part of the classes. To begin, we need to declare the class which will be used to construct the IndexedStateBlock, and that contains the methods associated with initializing the State Block. Whilst this might seem a little counter-intuitive (defining initialization before we define the actual model), it is necessary to do it this way as we need to have defined this class before we can define the StateBlockData class (for reasons that you will see later).

The StateBlock class¶

To begin, we declare a new class which inherits from the IDAES StateBlock base class, which we will call _HDAStateBlock here. Note that we need to be careful with naming here; later we will be automatically generating a class names HDAStateBlock which will be the actual state block that we use for our model (and incidentally the one we created a pointer to back in the Physical Parameter Block class). Thus, we add an underscore to the beginning of this class name to distinguish it from the class that will be generated later, and to indicate that this is for internal use only by the property package.

Next, the _HDAStateBlock class needs to have two methods added to it; initialize and release_state which are used to run the initialization routine for the State Block. Each of these have a set of arguments which must be part of the method definition (you can add more if you have a need for them, but the pre-defined set are a minimum).

The initialize method¶

As the name suggests, the initialize method is used to run the initialization routine for the State Block, and this is where we will use the prepare_state, initialize_state and restore_state methods we wrote previously. The initialize method requires the following arguments to be declared:

• blk: this will be a pointer to an instance of the State Block to be initialized.
• state_args: this is used to pass the ‘dict’ of initial guesses to the initialization routine. This should default to None if not provided. The fix_state_vars method will interpret a value of None as no guesses provided as use the current values instead.
• solver: this argument is used to allow tell the State Block to use a specific solver during initialization, and should be a string recognized by the Pyomo SolverFactory. We generally set this to None in order to signify that IDAES Should use the default solver (which is IPOPT).
• optarg: this argument is used to set any solver options the user desires. Again this is generally set to None to indicate that the default solver settings should be used.
• state_vars_fixed: argument to allow the unit model to inform the State Block that the state variables are already fixed. This should default to False.
• hold_state: argument used to allow the unit model to inform the State Block that the state variables should not be unfixed at the end of the initialization routine and that it should return the flags dict so that the unit model can do this later.
• outlvl: this argument is used to inform the State Block of the current logger output settings, and should default to idaeslog.NOTSET.

The code below shows the initialize method for our current example. The method begins by setting up the logger objects that will be used to record any messages generated by the initialization routine. Next, the IDAES get_solver is used to create the solver_obj object that will be used during initialization and sets the solver options. Then, the prepare_state, initialize_state and restore_state methods are called in order to initialize the State Block, providing any arguments as required. Finally, we log a message to inform the user that the initialization routine has successfully completed.

The release_state method¶

The release_state method is used to handle the unfixing of the state variables and reactivation of any constraints that may have remained deactivated at this point. The release_State method requires the following arguments:

• blk: a pointer to an instance of the State Block for which to unfix the state variables.
• flags: a dict of which state variables were fixed in the first stage of initialization, as returned from the fix_state_vars method. This argument should not have a default, as it needs to be provided.
• outlvl: as before, this argument is used to inform the State Block of the current logger output settings, and should default to idaeslog.NOTSET.

In the release_state method below we directly call the unfix_state method we wrote above with the same arguments as the release_State method.

In [17]:

class _HDAStateBlock(StateBlock):

    def initialize(blk, state_args=None, state_vars_fixed=False,
                   hold_state=False, outlvl=idaeslog.NOTSET,
                   solver=None, optarg=None):

        init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties")
        solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="properties")

        # Create solver
        solver_obj = get_solver(solver, optarg)
        
        flags = prepare_state(blk, state_args, state_vars_fixed)
        initialize_state(blk, solver_obj, init_log, solve_log)
        restore_state(blk, flags, hold_state)

        init_log.info("Initialization Complete")

    def release_state(blk, flags, outlvl=idaeslog.NOTSET):
        unfix_state(blk, flags, outlvl)

The StateBlockData class¶

Finally, we can build the StateBlockData class, which we will call HDAStateBlockData. First, we use the declare_process_block_class decorator but this time we provide two arguments. The first argument is the name of the class that will be automatically constructed for us (HDAStateBlock) whilst the second argument is a reference to the class we wish to use as the base when building the IndexedHDAStateBlock class – i.e. the _HDAStateBlock class we just declared. Then, we declare our new HDAStateBlockData class and inherit from the IDAES StateBlockData base class.

As usual, the first thing we need to define in our new class is a build method, where we will provide the instructions for constructing our property model, and once again the first thing we should do is call super().build() to construct all the underlying components defined by the parent class. After this, we can call the methods we wrote earlier to construct the state variables and the add the calculations for the properties of interest.

However, if you recall from when we defined the properties metadata at the beginning of the example, we decided that the specific molar enthalpy of the mixture would be a “build-on-demand” property (i.e., we provided a specific method in the properties metadata rather than None). Thus, we do not want to call the method to construct the specific molar enthalpy as part of the build method, meaning that we only call the add_state_variables, add_mole_fraction_constraint and add_molecular_weight_and_density as in the build method.

To add the specific molar enthalpy calculation as a “build-on-demand” property, we instead declare a separate method with the name we provided in the properties metadata. Whenever the specific molar enthalpy is required by a unit model it will check to see if the property already exists, and if not it will look up the properties metadata and call the method listed there; i.e. _enth_mol in this case. Thus, we declare another method on our HDAStateBlockData class named enth_mol which takes only the class instance as an argument (self), and then call the add_enth_mol method we created earlier to construct the required Expression.

That is all we need to do in order to construct the variables and constraints we need for the property calculations. However, there are a number of other things we need to define in our State Block Data class. In order to provide much of the flexibility present in the IDAES modeling framework, we defer making decisions on much of the form of the overall model for as long as possible. However, these decisions need to be made at some point, and the State Block Data is where this finally occurs.

The first thing we need to do is to define the state variables used by the property package using the return_state_var_dict dict we declared earlier. This dict is used for many important decisions, such as determining which variables should be fixed by the fix_state_vars method during initialization, and what information should be included in the Ports of the unit models. To do this, we need to define a method named define_state_vars in our State Block Data class, which returns the dict of state variables.

Next, we need to provide methods which tell the unit model (and associated control volumes) how the flow terms should be defined in the material and energy balances. This is a core aspect of the flexibility of the IDAES modeling framework, as different forms for the flow terms may be more or less amenable to different formulations of the property calculations (primarily the state variables). Thus, we let the property package decide on the best form to use, which is achieved by defining methods in the State Block Data which return the desired form of the flow terms. Whenever a unit model or control volume requires a flow term, it looks to these methods in the State Block in order to know what expression should be used for the flow terms. Thus, we need to define two methods in our property package:

• get_material_flow_terms(self, p, j): this method should return the desired expression for the material flow term, indexed by phase p and component j. Note that the p and j indexes are required, even if the property package only contains one phase and/or component. If you do have a case where you only have one phase or component (such as this case, where we only have 1 phase), you can just ignore the indices that are not required. For this example, the material flow term we need to use is total flow rate multiplied by component mole fractions, or self.flow_mol*self.mole_frac_comp[j]
• get_enthalpy_flow_terms(self, p): this method should return the desired enthalpy flow term, indexed only by phase p (enthalpy by component does not make a lot of sense). In this case, the enthalpy flow term we need to use is molar flow rate multiplied by specific enthalpy, i.e. self.flow_mol*self.enth_mol.

For dynamic flowsheets, it is also necessary to define methods for the material and internal energy density of the mixture (required to calculate holdup terms); however, we will not define these in this example.

Next, the IDAES modeling framework gives the user the ability to choose different forms for the material and energy balance equations; however, the tractability of each form is significantly determined by the form of the property calculations. The State Block should specify a default form for the material and energy balances that is expected to be tractable under most circumstance (the end-user can override this as necessary). To do this, we also need to define two methods which set the default balance type to use:

• default_material_balance_type should return an instance of the IDAES MaterialBalanceType Enum (imported from idaes.core).
• default_energy_balance_type should return an instance of the IDAES EnergyBalanceType Enum (imported from idaes.core).

Finally, we need to specify the basis of the material flow terms (mass, mole or other). This is used to automatically convert between different bases as required (e.g. a user can define a custom mass transfer term on a molar basis whilst using a mass basis for the actual material balance). Note that automatic conversion only works for mass and molar basis; the “other” basis is used to indicate forms which cannot be easily converted (i.e., the modeler needs to handle this manually). To define the material flow term basis we define a final method named get_material_flow_basis which returns an instance of the IDAES MaterialFlowBasis Enum (again imported from idaes.core).

In [18]:

@declare_process_block_class("HDAStateBlock",
                             block_class=_HDAStateBlock)
class HDAStateBlockData(StateBlockData):
    """
    Example property package for an ideal gas containing benzene, toluene
    hydrogen, methane and diphenyl.
    """

    def build(self):
        """Callable method for Block construction."""
        super(HDAStateBlockData, self).build()

        add_state_variables(self)
        add_mole_fraction_constraint(self)
        add_molecular_weight_and_density(self)
        
    def _enth_mol(self):
        add_enth_mol(self)

    def define_state_vars(self):
        return return_state_var_dict(self)
        
    def get_material_flow_terms(self, p, j):
        return self.flow_mol * self.mole_frac_comp[j]

    def get_enthalpy_flow_terms(self, p):
        """Create enthalpy flow terms."""
        return self.flow_mol * self.enth_mol

    def default_material_balance_type(self):
        return MaterialBalanceType.componentPhase

    def default_energy_balance_type(self):
        return EnergyBalanceType.enthalpyTotal

    def get_material_flow_basis(self):
        return MaterialFlowBasis.molar

Demonstration¶

We have now finished writing our example property package, so let us put it to use in a simple demonstration.

First, we need to import some basic components from pyomo.environ and idaes.core to build a flowsheet model.

In [19]:

from pyomo.environ import ConcreteModel
from idaes.core import FlowsheetBlock

Next, we create ConcreteModel and a steady-state FlowsheetBlock to contain our example, ot which we attach an instance of our new HDAPhysicalParameterBlock.

We can then create an instance of the HDAStateBlock directly from the parameter block using the build_state_block method. This uses the reference to the State Block class that we attached ot the Physical Parameter Block earlier in the example to automatically create an instance of the correct class. Note that we index the State Block by the time domain, so that it looks like a typical state block in a flowsheet, and we set defined_state = True so that we can set values for all the component mole fractions later.

In [20]:

m = ConcreteModel()

m.fs = FlowsheetBlock(dynamic=False)

m.fs.thermo_props = HDAParameterBlock()

m.fs.state = m.fs.thermo_props.build_state_block(
    m.fs.config.time, defined_state=True)

We now have an instance of our new State Block in our flowsheet, so let’s display it and see what it contains.

In [21]:

m.fs.state.display()

Block fs.state[0.0]

  Variables:
    flow_mol : Molar flow rate
        Size=1, Index=None, Units=mol/s
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None : 1e-08 :     1 :  1000 : False : False :  Reals
    mole_frac_comp : Component mole fractions
        Size=5, Index=fs.thermo_props.component_list, Units=dimensionless
        Key      : Lower : Value : Upper : Fixed : Stale : Domain
         benzene :     0 :   0.2 :  None : False : False :  Reals
        diphenyl :     0 :   0.2 :  None : False : False :  Reals
        hydrogen :     0 :   0.2 :  None : False : False :  Reals
         methane :     0 :   0.2 :  None : False : False :  Reals
         toluene :     0 :   0.2 :  None : False : False :  Reals
    pressure : State pressure
        Size=1, Index=None, Units=Pa
        Key  : Lower  : Value  : Upper  : Fixed : Stale : Domain
        None : 101325 : 101325 : 400000 : False : False :  Reals
    temperature : State temperature
        Size=1, Index=None, Units=K
        Key  : Lower  : Value  : Upper : Fixed : Stale : Domain
        None : 298.15 : 298.15 :  1500 : False : False :  Reals
    dens_mol : Mixture density
        Size=1, Index=None, Units=mol/m**3
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :  None :     1 :  None : False : False :  Reals

  Objectives:
    None

  Constraints:
    ideal_gas_eq : Size=1
        Key  : Lower : Body             : Upper
        None :   0.0 : 98846.0429704433 :   0.0

We can see that our State Block contains a single point in time, which in turn contains the five variables. These are our four state variables (molar flow rate, component mole fraction, temperature and pressure) as well as the mixture density. We also have a single constraint, which is the ideal gas equation used to calculate density. Note that we don’t see the component molecular weights as they are Params (References take on the appearance of the component being referenced) or the molar enthalpy as it is an Expression, not a variable (plus it hasn’t been constructed yet as we haven’t asked for it).

Next, let us check the degrees of freedom in our State Block.

In [22]:

print("Degrees of freedom: ", degrees_of_freedom(m))

Degrees of freedom:  2

This is unexpected: the degrees_of_freedom method is saying there are only 2 degrees of freedom in our State Block, but there are 8 state variables.

However, if we think about the constraints we have written, we are only actually using 2 of the state variables in any constraint (temperature and pressure appear in the ideal gas equation). The molar flowrate and component mole fractions are not actually used anywhere in our model, so they have been excluded from the degrees of freedom calculation. In Pyomo terminology, these variables are “Stale”, and they will not be sent to the solver when it is called. Thus, the two degrees of freedom is in fact correct.

Note that this is only the case because our property package is so simple. Also, the specific enthalpy calculation depends on the component mole fractions, so whilst we could solve the State Block by only specifying temperature and pressure, the value of the specific molar enthalpy would be meaningless.

So, lets set some values for all of the state variables as shown below.

In [23]:

m.fs.state[0].flow_mol.fix(100)
m.fs.state[0].temperature.fix(500)
m.fs.state[0].pressure.fix(350000)
m.fs.state[0].mole_frac_comp["benzene"].fix(0.1)
m.fs.state[0].mole_frac_comp["toluene"].fix(0.4)
m.fs.state[0].mole_frac_comp["hydrogen"].fix(0.4)
m.fs.state[0].mole_frac_comp["methane"].fix(0.1)
m.fs.state[0].mole_frac_comp["diphenyl"].fix(0.0)

Now that we have fixed the values for all the state variables, we would expect that the degrees of freedom should be zero (even though we fixed all 8 variables, only temperature and pressure actually contribute to the degrees of freedom). Let’s check this to be sure.

In [24]:

print("Degrees of freedom: ", degrees_of_freedom(m))

Degrees of freedom:  0

We have the expected zero degrees of freedom, so lets try to run the initialization routine we defined and see if we can solve our model.

In [25]:

m.fs.state.initialize()

m.fs.state[0].dens_mol.display()

2023-03-04 01:30:16 [INFO] idaes.init.fs.state: Properties Initialized optimal - Optimal Solution Found.
2023-03-04 01:30:16 [INFO] idaes.init.fs.state: Initialization Complete
dens_mol : Mixture density
    Size=1, Index=None, Units=mol/m**3
    Key  : Lower : Value             : Upper : Fixed : Stale : Domain
    None :  None : 84.19064853145991 :  None : False : False :  Reals

If all went according to plan, we should see that the initialization routine reports an optimal solution found, and that the calculated molar density is 84.19 $mol/m^3$. So far so good, but let's see if the calculation for specific enthalpy works by displaying its value as well.

In [26]:

m.fs.state[0].enth_mol.display()

enth_mol : Size=1
    Key  : Value
    None : -22169.95123146947

We should hopefully see a value of -22170 $J/mol$, meaning our property package appears to be working.

However, whilst the code ran smoothly and we managed to solve the model, there is still the possibility that an error slipped in somewhere. There are many things we could do to verify and validate our model, but a simple place to start is to check that all the units of measurement in the model make sense. When we declared all the variables in our model, we assigned them units of measurement so we would expect that the units should be consistent in all our model equations.

Pyomo contains a useful tool for automatically verifying that the units of measurement are consistent, so lets import it here and run it over our model.

In [27]:

from pyomo.util.check_units import assert_units_consistent

assert_units_consistent(m)

If all went according to plan, the assert_units_consistent method will return nothing, indicating that our units are consistent, and thus that the first step in verifying our model is complete. If our units were not consistent, we would have seen an error telling us that the units of our model were not consistent, and a pointer to the offending constraint and the inconsistent units (Pyomo only returns the first error encountered, so fixing the first error may reveal another).

Concluding Remarks¶

The above example has hopefully introduced you to the basic requirements for creating your own custom physical property packages. However, it is probably clear that it requires a significant amount of effort to write your own property packages, thus users are encouraged to look into the IDAES Modular Properties Framework if they are not already familiar with this.

The IDAES Modular Properties Framework is designed to automatically generate user-defined property packages for common applications based on a single configuration file. Users provide a list of phases and components of interest, and select from a library of common forms for properties of interest, and the Modular Properties Framework then does the hard work of assembling the necessary code to construct the desired model.

To demonstrate the application of the Modular Properties Framework, an equivalent property package to the one we wrote above can be generated using the following configuration dict.

In [28]:

# Import Modular Property Framework Libraries
from idaes.core import VaporPhase, Component

from idaes.models.properties.modular_properties.state_definitions import FTPx
from idaes.models.properties.modular_properties.eos.ideal import Ideal

import idaes.models.properties.modular_properties.pure.RPP3 as RPP

# Build configuration dictionary
configuration = {
    # Specifying components
    "components": {
        'benzene': {"type": Component,
                    "enth_mol_ig_comp": RPP,
                    "parameter_data": {
                        "mw": (78.1136E-3, pyunits.kg/pyunits.mol),
                        "cp_mol_ig_comp_coeff": {
                            'A': (-3.392E1, pyunits.J/pyunits.mol/pyunits.K),
                            'B': (4.739E-1, pyunits.J/pyunits.mol/pyunits.K**2),
                            'C': (-3.017E-4, pyunits.J/pyunits.mol/pyunits.K**3),
                            'D': (7.130E-8, pyunits.J/pyunits.mol/pyunits.K**4)},
                        "enth_mol_form_vap_comp_ref": (
                            82.9e3, pyunits.J/pyunits.mol)}},
        'toluene': {"type": Component,
                    "enth_mol_ig_comp": RPP,
                    "parameter_data": {
                        "mw": (92.1405E-3, pyunits.kg/pyunits.mol),
                        "cp_mol_ig_comp_coeff": {
                            'A': (-2.435E1, pyunits.J/pyunits.mol/pyunits.K),
                            'B': (5.125E-1, pyunits.J/pyunits.mol/pyunits.K**2),
                            'C': (-2.765E-4, pyunits.J/pyunits.mol/pyunits.K**3),
                            'D': (4.911E-8, pyunits.J/pyunits.mol/pyunits.K**4)},
                        "enth_mol_form_vap_comp_ref": (
                            50.1e3, pyunits.J/pyunits.mol)}},
        'hydrogen': {"type": Component,
                    "enth_mol_ig_comp": RPP,
                    "parameter_data": {
                        "mw": (2.016e-3, pyunits.kg/pyunits.mol),
                        "cp_mol_ig_comp_coeff": {
                            'A': (2.714e1, pyunits.J/pyunits.mol/pyunits.K),
                            'B': (9.274e-3, pyunits.J/pyunits.mol/pyunits.K**2),
                            'C': (-1.381e-5, pyunits.J/pyunits.mol/pyunits.K**3),
                            'D': (7.645e-9, pyunits.J/pyunits.mol/pyunits.K**4)},
                        "enth_mol_form_vap_comp_ref": (
                            0, pyunits.J/pyunits.mol)}},
        'methane': {"type": Component,
                    "enth_mol_ig_comp": RPP,
                    "parameter_data": {
                        "mw": (16.043e-3, pyunits.kg/pyunits.mol),
                        "cp_mol_ig_comp_coeff": {
                            'A': (1.925e1, pyunits.J/pyunits.mol/pyunits.K),
                            'B': (5.213e-2, pyunits.J/pyunits.mol/pyunits.K**2),
                            'C': (-8.855e-4, pyunits.J/pyunits.mol/pyunits.K**3),
                            'D': (-1.132e-8, pyunits.J/pyunits.mol/pyunits.K**4)},
                        "enth_mol_form_vap_comp_ref": (
                            -75e3, pyunits.J/pyunits.mol)}},
        'diphenyl': {"type": Component,
                    "enth_mol_ig_comp": RPP,
                    "parameter_data": {
                        "mw": (154.212e-4, pyunits.kg/pyunits.mol),
                        "cp_mol_ig_comp_coeff": {
                            'A': (-9.707e1, pyunits.J/pyunits.mol/pyunits.K),
                            'B': (1.106e0, pyunits.J/pyunits.mol/pyunits.K**2),
                            'C': (-8.855e-4, pyunits.J/pyunits.mol/pyunits.K**3),
                            'D': (2.790e-7, pyunits.J/pyunits.mol/pyunits.K**4)},
                        "enth_mol_form_vap_comp_ref": (
                            -180e3, pyunits.J/pyunits.mol)}}},

    # Specifying phases
    "phases":  {'Vap': {"type": VaporPhase,
                        "equation_of_state": Ideal}},

    # Set base units of measurement
    "base_units": {"time": pyunits.s,
                   "length": pyunits.m,
                   "mass": pyunits.kg,
                   "amount": pyunits.mol,
                   "temperature": pyunits.K},

    # Specifying state definition
    "state_definition": FTPx,
    "state_bounds": {"flow_mol": (1e-8, 1, 1000, pyunits.mol/pyunits.s),
                     "temperature": (298.15, 298.15, 1500, pyunits.K),
                     "pressure": (101325, 101325, 400000, pyunits.Pa)},
    "pressure_ref": (101325, pyunits.Pa),
    "temperature_ref": (298.15, pyunits.K)}