Merge pull request #16 from waikato-ahuora-smart-energy-systems/fix-fphx

Generic property package fixes

property_packages/base/state_block_constraints.py

@@ -0,0 +1,44 @@
+from pyomo.environ import Block, Constraint
+from pyomo.core.base.expression import ScalarExpression, Expression, _GeneralExpressionData, ExpressionData
+from pyomo.core.base.var import ScalarVar, _GeneralVarData, VarData, IndexedVar, Var
+
+from property_packages.utils.add_extra_expressions import add_extra_expressions
+
+
+class StateBlockConstraints:
+
+    """
+    Extended state blocks should inherit from this class.
+    Adds additional expressions and constraints to a state block.
+    """
+
+    def build(blk, *args):
+        add_extra_expressions(blk)
+        blk.constraints = Block()
+
+
+    def constrain(blk, name: str, value: float) -> Constraint | Var | None:
+        """constrain a component by name to a value"""
+        var = getattr(blk, name)
+        return blk.constrain_component(var, value)
+
+
+    def constrain_component(blk, component: Var | Expression, value: float) -> Constraint | Var | None:
+        """
+        Constrain a component to a value
+        """
+        if type(component) == ScalarExpression:
+            c = Constraint(expr=component == value)
+            c.defining_state_var = True
+            blk.constraints.add_component(component.local_name, c)
+            return c
+        elif type(component) in (ScalarVar, _GeneralVarData, VarData, IndexedVar):
+            component.fix(value)
+            return component
+        elif type(component) in (_GeneralExpressionData, ExpressionData):
+            # allowed, but we don't need to fix it (eg. mole_frac_comp in helmholtz)
+            return None
+        else:
+            raise Exception(
+                f"Component {component} is not a Var or Expression: {type(component)}"
+            )

property_packages/helmholtz/helmholtz_extended.py

@@ -1,101 +1,87 @@
-# extends the IDAES helmholtz property package to include additional properties and methods.
+from pyomo.environ import (
+    Expression,
+    SolverFactory,
+    check_optimal_termination,
+    units
+)
 from idaes.models.properties.general_helmholtz.helmholtz_state import HelmholtzStateBlockData, _StateBlock
 from idaes.models.properties.general_helmholtz.helmholtz_functions import HelmholtzParameterBlockData
 from idaes.core import declare_process_block_class
-from property_packages.utils.add_extra_expressions import add_extra_expressions
-from pyomo.environ import Constraint, Block
-from pyomo.core.base.expression import Expression, ScalarExpression, _GeneralExpressionData, ExpressionData
-from pyomo.core.base.var import IndexedVar, ScalarVar, Var, _GeneralVarData,VarData
+from idaes.core.util.initialization import solve_indexed_blocks, revert_state_vars
+from idaes.core.util.model_statistics import degrees_of_freedom
+from idaes.core.util.exceptions import InitializationError
 import idaes.logger as idaeslog
 
+from property_packages.base.state_block_constraints import StateBlockConstraints
+from property_packages.utils.fix_state_vars import fix_state_vars
+
 
 class _ExtendedStateBlock(_StateBlock):
-    """
-    This class contains methods which should be applied to Property Blocks as a
-    whole, rather than individual elements of indexed Property Blocks.
-    """
-    def __init__(self, *args, **kwargs):
-        super().__init__(*args, **kwargs) 
 
-    def initialize(self, *args, **kwargs):
-        hold_state = kwargs.pop("hold_state", False)
-        for i, v in self.items():
-            v.constraints.deactivate()
-        res = super().initialize(*args, **kwargs)
-        flags = {}
-        for i, v in self.items():
-            v.constraints.activate()
-            flags[i] = {}
-            if hold_state:
-                # Fix the required state variables for zero degrees of freedom, and return a dictionary of the flags.
-                if not hasattr(v.constraints,"flow_mass") and not v.flow_mol.is_fixed():
-                    # We need to fix the flow_mol variable
-                    flags[i]["flow_mol"] = True
-                    v.flow_mol.fix()
-                avaliable_constraints = ["enth_mass","temperature","total_energy_flow","entr_mass","entr_mol","smooth_temperature","custom_vapor_frac"]
-                if not v.enth_mol.is_fixed():
-                    # check if any of the constraints exist
-                    found_constraint = False
-                    for constraint in avaliable_constraints:
-                        if hasattr(v.constraints,constraint):
-                            # we don't need to fix the variable
-                            # but we need to remove this from the list of constraints (it can't be used to define pressure)
-                            avaliable_constraints.remove(constraint)
-                            found_constraint = True
-                            break
-                    if not found_constraint:
-                        # we need to fix the variable
-                        flags[i]["enth_mol"] = True
-                        v.enth_mol.fix()
-                if not v.pressure.is_fixed():
-                    # check if any of the constraints exist
-                    found_constraint = False
-                    for constraint in avaliable_constraints:
-                        if hasattr(v.constraints,constraint):
-                            # we don't need to fix the variable
-                            # but we need to remove this from the list of constraints (it can't be used to define pressure)
-                            avaliable_constraints.remove(constraint)
-                            found_constraint = True
-                            break
-                    if not found_constraint:
-                        # we need to fix the variable
-                        flags[i]["pressure"] = True
-                        v.pressure.fix() 
-        return flags
-
-    def release_state(self, flags, outlvl=idaeslog.NOTSET):
-        for i, v in self.items():
-            for key in flags[i]:
-                getattr(v,key).unfix()
+    def initialize(blk, *args, **kwargs):
+        outlvl = kwargs.get("outlvl", idaeslog.NOTSET)
+        init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties")
+        solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="properties")
+
+        flag_dict = fix_state_vars(blk, kwargs.get("state_args", None))
+
+        dof = degrees_of_freedom(blk)
+        if dof != 0:
+            raise InitializationError(
+                f"{blk.name} Unexpected degrees of freedom during "
+                f"initialization at property initialization step: {dof}."
+            )
+
+        res = None
+        opt = SolverFactory('ipopt')
+        with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+            try:
+                res = solve_indexed_blocks(opt, [blk], tee=slc.tee)
+            except ValueError as e:
+                if str(e).startswith("No variables appear"):
+                    # https://github.com/Pyomo/pyomo/pull/3445
+                    pass
+                else:
+                    raise e
+
+        if res is not None and not check_optimal_termination(res):
+            raise InitializationError(
+                f"{blk.name} failed to initialize successfully. Please check "
+                f"the output logs for more information."
+            )
+
+        if kwargs.get("hold_state") is True:
+            return flag_dict
+        else:
+            blk.release_state(flag_dict)
+
+        init_log.info(
+            "Property package initialization: {}.".format(idaeslog.condition(res))
+        )
 
+    def release_state(blk, flags, outlvl=idaeslog.NOTSET):
+        revert_state_vars(blk, flags)
+
 
 @declare_process_block_class("HelmholtzExtendedStateBlock", block_class=_ExtendedStateBlock)
-class HelmholtzExtendedStateBlockData(HelmholtzStateBlockData):
+class HelmholtzExtendedStateBlockData(HelmholtzStateBlockData, StateBlockConstraints):
+
+    def build(blk, *args):
+        HelmholtzStateBlockData.build(blk, *args)
+        StateBlockConstraints.build(blk, *args)
+
+        # rename temperature to old_temperature
+        old_temperature = blk.temperature
+        blk.del_component("temperature")
+        blk.add_component("old_temperature", old_temperature)
+        # create smooth temperature expression
+        blk.add_component("temperature", Expression(
+            expr=blk.old_temperature + (blk.enth_mol / (1 * units.J/units.mol))*0.000001 * units.K
+        ))
 
-    def build(self, *args):
-        super().build(*args)
-        # Add expressions for smooth_temperature, enthalpy in terms of mass, etc.
-        add_extra_expressions(self)
-        # Add a block for constraints, so we can disable or enable them in bulk
-        self.constraints = Block()
-
-    def constrain(self,name:str,value:float):
-        # Value must be a float. TODO: Handle unit conversion.
-        var = getattr(self,name)
-        if type(var) == ScalarExpression:
-            self.constraints.add_component(name, Constraint(expr=var == value))
-        elif type(var) in (ScalarVar, _GeneralVarData, VarData):
-            var.fix(value)
-        elif type(var) in ( _GeneralExpressionData, ExpressionData) :
-            # allowed, but we don't need to fix it (eg. mole_frac_comp in helmholtz)
-            print(f"Variable {self} {name} is an Expression: {type(var)}")
-            pass
-        else:
-            raise Exception(f"Variable {self} {name} is not a Var or Expression: {type(var)}")
 
 @declare_process_block_class("HelmholtzExtendedParameterBlock")
 class HelmholtzExtendedParameterBlockData(HelmholtzParameterBlockData):
     def build(self):
         super().build()
-        # set that we should use the extended state block
-        self._state_block_class = HelmholtzExtendedStateBlock # type: ignore because it'll get created with declare_process_block_class
+        self._state_block_class = HelmholtzExtendedStateBlock # noqa: F821

property_packages/tests/test_helmholtz.py → property_packages/helmholtz/tests/test_helmholtz.py

@@ -1,5 +1,5 @@
 # Build and solve a heater block.
-from ..build_package import build_package
+from property_packages.build_package import build_package
 from pytest import approx
 
 # Import objects from pyomo package 
@@ -9,8 +9,6 @@
 from idaes.core import FlowsheetBlock
 from idaes.core.util.model_statistics import degrees_of_freedom
 from idaes.models.unit_models.heater import Heater
-from idaes.core.util.tables import _get_state_from_port
-
 
 
 def test_helmholtz():
@@ -23,12 +21,12 @@ def test_helmholtz():
     m.fs.heater.inlet.flow_mol.fix(1)
     #m.fs.heater.inlet.vapor_frac.fix(0)
     #m.fs.heater.inlet.temperature.fix(298) # room temperature in K
-    m.fs.heater.inlet.enth_mol.fix(1878.87)
+    m.fs.heater.inlet.enth_mol.fix(1878.71)
     m.fs.heater.inlet.pressure.fix(101325)
     assert degrees_of_freedom(m) == 0
     m.fs.heater.initialize()
     solver = SolverFactory('ipopt')
     solver.solve(m, tee=True)
-    assert value(_get_state_from_port(m.fs.heater.outlet,0).temperature) == approx(298)
+    assert value(m.fs.heater.control_volume.properties_out[0].temperature) == approx(298)
     assert value(m.fs.heater.outlet.pressure[0]) == approx(101325)
     assert value(m.fs.heater.outlet.flow_mol[0]) == approx(1)

property_packages/helmholtz/tests/test_state_definitions.py

@@ -1,39 +1,68 @@
 from ..helmholtz_builder import build_helmholtz_package
 from pytest import approx
 from pyomo.environ import ConcreteModel, SolverFactory, value, units
+from pyomo.core.base.constraint import Constraint, ScalarConstraint
 from idaes.core import FlowsheetBlock
 from idaes.models.unit_models import Compressor
+from idaes.core.util.model_statistics import degrees_of_freedom
 
 def flowsheet():
     m = ConcreteModel()
     m.fs = FlowsheetBlock(dynamic=False)
     m.fs.properties = build_helmholtz_package(["h2o"])
     return m
 
+
+def build_state(m):
+    m.fs.state = m.fs.properties.build_state_block([0], defined_state=True)
+    return m.fs.state[0]
+
+
+def initialize(m):
+    assert degrees_of_freedom(m) == 0
+    m.fs.state.initialize()
+    assert degrees_of_freedom(m) == 0
+
+
 def solve(m):
     solver = SolverFactory('ipopt')
     solver.solve(m, tee=True)
 
+
+def test_constrain():
+    m = flowsheet()
+    sb = build_state(m)
+    # fix flow_mol directly
+    c = sb.constrain("flow_mol", 1)
+    assert c == sb.flow_mol
+    assert c.value == 1
+    assert c.is_fixed()
+
+    # add a constraint for flow_mass
+    c = sb.constrain("flow_mass", 1)
+    assert type(c) == ScalarConstraint
+    assert c in sb.component_data_objects(Constraint)
+    assert getattr(sb.constraints, "flow_mass") == c
+
+
 def test_state_definition():
     m = flowsheet()
-    m.fs.sb = m.fs.properties.build_state_block(m.fs.time)
-    sb = m.fs.sb[0]
+    sb = build_state(m)
     sb.constrain("enth_mass", 1878.87)
     sb.constrain("pressure", 101325)
     sb.constrain("flow_mass", 1)
-    m.fs.sb.initialize()
+    initialize(m)
     solve(m)
     assert value(sb.temperature) == approx(273.5809)
 
 
 def test_state_definition_temp():
     m = flowsheet()
-    m.fs.sb = m.fs.properties.build_state_block(m.fs.time)
-    sb = m.fs.sb[0]
-    sb.constrain("smooth_temperature", 273.5809)
+    sb = build_state(m)
+    sb.constrain("temperature", 273.5809)
     sb.constrain("pressure", 101325)
     sb.constrain("flow_mass", 1)
-    m.fs.sb.initialize()
+    initialize(m)
     solve(m)
     assert value(sb.enth_mass) == approx(1878.712)
 
@@ -48,12 +77,13 @@ def test_initialise_compressor():
     m.fs.compressor = Compressor(property_package=m.fs.properties)
     inlet = m.fs.compressor.control_volume.properties_in[0]
     outlet = m.fs.compressor.control_volume.properties_out[0]
-    inlet.constrain("smooth_temperature", 273.5809)
+    inlet.constrain("temperature", 273.5809)
     inlet.constrain("pressure", 101325)
     inlet.constrain("flow_mass", 1)
+    m.fs.compressor.efficiency_isentropic.fix(1)
     m.fs.compressor.deltaP.fix(100000)
-    m.fs.compressor.initialize()
+    m.fs.compressor.initialize(outlvl=1)
     solve(m)
-    assert value(outlet.temperature) == approx(393.5689573)
+    assert value(outlet.temperature) == approx(273.58047830928655)
 
 

property_packages/modular/builder/common_parsers.py

@@ -3,7 +3,7 @@
 from .base_parser import BuildBase
 from compounds.Compound import Compound
 from pyomo.environ import units as pyunits
-from idaes.models.properties.modular_properties.state_definitions import FTPx
+from idaes.models.properties.modular_properties.state_definitions import FTPx, FPhx
 from idaes.models.properties.modular_properties.phase_equil.bubble_dew import (LogBubbleDew)
 from idaes.core import LiquidPhase, VaporPhase, Component, PhaseType as PT
 from idaes.models.properties.modular_properties.phase_equil import (SmoothVLE)
@@ -246,30 +246,12 @@ def serialise(compounds: List[Compound], valid_states: List[States]) -> Dict[str
         min_critical_temperature = min([compound["CriticalTemperature"].value for compound in compounds])
 
         # TODO: Refactor this logic, need a more versatile approach
-        if (min_critical_temperature > 500):
-            return {
-                "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s),
-                "temperature": (min_melting_point, 300, 500, pyunits.K),
-                "pressure": (5e4, 1e5, 1e6, pyunits.Pa),
-            }
-        elif (min_critical_temperature > 300):
-            return {
-                "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s),
-                "temperature": (min_melting_point, 200, 400, pyunits.K),
-                "pressure": (5e4, 1e5, 1e6, pyunits.Pa),
-            }
-        elif (min_critical_temperature > 120):
-            return {
-                "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s),
-                "temperature": (min_melting_point, 150, 500, pyunits.K),
-                "pressure": (5e4, 1e5, 1e6, pyunits.Pa),
-            }
-        else:
-            return {
-                "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s),
-                "temperature": (min_melting_point, 150, 350, pyunits.K),
-                "pressure": (5e4, 1e5, 1e6, pyunits.Pa),
-            }
+        return {
+            "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s),
+            "temperature": (max(min_melting_point-50,1), 300, 3000, pyunits.K),
+            "pressure": (5e4, 1e5, 1e6, pyunits.Pa),
+        }
+
 
 class state_definition_parser(BuildBase):
     @staticmethod