Scaler for equilibrium reactor and saponification properties

docs/explanations/scaling_toolbox/scaling_theory.rst

@@ -14,7 +14,7 @@ Types of Scaling
 * Constraint residual scaling refers to the magnitude of the residual of all constraints in the model. This is important as this is required to determine whether or not a model has converged, thus it is important that a residual equal to the solver tolerance does not significantly alter the solution to the model.
 
   * E.g. consider a constraint A=B. If the magnitude of A and B are 1e6 and the solver tolerance is 1e-6, this means that A and B need to be solved to 12 significant figures of precision to converge the constraint (which may be unnecessarily onerous). Similarly, if A and B were of order 1e-6, then they will be considered equal if A=1e-6 and B=0; i.e. ±100% error.
-  
+
 * Jacobian scaling refers to the overall scaling and conditioning of the problem Jacobian. This is important as this determines the ability of the solver to find a path from the initial state to the final solution. It is important to ensure that the Jacobian is well-conditioned both in order to get good numerical behavior from the solver and also to ensure that floating point round-off error does not result in large changes to the variables at the solution.
 
 These aspects are not always complimentary, and there are often cases where improving one aspect of the model scaling can negatively affect another aspect (e.g., focusing too much on Jacobian condition number can often result in poor constraint residual tolerances). Additionally, each aspect affects different parts of the solver routine; Jacobian scaling is important for determining the step direction and size at each iteration, whilst constraint residual scaling is important for determining if and when a solution is found. The relative importance of each of these depends on the solver being used (e.g., for IPOPT constraint residual scaling is more important than Jacobian scaling as IPOPT does its own internal scaling of the Jacobian), however in general all of these will have an impact upon the solver behavior and users should endeavor to have good scaling for all aspects of the model. 

docs/reference_guides/model_libraries/generic/unit_models/equilibrium.rst

@@ -51,3 +51,9 @@ EquilibriumReactorData Class
 
 .. autoclass:: EquilibriumReactorData
   :members:
+
+EquilibriumReactorScaler Class
+------------------------------
+
+.. autoclass:: EquilibriumReactorScaler
+  :members:

idaes/core/base/property_base.py

@@ -593,6 +593,14 @@ def include_inherent_reactions(self):
         # pylint: disable-next=protected-access
         return self.parent_component()._include_inherent_reactions()
 
+    @property
+    def default_initializer(self):
+        return self.parent_component().default_initializer
+
+    @property
+    def default_scaler(self):
+        return self.parent_component().default_scaler
+
     def build(self):
         """
         General build method for StateBlockDatas.

idaes/core/base/reaction_base.py

@@ -267,6 +267,14 @@ def phase_list(self):
     def phase_component_set(self):
         return self.state_ref.phase_component_set
 
+    @property
+    def default_initializer(self):
+        return self.parent_component().default_initializer
+
+    @property
+    def default_scaler(self):
+        return self.parent_component().default_scaler
+
     def lock_attribute_creation_context(self):
         """Returns a context manager that does not allow attributes to be created
         while in the context and allows attributes to be created normally outside

idaes/core/scaling/scaler_profiling.py

@@ -19,7 +19,11 @@
 from pyomo.common.tempfiles import TempfileManager
 
 from idaes.core.util.scaling import jacobian_cond
-from idaes.core.scaling import AutoScaler, CustomScalerBase
+from idaes.core.scaling.autoscaling import AutoScaler
+from idaes.core.scaling.custom_scaler_base import (
+    CustomScalerBase,
+    ConstraintScalingScheme,
+)
 from idaes.core.solvers import get_solver
 
 
@@ -83,23 +87,23 @@ def __init__(
                 "Vars Only": (None, {}),
                 "Harmonic": (
                     cscaler.scale_constraint_by_nominal_value,
-                    {"scheme": "harmonic_mean"},
+                    {"scheme": ConstraintScalingScheme.harmonicMean},
                 ),
                 "Inverse Sum": (
                     cscaler.scale_constraint_by_nominal_value,
-                    {"scheme": "inverse_sum"},
+                    {"scheme": ConstraintScalingScheme.inverseSum},
                 ),
                 "Inverse Root Sum Squares": (
                     cscaler.scale_constraint_by_nominal_value,
-                    {"scheme": "inverse_root_sum_squared"},
+                    {"scheme": ConstraintScalingScheme.inverseRSS},
                 ),
                 "Inverse Maximum": (
                     cscaler.scale_constraint_by_nominal_value,
-                    {"scheme": "inverse_maximum"},
+                    {"scheme": ConstraintScalingScheme.inverseMaximum},
                 ),
                 "Inverse Minimum": (
                     cscaler.scale_constraint_by_nominal_value,
-                    {"scheme": "inverse_minimum"},
+                    {"scheme": ConstraintScalingScheme.inverseMinimum},
                 ),
                 "Nominal L1 Norm": (
                     cscaler.scale_constraint_by_nominal_derivative_norm,

idaes/core/scaling/tests/test_custom_scaling_integration.py

@@ -38,8 +38,6 @@
 from idaes.models.properties.activity_coeff_models.methane_combustion_ideal import (
     MethaneParameterBlock as MethaneCombustionParameterBlock,
 )
-from idaes.core.util.testing import PhysicalParameterTestBlock, initialization_tester
-from idaes.core.solvers import get_solver
 from idaes.core.util import to_json, from_json, StoreSpec
 from idaes.core.util.scaling import jacobian_cond
 from idaes.core.scaling import AutoScaler, CustomScalerBase, set_scaling_factor

idaes/models/properties/examples/saponification_reactions.py

@@ -31,11 +31,13 @@
 from idaes.core.util.misc import add_object_reference
 from idaes.core.util.constants import Constants as const
 import idaes.logger as idaeslog
+from idaes.core.scaling import CustomScalerBase
 
 
 # Some more information about this module
 __author__ = "Andrew Lee"
 
+from idaes.core.util.scaling import get_scaling_factor
 
 # Set up logger
 _log = idaeslog.getLogger(__name__)
@@ -111,12 +113,77 @@
         )
 
 
+class SaponificationReactionScaler(CustomScalerBase):
+    DEFAULT_SCALING_FACTORS = {"reaction_rate": 1e2}
+
+    def variable_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+        if model.is_property_constructed("k_rxn"):
+            # First check to see if k_rxn is already scaled
+            sf = get_scaling_factor(model.k_rxn)
+
+            if sf is not None and not overwrite:
+                # k_rxn already has a scaling factor and we are not set to overwrite - move on
+                pass
+            else:
+                # Hopefully temperature has been scaled, so we can get the nominal value of k_rxn
+                # by walking the expression in the constraint.
+                nominals = self.get_expression_nominal_values(model.arrhenius_eqn)
+
+                # We should get two values, k_rxn (LHS) and the Arrhenius equation (RHS)
+                # As of 10/3/2024, the LHS will be the 0-th element of the list, and the RHS the 1st
+                # However, we cannot assume this will always be the case
+
+                # If LHS has been scaled, nominal will be 1/sf, otherwise it will be k_rxn.value
+                # Find the value which does NOT match this - guess that this is the 1st element
+                if nominals[1] != model.k_rxn.value and sf is None:
+                    # This is the most likely case, so check it first
+                    nominal = nominals[1]
+                elif sf is not None and nominals[1] != 1 / sf:
+                    # Next, check for case where k_rxn was already scaled
+                    nominal = nominals[1]
+                else:
+                    # Otherwise we have the case where something changed in Pyomo since 10/3/2024
+                    nominal = nominals[0]
+
+                self.set_variable_scaling_factor(
+                    model.k_rxn,
+                    1 / nominal,
+                    overwrite=overwrite,
+                )
+
+        if model.is_property_constructed("reaction_rate"):
+            for j in model.reaction_rate.values():
+                self.scale_variable_by_default(j, overwrite=overwrite)
+
+    def constraint_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+        if model.is_property_constructed("arrhenius_eqn"):
+            self.scale_constraint_by_nominal_value(
+                model.arrhenius_eqn,
+                scheme="inverse_maximum",
+                overwrite=overwrite,
+            )
+
+        if model.is_property_constructed("rate_expression"):
+            for j in model.rate_expression.values():
+                self.scale_constraint_by_nominal_value(
+                    j,
+                    scheme="inverse_maximum",
+                    overwrite=overwrite,
+                )
+
+
 class _ReactionBlock(ReactionBlockBase):
     """
     This Class contains methods which should be applied to Reaction Blocks as a
     whole, rather than individual elements of indexed Reaction Blocks.
     """
 
+    default_scaler = SaponificationReactionScaler
+
     def initialize(blk, outlvl=idaeslog.NOTSET, **kwargs):
         """
         Initialization routine for reaction package.

idaes/models/properties/examples/saponification_thermo.py

@@ -44,6 +44,7 @@
 from idaes.core.util.model_statistics import degrees_of_freedom
 from idaes.core.util.initialization import fix_state_vars, revert_state_vars
 import idaes.logger as idaeslog
+from idaes.core.scaling import CustomScalerBase
 
 # Some more information about this module
 __author__ = "Andrew Lee"
@@ -135,12 +136,48 @@
         )
 
 
+class SaponificationPropertiesScaler(CustomScalerBase):
+    UNIT_SCALING_FACTORS = {
+        # "QuantityName: (reference units, scaling factor)
+        "Pressure": (units.Pa, 1e-5),
+    }
+
+    DEFAULT_SCALING_FACTORS = {
+        "flow_vol": 1e2,
+        "conc_mol_comp": 1e-2,
+    }
+
+    def variable_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+        self.scale_variable_by_default(model.flow_vol, overwrite=overwrite)
+        self.scale_variable_by_units(model.pressure, overwrite=overwrite)
+        self.scale_variable_by_bounds(model.temperature, overwrite=overwrite)
+        for k, v in model.conc_mol_comp.items():
+            if k == "H2O":
+                self.set_variable_scaling_factor(v, 1e-4, overwrite=overwrite)
+            else:
+                self.scale_variable_by_default(v, overwrite=overwrite)
+
+    def constraint_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+        if model.is_property_constructed("conc_water_eqn"):
+            self.set_constraint_scaling_factor(
+                model.conc_water_eqn,
+                1e-4,
+                overwrite=overwrite,
+            )
+
+
 class _StateBlock(StateBlock):
     """
     This Class contains methods which should be applied to Property Blocks as a
     whole, rather than individual elements of indexed Property Blocks.
     """
 
+    default_scaler = SaponificationPropertiesScaler
+
     def fix_initialization_states(self):
         """
         Fixes state variables for state blocks.

idaes/models/properties/examples/tests/test_saponification_reaction.py

@@ -14,18 +14,20 @@
 Tests for saponification property package example.
 Authors: Andrew Lee
 """
+from math import exp
 import pytest
-from pyomo.environ import ConcreteModel, Constraint, Param, units, value, Var
-from idaes.core import MaterialFlowBasis
 
+from pyomo.environ import ConcreteModel, Constraint, Param, Suffix, units, value, Var
+
+from idaes.core import MaterialFlowBasis
 from idaes.models.properties.examples.saponification_reactions import (
     SaponificationReactionParameterBlock,
     ReactionBlock,
+    SaponificationReactionScaler,
 )
 from idaes.models.properties.examples.saponification_thermo import (
     SaponificationParameterBlock,
 )
-
 from idaes.core.solvers import get_solver
 
 
@@ -100,6 +102,8 @@ def test_build(self, model):
         assert model.rxns[1].temperature_ref is model.props[1].temperature
         assert model.rxns[1].dh_rxn is model.rparams.dh_rxn
 
+        assert model.rxns[1].default_scaler is SaponificationReactionScaler
+
     @pytest.mark.unit
     def test_rate_constant(self, model):
         assert isinstance(model.rxns[1].k_rxn, Var)
@@ -124,3 +128,90 @@ def test_initialize(self, model):
 
     def check_units(self, model):
         units.assert_units_consistent(model)
+
+
+class TestSaponificationReactionScaler(object):
+    @pytest.mark.unit
+    def test_variable_scaling_routine(self):
+        model = ConcreteModel()
+        model.pparams = SaponificationParameterBlock()
+        model.rparams = SaponificationReactionParameterBlock(
+            property_package=model.pparams
+        )
+
+        model.props = model.pparams.build_state_block([1])
+        model.rxns = model.rparams.build_reaction_block([1], state_block=model.props)
+
+        # Trigger build of reaction properties
+        model.rxns[1].reaction_rate
+
+        scaler = model.rxns[1].default_scaler()
+        assert isinstance(scaler, SaponificationReactionScaler)
+
+        scaler.variable_scaling_routine(model.rxns[1])
+
+        assert isinstance(model.rxns[1].scaling_factor, Suffix)
+
+        sfx = model.rxns[1].scaling_factor
+        assert len(sfx) == 2
+        assert sfx[model.rxns[1].k_rxn] == pytest.approx(
+            1 / (3.132e6 * exp(-43000 / (8.31446262 * 298.15))), rel=1e-8
+        )
+        assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8)
+
+    @pytest.mark.unit
+    def test_constraint_scaling_routine(self):
+        model = ConcreteModel()
+        model.pparams = SaponificationParameterBlock()
+        model.rparams = SaponificationReactionParameterBlock(
+            property_package=model.pparams
+        )
+
+        model.props = model.pparams.build_state_block([1])
+        model.rxns = model.rparams.build_reaction_block([1], state_block=model.props)
+
+        # Trigger build of reaction properties
+        model.rxns[1].reaction_rate
+
+        scaler = model.rxns[1].default_scaler()
+        assert isinstance(scaler, SaponificationReactionScaler)
+
+        scaler.constraint_scaling_routine(model.rxns[1])
+
+        assert isinstance(model.rxns[1].scaling_factor, Suffix)
+
+        sfx = model.rxns[1].scaling_factor
+        assert len(sfx) == 2
+        assert sfx[model.rxns[1].arrhenius_eqn] == pytest.approx(1, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx(1e-4, rel=1e-8)
+
+    @pytest.mark.unit
+    def test_scale_model(self):
+        model = ConcreteModel()
+        model.pparams = SaponificationParameterBlock()
+        model.rparams = SaponificationReactionParameterBlock(
+            property_package=model.pparams
+        )
+
+        model.props = model.pparams.build_state_block([1])
+        model.rxns = model.rparams.build_reaction_block([1], state_block=model.props)
+
+        # Trigger build of reaction properties
+        model.rxns[1].reaction_rate
+
+        scaler = model.rxns[1].default_scaler()
+        assert isinstance(scaler, SaponificationReactionScaler)
+
+        scaler.scale_model(model.rxns[1])
+
+        assert isinstance(model.rxns[1].scaling_factor, Suffix)
+
+        sfx = model.rxns[1].scaling_factor
+        k_rxn_sf = 1 / (3.132e6 * exp(-43000 / (8.31446262 * 298.15)))
+        assert len(sfx) == 4
+        assert sfx[model.rxns[1].k_rxn] == pytest.approx(k_rxn_sf, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].arrhenius_eqn] == pytest.approx(k_rxn_sf, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx(
+            1e-4 * k_rxn_sf, rel=1e-8
+        )