Merge pull request #14 from waikato-ahuora-smart-energy-systems/state…

…_initialisation

State initialisation

property_packages/helmholtz/helmholtz_builder.py

@@ -6,6 +6,7 @@
     StateVars,
     AmountBasis
 )
+from .helmholtz_extended import HelmholtzExtendedParameterBlock
 from .parameters import register_compounds
 
 # Add the "parameters" directory to the path so that the Helmholtz EOS can find the parameter files.
@@ -22,7 +23,7 @@ def build_helmholtz_package(compound_list: List[str]):
         component = compound_list[0]
         if component in registered_components():
             # Build and return a helmholtz property package for this compound
-            return HelmholtzParameterBlock(pure_component=component,
+            return HelmholtzExtendedParameterBlock(pure_component=component,
                                             phase_presentation=PhaseType.MIX,
                                             state_vars=StateVars.PH,
                                             amount_basis=AmountBasis.MOLE)

property_packages/helmholtz/helmholtz_extended.py

@@ -0,0 +1,101 @@
+# extends the IDAES helmholtz property package to include additional properties and methods.
+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
+import idaes.logger as idaeslog
+
+
+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","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()
+
+
+@declare_process_block_class("HelmholtzExtendedStateBlock", block_class=_ExtendedStateBlock)
+class HelmholtzExtendedStateBlockData(HelmholtzStateBlockData):
+
+    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

property_packages/helmholtz/tests/test_state_definitions.py

@@ -0,0 +1,59 @@
+from ..helmholtz_builder import build_helmholtz_package
+from pytest import approx
+from pyomo.environ import ConcreteModel, SolverFactory, value, units
+from idaes.core import FlowsheetBlock
+from idaes.models.unit_models import Compressor
+
+def flowsheet():
+    m = ConcreteModel()
+    m.fs = FlowsheetBlock(dynamic=False)
+    m.fs.properties = build_helmholtz_package(["h2o"])
+    return m
+
+def solve(m):
+    solver = SolverFactory('ipopt')
+    solver.solve(m, tee=True)
+
+def test_state_definition():
+    m = flowsheet()
+    m.fs.sb = m.fs.properties.build_state_block(m.fs.time)
+    sb = m.fs.sb[0]
+    sb.constrain("enth_mass", 1878.87)
+    sb.constrain("pressure", 101325)
+    sb.constrain("flow_mass", 1)
+    m.fs.sb.initialize()
+    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.constrain("pressure", 101325)
+    sb.constrain("flow_mass", 1)
+    m.fs.sb.initialize()
+    solve(m)
+    assert value(sb.enth_mass) == approx(1878.712)
+
+
+def test_initialise_compressor():
+    # The purpose of this test is to make sure the compressor can initialise.
+    # We have had some errors with Too Few Degrees of Freedom
+    # being thrown in initialisation. This is because constraints are 
+    # being fixed in the inlet state block instead of variables.
+    # This asserts that the appropriate variables are unfixed.
+    m = flowsheet()
+    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("pressure", 101325)
+    inlet.constrain("flow_mass", 1)
+    m.fs.compressor.deltaP.fix(100000)
+    m.fs.compressor.initialize()
+    solve(m)
+    assert value(outlet.temperature) == approx(393.5689573)
+
+

property_packages/modular/modular_extended.py

@@ -0,0 +1,115 @@
+# extends the IDAES helmholtz property package to include additional properties and methods.
+from idaes.models.properties.general_helmholtz.helmholtz_state import HelmholtzStateBlockData, _StateBlock
+from idaes.models.properties.general_helmholtz.helmholtz_functions import HelmholtzParameterBlockData
+from idaes.models.properties.modular_properties.base.generic_property import GenericParameterBlock, _GenericStateBlock, GenericParameterData, GenericStateBlockData
+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
+import idaes.logger as idaeslog
+
+# NOTE:
+# THis only works for FTPx formulation right now.
+
+
+class _ExtendedGenericStateBlock(_GenericStateBlock):
+    """
+    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)  # Missing argument
+
+    def initialize(self, *args, **kwargs):
+        print("GETTING STATE")
+        hold_state = kwargs.pop("hold_state", False)
+        for i, v in self.items():
+            print(f"State block {i}")
+            print(f"block data: {v}")
+            v.constraints.deactivate()
+        print("ACTIVATING CONSTRAINTS")
+        res = super().initialize(*args, **kwargs)
+        flags = {}
+        print("STATE BLOCKS")
+        print(self.items())
+        for i, v in self.items():
+            print(f"STATE {i}")
+            print(f"BLOCK DATA {v}")
+            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","entr_mass","entr_mol","smooth_temperature","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() 
+        print("WE DONE")
+        return flags
+
+    def release_state(self, flags, outlvl=idaeslog.NOTSET):
+        print("RELEASING STATE")
+        for i, v in self.items():
+            for key in flags[i]:
+                getattr(v,key).unfix()
+
+
+@declare_process_block_class("GenericExtendedStateBlock", block_class=_ExtendedGenericStateBlock)
+class GenericExtendedStateBlockData(GenericStateBlockData):
+
+    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("GenericExtendedParameterBlock")
+class GenericExtendedParameterData(GenericParameterData):
+    def build(self):
+        super().build()
+        # set that we should use the extended state block
+        self._state_block_class = GenericExtendedStateBlock # type: ignore because it'll get created with declare_process_block_class

property_packages/modular/tests/test_state_definitions.py

@@ -0,0 +1,65 @@
+from ..template_builder import build_config
+from pytest import approx
+from pyomo.environ import ConcreteModel, SolverFactory, value, units
+from idaes.core import FlowsheetBlock
+from idaes.models.unit_models import Compressor
+
+def flowsheet():
+    m = ConcreteModel()
+    m.fs = FlowsheetBlock(dynamic=False)
+    m.fs.properties = build_config("peng-robinson",["benzene","toluene"],["Liq","Vap"])
+    return m
+
+def solve(m):
+    solver = SolverFactory('ipopt')
+    solver.solve(m, tee=True)
+
+def test_state_definition2():
+    m = flowsheet()
+    m.fs.sb = m.fs.properties.build_state_block(m.fs.time)
+    sb = m.fs.sb[0]
+    sb.constrain("temperature", 280)
+    sb.constrain("pressure", 101325)
+    sb.constrain("flow_mass", 1)
+    sb.mole_frac_comp["benzene"].fix(0.5)
+    sb.mole_frac_comp["toluene"].fix(0.5)
+    m.fs.sb.initialize()
+    solve(m)
+    assert value(sb.temperature) == approx(280)
+
+
+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.constrain("pressure", 101325)
+    sb.constrain("flow_mass", 1)
+    sb.mole_frac_comp["benzene"].fix(0.5)
+    sb.mole_frac_comp["toluene"].fix(0.5)
+    m.fs.sb.initialize()
+    solve(m)
+    assert value(sb.enth_mass) == approx(1878.712)
+
+
+def test_initialise_compressor():
+    # The purpose of this test is to make sure the compressor can initialise.
+    # We have had some errors with Too Few Degrees of Freedom
+    # being thrown in initialisation. This is because constraints are 
+    # being fixed in the inlet state block instead of variables.
+    # This asserts that the appropriate variables are unfixed.
+    m = flowsheet()
+    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("pressure", 101325)
+    inlet.constrain("flow_mass", 1)
+    inlet.mole_frac_comp["benzene"].fix(0.5)
+    inlet.mole_frac_comp["toluene"].fix(0.5)
+    m.fs.compressor.deltaP.fix(100000)
+    m.fs.compressor.initialize()
+    solve(m)
+    assert value(outlet.temperature) == approx(393.5689573)
+
+

property_packages/utils/add_extra_expressions.py

@@ -0,0 +1,18 @@
+from pyomo.environ import Block, Expression, value, units
+
+def add_extra_expressions(sb: Block):
+    """
+    idaes state blocks don't support all the properties we need, so we add some extra expressions here
+    """
+    #if not hasattr(sb, "smooth_pressure"):
+    #    sb.add_component("smooth_pressure", Expression(expr=sb.pressure + value(sb.enth_mol)*0.0001 * units.Pa))
+    if not hasattr(sb, "smooth_temperature"):
+        sb.add_component("smooth_temperature", Expression(expr=sb.temperature + (sb.enth_mol / (1 * units.J/units.mol))*0.000001 * units.K))
+    if not hasattr(sb, "enth_mass"):
+        sb.add_component("enth_mass", Expression(expr=(sb.flow_mol * sb.enth_mol) / sb.flow_mass))
+    if not hasattr(sb, "entr_mass"):
+        sb.add_component("entr_mass", Expression(expr=(sb.flow_mol * sb.entr_mol) / sb.flow_mass))
+    if not hasattr(sb, "entr_mol"):
+        sb.add_component("entr_mol", Expression(expr=(sb.flow_mol * sb.entr_mass) / sb.flow_mass))
+    if not hasattr(sb, "total_energy_flow"):
+        sb.add_component("total_energy_flow", Expression(expr=sb.flow_mass * sb.enth_mass))

pyproject.toml

@@ -6,7 +6,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "ahuora-compounds"
-version = "0.0.3"
+version = "0.0.4"
 
 authors = [
   { name="Example Author", email="[email protected]" },