Anaerobic Digester Re-Scaling

watertap/flowsheets/full_water_resource_recovery_facility/BSM2.py

@@ -10,7 +10,7 @@
 # "https://github.com/watertap-org/watertap/"
 #################################################################################
 """
-Flowsheet example full Water Resource Recovery Facility 
+Flowsheet example full Water Resource Recovery Facility
 (WRRF; a.k.a., wastewater treatment plant) with ASM1 and ADM1.
 
 The flowsheet follows the same formulation as benchmark simulation model no.2 (BSM2)

watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension.py

@@ -88,6 +88,11 @@
     cost_circular_clarifier,
     cost_primary_clarifier,
 )
+from idaes.core.scaling.custom_scaler_base import (
+    CustomScalerBase,
+    ConstraintScalingScheme,
+)
+
 
 # Set up logger
 _log = idaeslog.getLogger(__name__)
@@ -373,11 +378,6 @@ def build(bio_P=False):
         doc="Dissolved oxygen concentration at equilibrium",
     )
 
-    m.fs.aerobic_reactors = (m.fs.R5, m.fs.R6, m.fs.R7)
-    for R in m.fs.aerobic_reactors:
-        iscale.set_scaling_factor(R.KLa, 1e-2)
-        iscale.set_scaling_factor(R.hydraulic_retention_time[0], 1e-3)
-
     @m.fs.R5.Constraint(m.fs.time, doc="Mass transfer constraint for R3")
     def mass_transfer_R5(self, t):
         return pyo.units.convert(
@@ -526,53 +526,55 @@ def set_operating_conditions(m, bio_P=False):
     m.fs.thickener.hydraulic_retention_time.fix(86400 * pyo.units.s)
     m.fs.thickener.diameter.fix(10 * pyo.units.m)
 
+    scaler = CustomScalerBase()
+
     def scale_variables(m):
         for var in m.fs.component_data_objects(pyo.Var, descend_into=True):
             if "flow_vol" in var.name:
-                iscale.set_scaling_factor(var, 1e0)
+                iscale.set_scaling_factor(var, 1e3)
             if "temperature" in var.name:
                 iscale.set_scaling_factor(var, 1e-2)
             if "pressure" in var.name:
                 iscale.set_scaling_factor(var, 1e-5)
             if "conc_mass_comp" in var.name:
-                iscale.set_scaling_factor(var, 1e1)
+                if bio_P:
+                    iscale.set_scaling_factor(var, 1e1)
+                else:
+                    iscale.set_scaling_factor(var, 1e2)
             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)
-        iscale.set_scaling_factor(block.hydraulic_retention_time, 1e-3)
+    def scale_constraints(m):
+        for c in m.fs.component_data_objects(pyo.Constraint, descend_into=True):
+            if "flow_vol_equality" in c.name:
+                scaler.scale_constraint_by_nominal_value(
+                    c,
+                    scheme=ConstraintScalingScheme.inverseMaximum,
+                    overwrite=True,
+                )
 
-    for unit in ("R1", "R2", "R3", "R4", "R5", "R6", "R7"):
-        block = getattr(m.fs, unit)
-        iscale.set_scaling_factor(
-            block.control_volume.reactions[0.0].rate_expression, 1e3
-        )
-        iscale.set_scaling_factor(block.cstr_performance_eqn, 1e3)
-        iscale.set_scaling_factor(
-            block.control_volume.rate_reaction_stoichiometry_constraint, 1e3
-        )
-        iscale.set_scaling_factor(block.control_volume.material_balances, 1e3)
-
-    iscale.set_scaling_factor(m.fs.AD.KH_co2, 1e1)
-    iscale.set_scaling_factor(m.fs.AD.KH_ch4, 1e1)
-    iscale.set_scaling_factor(m.fs.AD.KH_h2, 1e2)
+    m.fs.aerobic_reactors = (m.fs.R5, m.fs.R6, m.fs.R7)
+    for R in m.fs.aerobic_reactors:
+        iscale.set_scaling_factor(R.KLa, 1e-2)
+        iscale.set_scaling_factor(R.hydraulic_retention_time[0], 1e-3)
 
     if bio_P:
         iscale.set_scaling_factor(m.fs.AD.liquid_phase.heat, 1e3)
-        iscale.constraint_scaling_transform(
-            m.fs.AD.liquid_phase.enthalpy_balances[0], 1e-6
-        )
     else:
-        iscale.set_scaling_factor(m.fs.AD.liquid_phase.heat, 1e2)
-        iscale.constraint_scaling_transform(
-            m.fs.AD.liquid_phase.enthalpy_balances[0], 1e-3
+        iscale.set_scaling_factor(m.fs.AD.liquid_phase.heat, 1e4)
+        scaler.scale_constraint_by_nominal_value(
+            m.fs.AD.liquid_phase.enthalpy_balances[0],
+            scheme=ConstraintScalingScheme.inverseMinimum,
+            overwrite=True,
         )
+        iscale.set_scaling_factor(m.fs.AD.hydraulic_retention_time[0], 1e-6)
+        iscale.constraint_scaling_transform(m.fs.AD.AD_retention_time[0], 1e-4)
 
     # Apply scaling
     scale_variables(m)
+    scale_constraints(m)
     iscale.calculate_scaling_factors(m)
 
 

watertap/flowsheets/full_water_resource_recovery_facility/test/test_BSM2_P_extension.py

@@ -56,71 +56,71 @@ def test_solve(self, system_frame):
             0.24219, rel=1e-3
         )
         assert value(m.fs.Treated.properties[0].conc_mass_comp["S_A"]) == pytest.approx(
-            6.4300e-07, abs=1e-6
+            6.43000895e-07, abs=1e-6
         )
         assert value(m.fs.Treated.properties[0].conc_mass_comp["S_F"]) == pytest.approx(
-            0.00027610, rel=1e-3
+            0.0002760866, rel=1e-3
         )
         assert value(m.fs.Treated.properties[0].conc_mass_comp["S_I"]) == pytest.approx(
-            0.057450, rel=1e-3
+            0.057450006, rel=1e-3
         )
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["S_N2"]
-        ) == pytest.approx(0.052060, rel=1e-3)
+        ) == pytest.approx(0.0520459, rel=1e-3)
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["S_NH4"]
-        ) == pytest.approx(0.00017686, rel=1e-3)
+        ) == pytest.approx(0.0001769519, rel=1e-3)
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["S_NO3"]
-        ) == pytest.approx(0.0060483, rel=1e-3)
+        ) == pytest.approx(0.006044579, rel=1e-3)
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["S_O2"]
-        ) == pytest.approx(0.0076767, rel=1e-3)
+        ) == pytest.approx(0.00767666, rel=1e-3)
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["S_PO4"]
-        ) == pytest.approx(0.65267, rel=1e-3)
+        ) == pytest.approx(0.6527876, rel=1e-3)
         assert value(m.fs.Treated.properties[0].conc_mass_comp["S_K"]) == pytest.approx(
-            0.36810, rel=1e-3
+            0.36809549, rel=1e-3
         )
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["S_Mg"]
-        ) == pytest.approx(0.018654, rel=1e-3)
+        ) == pytest.approx(0.01865257, rel=1e-3)
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["S_IC"]
-        ) == pytest.approx(0.15148, rel=1e-3)
+        ) == pytest.approx(0.1507046, rel=1e-3)
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["X_AUT"]
-        ) == pytest.approx(0.00041597, rel=1e-3)
+        ) == pytest.approx(0.000415746, rel=1e-3)
         assert value(m.fs.Treated.properties[0].conc_mass_comp["X_H"]) == pytest.approx(
-            0.013371, rel=1e-3
+            0.0133707, rel=1e-3
         )
         assert value(m.fs.Treated.properties[0].conc_mass_comp["X_I"]) == pytest.approx(
-            0.012360, rel=1e-3
+            0.01235979, rel=1e-3
         )
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["X_PAO"]
-        ) == pytest.approx(0.011077, rel=1e-3)
+        ) == pytest.approx(0.0110803, rel=1e-3)
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["X_PHA"]
-        ) == pytest.approx(5.0389e-06, abs=1e-6)
+        ) == pytest.approx(5.0389213e-06, abs=1e-6)
         assert value(
             m.fs.Treated.properties[0].conc_mass_comp["X_PP"]
-        ) == pytest.approx(0.0036998, rel=1e-3)
+        ) == pytest.approx(0.00370076, rel=1e-3)
         assert value(m.fs.Treated.properties[0].conc_mass_comp["X_S"]) == pytest.approx(
-            0.00021570, rel=1e-3
+            0.0002157186, rel=1e-3
         )
 
     @pytest.mark.component
     def test_costing(self, system_frame):
         m = system_frame
 
         # check costing
-        assert value(m.fs.costing.LCOW) == pytest.approx(0.468069, rel=1e-3)
+        assert value(m.fs.costing.LCOW) == pytest.approx(0.468069088, rel=1e-3)
         assert value(m.fs.costing.total_capital_cost) == pytest.approx(
-            23935727.742, rel=1e-3
+            23935727.743, rel=1e-3
         )
         assert value(m.fs.costing.total_operating_cost) == pytest.approx(
-            827635.25, rel=1e-3
+            827635.247, rel=1e-3
         )
 
 

watertap/property_models/unit_specific/anaerobic_digestion/tests/test_modified_adm1_reaction.py

@@ -940,6 +940,9 @@ def model(self):
         m.fs.unit.liquid_phase.properties_in[0].TSS
         m.fs.unit.liquid_phase.properties_out[0].TSS
 
+        iscale.calculate_scaling_factors(m)
+        iscale.set_scaling_factor(m.fs.unit.liquid_phase.heat, 1e-6)
+
         return m
 
     @pytest.mark.unit

watertap/unit_models/anaerobic_digester.py

@@ -67,18 +67,125 @@
 from idaes.core.util.constants import Constants
 from idaes.core.util.exceptions import ConfigurationError, InitializationError
 from idaes.core.util.tables import create_stream_table_dataframe
+from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme
 
 from watertap.costing.unit_models.anaerobic_digester import cost_anaerobic_digester
 
 __author__ = "Alejandro Garciadiego, Andrew Lee, Xinhong Liu"
 
 
+class ADScaler(CustomScalerBase):
+    """
+    Default modular scaler for anaerobic digester.
+
+    This Scaler relies on the associated property and reaction packages,
+    either through user provided options (submodel_scalers argument) or by default
+    Scalers assigned to the packages.
+    """
+
+    DEFAULT_SCALING_FACTORS = {
+        "volume": 1e-2,
+    }
+
+    def variable_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+        """
+        Routine to apply scaling factors to variables in model.
+
+        Args:
+            model: model to be scaled
+            overwrite: whether to overwrite existing scaling factors
+            submodel_scalers: dict of Scalers to use for sub-models, keyed by submodel local name
+
+        Returns:
+            None
+        """
+        # Call scaling methods for sub-models
+        self.call_submodel_scaler_method(
+            submodel=model.liquid_phase.properties_in,
+            method="variable_scaling_routine",
+            submodel_scalers=submodel_scalers,
+            overwrite=overwrite,
+        )
+        self.propagate_state_scaling(
+            target_state=model.liquid_phase.properties_out,
+            source_state=model.liquid_phase.properties_in,
+            overwrite=overwrite,
+        )
+
+        self.call_submodel_scaler_method(
+            submodel=model.liquid_phase.properties_out,
+            method="variable_scaling_routine",
+            submodel_scalers=submodel_scalers,
+            overwrite=overwrite,
+        )
+        self.call_submodel_scaler_method(
+            submodel=model.liquid_phase.reactions,
+            method="variable_scaling_routine",
+            submodel_scalers=submodel_scalers,
+            overwrite=overwrite,
+        )
+
+        # Scaling control volume variables
+        self.scale_variable_by_default(
+            model.liquid_phase.volume[0], overwrite=overwrite
+        )
+
+    def constraint_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+        """
+        Routine to apply scaling factors to constraints in model.
+
+        Submodel Scalers are called for the property and reaction blocks. All other constraints
+        are scaled using the inverse maximum scheme.
+
+        Args:
+            model: model to be scaled
+            overwrite: whether to overwrite existing scaling factors
+            submodel_scalers: dict of Scalers to use for sub-models, keyed by submodel local name
+
+        Returns:
+            None
+        """
+        # Call scaling methods for sub-models
+        self.call_submodel_scaler_method(
+            submodel=model.liquid_phase.properties_in,
+            method="constraint_scaling_routine",
+            submodel_scalers=submodel_scalers,
+            overwrite=overwrite,
+        )
+        self.call_submodel_scaler_method(
+            submodel=model.liquid_phase.properties_out,
+            method="constraint_scaling_routine",
+            submodel_scalers=submodel_scalers,
+            overwrite=overwrite,
+        )
+        self.call_submodel_scaler_method(
+            submodel=model.liquid_phase.reactions,
+            method="constraint_scaling_routine",
+            submodel_scalers=submodel_scalers,
+            overwrite=overwrite,
+        )
+
+        # Scale unit level constraints
+        if hasattr(model, "AD_retention_time"):
+            self.scale_constraint_by_nominal_value(
+                model.AD_retention_time[0],
+                scheme=ConstraintScalingScheme.inverseMaximum,
+                overwrite=overwrite,
+            )
+
+
 @declare_process_block_class("AD")
 class ADData(UnitModelBlockData):
     """
     AD Unit Model Class
     """
 
+    default_scaler = ADScaler
+
     CONFIG = UnitModelBlockData.CONFIG()
 
     CONFIG.declare(
@@ -366,7 +473,7 @@ def build(self):
         )
 
         # ---------------------------------------------------------------------
-        # Check flow basis is compatable
+        # Check flow basis is compatible
         # TODO : Could add code to convert flow bases, but not now
         t_init = self.flowsheet().time.first()
         if (
@@ -792,22 +899,6 @@ def rule_electricity_consumption(self, t):
         # Set references to balance terms at unit level
         self.heat_duty = Reference(self.liquid_phase.heat[:])
 
-        iscale.set_scaling_factor(self.KH_co2, 1e2)
-        iscale.set_scaling_factor(self.KH_ch4, 1e2)
-        iscale.set_scaling_factor(self.KH_h2, 1e2)
-        iscale.set_scaling_factor(self.hydraulic_retention_time, 1e-6)
-        iscale.set_scaling_factor(self.volume_AD, 1e-2)
-        iscale.set_scaling_factor(self.volume_vapor, 1e-2)
-        iscale.set_scaling_factor(self.liquid_phase.rate_reaction_generation, 1e4)
-        iscale.set_scaling_factor(self.liquid_phase.mass_transfer_term, 1e2)
-        iscale.set_scaling_factor(self.liquid_phase.heat, 1e0)
-        iscale.set_scaling_factor(self.liquid_phase.rate_reaction_extent, 1e4)
-        iscale.set_scaling_factor(self.liquid_phase.enthalpy_transfer, 1e0)
-        iscale.set_scaling_factor(self.liquid_phase.volume, 1e-2)
-        iscale.set_scaling_factor(self.electricity_consumption, 1e0)
-        for i, c in self.ad_performance_eqn.items():
-            iscale.constraint_scaling_transform(c, 1e2)
-
     def _get_stream_table_contents(self, time_point=0):
         return create_stream_table_dataframe(
             {
@@ -836,6 +927,21 @@ def calculate_scaling_factors(self):
             & self.liquid_phase.properties_out.component_list
         )
 
+        iscale.set_scaling_factor(self.KH_co2, 1e2)
+        iscale.set_scaling_factor(self.KH_ch4, 1e2)
+        iscale.set_scaling_factor(self.KH_h2, 1e2)
+        iscale.set_scaling_factor(self.hydraulic_retention_time, 1e-6)
+        iscale.set_scaling_factor(self.volume_AD, 1e-2)
+        iscale.set_scaling_factor(self.volume_vapor, 1e-2)
+        iscale.set_scaling_factor(self.liquid_phase.rate_reaction_generation, 1e4)
+        iscale.set_scaling_factor(self.liquid_phase.mass_transfer_term, 1e2)
+        iscale.set_scaling_factor(self.liquid_phase.rate_reaction_extent, 1e4)
+        iscale.set_scaling_factor(self.liquid_phase.enthalpy_transfer, 1e0)
+        iscale.set_scaling_factor(self.liquid_phase.volume, 1e-2)
+        iscale.set_scaling_factor(self.electricity_consumption, 1e0)
+        for i, c in self.ad_performance_eqn.items():
+            iscale.constraint_scaling_transform(c, 1e2)
+
         # TODO: improve this later; for now, this resolved some scaling issues for modified adm1 test file
         if "S_IP" in self.config.liquid_property_package.component_list:
             sf = iscale.get_scaling_factor(