ADM1 Property Package Re-Scaling (#1530)

* Create initial files for ADM1 scalers

* Add new scaler objects

* Update tests

* Isolate iscale functionality to calculate_scaling_factors

* Update testing

* Undo all changes made to adm1 vapor properties

* Try to fix MAC failure

* Add additional testing for adm1 reactions

watertap/flowsheets/anaerobic_digester/modified_ADM1_flowsheet_with_translators.py

@@ -293,6 +293,8 @@ def set_operating_conditions(m, bio_P=False):
     m.fs.AD.volume_vapor.fix(300)
     m.fs.AD.liquid_outlet.temperature.fix(308.15)
 
+    iscale.calculate_scaling_factors(m.fs)
+
     def scale_variables(m):
         for var in m.fs.component_data_objects(pyo.Var, descend_into=True):
             if "flow_vol" in var.name:
@@ -306,7 +308,6 @@ def scale_variables(m):
 
     # Apply scaling
     scale_variables(m)
-    iscale.calculate_scaling_factors(m.fs)
 
 
 def initialize_system(m):

watertap/property_models/unit_specific/anaerobic_digestion/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__ = "Alejandro Garciadiego, Adam Atia, Xinhong Liu"
@@ -247,13 +248,46 @@ def release_state(self, flags, outlvl=idaeslog.NOTSET):
         init_log.info("State Released.")
 
 
+class ADM1PropertiesScaler(CustomScalerBase):
+    """
+    Scaler for the 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)
+
+    # There are currently no constraints in this model
+    def constraint_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+        pass
+
+
 @declare_process_block_class("ADM1StateBlock", block_class=_ADM1StateBlock)
 class ADM1StateBlockData(StateBlockData):
     """
     StateBlock for calculating thermophysical properties associated with the ADM1
     reaction system.
     """
 
+    default_scaler = ADM1PropertiesScaler
+
     def build(self):
         """
         Callable method for Block construction

watertap/property_models/unit_specific/anaerobic_digestion/adm1_reactions.py

@@ -41,6 +41,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__ = "Adam Atia, Alejandro Garciadiego, Xinhong Liu"
@@ -1159,12 +1160,137 @@ def define_metadata(cls, obj):
         )
 
 
+class ADM1ReactionScaler(CustomScalerBase):
+    """
+    Scaler for the 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("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 _ADM1ReactionBlock(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 = ADM1ReactionScaler
+
     def initialize(self, outlvl=idaeslog.NOTSET, **kwargs):
         """
         Initialization routine for reaction package.
@@ -1769,8 +1895,11 @@ def rate_expression_rule(b, r):
             self.del_component(self.rate_expression)
             raise
 
-        for i, c in self.rates.items():
-            iscale.set_scaling_factor(self.reaction_rate[i], 1 / c)
+    def get_reaction_rate_basis(self):
+        return MaterialFlowBasis.mass
+
+    def calculate_scaling_factors(self):
+        super().calculate_scaling_factors()
 
         iscale.set_scaling_factor(self.I, 1e1)
         iscale.set_scaling_factor(self.conc_mass_va, 1e2)
@@ -1787,11 +1916,8 @@ 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.rates.items():
+            iscale.set_scaling_factor(self.reaction_rate[i], 1 / c)
 
         for i, c in self.rate_expression.items():
             iscale.constraint_scaling_transform(