Scaler for equilibrium reactor and saponification properties (#1500)

* Adding infrastructure to support ipopt_v2 * Moving core/util to ipopt_v2 * Moving MH initializer to ipopt_v2 * Fixing pint version issue * Set TSA to use old IPOPT interface * Trying to resolve Windows failures * Working on platofrm dependent failure * BTInitializer with presolve * Moving last bits of core code to ipopt_v2 * Starting on idaes/models * Removing ma57_automatic_scaling default and updating idaes/models/control * idaes/model/properties part 1 * Remaining parts of idaes/models/proeprties * Fixing typo * Switching idaes/models/unit_models to ipopt_v2 * Attempt to work around HXLC issues for now * Some clean up * Switching modular properties initializer to solver indexed blocks * Addressing comments * Recovering from previous branch * Some clean up * Adding ScalerBase class and tests * Working on CustomScalerBase * Nominal value constraint scaling * Adding some initial integration tests for scaling * Some more nominal magnitude scaling approaches * Prototyping pseudojacobian scaler * Trying to debug pseudojacobian * Removing unnecessary import * Addressing pylint issues * Cleaning up nominal jacobian * More methods for CustomScalerBase * Prototyping Gibbs reactor scaler * Gibbs reactor constraint scaling * Working on testing and profiling * Refining Gibbs scaler * Refining nominal value walker * Fixing walker tests * Testing GibbsScaler with initialization * Fixing auto norm scaling on indexed blocks * Testing scaling profiler * Fixing typos * Fixing pylint issue * Imrpoving some doc strings * Apply suggestions from code review Co-authored-by: MarcusHolly <[email protected]> * Fixing issue with autoscaling vars with value None * Adding profiler to __init__ * Fixing name for RSS method * Fixing import of pyomo.environ * Allowing default scaling for indexed components * Adding catch for critical solver failure in profiler * Starting on docs * Finishing docs and unifiying method names * Profiler report methods and docs * Fixing typos * Pylint: fix unnecessary f-string * Working on equilibrium reactor scaler * Finishing tests for equilibrium reactor scaler * Updaing requil scaler * Fixing typo * Addressing pylint issue * Improving guess for heat duty scaling * Updating test value for heat duty scaling * Minor fix * Adding Enum to __init__ * Adding doc strings for saponification scalers * Fixing typo --------- Co-authored-by: MarcusHolly <[email protected]>
IDAES · Oct 22, 2024 · 96237b9 · 96237b9
1 parent 7369728
commit 96237b9
Show file tree

Hide file tree

Showing 14 changed files with 900 additions and 18 deletions.
diff --git a/docs/explanations/scaling_toolbox/scaling_theory.rst b/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. 

diff --git a/docs/reference_guides/model_libraries/generic/unit_models/equilibrium.rst b/docs/reference_guides/model_libraries/generic/unit_models/equilibrium.rst
@@ -51,3 +51,9 @@ EquilibriumReactorData Class
 
 .. autoclass:: EquilibriumReactorData
   :members:
+
+EquilibriumReactorScaler Class
+------------------------------
+
+.. autoclass:: EquilibriumReactorScaler
+  :members:
diff --git a/idaes/core/base/property_base.py b/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.

diff --git a/idaes/core/base/reaction_base.py b/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

diff --git a/idaes/core/scaling/__init__.py b/idaes/core/scaling/__init__.py
@@ -11,7 +11,7 @@
 # for full copyright and license information.
 #################################################################################
 from .autoscaling import AutoScaler
-from .custom_scaler_base import CustomScalerBase
+from .custom_scaler_base import CustomScalerBase, ConstraintScalingScheme
 from .scaler_profiling import ScalingProfiler
 from .util import (
     scaling_factors_from_json_file,

diff --git a/idaes/core/scaling/custom_scaler_base.py b/idaes/core/scaling/custom_scaler_base.py
@@ -571,6 +571,7 @@ def call_submodel_scaler_method(
             submodel: local name of submodel to be scaled as str
             submodel_scalers: user provided ComponentMap of Scalers to use for submodels
             method: name of method to call from submodel (as string)
+            overwrite: whether to overwrite existing scaling factors
 
         Returns:
             None

diff --git a/idaes/core/scaling/scaler_profiling.py b/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,

diff --git a/idaes/core/scaling/tests/test_custom_scaling_integration.py b/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

diff --git a/idaes/models/properties/examples/saponification_reactions.py b/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,84 @@ def define_metadata(cls, obj):
         )
 
 
+class SaponificationReactionScaler(CustomScalerBase):
+    """
+    Scaler for saponification reaction package.
+
+    Variables are scaled by nominal order of magnitude, and constraints
+    using the inverse maximum scheme.
+    """
+
+    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.

diff --git a/idaes/models/properties/examples/saponification_thermo.py b/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,57 @@ def define_metadata(cls, obj):
         )
 
 
+class SaponificationPropertiesScaler(CustomScalerBase):
+    """
+    Scaler for saponification properties package.
+
+    Flow and concentration are scaled by default value (if no user input provided),
+    pressure is scaled assuming order of magnitude of 1e5 Pa, and temperature is
+    scaled using the average of the bounds. Constraints using the inverse maximum
+    scheme.
+    """
+
+    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.