Mixer Scaler (#1676)

* additional test for propagating scaling factors within a block

* Control Volume 1D scaler

* phase material balance cv0d

* test for submodel scaler for blockdata

* Tests for flow direction attribute

* black and typo

* remove unused variables

* protected property error message depends on python version

* test no override

* test propagate scaling factors

* additional test coverage

* Mixer scaler

* test with inherent reactions

* test inlet_blocks

* Pylint

* Apply suggestion from @bpaul4

Co-authored-by: Brandon Paul <[email protected]>

* default_scaler attribute should live on indexed block

* Apply suggestion from @dallan-keylogic

* Apply suggestion from @dallan-keylogic

---------

Co-authored-by: Brandon Paul <[email protected]>

Author: dallan-keylogic

idaes/core/base/control_volume_base.py

@@ -120,7 +120,7 @@ class ControlVolumeScalerBase(CustomScalerBase):
     # terms. Presently (9/25/25), this attribute exists to take into
     # account the fact that all the material and energy terms in the
     # ControlVolume1D are given on the basis of material or energy per
-    # unit length, so we want to weigh them accordingly.
+    # unit length, so we want to weight them accordingly.
     _weight_attr_name = None
 
     def variable_scaling_routine(
@@ -517,6 +517,7 @@ def constraint_scaling_routine(
         else:
             props = getattr(model, self._state_block_ref)
             phase_list = props.phase_list
+            pc_set = props.phase_component_set
 
         if hasattr(model, "reactions"):
             self.call_submodel_scaler_method(
@@ -525,7 +526,7 @@ def constraint_scaling_routine(
                 method="constraint_scaling_routine",
                 overwrite=overwrite,
             )
-        # Transform constraints in order of appearance
+
         if hasattr(model, "material_holdup_calculation"):
             for idx in model.material_holdup_calculation:
                 self.scale_constraint_by_component(
@@ -550,18 +551,34 @@ def constraint_scaling_routine(
                     overwrite=overwrite,
                 )
 
+        inh_rxn_con = None
         if hasattr(model, "inherent_reaction_stoichiometry_constraint"):
-            for idx in model.inherent_reaction_stoichiometry_constraint:
+            # ControlVolume0D and ControlVolume1D
+            inh_rxn_con = model.inherent_reaction_stoichiometry_constraint
+        elif hasattr(model, "inherent_reaction_constraint"):
+            # Mixer
+            inh_rxn_con = model.inherent_reaction_constraint
+
+        if inh_rxn_con is not None:
+            for idx in inh_rxn_con:
                 self.scale_constraint_by_component(
-                    model.inherent_reaction_stoichiometry_constraint[idx],
+                    inh_rxn_con[idx],
                     model.inherent_reaction_generation[idx],
                     overwrite=overwrite,
                 )
 
+        mb_eqn = None
         if hasattr(model, "material_balances"):
+            # ControlVolume0D and ControlVolume1D
+            mb_eqn = model.material_balances
+        elif hasattr(model, "material_mixing_equations"):
+            # Mixer
+            mb_eqn = model.material_mixing_equations
+
+        if mb_eqn is not None:
             mb_type = model._constructed_material_balance_type  # pylint: disable=W0212
             if mb_type == MaterialBalanceType.componentTotal:
-                for idx in model.material_balances:
+                for idx in mb_eqn:
                     c = idx[-1]
                     nom_list = []
                     for p in phase_list:
@@ -572,22 +589,38 @@ def constraint_scaling_routine(
                         )
                     nom = max(nom_list)
                     self.set_component_scaling_factor(
-                        model.material_balances[idx], 1 / nom, overwrite=overwrite
+                        mb_eqn[idx], 1 / nom, overwrite=overwrite
                     )
             elif mb_type == MaterialBalanceType.componentPhase:
-                for idx in model.material_balances:
+                for idx in mb_eqn:
                     p = idx[-2]
                     c = idx[-1]
                     nom = self.get_expression_nominal_value(
                         props[idx[:-2]].get_material_flow_terms(p, c)
                     )
                     self.set_component_scaling_factor(
-                        model.material_balances[idx], 1 / nom, overwrite=overwrite
+                        mb_eqn[idx], 1 / nom, overwrite=overwrite
+                    )
+            elif mb_type == MaterialBalanceType.total:
+                for idx in mb_eqn:
+                    nom_list = []
+                    for p, c in pc_set:
+                        nom_list.append(
+                            self.get_expression_nominal_value(
+                                props[idx[:-1]].get_material_flow_terms(p, c)
+                            )
+                        )
+                    nom = max(nom_list)
+                    self.set_component_scaling_factor(
+                        mb_eqn[idx], 1 / nom, overwrite=overwrite
                     )
             else:
                 # There are some other material balance types but they create
                 # constraints with different names.
-                _log.warning(f"Unknown material balance type {mb_type}")
+                _log.warning(
+                    f"Unknown material balance type {mb_type}. It cannot be "
+                    "automatically scaled."
+                )
 
         # TODO element balances
         # if hasattr(self, "element_balances"):
@@ -602,10 +635,16 @@ def constraint_scaling_routine(
         #         sf = iscale.get_scaling_factor(self.element_holdup[t, e])
         #         iscale.constraint_scaling_transform(c, sf, overwrite=False)
 
+        eb_eqn = None
         if hasattr(model, "enthalpy_balances"):
+            eb_eqn = model.enthalpy_balances
+        elif hasattr(model, "enthalpy_mixing_equations"):
+            eb_eqn = model.enthalpy_mixing_equations
+
+        if eb_eqn is not None:
             # Phase enthalpy balances are not implemented
             # as of 9/26/25
-            for idx in model.enthalpy_balances:
+            for idx in eb_eqn:
                 nom_list = []
                 for p in phase_list:
                     nom_list.append(
@@ -615,7 +654,7 @@ def constraint_scaling_routine(
                     )
                 nom = max(nom_list)
                 self.set_component_scaling_factor(
-                    model.enthalpy_balances[idx], 1 / nom, overwrite=overwrite
+                    eb_eqn[idx], 1 / nom, overwrite=overwrite
                 )
 
         if hasattr(model, "energy_holdup_calculation"):
@@ -626,11 +665,16 @@ def constraint_scaling_routine(
                     overwrite=overwrite,
                 )
 
+        pb_eqn = None
         if hasattr(model, "pressure_balance"):
-            for con in model.pressure_balance.values():
-                self.scale_constraint_by_nominal_value(
+            # ControlVolume0D and ControlVolume1D
+            pb_eqn = model.pressure_balance
+
+        if pb_eqn is not None:
+            for idx, con in pb_eqn.items():
+                self.scale_constraint_by_component(
                     con,
-                    scheme=ConstraintScalingScheme.inverseMaximum,
+                    props[idx].pressure,
                     overwrite=overwrite,
                 )
 

idaes/core/scaling/custom_scaler_base.py

@@ -299,6 +299,7 @@ def scale_variable_by_component(
                 variable=target_variable, scaling_factor=sf, overwrite=overwrite
             )
         else:
+            # TODO add infrastructure to log a warning
             _log.debug(
                 f"Could not set scaling factor for {target_variable.name}, "
                 f"no scaling factor set for {scaling_component.name}"

idaes/models/properties/examples/saponification_thermo.py

@@ -60,7 +60,8 @@ class PhysicalParameterData(PhysicalParameterBlock):
     Property Parameter Block Class
 
     Contains parameters and indexing sets associated with properties for
-    superheated steam.
+    a dilute solution of NaOH, Ethyl Acetate, Sodium Acetate, and Ethanol
+    in water.
 
     """
 

idaes/models/unit_models/mixer.py

@@ -18,6 +18,7 @@
 from pyomo.environ import (
     Block,
     check_optimal_termination,
+    ComponentMap,
     Constraint,
     Param,
     PositiveReals,
@@ -34,6 +35,7 @@
     MaterialBalanceType,
     MaterialFlowBasis,
 )
+from idaes.core.base.control_volume_base import ControlVolumeScalerBase
 from idaes.core.util.config import (
     is_physical_parameter_block,
     is_state_block,
@@ -52,7 +54,7 @@
 
 import idaes.logger as idaeslog
 
-__author__ = "Andrew Lee"
+__author__ = "Andrew Lee, Douglas Allan"
 
 
 # Set up logger
@@ -80,6 +82,96 @@ class MomentumMixingType(Enum):
     minimize_and_equality = 3
 
 
+class MixerScaler(ControlVolumeScalerBase):
+    """
+    Scaler object for the Mixer unit model
+    """
+
+    # This attribute gives the parent ControlVolumeScalerBase
+    # methods a state block with the same index as the material
+    # and energy balances to get scaling information from
+    _state_block_ref = "mixed_state"
+
+    def variable_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None
+    ):
+        for _, iblock in model.inlet_blocks:
+            self.call_submodel_scaler_method(
+                submodel=iblock,
+                submodel_scalers=submodel_scalers,
+                method="variable_scaling_routine",
+                overwrite=overwrite,
+            )
+
+        if model.config.mixed_state_block is None:
+            mblock = model.mixed_state
+        else:
+            mblock = model.config.mixed_state_block
+
+        self.call_submodel_scaler_method(
+            submodel=mblock,
+            submodel_scalers=submodel_scalers,
+            method="variable_scaling_routine",
+            overwrite=overwrite,
+        )
+
+        # Scale inherent reaction and phase equilibrium
+        # variables, if they exist
+        super().variable_scaling_routine(
+            model, overwrite=overwrite, submodel_scalers=submodel_scalers
+        )
+
+        if hasattr(model, "minimum_pressure"):
+            for (t, _), v in model.minimum_pressure.items():
+                self.scale_variable_by_component(
+                    v, model.mixed_state[t].pressure, overwrite=overwrite
+                )
+
+    def constraint_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None
+    ):
+        for _, iblock in model.inlet_blocks:
+            self.call_submodel_scaler_method(
+                submodel=iblock,
+                submodel_scalers=submodel_scalers,
+                method="constraint_scaling_routine",
+                overwrite=overwrite,
+            )
+
+        if model.config.mixed_state_block is None:
+            mblock = model.mixed_state
+        else:
+            mblock = model.config.mixed_state_block
+
+        self.call_submodel_scaler_method(
+            submodel=mblock,
+            submodel_scalers=submodel_scalers,
+            method="constraint_scaling_routine",
+            overwrite=overwrite,
+        )
+
+        # Scale material and energy balance equations
+        super().constraint_scaling_routine(
+            model, overwrite=overwrite, submodel_scalers=submodel_scalers
+        )
+
+        if hasattr(model, "pressure_equality_constraints"):
+            for (t, _), con in model.pressure_equality_constraints.items():
+                self.scale_constraint_by_component(
+                    con, model.mixed_state[t].pressure, overwrite=overwrite
+                )
+        if hasattr(model, "minimum_pressure_constraint"):
+            for (t, _), con in model.minimum_pressure_constraint.items():
+                self.scale_constraint_by_component(
+                    con, model.mixed_state[t].pressure, overwrite=overwrite
+                )
+        if hasattr(model, "mixture_pressure"):
+            for t, con in model.mixture_pressure.items():
+                self.scale_constraint_by_component(
+                    con, model.mixed_state[t].pressure, overwrite=overwrite
+                )
+
+
 class MixerInitializer(ModularInitializerBase):
     """
     Hierarchical Initializer for Mixer blocks.
@@ -113,15 +205,11 @@ def initialization_routine(
         solver = self._get_solver()
 
         # Initialize inlet state blocks
-        inlet_list = model.create_inlet_list()
         i_block_list = []
-        for i in inlet_list:
-            i_block = getattr(model, i + "_state")
-            i_block_list.append(i_block)
 
-            # Get initializer for inlet
+        for _, i_block in model.inlet_blocks:
+            i_block_list.append(i_block)
             iinit = self.get_submodel_initializer(i_block)
-
             iinit.initialize(i_block)
 
         # Initialize mixed state block
@@ -225,6 +313,9 @@ class MixerData(UnitModelBlockData):
     """
 
     default_initializer = MixerInitializer
+    default_scaler = MixerScaler
+
+    _inlet_dict = None
 
     CONFIG = ConfigBlock()
     CONFIG.declare(
@@ -397,6 +488,19 @@ class MixerData(UnitModelBlockData):
         ),
     )
 
+    @property
+    def inlet_blocks(self):
+        """
+        Allows the user to iterate over the inlet stream names and state blocks.
+
+        Returns:
+            dict_items with the inlet stream names as keys and the
+                inlet state blocks as values
+        """
+        # Return an iterator so the user cannot
+        # mutate _inlet_dict
+        return self._inlet_dict.items()
+
     def build(self):
         """
         General build method for MixerData. This method calls a number
@@ -424,6 +528,8 @@ def build(self):
         # Build StateBlocks
         inlet_blocks = self.add_inlet_state_blocks(inlet_list)
 
+        self._inlet_dict = {key: blk for key, blk in zip(inlet_list, inlet_blocks)}
+
         if self.config.mixed_state_block is None:
             mixed_block = self.add_mixed_state_block()
         else:
@@ -635,6 +741,7 @@ def add_material_mixing_equations(self, inlet_blocks, mixed_block, mb_type):
         else:
             # Let this pass for now with no units
             flow_units = None
+        self._constructed_material_balance_type = mb_type
 
         if mixed_block.include_inherent_reactions:
             if mb_type == MaterialBalanceType.total: