ASM1 Property Package Re-Scaling (#1519)

* Update property packages

* Add scaler testing to asm1_thermo

* Undo changes to adm1 properties

* Add iscale functionality back to ASM1 props

* Add scalers for asm1 properties and reactions

* Minor updates to BSM2 flowsheet

* Update testing for BSM2 flowsheet

* Address pylint issue

* Update BSM2 Jupyter Notebook

* Update test file

* Minor updates to BSM2 Jupyter Notebook

* Revert small change in asm1_properties

* Update asm1_thermo test

* Correct typo

* Have new scaling factors match old scaling factors

* Isolate iscale functionality to calculate_scaling_factors

* Update iscale in unit model testing

* Try resolving BSM2 issue

* BSM2 fix attempt #2

* Update constraint scaling in BSM2

* UI check passing locally - trying to unfix scaling factors

* Call calculate scaling factors in BSM2 costing as well

* Clean up flowsheet

* Try reverting some BSM2 changes

* Add manual scaling to BSM2 flowsheet

* Minor cleanup

* Clean up property packages

* Update scaling in unit model tests

* Update scaling in tests

* Revert changes to BSM2 tutorial

* Add more costing scaling to BSM2

* Update BSM2 tutorial

* Address pylint

* Update costing scaling factors again

* Update costing scaling factors

* Try to resolve MAC failure

---------

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

watertap/flowsheets/activated_sludge/ASM1_flowsheet.py

@@ -246,6 +246,15 @@ def mass_transfer_R4(self, t):
     assert degrees_of_freedom(m) == 0
 
     # Apply scaling
+    for var in m.fs.component_data_objects(pyo.Var, descend_into=True):
+        if "flow_vol" in var.name:
+            iscale.set_scaling_factor(var, 1e1)
+        if "temperature" in var.name:
+            iscale.set_scaling_factor(var, 1e-1)
+        if "pressure" in var.name:
+            iscale.set_scaling_factor(var, 1e-6)
+        if "conc_mass_comp" in var.name:
+            iscale.set_scaling_factor(var, 1e1)
     iscale.calculate_scaling_factors(m.fs)
 
     # Initialize flowsheet

watertap/flowsheets/anaerobic_digester/ADM1_flowsheet.py

@@ -16,6 +16,7 @@
     units,
 )
 from idaes.core import FlowsheetBlock
+import idaes.core.util.scaling as iscale
 from watertap.core.solvers import get_solver
 import idaes.logger as idaeslog
 from watertap.unit_models.anaerobic_digester import AD
@@ -89,6 +90,8 @@ def build_flowsheet():
 
     m.fs.R1.liquid_outlet.temperature.fix(308.15 * pyo.units.K)
 
+    iscale.calculate_scaling_factors(m)
+
     # TO DO: Fix initialization
     m.fs.R1.initialize(outlvl=idaeslog.INFO_HIGH)
 

watertap/flowsheets/full_water_resource_recovery_facility/BSM2.py

@@ -46,7 +46,7 @@
     ADM1_vaporParameterBlock,
 )
 
-from idaes.core import FlowsheetBlock, UnitModelCostingBlock
+from idaes.core import FlowsheetBlock, UnitModelCostingBlock, UnitModelBlockData
 from idaes.models.unit_models import (
     Feed,
     Mixer,
@@ -277,8 +277,6 @@ def build():
 
     pyo.TransformationFactory("network.expand_arcs").apply_to(m)
 
-    iscale.calculate_scaling_factors(m.fs)
-
     # keep handy all the mixers
     m.mixers = (m.fs.MX1, m.fs.MX2, m.fs.MX3, m.fs.MX4, m.fs.MX6)
 
@@ -387,6 +385,18 @@ def set_operating_conditions(m):
     for mx in m.mixers:
         mx.pressure_equality_constraints[0.0, 2].deactivate()
 
+    for var in m.fs.component_data_objects(pyo.Var, descend_into=True):
+        if "flow_vol" in var.name:
+            iscale.set_scaling_factor(var, 1e1)
+        if "temperature" in var.name:
+            iscale.set_scaling_factor(var, 1e-1)
+        if "pressure" in var.name:
+            iscale.set_scaling_factor(var, 1e-6)
+        if "conc_mass_comp" in var.name:
+            iscale.set_scaling_factor(var, 1e1)
+
+    iscale.calculate_scaling_factors(m)
+
 
 def initialize_system(m):
     # Initialize flowsheet
@@ -494,10 +504,14 @@ def add_costing(m):
 
     m.fs.objective = pyo.Objective(expr=m.fs.costing.LCOW)
     iscale.set_scaling_factor(m.fs.costing.LCOW, 1e3)
-    iscale.set_scaling_factor(m.fs.costing.total_capital_cost, 1e-7)
     iscale.set_scaling_factor(m.fs.costing.total_capital_cost, 1e-5)
 
-    iscale.calculate_scaling_factors(m.fs)
+    for block in m.fs.component_objects(pyo.Block, descend_into=True):
+        if isinstance(block, UnitModelBlockData) and hasattr(block, "costing"):
+            iscale.set_scaling_factor(block.costing.capital_cost, 1e-6)
+
+    iscale.constraint_scaling_transform(m.fs.DU.costing.capital_cost_constraint, 1e-6)
+    iscale.constraint_scaling_transform(m.fs.RADM.costing.capital_cost_constraint, 1e-6)
 
 
 def setup_optimization(m, reactor_volume_equalities=False):

watertap/flowsheets/full_water_resource_recovery_facility/test/test_full_WRRF_with_ASM1_ADM1.py

@@ -143,12 +143,12 @@ def test_optimization(self, system_frame):
         assert degrees_of_freedom(system_frame) == 10
 
         # check costing
-        assert value(m.fs.costing.LCOW) == pytest.approx(0.3497531, rel=1e-5)
+        assert value(m.fs.costing.LCOW) == pytest.approx(0.349772203, rel=1e-5)
         assert value(m.fs.costing.total_capital_cost) == pytest.approx(
-            17441736.89749642, rel=1e-5
+            17379540.339857, rel=1e-5
         )
         assert value(m.fs.costing.total_operating_cost) == pytest.approx(
-            629780.1104274583, rel=1e-5
+            636129.6209807, rel=1e-5
         )
 
 
@@ -258,10 +258,10 @@ def test_optimization(self, system_frame):
         assert degrees_of_freedom(system_frame) == 8
 
         # check costing
-        assert value(m.fs.costing.LCOW) == pytest.approx(0.3497531, rel=1e-5)
+        assert value(m.fs.costing.LCOW) == pytest.approx(0.349560273, rel=1e-5)
         assert value(m.fs.costing.total_capital_cost) == pytest.approx(
-            17441740.61915915, rel=1e-5
+            17370674.42102, rel=1e-5
         )
         assert value(m.fs.costing.total_operating_cost) == pytest.approx(
-            629779.9546967598, rel=1e-5
+            635577.7320509, rel=1e-5
         )

watertap/property_models/unit_specific/activated_sludge/asm1_properties.py

@@ -20,7 +20,6 @@
 from idaes.core import (
     declare_process_block_class,
     MaterialFlowBasis,
-    PhysicalParameterBlock,
     StateBlockData,
     StateBlock,
     MaterialBalanceType,
@@ -32,7 +31,9 @@
 )
 from idaes.core.util.model_statistics import degrees_of_freedom
 from idaes.core.util.initialization import fix_state_vars, revert_state_vars
+from idaes.core.scaling import CustomScalerBase
 import idaes.logger as idaeslog
+from idaes.core.base.property_base import PhysicalParameterBlock
 import idaes.core.util.scaling as iscale
 
 # Some more information about this module
@@ -191,12 +192,46 @@ def define_metadata(cls, obj):
         )
 
 
+class ASM1PropertiesScaler(CustomScalerBase):
+    """
+    Scaler for the Activated Sludge 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": 1e1,
+        "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
+
+
 class _ASM1StateBlock(StateBlock):
     """
     This Class contains methods which should be applied to Property Blocks as a
     whole, rather than individual elements of indexed Property Blocks.
     """
 
+    default_scaler = ASM1PropertiesScaler
+
     def initialize(
         self,
         state_args=None,
@@ -292,7 +327,7 @@ def release_state(self, flags, outlvl=idaeslog.NOTSET):
 @declare_process_block_class("ASM1StateBlock", block_class=_ASM1StateBlock)
 class ASM1StateBlockData(StateBlockData):
     """
-    StateBlock for calculating thermophysical proeprties associated with the ASM1
+    StateBlock for calculating thermophysical properties associated with the ASM1
     reaction system.
     """
 
@@ -306,7 +341,7 @@ def build(self):
         self.flow_vol = pyo.Var(
             initialize=1.0,
             domain=pyo.NonNegativeReals,
-            doc="Total volumentric flowrate",
+            doc="Total volumetric flowrate",
             units=pyo.units.m**3 / pyo.units.s,
         )
         self.pressure = pyo.Var(
@@ -470,11 +505,6 @@ def _Total_N(self):
             doc="Total Nitrogen",
         )
 
-        iscale.set_scaling_factor(self.flow_vol, 1e1)
-        iscale.set_scaling_factor(self.temperature, 1e-1)
-        iscale.set_scaling_factor(self.pressure, 1e-6)
-        iscale.set_scaling_factor(self.conc_mass_comp, 1e1)
-
     def get_material_flow_terms(self, p, j):
         return self.material_flow_expression[j]
 

watertap/property_models/unit_specific/activated_sludge/asm1_reactions.py

@@ -34,9 +34,9 @@
 )
 from idaes.core.util.misc import add_object_reference
 from idaes.core.util.exceptions import BurntToast
-import idaes.core.util.scaling as iscale
 import idaes.logger as idaeslog
-
+from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme
+import idaes.core.util.scaling as iscale
 
 # Some more information about this module
 __author__ = "Andrew Lee, Xinhong Liu, Adam Atia"
@@ -330,12 +330,47 @@ def define_metadata(cls, obj):
         )
 
 
+class ASM1ReactionScaler(CustomScalerBase):
+    """
+    Scaler for the Activated Sludge 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}
+
+    def variable_scaling_routine(
+        self, model, overwrite: bool = False, submodel_scalers: dict = None
+    ):
+
+        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 this scaling methodology
+        # Consider 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,
+                )
+
+
 class _ASM1ReactionBlock(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 = ASM1ReactionScaler
+
     def initialize(self, outlvl=idaeslog.NOTSET, **kwargs):
         """
         Initialization routine for reaction package.