Fixing spurious DoF in Gibbs reactor model (#1393)

* Catching unneccessary terms in lagrange multiplier expression

* Better implementation

* Adding tests for filtering inert elements

idaes/models/unit_models/gibbs_reactor.py

@@ -14,7 +14,7 @@
 Standard IDAES Gibbs reactor model.
 """
 # Import Pyomo libraries
-from pyomo.environ import Constraint, Param, Reals, Reference, Var
+from pyomo.environ import Constraint, Param, Reals, Reference, Set, Var
 from pyomo.common.config import ConfigBlock, ConfigValue, In, ListOf, Bool
 
 # Import IDAES cores
@@ -212,19 +212,33 @@ def build(self):
 
         # Add performance equations
         # Add Lagrangian multiplier variables
+        # Only need multipliers for species which participate in reactions
+        l_set = []
+        for e in self.control_volume.config.property_package.element_list:
+            skip = True
+            for j in self.control_volume.properties_in.component_list:
+                if j in self.config.inert_species:
+                    continue
+                elif self.control_volume.properties_out.params.element_comp[j][e] != 0:
+                    skip = False
+
+            if not skip:
+                l_set.append(e)
+        self.lagrange_set = Set(initialize=l_set)
+
         e_units = self.control_volume.config.property_package.get_metadata().get_derived_units(
             "energy_mole"
         )
         self.lagrange_mult = Var(
             self.flowsheet().time,
-            self.control_volume.config.property_package.element_list,
+            self.lagrange_set,
             domain=Reals,
             initialize=100,
             doc="Lagrangian multipliers",
             units=e_units,
         )
 
-        # TODO : Remove this once sacling is properly implemented
+        # TODO : Remove this once scaling is properly implemented
         self.gibbs_scaling = Param(default=1, mutable=True)
 
         # Use Lagrangian multiple method to derive equations for Out_Fi
@@ -245,10 +259,8 @@ def gibbs_minimization(b, t, p, j):
                 b.control_volume.properties_out[t].gibbs_mol_phase_comp[p, j]
                 + sum(
                     b.lagrange_mult[t, e]
-                    * b.control_volume.properties_out[t].config.parameters.element_comp[
-                        j
-                    ][e]
-                    for e in b.control_volume.config.property_package.element_list
+                    * b.control_volume.properties_out[t].params.element_comp[j][e]
+                    for e in self.lagrange_set
                 )
             )
 

idaes/models/unit_models/tests/test_gibbs.py

@@ -50,6 +50,12 @@
 )
 from idaes.core.util import DiagnosticsToolbox
 
+# Natural gas property package for integration testing
+from idaes.models.properties.modular_properties.base.generic_property import (
+    GenericParameterBlock,
+)
+from idaes.models_extra.power_generation.properties.natural_gas_PR import get_prop
+
 
 # -----------------------------------------------------------------------------
 # Get default solver for testing
@@ -84,88 +90,139 @@ def test_config():
     assert not hasattr(m.fs.unit, "inert_species_balance")
 
 
-@pytest.mark.unit
-def test_inerts():
-    m = ConcreteModel()
-    m.fs = FlowsheetBlock(dynamic=False)
+class TestGibbsInerts:
+    @pytest.mark.unit
+    def test_no_inerts(self):
+        m = ConcreteModel()
+        m.fs = FlowsheetBlock(dynamic=False)
 
-    m.fs.properties = PhysicalParameterTestBlock()
+        m.fs.properties = PhysicalParameterTestBlock()
 
-    m.fs.unit = GibbsReactor(property_package=m.fs.properties, inert_species=["c1"])
+        m.fs.unit = GibbsReactor(property_package=m.fs.properties)
 
-    assert isinstance(m.fs.unit.inert_species_balance, Constraint)
-    assert len(m.fs.unit.inert_species_balance) == 2
-    assert m.fs.unit.inert_species_balance[0, "p1", "c1"] != Constraint.Skip
-    assert m.fs.unit.inert_species_balance[0, "p2", "c1"] != Constraint.Skip
+        for e in m.fs.unit.lagrange_set:
+            assert e in ["H", "He", "Li"]
 
-    assert isinstance(m.fs.unit.gibbs_minimization, Constraint)
-    assert len(m.fs.unit.gibbs_minimization) == 2
+        assert not hasattr(m.fs.unit, "inert_species_balance")
 
+    @pytest.mark.unit
+    def test_inerts(self):
+        m = ConcreteModel()
+        m.fs = FlowsheetBlock(dynamic=False)
 
-@pytest.mark.unit
-def test_inerts_dependent_w_multi_phase():
-    m = ConcreteModel()
-    m.fs = FlowsheetBlock(dynamic=False)
+        m.fs.properties = PhysicalParameterTestBlock()
 
-    m.fs.properties = PhysicalParameterTestBlock()
-    # Change elemental composition to introduce dependency
-    m.fs.properties.element_comp = {
-        "c1": {"H": 0, "He": 0, "Li": 3},
-        "c2": {"H": 4, "He": 5, "Li": 0},
-    }
+        m.fs.unit = GibbsReactor(property_package=m.fs.properties, inert_species=["c1"])
 
-    m.fs.unit = GibbsReactor(property_package=m.fs.properties, inert_species=["c1"])
+        for e in m.fs.unit.lagrange_set:
+            assert e in ["H", "He", "Li"]
 
-    assert isinstance(m.fs.unit.inert_species_balance, Constraint)
-    assert len(m.fs.unit.inert_species_balance) == 2
-    assert m.fs.unit.inert_species_balance[0, "p1", "c1"] != Constraint.Skip
-    assert m.fs.unit.inert_species_balance[0, "p2", "c1"] != Constraint.Skip
+        assert isinstance(m.fs.unit.inert_species_balance, Constraint)
+        assert len(m.fs.unit.inert_species_balance) == 2
+        assert m.fs.unit.inert_species_balance[0, "p1", "c1"] != Constraint.Skip
+        assert m.fs.unit.inert_species_balance[0, "p2", "c1"] != Constraint.Skip
 
-    assert isinstance(m.fs.unit.gibbs_minimization, Constraint)
-    assert len(m.fs.unit.gibbs_minimization) == 2
+        assert isinstance(m.fs.unit.gibbs_minimization, Constraint)
+        assert len(m.fs.unit.gibbs_minimization) == 2
 
+    @pytest.mark.unit
+    def test_inerts_dependent_w_multi_phase(self):
+        m = ConcreteModel()
+        m.fs = FlowsheetBlock(dynamic=False)
 
-@pytest.mark.unit
-def test_inerts_dependent_w_single_phase():
-    m = ConcreteModel()
-    m.fs = FlowsheetBlock(dynamic=False)
+        m.fs.properties = PhysicalParameterTestBlock()
+        # Change elemental composition to introduce dependency
+        m.fs.properties.element_comp = {
+            "c1": {"H": 0, "He": 0, "Li": 3},
+            "c2": {"H": 4, "He": 5, "Li": 0},
+        }
 
-    m.fs.properties = PhysicalParameterTestBlock()
-    # Set phase list to only have 1 phase
-    m.fs.properties.phase_list = ["p1"]
-    # Change elemental composition to introduce dependency
-    m.fs.properties.element_comp = {
-        "c1": {"H": 0, "He": 0, "Li": 3},
-        "c2": {"H": 4, "He": 5, "Li": 0},
-    }
+        m.fs.unit = GibbsReactor(property_package=m.fs.properties, inert_species=["c1"])
 
-    m.fs.unit = GibbsReactor(property_package=m.fs.properties, inert_species=["c1"])
+        for e in m.fs.unit.lagrange_set:
+            assert e in ["H", "He"]
 
-    assert isinstance(m.fs.unit.inert_species_balance, Constraint)
-    assert len(m.fs.unit.inert_species_balance) == 0
-    assert (0, "p1", "c1") not in m.fs.unit.inert_species_balance
+        assert isinstance(m.fs.unit.inert_species_balance, Constraint)
+        assert len(m.fs.unit.inert_species_balance) == 2
+        assert m.fs.unit.inert_species_balance[0, "p1", "c1"] != Constraint.Skip
+        assert m.fs.unit.inert_species_balance[0, "p2", "c1"] != Constraint.Skip
 
-    assert isinstance(m.fs.unit.gibbs_minimization, Constraint)
-    assert len(m.fs.unit.gibbs_minimization) == 1
+        assert isinstance(m.fs.unit.gibbs_minimization, Constraint)
+        assert len(m.fs.unit.gibbs_minimization) == 2
 
+    @pytest.mark.unit
+    def test_inerts_dependent_w_single_phase(self):
+        m = ConcreteModel()
+        m.fs = FlowsheetBlock(dynamic=False)
 
-@pytest.mark.unit
-def test_invalid_inert():
-    m = ConcreteModel()
-    m.fs = FlowsheetBlock(dynamic=False)
+        m.fs.properties = PhysicalParameterTestBlock()
+        # Set phase list to only have 1 phase
+        m.fs.properties.phase_list = ["p1"]
+        # Change elemental composition to introduce dependency
+        m.fs.properties.element_comp = {
+            "c1": {"H": 0, "He": 0, "Li": 3},
+            "c2": {"H": 4, "He": 5, "Li": 0},
+        }
 
-    m.fs.properties = PhysicalParameterTestBlock()
+        m.fs.unit = GibbsReactor(property_package=m.fs.properties, inert_species=["c1"])
+
+        for e in m.fs.unit.lagrange_set:
+            assert e in ["H", "He"]
+
+        assert isinstance(m.fs.unit.inert_species_balance, Constraint)
+        assert len(m.fs.unit.inert_species_balance) == 0
+        assert (0, "p1", "c1") not in m.fs.unit.inert_species_balance
+
+        assert isinstance(m.fs.unit.gibbs_minimization, Constraint)
+        assert len(m.fs.unit.gibbs_minimization) == 1
 
-    with pytest.raises(
-        ConfigurationError,
-        match="fs.unit invalid component in inert_species "
-        "argument. foo is not in the property package "
-        "component list.",
-    ):
+    @pytest.mark.unit
+    def test_invalid_inert(self):
+        m = ConcreteModel()
+        m.fs = FlowsheetBlock(dynamic=False)
+
+        m.fs.properties = PhysicalParameterTestBlock()
+
+        with pytest.raises(
+            ConfigurationError,
+            match="fs.unit invalid component in inert_species "
+            "argument. foo is not in the property package "
+            "component list.",
+        ):
+            m.fs.unit = GibbsReactor(
+                property_package=m.fs.properties, inert_species=["foo"]
+            )
+
+    @pytest.mark.integration
+    def test_natural_gas_w_inerts(self):
+        m = ConcreteModel()
+        m.fs = FlowsheetBlock(dynamic=False)
+        ng_config = get_prop(
+            components=[
+                "H2",
+                "CO",
+                "H2O",
+                "CO2",
+                "CH4",
+                "C2H6",
+                "C3H8",
+                "C4H10",
+                "N2",
+                "O2",
+                "Ar",
+            ]
+        )
+        m.fs.ng_props = GenericParameterBlock(**ng_config)
         m.fs.unit = GibbsReactor(
-            property_package=m.fs.properties, inert_species=["foo"]
+            has_heat_transfer=True,
+            has_pressure_change=True,
+            inert_species=["N2", "Ar"],
+            property_package=m.fs.ng_props,
         )
 
+        for e in m.fs.unit.lagrange_set:
+            assert e in ["H", "C", "O"]
+
 
 # -----------------------------------------------------------------------------
 class TestMethane(object):