Modified ASM2d model with new scaling tool (#1533)

* add

* delete redundant files

* add new scaling tool to modified ASM2d

---------

Co-authored-by: Ludovico Bianchi <lbianchi@lbl.gov>

watertap/property_models/unit_specific/activated_sludge/modified_asm2d_properties.py

@@ -31,7 +31,6 @@
 from idaes.core import (
     declare_process_block_class,
     MaterialFlowBasis,
-    PhysicalParameterBlock,
     StateBlockData,
     StateBlock,
     MaterialBalanceType,
@@ -43,6 +42,8 @@
 from idaes.core.util.model_statistics import degrees_of_freedom
 from idaes.core.util.initialization import fix_state_vars, revert_state_vars
 import idaes.logger as idaeslog
+from idaes.core.scaling import CustomScalerBase
+from idaes.core.base.property_base import PhysicalParameterBlock
 
 # Some more information about this module
 __author__ = "Marcus Holly, Adam Atia, Xinhong Liu"
@@ -345,12 +346,45 @@ def define_metadata(cls, obj):
         )
 
 
+class ModifiedASM2dPropertiesScaler(CustomScalerBase):
+    """
+    Scaler for the Activated Sludge Model No.2d 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-5),
+    }
+
+    DEFAULT_SCALING_FACTORS = {
+        "flow_vol": 1e1,
+        "temperature": 1e-2,
+    }
+
+    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 _ModifiedASM2dStateBlock(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 = ModifiedASM2dPropertiesScaler
+
     def initialize(
         self,
         state_args=None,

watertap/property_models/unit_specific/activated_sludge/modified_asm2d_reactions.py

@@ -36,6 +36,7 @@
 
 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
@@ -986,12 +987,46 @@ def define_metadata(cls, obj):
         )
 
 
+class ModifiedASM2dReactionScaler(CustomScalerBase):
+    """
+    Scaler for the Activated Sludge Model No.2d 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 _ModifiedASM2dReactionBlock(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 = ModifiedASM2dReactionScaler
+
     def initialize(self, outlvl=idaeslog.NOTSET, **kwargs):
         """
         Initialization routine for reaction package.

watertap/property_models/unit_specific/activated_sludge/tests/test_modified_asm2d_reaction.py

@@ -25,6 +25,7 @@
     check_optimal_termination,
     ConcreteModel,
     Constraint,
+    Suffix,
     units,
     value,
     Var,
@@ -43,6 +44,7 @@
 from watertap.property_models.unit_specific.activated_sludge.modified_asm2d_reactions import (
     ModifiedASM2dReactionParameterBlock,
     ModifiedASM2dReactionBlock,
+    ModifiedASM2dReactionScaler,
 )
 import idaes.core.util.scaling as iscale
 
@@ -308,6 +310,191 @@ def check_units(self, model):
         assert_units_consistent(model)
 
 
+class TestASM1ReactionScaler(object):
+    @pytest.mark.unit
+    def test_variable_scaling_routine(self):
+        model = ConcreteModel()
+        model.pparams = ModifiedASM2dParameterBlock()
+        model.rparams = ModifiedASM2dReactionParameterBlock(
+            property_package=model.pparams
+        )
+
+        model.props = model.pparams.build_state_block([1])
+        model.rxns = model.rparams.build_reaction_block([1], state_block=model.props)
+
+        # Trigger build of reaction properties
+        model.rxns[1].reaction_rate
+
+        scaler = model.rxns[1].default_scaler()
+        assert isinstance(scaler, ModifiedASM2dReactionScaler)
+
+        scaler.variable_scaling_routine(model.rxns[1])
+
+        assert isinstance(model.rxns[1].scaling_factor, Suffix)
+
+        sfx = model.rxns[1].scaling_factor
+        assert len(sfx) == 19
+        assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R2"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R3"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R4"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R5"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R6"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R7"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R8"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R9"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R10"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R11"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R12"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R13"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R14"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R15"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R16"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R17"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R18"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R19"]] == pytest.approx(1e2, rel=1e-8)
+
+    @pytest.mark.unit
+    def test_constraint_scaling_routine(self):
+        model = ConcreteModel()
+        model.pparams = ModifiedASM2dParameterBlock()
+        model.rparams = ModifiedASM2dReactionParameterBlock(
+            property_package=model.pparams
+        )
+
+        model.props = model.pparams.build_state_block([1])
+        model.rxns = model.rparams.build_reaction_block([1], state_block=model.props)
+
+        # Trigger build of reaction properties
+        model.rxns[1].reaction_rate
+
+        scaler = model.rxns[1].default_scaler()
+        assert isinstance(scaler, ModifiedASM2dReactionScaler)
+
+        scaler.constraint_scaling_routine(model.rxns[1])
+
+        assert isinstance(model.rxns[1].scaling_factor, Suffix)
+
+        sfx = model.rxns[1].scaling_factor
+        assert len(sfx) == 19
+        assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx(
+            387114.1, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R2"]] == pytest.approx(
+            324208097.5, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R3"]] == pytest.approx(1e10, rel=1e-5)
+        assert sfx[model.rxns[1].rate_expression["R4"]] == pytest.approx(
+            425956.2, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R5"]] == pytest.approx(
+            425956.2, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R6"]] == pytest.approx(
+            267553743, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R7"]] == pytest.approx(
+            267553743, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R8"]] == pytest.approx(1e10, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R9"]] == pytest.approx(
+            3091885.7, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R10"]] == pytest.approx(
+            368921.0, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R11"]] == pytest.approx(
+            690719.1, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R12"]] == pytest.approx(
+            578477274, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R13"]] == pytest.approx(
+            1066963, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R14"]] == pytest.approx(
+            893581814, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R15"]] == pytest.approx(
+            6183771.4, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R16"]] == pytest.approx(
+            6183771.4, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R17"]] == pytest.approx(
+            6183771.4, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R18"]] == pytest.approx(
+            1437853.9, rel=1e-5
+        )
+        assert sfx[model.rxns[1].rate_expression["R19"]] == pytest.approx(
+            9648000, rel=1e-5
+        )
+
+    @pytest.mark.unit
+    def test_scale_model(self):
+        model = ConcreteModel()
+        model.pparams = ModifiedASM2dParameterBlock()
+        model.rparams = ModifiedASM2dReactionParameterBlock(
+            property_package=model.pparams
+        )
+
+        model.props = model.pparams.build_state_block([1])
+        model.rxns = model.rparams.build_reaction_block([1], state_block=model.props)
+
+        # Trigger build of reaction properties
+        model.rxns[1].reaction_rate
+
+        scaler = model.rxns[1].default_scaler()
+        assert isinstance(scaler, ModifiedASM2dReactionScaler)
+
+        scaler.scale_model(model.rxns[1])
+
+        assert isinstance(model.rxns[1].scaling_factor, Suffix)
+
+        sfx = model.rxns[1].scaling_factor
+        assert len(sfx) == 38
+        assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R2"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R3"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R4"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R5"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R6"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R7"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R8"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R9"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R10"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R11"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R12"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R13"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R14"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R15"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R16"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R17"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R18"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].reaction_rate["R19"]] == pytest.approx(1e2, rel=1e-8)
+
+        assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R2"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R3"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R4"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R5"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R6"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R7"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R8"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R9"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R10"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R11"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R12"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R13"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R14"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R15"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R16"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R17"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R18"]] == pytest.approx(1e2, rel=1e-8)
+        assert sfx[model.rxns[1].rate_expression["R19"]] == pytest.approx(1e2, rel=1e-8)
+
+
 class TestAerobic:
     @pytest.fixture(scope="class")
     def model(self):