Modified ADM1 Property Package Re-Scaling (#1537)

* Creat modified ADM1 scalers and tests

* Address pylint issue

* Add scaling factors to BSM2_P flowsheet

* Update modified adm1 reaction test

---------

Co-authored-by: Keith Beattie <[email protected]>

watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension.py

@@ -536,6 +536,10 @@ def scale_variables(m):
                 iscale.set_scaling_factor(var, 1e-5)
             if "conc_mass_comp" in var.name:
                 iscale.set_scaling_factor(var, 1e1)
+            if "anions" in var.name:
+                iscale.set_scaling_factor(var, 1e2)
+            if "cations" in var.name:
+                iscale.set_scaling_factor(var, 1e2)
 
     for unit in ("R1", "R2", "R3", "R4"):
         block = getattr(m.fs, unit)
@@ -972,7 +976,7 @@ def display_performance_metrics(m):
 
 
 if __name__ == "__main__":
-    m, results = main(bio_P=True)
+    m, results = main(bio_P=False)
 
     stream_table = create_stream_table_dataframe(
         {

watertap/property_models/unit_specific/anaerobic_digestion/modified_adm1_properties.py

@@ -34,6 +34,7 @@
 from idaes.core.util.initialization import fix_state_vars, revert_state_vars
 import idaes.logger as idaeslog
 import idaes.core.util.scaling as iscale
+from idaes.core.scaling import CustomScalerBase
 
 # Some more information about this module
 __author__ = "Chenyu Wang, Marcus Holly, Adam Atia, Xinhong Liu"
@@ -309,6 +310,36 @@ def release_state(self, flags, outlvl=idaeslog.NOTSET):
         init_log.info("State Released.")
 
 
+class ModifiedADM1PropertiesScaler(CustomScalerBase):
+    """
+    Scaler for the Modified Anaerobic Digestion Model No.1 property package.
+    Flow and temperature are scaled by the default value (if no user input provided), and
+    pressure is scaled assuming an order of magnitude of 1e5 Pa.
+    """
+
+    UNIT_SCALING_FACTORS = {
+        # "QuantityName: (reference units, scaling factor)
+        "Pressure": (pyo.units.Pa, 1e-6),
+    }
+
+    DEFAULT_SCALING_FACTORS = {
+        "flow_vol": 1e5,
+        "temperature": 1e-1,
+    }
+
+    def variable_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+        self.scale_variable_by_default(model.temperature, overwrite=overwrite)
+        self.scale_variable_by_default(model.flow_vol, overwrite=overwrite)
+        self.scale_variable_by_units(model.pressure, overwrite=overwrite)
+
+    def constraint_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+        pass
+
+
 @declare_process_block_class(
     "ModifiedADM1StateBlock", block_class=_ModifiedADM1StateBlock
 )
@@ -318,6 +349,8 @@ class ModifiedADM1StateBlockData(StateBlockData):
     reaction system.
     """
 
+    default_scaler = ModifiedADM1PropertiesScaler
+
     def build(self):
         """
         Callable method for Block construction
@@ -459,13 +492,6 @@ def energy_density_expression(self):
             rule=energy_density_expression, doc="Energy density term"
         )
 
-        iscale.set_scaling_factor(self.flow_vol, 1e5)
-        iscale.set_scaling_factor(self.temperature, 1e-1)
-        iscale.set_scaling_factor(self.pressure, 1e-6)
-        iscale.set_scaling_factor(self.conc_mass_comp, 1e2)
-        iscale.set_scaling_factor(self.anions, 1e2)
-        iscale.set_scaling_factor(self.cations, 1e2)
-
     # On-demand properties
     def _VSS(self):
         self.VSS = pyo.Var(

watertap/property_models/unit_specific/anaerobic_digestion/modified_adm1_reactions.py

@@ -38,6 +38,7 @@
 import idaes.logger as idaeslog
 import idaes.core.util.scaling as iscale
 from idaes.core.util.math import smooth_max
+from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme
 
 # Some more information about this module
 __author__ = "Chenyu Wang, Marcus Holly, Xinhong Liu"
@@ -1754,12 +1755,148 @@ def define_metadata(cls, obj):
         )
 
 
+class ModifiedADM1ReactionScaler(CustomScalerBase):
+    """
+    Scaler for the Modified Anaerobic Digestion Model No.1 reaction package.
+    Variables are scaled by their default scaling factor (if no user input provided), and constraints
+    are scaled using the inverse maximum scheme.
+    """
+
+    # TODO: Revisit this scaling factor
+    DEFAULT_SCALING_FACTORS = {
+        "reaction_rate": 1e2,
+        "I": 1e1,
+    }
+
+    def variable_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+        for r in model.params.rate_reaction_idx:
+            self.scale_variable_by_default(model.I[r], 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
+    ):
+        # TODO: Revisit these scaling methodologies
+        # Consider other schemes, scale_constraint_by_default, or scale_constraints_by_jacobian_norm
+        if model.is_property_constructed("rate_expression"):
+            for j in model.rate_expression.values():
+                self.scale_constraint_by_nominal_value(
+                    j,
+                    scheme=ConstraintScalingScheme.inverseMaximum,
+                    overwrite=overwrite,
+                )
+        if model.is_property_constructed("Dissociation"):
+            self.scale_constraint_by_nominal_value(
+                model.Dissociation,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("CO2_acid_base_equilibrium"):
+            self.scale_constraint_by_nominal_value(
+                model.CO2_acid_base_equilibrium,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("IN_acid_base_equilibrium"):
+            self.scale_constraint_by_nominal_value(
+                model.IN_acid_base_equilibrium,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("pH_calc"):
+            self.scale_constraint_by_nominal_value(
+                model.pH_calc,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("concentration_of_va"):
+            self.scale_constraint_by_nominal_value(
+                model.concentration_of_va,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("concentration_of_bu"):
+            self.scale_constraint_by_nominal_value(
+                model.concentration_of_bu,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("concentration_of_pro"):
+            self.scale_constraint_by_nominal_value(
+                model.concentration_of_pro,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("concentration_of_ac"):
+            self.scale_constraint_by_nominal_value(
+                model.concentration_of_ac,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("concentration_of_hco3"):
+            self.scale_constraint_by_nominal_value(
+                model.concentration_of_hco3,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("concentration_of_nh3"):
+            self.scale_constraint_by_nominal_value(
+                model.concentration_of_nh3,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("concentration_of_co2"):
+            self.scale_constraint_by_nominal_value(
+                model.concentration_of_co2,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("concentration_of_nh4"):
+            self.scale_constraint_by_nominal_value(
+                model.concentration_of_nh4,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("concentration_of_Mg"):
+            self.scale_constraint_by_nominal_value(
+                model.concentration_of_Mg,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("concentration_of_K"):
+            self.scale_constraint_by_nominal_value(
+                model.concentration_of_K,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("S_H_cons"):
+            self.scale_constraint_by_nominal_value(
+                model.S_H_cons,
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+        if model.is_property_constructed("I_fun"):
+            for r in model.params.rate_reaction_idx:
+                self.scale_constraint_by_nominal_value(
+                    model.I_fun[r],
+                    scheme=ConstraintScalingScheme.inverseMaximum,
+                    overwrite=overwrite,
+                )
+
+
 class _ModifiedADM1ReactionBlock(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 = ModifiedADM1ReactionScaler
+
     def initialize(self, outlvl=idaeslog.NOTSET, **kwargs):
         """
         Initialization routine for reaction package.
@@ -2577,6 +2714,12 @@ def rate_expression_rule(b, r):
             self.del_component(self.rate_expression)
             raise
 
+    def get_reaction_rate_basis(self):
+        return MaterialFlowBasis.mass
+
+    def calculate_scaling_factors(self):
+        super().calculate_scaling_factors()
+
         for i, c in self.rates.items():
             iscale.set_scaling_factor(self.reaction_rate[i], 1 / c)
 
@@ -2597,12 +2740,6 @@ def rate_expression_rule(b, r):
         iscale.set_scaling_factor(self.pK_a_IN, 1e0)
         iscale.set_scaling_factor(self.pH, 1e0)
 
-    def get_reaction_rate_basis(self):
-        return MaterialFlowBasis.mass
-
-    def calculate_scaling_factors(self):
-        super().calculate_scaling_factors()
-
         for i, c in self.rate_expression.items():
             iscale.constraint_scaling_transform(
                 c,