ADM1 Property Package Re-Scaling

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_properties_vapor.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, ConstraintScalingScheme
 
 # Some more information about this module
 __author__ = "Alejandro Garciadiego, Xinhong Liu"
@@ -217,13 +218,53 @@ def release_state(self, flags, outlvl=idaeslog.NOTSET):
         init_log.info("State Released.")
 
 
+class ADM1_vaporPropertiesScaler(CustomScalerBase):
+    """
+    Scaler for the vapor 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-3),
+    }
+
+    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
+    ):
+        # TODO: Revisit this scaling methodologies
+        # Consider other schemes, scale_constraint_by_default, or scale_constraints_by_jacobian_norm
+        if model.is_property_constructed("pressure_sat"):
+            for j in model._pressure_sat.values():
+                self.scale_constraint_by_nominal_value(
+                    j,
+                    scheme=ConstraintScalingScheme.inverseMaximum,
+                    overwrite=overwrite,
+                )
+
+
 @declare_process_block_class("ADM1_vaporStateBlock", block_class=_ADM1_vaporStateBlock)
 class ADM1_vaporStateBlockData(StateBlockData):
     """
     StateBlock for calculating thermophysical properties associated with the ADM1
     reaction system.
     """
 
+    default_scaler = ADM1_vaporPropertiesScaler
+
     def build(self):
         """
         Callable method for Block construction
@@ -413,8 +454,6 @@ def calculate_scaling_factors(self):
         # Get default scale factors and do calculations from base classes
         super().calculate_scaling_factors()
 
-        # No constraints in this model as yet, just need to set scaling factors
-        # for expressions
         sf_F = iscale.get_scaling_factor(self.flow_vol, default=1e2, warning=True)
         sf_T = iscale.get_scaling_factor(self.temperature, default=1e-2, warning=True)
 

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,9 +1895,6 @@ 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)
-
         iscale.set_scaling_factor(self.I, 1e1)
         iscale.set_scaling_factor(self.conc_mass_va, 1e2)
         iscale.set_scaling_factor(self.conc_mass_bu, 1e2)
@@ -1793,6 +1916,9 @@ def get_reaction_rate_basis(self):
     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(
                 c,