diff --git a/watertap/unit_models/aeration_tank.py b/watertap/unit_models/aeration_tank.py index f806a98b8d..5531eaba9b 100644 --- a/watertap/unit_models/aeration_tank.py +++ b/watertap/unit_models/aeration_tank.py @@ -20,6 +20,7 @@ from idaes.core import ( declare_process_block_class, ) +from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme from watertap.unit_models.cstr_injection import ( @@ -30,12 +31,132 @@ __author__ = "Adam Atia" +class AerationTankScaler(CustomScalerBase): + """ + Default modular scaler for the aeration tank unit model. + + 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-3, + } + + 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.control_volume.properties_in, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.propagate_state_scaling( + target_state=model.control_volume.properties_out, + source_state=model.control_volume.properties_in, + overwrite=overwrite, + ) + + self.call_submodel_scaler_method( + submodel=model.control_volume.properties_out, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.control_volume.reactions, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scaling control volume variables + self.scale_variable_by_default( + model.control_volume.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.control_volume.properties_in, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.control_volume.properties_out, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.control_volume.reactions, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level constraints + if hasattr(model, "eq_hydraulic_retention_time"): + self.scale_constraint_by_nominal_value( + model.eq_hydraulic_retention_time[0], + scheme=ConstraintScalingScheme.inverseMaximum, + overwrite=overwrite, + ) + + if hasattr(model, "eq_mass_transfer"): + self.scale_constraint_by_nominal_value( + model.eq_mass_transfer[0], + scheme=ConstraintScalingScheme.inverseMaximum, + overwrite=overwrite, + ) + + if hasattr(model, "eq_electricity_consumption"): + self.scale_constraint_by_nominal_value( + model.eq_electricity_consumption[0], + scheme=ConstraintScalingScheme.inverseMaximum, + overwrite=overwrite, + ) + + @declare_process_block_class("AerationTank") class AerationTankData(CSTR_InjectionData): """ CSTR Unit Model with Injection Class """ + default_scaler = AerationTankScaler + CONFIG = CSTR_InjectionData.CONFIG() CONFIG.electricity_consumption = ElectricityConsumption.calculated CONFIG.get("has_aeration")._domain = In([True]) diff --git a/watertap/unit_models/cstr.py b/watertap/unit_models/cstr.py index d77ebb8b7d..e5da6b7a7a 100644 --- a/watertap/unit_models/cstr.py +++ b/watertap/unit_models/cstr.py @@ -21,6 +21,7 @@ from idaes.models.unit_models.cstr import CSTRData as CSTRIDAESData import idaes.logger as idaeslog +from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme from pyomo.environ import ( Constraint, @@ -38,12 +39,118 @@ _log = idaeslog.getLogger(__name__) +class CSTRScaler(CustomScalerBase): + """ + Default modular scaler for CSTR. + + 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-3, + } + + 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.control_volume.properties_in, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.propagate_state_scaling( + target_state=model.control_volume.properties_out, + source_state=model.control_volume.properties_in, + overwrite=overwrite, + ) + + self.call_submodel_scaler_method( + submodel=model.control_volume.properties_out, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.control_volume.reactions, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scaling control volume variables + self.scale_variable_by_default( + model.control_volume.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.control_volume.properties_in, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.control_volume.properties_out, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.control_volume.reactions, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level constraints + if hasattr(model, "CSTR_retention_time"): + self.scale_constraint_by_nominal_value( + model.CSTR_retention_time[0], + scheme=ConstraintScalingScheme.inverseMaximum, + overwrite=overwrite, + ) + + @declare_process_block_class("CSTR") class CSTRData(CSTRIDAESData): """ CSTR unit block for BSM2 """ + default_scaler = CSTRScaler + CONFIG = CSTRIDAESData.CONFIG() def build(self): diff --git a/watertap/unit_models/cstr_injection.py b/watertap/unit_models/cstr_injection.py index 05e751e1a2..dea4fa6b35 100644 --- a/watertap/unit_models/cstr_injection.py +++ b/watertap/unit_models/cstr_injection.py @@ -36,6 +36,7 @@ UnitModelBlockData, useDefault, ) +from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme from idaes.core.util.config import ( is_physical_parameter_block, is_reaction_parameter_block, @@ -49,8 +50,6 @@ __author__ = "Andrew Lee, Adam Atia, Vibhav Dabadghao" -from enum import Enum, auto - class ElectricityConsumption(Enum): """ @@ -64,12 +63,132 @@ class ElectricityConsumption(Enum): calculated = auto() +class CSTR_InjectionScaler(CustomScalerBase): + """ + Default modular scaler for CSTR with injection. + + 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-3, + } + + 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.control_volume.properties_in, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.propagate_state_scaling( + target_state=model.control_volume.properties_out, + source_state=model.control_volume.properties_in, + overwrite=overwrite, + ) + + self.call_submodel_scaler_method( + submodel=model.control_volume.properties_out, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.control_volume.reactions, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scaling control volume variables + self.scale_variable_by_default( + model.control_volume.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.control_volume.properties_in, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.control_volume.properties_out, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.control_volume.reactions, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level constraints + if hasattr(model, "eq_hydraulic_retention_time"): + self.scale_constraint_by_nominal_value( + model.eq_hydraulic_retention_time[0], + scheme=ConstraintScalingScheme.inverseMaximum, + overwrite=overwrite, + ) + + if hasattr(model, "eq_mass_transfer"): + self.scale_constraint_by_nominal_value( + model.eq_mass_transfer[0], + scheme=ConstraintScalingScheme.inverseMaximum, + overwrite=overwrite, + ) + + if hasattr(model, "eq_electricity_consumption"): + self.scale_constraint_by_nominal_value( + model.eq_electricity_consumption[0], + scheme=ConstraintScalingScheme.inverseMaximum, + overwrite=overwrite, + ) + + @declare_process_block_class("CSTR_Injection") class CSTR_InjectionData(InitializationMixin, UnitModelBlockData): """ CSTR Unit Model with Injection Class """ + default_scaler = CSTR_InjectionScaler + CONFIG = UnitModelBlockData.CONFIG() CONFIG.declare( diff --git a/watertap/unit_models/tests/test_aeration_tank.py b/watertap/unit_models/tests/test_aeration_tank.py index 343a5d961e..36d91c57cc 100644 --- a/watertap/unit_models/tests/test_aeration_tank.py +++ b/watertap/unit_models/tests/test_aeration_tank.py @@ -14,12 +14,15 @@ Authors: Andrew Lee, Vibhav Dabadghao """ +from io import StringIO import pytest from pyomo.environ import ( ConcreteModel, units, value, assert_optimal_termination, + Suffix, + TransformationFactory, ) from idaes.core import ( FlowsheetBlock, @@ -34,6 +37,13 @@ number_total_constraints, number_unused_variables, ) +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) +from idaes.core.scaling.scaler_profiling import ScalingProfiler +import idaes.core.util.scaling as iscale +from idaes.core.scaling.scaling_base import ScalerBase from idaes.core.util.testing import ( initialization_tester, ) @@ -41,15 +51,21 @@ from watertap.core.solvers import get_solver from pyomo.util.check_units import assert_units_consistent, assert_units_equivalent -from watertap.unit_models.aeration_tank import AerationTank, ElectricityConsumption +from watertap.unit_models.aeration_tank import ( + AerationTank, + ElectricityConsumption, + AerationTankScaler, +) from idaes.core import UnitModelCostingBlock from watertap.costing import WaterTAPCosting from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( ASM1ParameterBlock, + ASM1PropertiesScaler, ) from watertap.property_models.unit_specific.activated_sludge.asm1_reactions import ( ASM1ReactionParameterBlock, + ASM1ReactionScaler, ) from watertap.property_models.unit_specific.activated_sludge.asm2d_properties import ( ASM2dParameterBlock, @@ -333,3 +349,639 @@ def test_error_without_oxygen(): property_package=m.fs.properties, reaction_package=m.fs.reactions, ) + + +class TestAerationTankScaler: + @pytest.fixture + def model(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = ASM1ParameterBlock() + m.fs.reactions = ASM1ReactionParameterBlock(property_package=m.fs.properties) + + m.fs.unit = AerationTank( + property_package=m.fs.properties, + reaction_package=m.fs.reactions, + electricity_consumption=ElectricityConsumption.calculated, + ) + + m.fs.unit.inlet.flow_vol.fix(20648 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(58 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(92 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(363 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(50 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(23 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(5 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(16 * units.g / units.m**3) + m.fs.unit.inlet.alkalinity.fix(7 * units.mol / units.m**3) + + m.fs.unit.volume.fix(500) + m.fs.unit.injection.fix(0) + m.fs.unit.injection[0, "Liq", "S_O"].fix(2e-3) + + return m + + @pytest.mark.component + def test_variable_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, AerationTankScaler) + + scaler.variable_scaling_routine(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.control_volume.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Outlet state - should be the same as the inlet + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 3 + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Reaction block + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 8 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R1"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R2"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R3"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R4"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R5"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R6"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R7"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R8"] + ] == pytest.approx(1e2, rel=1e-8) + + # Check that unit model has scaling factors + sfx_cv = model.fs.unit.control_volume.scaling_factor + assert isinstance(sfx_cv, Suffix) + assert len(sfx_cv) == 1 + assert sfx_cv[model.fs.unit.control_volume.volume[0]] == pytest.approx( + 1e-3, rel=1e-3 + ) + + @pytest.mark.component + def test_constraint_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, AerationTankScaler) + + scaler.constraint_scaling_routine(model.fs.unit) + + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 0 + + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 8 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R1"] + ] == pytest.approx(2.380752e5, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R2"] + ] == pytest.approx(1.49540985e8, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R3"] + ] == pytest.approx(1.75226112e6, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R4"] + ] == pytest.approx(2.88e6, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R5"] + ] == pytest.approx(1.728e7, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R6"] + ] == pytest.approx(1.728e5, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R7"] + ] == pytest.approx(3.174336e5, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R8"] + ] == pytest.approx(1, rel=1e-8) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 3 + assert sfx_unit[model.fs.unit.eq_hydraulic_retention_time[0]] == pytest.approx( + 4.77962962e-4, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_mass_transfer[0]] == pytest.approx( + 0.004, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_electricity_consumption[0]] == pytest.approx( + 0.09, rel=1e-8 + ) + + @pytest.mark.component + def test_scale_model(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, AerationTankScaler) + + scaler.scale_model(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.control_volume.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Outlet state - should be the same as the inlet + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 3 + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Reaction block + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 16 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R1"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R2"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R3"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R4"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R5"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R6"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R7"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R8"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R1"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R2"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R3"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R4"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R5"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R6"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R7"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R8"] + ] == pytest.approx(1e2, rel=1e-8) + + # Check that unit model has scaling factors + sfx_cv = model.fs.unit.control_volume.scaling_factor + assert isinstance(sfx_cv, Suffix) + assert len(sfx_cv) == 1 + assert sfx_cv[model.fs.unit.control_volume.volume[0]] == pytest.approx( + 1e-3, rel=1e-3 + ) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 3 + assert sfx_unit[model.fs.unit.eq_hydraulic_retention_time[0]] == pytest.approx( + 0.00047796296, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_mass_transfer[0]] == pytest.approx( + 0.004, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_electricity_consumption[0]] == pytest.approx( + 0.09, rel=1e-8 + ) + + # TODO: Remove test once iscale is deprecated + @pytest.mark.integration + def test_example_case_iscale(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = ASM1ParameterBlock() + m.fs.reactions = ASM1ReactionParameterBlock(property_package=m.fs.properties) + + m.fs.unit = AerationTank( + property_package=m.fs.properties, + reaction_package=m.fs.reactions, + electricity_consumption=ElectricityConsumption.calculated, + ) + + m.fs.unit.inlet.flow_vol.fix(20648 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(58 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(92 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(363 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(50 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(23 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(5 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(16 * units.g / units.m**3) + m.fs.unit.inlet.alkalinity.fix(7 * units.mol / units.m**3) + + m.fs.unit.volume.fix(500) + m.fs.unit.injection.fix(0) + m.fs.unit.injection[0, "Liq", "S_O"].fix(2e-3) + + # Set scaling factors for badly scaled variables + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0.0].pressure, 1e-5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0.0].conc_mass_comp["X_P"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_S"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_S"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_BH"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_P"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_O"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_NH"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_ND"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_ND"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_ALK"], 1e3 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R4"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R6"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R7"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R8"], 1e5 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R1"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R4"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R6"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R7"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R8"], 1e5 + ) + + iscale.calculate_scaling_factors(m.fs.unit) + + # Check condition number to confirm scaling + sm = TransformationFactory("core.scale_model").create_using(m, rename=False) + jac, _ = get_jacobian(sm, scaled=False) + assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx( + 1.0843930927394e13, rel=1e-3 + ) + + @pytest.mark.integration + def test_example_case_scaler(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = ASM1ParameterBlock() + m.fs.reactions = ASM1ReactionParameterBlock(property_package=m.fs.properties) + + m.fs.unit = AerationTank( + property_package=m.fs.properties, + reaction_package=m.fs.reactions, + electricity_consumption=ElectricityConsumption.calculated, + ) + + m.fs.unit.inlet.flow_vol.fix(20648 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(58 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(92 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(363 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(50 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(23 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(5 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(16 * units.g / units.m**3) + m.fs.unit.inlet.alkalinity.fix(7 * units.mol / units.m**3) + + m.fs.unit.volume.fix(500) + m.fs.unit.injection.fix(0) + m.fs.unit.injection[0, "Liq", "S_O"].fix(2e-3) + + # Set scaling factors for badly scaled variables + sb = ScalerBase() + sb.set_variable_scaling_factor(m.fs.unit.hydraulic_retention_time[0], 1e-3) + + scaler = AerationTankScaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.control_volume.properties_in: ASM1PropertiesScaler, + m.fs.unit.control_volume.properties_out: ASM1PropertiesScaler, + m.fs.unit.control_volume.reactions: ASM1ReactionScaler, + }, + ) + + # Check condition number to confirm scaling + sm = TransformationFactory("core.scale_model").create_using(m, rename=False) + jac, _ = get_jacobian(sm, scaled=False) + assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx( + 5.09707868988e11, rel=1e-3 + ) + + +def build_model(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = ASM1ParameterBlock() + m.fs.reactions = ASM1ReactionParameterBlock(property_package=m.fs.properties) + + m.fs.unit = AerationTank( + property_package=m.fs.properties, + reaction_package=m.fs.reactions, + electricity_consumption=ElectricityConsumption.calculated, + ) + + m.fs.unit.inlet.flow_vol.fix(20648 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(58 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(92 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(363 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(50 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(23 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(5 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(16 * units.g / units.m**3) + m.fs.unit.inlet.alkalinity.fix(7 * units.mol / units.m**3) + + m.fs.unit.volume.fix(500) + m.fs.unit.injection.fix(0) + m.fs.unit.injection[0, "Liq", "S_O"].fix(2e-3) + + solver = get_solver() + solver.solve(m) + + return m + + +def scale_vars_with_scalers(m): + scaler = AerationTankScaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.control_volume.properties_in: ASM1PropertiesScaler, + m.fs.unit.control_volume.properties_out: ASM1PropertiesScaler, + m.fs.unit.control_volume.reactions: ASM1ReactionScaler, + }, + ) + + +def scale_vars_with_iscale(m): + # Set scaling factors for badly scaled variables + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0.0].pressure, 1e-5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0.0].conc_mass_comp["X_P"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_S"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_S"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_BH"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_P"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_O"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_NH"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_ND"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_ND"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_ALK"], 1e3 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R4"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R6"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R7"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R8"], 1e5 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R1"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R4"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R6"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R7"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R8"], 1e5 + ) + + iscale.calculate_scaling_factors(m.fs.unit) + + +def perturb_solution(m): + m.fs.unit.inlet.flow_vol.fix(20648 * 0.9 * units.m**3 / units.day) + m.fs.unit.volume.fix(500 * 0.85) + m.fs.unit.injection[0, "Liq", "S_O"].fix(2e-3 * 0.7) + + +@pytest.mark.requires_idaes_solver +@pytest.mark.unit +def test_scaling_profiler_with_scalers(): + sp = ScalingProfiler( + build_model=build_model, + user_scaling=scale_vars_with_scalers, + perturb_state=perturb_solution, + ) + + stream = StringIO() + + sp.report_scaling_profiles(stream=stream) + + expected = """ +============================================================================ +Scaling Profile Report +---------------------------------------------------------------------------- +Scaling Method || User Scaling || Perfect Scaling +Unscaled || 1.826E+16 | Solved 14 || +Vars Only || 4.843E+13 | Solved 13 || 2.014E+21 | Solved 4 +Harmonic || 9.974E+17 | Solved 50 || 4.443E+22 | Solved 27 +Inverse Sum || 3.001E+17 | Solved 4 || 2.399E+14 | Solved 4 +Inverse Root Sum Squares || 3.001E+17 | Solved 4 || 3.412E+14 | Solved 4 +Inverse Maximum || 3.001E+17 | Solved 4 || 4.809E+14 | Solved 4 +Inverse Minimum || 9.974E+17 | Solved 29 || 4.455E+22 | Solved 24 +Nominal L1 Norm || 2.365E+09 | Solved 14 || 2.842E+14 | Solved 3 +Nominal L2 Norm || 1.648E+09 | Solved 14 || 3.755E+14 | Solved 3 +Actual L1 Norm || 8.636E+08 | Solved 14 || 5.461E+13 | Solved 4 +Actual L2 Norm || 7.902E+08 | Solved 11 || 6.491E+13 | Solved 4 +============================================================================ +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +@pytest.mark.requires_idaes_solver +def test_scaling_profiler_with_iscale(): + sp = ScalingProfiler( + build_model=build_model, + user_scaling=scale_vars_with_iscale, + perturb_state=perturb_solution, + ) + + stream = StringIO() + + sp.report_scaling_profiles(stream=stream) + + expected = """ +============================================================================ +Scaling Profile Report +---------------------------------------------------------------------------- +Scaling Method || User Scaling || Perfect Scaling +Unscaled || 1.826E+16 | Solved 14 || +Vars Only || 8.948E+12 | Solved 4 || 2.014E+21 | Solved 4 +Harmonic || 1.044E+17 | Solved 44 || 4.443E+22 | Solved 27 +Inverse Sum || 5.247E+17 | Solved 66 || 2.399E+14 | Solved 4 +Inverse Root Sum Squares || 5.220E+17 | Solved 73 || 3.412E+14 | Solved 4 +Inverse Maximum || 5.208E+17 | Solved 66 || 4.809E+14 | Solved 4 +Inverse Minimum || 2.103E+17 | Solved 85 || 4.455E+22 | Solved 24 +Nominal L1 Norm || 7.817E+09 | Solved 6 || 2.842E+14 | Solved 3 +Nominal L2 Norm || 1.278E+10 | Solved 6 || 3.755E+14 | Solved 3 +Actual L1 Norm || 3.950E+09 | Solved 3 || 5.461E+13 | Solved 4 +Actual L2 Norm || 4.339E+09 | Solved 3 || 6.491E+13 | Solved 4 +============================================================================ +""" + + assert stream.getvalue() == expected diff --git a/watertap/unit_models/tests/test_cstr.py b/watertap/unit_models/tests/test_cstr.py index 346ea4ecf1..41a7f34ef1 100644 --- a/watertap/unit_models/tests/test_cstr.py +++ b/watertap/unit_models/tests/test_cstr.py @@ -13,7 +13,7 @@ Tests for CSTR unit model. Authors: Marcus Holly """ - +from io import StringIO import pytest from pyomo.environ import ( @@ -21,23 +21,31 @@ units, value, Objective, + Suffix, + TransformationFactory, ) from watertap.unit_models.tests.unit_test_harness import UnitTestHarness import idaes.core.util.scaling as iscale - from idaes.core import ( FlowsheetBlock, UnitModelCostingBlock, ) -from watertap.unit_models.cstr import CSTR +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) +from idaes.core.scaling.scaler_profiling import ScalingProfiler +from watertap.unit_models.cstr import CSTR, CSTRScaler from watertap.costing import WaterTAPCosting from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( ASM1ParameterBlock, + ASM1PropertiesScaler, ) from watertap.property_models.unit_specific.activated_sludge.asm1_reactions import ( ASM1ReactionParameterBlock, + ASM1ReactionScaler, ) from idaes.models.properties.examples.saponification_thermo import ( @@ -391,21 +399,6 @@ def configure(self): assert pytest.approx(0.002127, rel=1e-3) == value(m.fs.costing.LCOW) - component_list = [ - "S_I", - "S_S", - "X_I", - "X_S", - "X_BH", - "X_BA", - "X_P", - "S_O", - "S_NO", - "S_NH", - "S_ND", - "X_ND", - ] - self.conservation_equality = { "Check 1": { "in": m.fs.unit.inlet.flow_vol[0], @@ -414,3 +407,608 @@ def configure(self): } return m + + +class TestCSTRScaler: + @pytest.fixture + def model(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM1 = ASM1ParameterBlock() + m.fs.ASM1_rxn_props = ASM1ReactionParameterBlock( + property_package=m.fs.props_ASM1 + ) + + m.fs.unit = CSTR( + property_package=m.fs.props_ASM1, reaction_package=m.fs.ASM1_rxn_props + ) + + m.fs.unit.inlet.flow_vol[0].fix(1.2199 * units.m**3 / units.s) + m.fs.unit.inlet.alkalinity[0].fix(4.5102 * units.mole / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(0.061909 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(0.012366 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(1.4258 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(0.090508 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(2.8404 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0.20512 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0.58681 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0.00036092 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0.012424 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(0.0076936 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.0019068 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(0.0053166 * units.kg / units.m**3) + + m.fs.unit.inlet.temperature[0].fix(308.15 * units.K) + m.fs.unit.inlet.pressure[0].fix(84790.0 * units.Pa) + + m.fs.unit.volume[0].fix(1000 * units.m**3) + + return m + + @pytest.mark.component + def test_variable_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, CSTRScaler) + + scaler.variable_scaling_routine(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.control_volume.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Outlet state - should be the same as the inlet + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 3 + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Reaction block + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 8 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R1"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R2"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R3"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R4"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R5"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R6"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R7"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R8"] + ] == pytest.approx(1e2, rel=1e-8) + + # Check that unit model has scaling factors + sfx_cv = model.fs.unit.control_volume.scaling_factor + assert isinstance(sfx_cv, Suffix) + assert len(sfx_cv) == 1 + assert sfx_cv[model.fs.unit.control_volume.volume[0]] == pytest.approx( + 1e-3, rel=1e-3 + ) + + @pytest.mark.component + def test_constraint_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, CSTRScaler) + + scaler.constraint_scaling_routine(model.fs.unit) + + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 0 + + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 8 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R1"] + ] == pytest.approx(2.380752e5, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R2"] + ] == pytest.approx(1.49540985e8, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R3"] + ] == pytest.approx(1.75226112e6, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R4"] + ] == pytest.approx(2.88e6, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R5"] + ] == pytest.approx(1.728e7, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R6"] + ] == pytest.approx(1.728e5, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R7"] + ] == pytest.approx(3.174336e5, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R8"] + ] == pytest.approx(1, rel=1e-8) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 1 + assert sfx_unit[model.fs.unit.CSTR_retention_time[0]] == pytest.approx( + 1.2199e-3, rel=1e-8 + ) + + @pytest.mark.component + def test_scale_model(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, CSTRScaler) + + scaler.scale_model(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.control_volume.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Outlet state - should be the same as the inlet + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 3 + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Reaction block + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 16 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R1"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R2"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R3"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R4"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R5"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R6"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R7"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R8"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R1"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R2"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R3"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R4"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R5"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R6"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R7"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R8"] + ] == pytest.approx(1e2, rel=1e-8) + + # Check that unit model has scaling factors + sfx_cv = model.fs.unit.control_volume.scaling_factor + assert isinstance(sfx_cv, Suffix) + assert len(sfx_cv) == 1 + assert sfx_cv[model.fs.unit.control_volume.volume[0]] == pytest.approx( + 1e-3, rel=1e-3 + ) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 1 + assert sfx_unit[model.fs.unit.CSTR_retention_time[0]] == pytest.approx( + 0.0012199, rel=1e-8 + ) + + # TODO: Remove test once iscale is deprecated + @pytest.mark.integration + def test_example_case_iscale(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM1 = ASM1ParameterBlock() + m.fs.ASM1_rxn_props = ASM1ReactionParameterBlock( + property_package=m.fs.props_ASM1 + ) + + m.fs.unit = CSTR( + property_package=m.fs.props_ASM1, reaction_package=m.fs.ASM1_rxn_props + ) + + m.fs.unit.inlet.flow_vol[0].fix(1.2199 * units.m**3 / units.s) + m.fs.unit.inlet.alkalinity[0].fix(4.5102 * units.mole / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(0.061909 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(0.012366 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(1.4258 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(0.090508 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(2.8404 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0.20512 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0.58681 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0.00036092 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0.012424 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(0.0076936 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.0019068 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(0.0053166 * units.kg / units.m**3) + + m.fs.unit.inlet.temperature[0].fix(308.15 * units.K) + m.fs.unit.inlet.pressure[0].fix(84790.0 * units.Pa) + + m.fs.unit.volume[0].fix(1000 * units.m**3) + + # Set scaling factors for badly scaled variables + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0].pressure, 1e-5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0.0].conc_mass_comp["S_O"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_BA"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_P"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_O"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_NH"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_ND"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_ND"], 1e5 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R1"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R3"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R5"], 1e5 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R1"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R2"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R3"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R4"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R5"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R6"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R7"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R8"], 1e7 + ) + + iscale.calculate_scaling_factors(m.fs.unit) + + # Check condition number to confirm scaling + sm = TransformationFactory("core.scale_model").create_using(m, rename=False) + jac, _ = get_jacobian(sm, scaled=False) + assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx( + 6.1277317e13, rel=1e-3 + ) + + @pytest.mark.integration + def test_example_case_scaler(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM1 = ASM1ParameterBlock() + m.fs.ASM1_rxn_props = ASM1ReactionParameterBlock( + property_package=m.fs.props_ASM1 + ) + + m.fs.unit = CSTR( + property_package=m.fs.props_ASM1, reaction_package=m.fs.ASM1_rxn_props + ) + + m.fs.unit.inlet.flow_vol[0].fix(1.2199 * units.m**3 / units.s) + m.fs.unit.inlet.alkalinity[0].fix(4.5102 * units.mole / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(0.061909 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(0.012366 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(1.4258 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(0.090508 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(2.8404 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0.20512 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0.58681 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0.00036092 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0.012424 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(0.0076936 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.0019068 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(0.0053166 * units.kg / units.m**3) + + m.fs.unit.inlet.temperature[0].fix(308.15 * units.K) + m.fs.unit.inlet.pressure[0].fix(84790.0 * units.Pa) + + m.fs.unit.volume[0].fix(1000 * units.m**3) + + scaler = CSTRScaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.control_volume.properties_in: ASM1PropertiesScaler, + m.fs.unit.control_volume.properties_out: ASM1PropertiesScaler, + m.fs.unit.control_volume.reactions: ASM1ReactionScaler, + }, + ) + + # Check condition number to confirm scaling + sm = TransformationFactory("core.scale_model").create_using(m, rename=False) + jac, _ = get_jacobian(sm, scaled=False) + assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx( + 3.626085e10, rel=1e-3 + ) + + +def build_model(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM1 = ASM1ParameterBlock() + m.fs.ASM1_rxn_props = ASM1ReactionParameterBlock(property_package=m.fs.props_ASM1) + + m.fs.unit = CSTR( + property_package=m.fs.props_ASM1, reaction_package=m.fs.ASM1_rxn_props + ) + + m.fs.unit.inlet.flow_vol[0].fix(1.2199 * units.m**3 / units.s) + m.fs.unit.inlet.alkalinity[0].fix(4.5102 * units.mole / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(0.061909 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(0.012366 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(1.4258 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(0.090508 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(2.8404 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0.20512 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0.58681 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0.00036092 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0.012424 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(0.0076936 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.0019068 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(0.0053166 * units.kg / units.m**3) + + m.fs.unit.inlet.temperature[0].fix(308.15 * units.K) + m.fs.unit.inlet.pressure[0].fix(84790.0 * units.Pa) + + m.fs.unit.volume[0].fix(1000 * units.m**3) + + solver = get_solver() + solver.solve(m) + + return m + + +def scale_vars_with_scalers(m): + scaler = CSTRScaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.control_volume.properties_in: ASM1PropertiesScaler, + m.fs.unit.control_volume.properties_out: ASM1PropertiesScaler, + m.fs.unit.control_volume.reactions: ASM1ReactionScaler, + }, + ) + + +def scale_vars_with_iscale(m): + # Set scaling factors for badly scaled variables + iscale.set_scaling_factor(m.fs.unit.control_volume.properties_out[0].pressure, 1e-5) + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0.0].conc_mass_comp["S_O"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_BA"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_P"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_O"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_NH"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_ND"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_ND"], 1e5 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R1"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R3"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R5"], 1e5 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R1"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R2"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R3"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R4"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R5"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R6"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R7"], 1e7 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R8"], 1e7 + ) + + iscale.calculate_scaling_factors(m.fs.unit) + + +def perturb_solution(m): + m.fs.unit.inlet.flow_vol[0].fix(1.2199 * 0.5 * units.m**3 / units.s) + m.fs.unit.volume[0].fix(1000 * 0.75 * units.m**3) + + +@pytest.mark.requires_idaes_solver +@pytest.mark.unit +def test_scaling_profiler_with_scalers(): + sp = ScalingProfiler( + build_model=build_model, + user_scaling=scale_vars_with_scalers, + perturb_state=perturb_solution, + ) + + stream = StringIO() + + sp.report_scaling_profiles(stream=stream) + + expected = """ +============================================================================ +Scaling Profile Report +---------------------------------------------------------------------------- +Scaling Method || User Scaling || Perfect Scaling +Unscaled || 1.196E+12 | Solved 4 || +Vars Only || 4.434E+10 | Solved 5 || 1.674E+17 | Solved 1 +Harmonic || 1.124E+10 | Solved 5 || 1.161E+05 | Solved 3 +Inverse Sum || 2.031E+07 | Solved 5 || 2.054E+02 | Solved 3 +Inverse Root Sum Squares || 3.266E+07 | Solved 5 || 2.131E+02 | Solved 3 +Inverse Maximum || 4.057E+07 | Solved 5 || 2.427E+02 | Solved 3 +Inverse Minimum || 1.089E+10 | Solved 5 || 1.317E+05 | Solved 3 +Nominal L1 Norm || 1.799E+04 | Solved 5 || 2.183E+02 | Solved 3 +Nominal L2 Norm || 1.789E+04 | Solved 5 || 2.051E+02 | Solved 3 +Actual L1 Norm || 1.548E+04 | Solved 5 || 2.146E+02 | Solved 3 +Actual L2 Norm || 1.569E+04 | Solved 5 || 2.114E+02 | Solved 3 +============================================================================ +""" + + assert stream.getvalue() == expected + + +@pytest.mark.requires_idaes_solver +@pytest.mark.unit +def test_scaling_profiler_with_iscale(): + sp = ScalingProfiler( + build_model=build_model, + user_scaling=scale_vars_with_iscale, + perturb_state=perturb_solution, + ) + + stream = StringIO() + + sp.report_scaling_profiles(stream=stream) + + expected = """ +============================================================================ +Scaling Profile Report +---------------------------------------------------------------------------- +Scaling Method || User Scaling || Perfect Scaling +Unscaled || 1.196E+12 | Solved 4 || +Vars Only || 3.003E+13 | Solved 3 || 1.674E+17 | Solved 1 +Harmonic || 9.790E+09 | Solved 5 || 1.161E+05 | Solved 3 +Inverse Sum || 2.885E+07 | Solved 5 || 2.054E+02 | Solved 3 +Inverse Root Sum Squares || 3.046E+07 | Solved 5 || 2.131E+02 | Solved 3 +Inverse Maximum || 3.164E+07 | Solved 5 || 2.427E+02 | Solved 3 +Inverse Minimum || 1.879E+10 | Solved 5 || 1.317E+05 | Solved 3 +Nominal L1 Norm || 1.261E+09 | Solved 5 || 2.183E+02 | Solved 3 +Nominal L2 Norm || 9.543E+08 | Solved 5 || 2.051E+02 | Solved 3 +Actual L1 Norm || 1.665E+06 | Solved 4 || 2.146E+02 | Solved 3 +Actual L2 Norm || 1.736E+06 | Solved 4 || 2.114E+02 | Solved 3 +============================================================================ +""" + + assert stream.getvalue() == expected diff --git a/watertap/unit_models/tests/test_cstr_injection.py b/watertap/unit_models/tests/test_cstr_injection.py index 31dcfb8bfc..915ff1b96e 100644 --- a/watertap/unit_models/tests/test_cstr_injection.py +++ b/watertap/unit_models/tests/test_cstr_injection.py @@ -14,12 +14,15 @@ Authors: Andrew Lee, Adam Atia, Vibhav Dabadghao """ +from io import StringIO import pytest from pyomo.environ import ( ConcreteModel, units, value, Objective, + Suffix, + TransformationFactory, ) from idaes.core import ( FlowsheetBlock, @@ -27,6 +30,12 @@ from watertap.unit_models.tests.unit_test_harness import UnitTestHarness import idaes.core.util.scaling as iscale +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) +from idaes.core.scaling.scaler_profiling import ScalingProfiler +from idaes.core.scaling.scaling_base import ScalerBase from idaes.models.properties.examples.saponification_thermo import ( SaponificationParameterBlock, ) @@ -36,14 +45,20 @@ from idaes.core.util.exceptions import ConfigurationError from watertap.core.solvers import get_solver -from watertap.unit_models.cstr_injection import CSTR_Injection, ElectricityConsumption +from watertap.unit_models.cstr_injection import ( + CSTR_Injection, + ElectricityConsumption, + CSTR_InjectionScaler, +) from idaes.core import UnitModelCostingBlock from watertap.costing import WaterTAPCosting from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( ASM1ParameterBlock, + ASM1PropertiesScaler, ) from watertap.property_models.unit_specific.activated_sludge.asm1_reactions import ( ASM1ReactionParameterBlock, + ASM1ReactionScaler, ) from watertap.property_models.unit_specific.activated_sludge.asm2d_properties import ( ASM2dParameterBlock, @@ -467,3 +482,642 @@ def test_error_without_oxygen(): reaction_package=m.fs.reactions, has_aeration=True, ) + + +class TestCSTR_InjectionScaler: + @pytest.fixture + def model(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = ASM1ParameterBlock() + m.fs.reactions = ASM1ReactionParameterBlock(property_package=m.fs.properties) + + m.fs.unit = CSTR_Injection( + property_package=m.fs.properties, + reaction_package=m.fs.reactions, + has_aeration=True, + electricity_consumption=ElectricityConsumption.calculated, + ) + + m.fs.unit.inlet.flow_vol.fix(20648 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(58 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(92 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(363 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(50 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(23 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(5 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(16 * units.g / units.m**3) + m.fs.unit.inlet.alkalinity.fix(7 * units.mol / units.m**3) + + m.fs.unit.volume.fix(500) + m.fs.unit.injection.fix(0) + m.fs.unit.injection[0, "Liq", "S_O"].fix(2e-3) + + return m + + @pytest.mark.component + def test_variable_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, CSTR_InjectionScaler) + + scaler.variable_scaling_routine(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.control_volume.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Outlet state - should be the same as the inlet + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 3 + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Reaction block + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 8 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R1"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R2"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R3"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R4"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R5"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R6"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R7"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R8"] + ] == pytest.approx(1e2, rel=1e-8) + + # Check that unit model has scaling factors + sfx_cv = model.fs.unit.control_volume.scaling_factor + assert isinstance(sfx_cv, Suffix) + assert len(sfx_cv) == 1 + assert sfx_cv[model.fs.unit.control_volume.volume[0]] == pytest.approx( + 1e-3, rel=1e-3 + ) + + @pytest.mark.component + def test_constraint_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, CSTR_InjectionScaler) + + scaler.constraint_scaling_routine(model.fs.unit) + + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 0 + + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 8 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R1"] + ] == pytest.approx(2.380752e5, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R2"] + ] == pytest.approx(1.49540985e8, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R3"] + ] == pytest.approx(1.75226112e6, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R4"] + ] == pytest.approx(2.88e6, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R5"] + ] == pytest.approx(1.728e7, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R6"] + ] == pytest.approx(1.728e5, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R7"] + ] == pytest.approx(3.174336e5, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R8"] + ] == pytest.approx(1, rel=1e-8) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 3 + assert sfx_unit[model.fs.unit.eq_hydraulic_retention_time[0]] == pytest.approx( + 4.77962962e-4, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_mass_transfer[0]] == pytest.approx( + 0.004, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_electricity_consumption[0]] == pytest.approx( + 0.09, rel=1e-8 + ) + + @pytest.mark.component + def test_scale_model(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, CSTR_InjectionScaler) + + scaler.scale_model(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.control_volume.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Outlet state - should be the same as the inlet + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 3 + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Reaction block + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 16 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R1"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R2"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R3"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R4"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R5"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R6"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R7"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R8"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R1"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R2"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R3"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R4"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R5"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R6"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R7"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R8"] + ] == pytest.approx(1e2, rel=1e-8) + + # Check that unit model has scaling factors + sfx_cv = model.fs.unit.control_volume.scaling_factor + assert isinstance(sfx_cv, Suffix) + assert len(sfx_cv) == 1 + assert sfx_cv[model.fs.unit.control_volume.volume[0]] == pytest.approx( + 1e-3, rel=1e-3 + ) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 3 + assert sfx_unit[model.fs.unit.eq_hydraulic_retention_time[0]] == pytest.approx( + 0.00047796296, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_mass_transfer[0]] == pytest.approx( + 0.004, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_electricity_consumption[0]] == pytest.approx( + 0.09, rel=1e-8 + ) + + # TODO: Remove test once iscale is deprecated + @pytest.mark.integration + def test_example_case_iscale(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = ASM1ParameterBlock() + m.fs.reactions = ASM1ReactionParameterBlock(property_package=m.fs.properties) + + m.fs.unit = CSTR_Injection( + property_package=m.fs.properties, + reaction_package=m.fs.reactions, + has_aeration=True, + electricity_consumption=ElectricityConsumption.calculated, + ) + + m.fs.unit.inlet.flow_vol.fix(20648 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(58 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(92 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(363 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(50 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(23 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(5 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(16 * units.g / units.m**3) + m.fs.unit.inlet.alkalinity.fix(7 * units.mol / units.m**3) + + m.fs.unit.volume.fix(500) + m.fs.unit.injection.fix(0) + m.fs.unit.injection[0, "Liq", "S_O"].fix(2e-3) + + # Set scaling factors for badly scaled variables + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0.0].pressure, 1e-5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0.0].conc_mass_comp["X_P"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_S"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_S"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_BH"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_P"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_O"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_NH"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_ND"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_ND"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_ALK"], 1e3 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R4"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R6"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R7"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R8"], 1e5 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R1"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R4"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R6"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R7"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R8"], 1e5 + ) + + iscale.calculate_scaling_factors(m.fs.unit) + + # Check condition number to confirm scaling + sm = TransformationFactory("core.scale_model").create_using(m, rename=False) + jac, _ = get_jacobian(sm, scaled=False) + assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx( + 1.0843930927394e13, rel=1e-3 + ) + + @pytest.mark.integration + def test_example_case_scaler(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = ASM1ParameterBlock() + m.fs.reactions = ASM1ReactionParameterBlock(property_package=m.fs.properties) + + m.fs.unit = CSTR_Injection( + property_package=m.fs.properties, + reaction_package=m.fs.reactions, + has_aeration=True, + electricity_consumption=ElectricityConsumption.calculated, + ) + + m.fs.unit.inlet.flow_vol.fix(20648 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(58 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(92 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(363 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(50 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(23 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(5 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(16 * units.g / units.m**3) + m.fs.unit.inlet.alkalinity.fix(7 * units.mol / units.m**3) + + m.fs.unit.volume.fix(500) + m.fs.unit.injection.fix(0) + m.fs.unit.injection[0, "Liq", "S_O"].fix(2e-3) + + # Set scaling factors for badly scaled variables + sb = ScalerBase() + sb.set_variable_scaling_factor(m.fs.unit.hydraulic_retention_time[0], 1e-3) + + scaler = CSTR_InjectionScaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.control_volume.properties_in: ASM1PropertiesScaler, + m.fs.unit.control_volume.properties_out: ASM1PropertiesScaler, + m.fs.unit.control_volume.reactions: ASM1ReactionScaler, + }, + ) + + # Check condition number to confirm scaling + sm = TransformationFactory("core.scale_model").create_using(m, rename=False) + jac, _ = get_jacobian(sm, scaled=False) + assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx( + 5.09707868988e11, rel=1e-3 + ) + + +def build_model(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = ASM1ParameterBlock() + m.fs.reactions = ASM1ReactionParameterBlock(property_package=m.fs.properties) + + m.fs.unit = CSTR_Injection( + property_package=m.fs.properties, + reaction_package=m.fs.reactions, + has_aeration=True, + electricity_consumption=ElectricityConsumption.calculated, + ) + + m.fs.unit.inlet.flow_vol.fix(20648 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(58 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(92 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(363 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(50 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(0 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(23 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(5 * units.g / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(16 * units.g / units.m**3) + m.fs.unit.inlet.alkalinity.fix(7 * units.mol / units.m**3) + + m.fs.unit.volume.fix(500) + m.fs.unit.injection.fix(0) + m.fs.unit.injection[0, "Liq", "S_O"].fix(2e-3) + + solver = get_solver() + solver.solve(m) + + return m + + +def scale_vars_with_scalers(m): + scaler = CSTR_InjectionScaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.control_volume.properties_in: ASM1PropertiesScaler, + m.fs.unit.control_volume.properties_out: ASM1PropertiesScaler, + m.fs.unit.control_volume.reactions: ASM1ReactionScaler, + }, + ) + + +def scale_vars_with_iscale(m): + # Set scaling factors for badly scaled variables + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0.0].pressure, 1e-5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.properties_out[0.0].conc_mass_comp["X_P"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_S"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_S"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_BH"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_P"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_O"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_NH"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_ND"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "X_ND"], 1e3 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_generation[0.0, "Liq", "S_ALK"], 1e3 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R4"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R6"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R7"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.rate_reaction_extent[0.0, "R8"], 1e5 + ) + + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R1"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R4"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R6"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R7"], 1e5 + ) + iscale.set_scaling_factor( + m.fs.unit.control_volume.reactions[0.0].reaction_rate["R8"], 1e5 + ) + + iscale.calculate_scaling_factors(m.fs.unit) + + +def perturb_solution(m): + m.fs.unit.inlet.flow_vol.fix(20648 * 0.9 * units.m**3 / units.day) + m.fs.unit.volume.fix(500 * 0.85) + + +@pytest.mark.requires_idaes_solver +@pytest.mark.unit +def test_scaling_profiler_with_scalers(): + sp = ScalingProfiler( + build_model=build_model, + user_scaling=scale_vars_with_scalers, + perturb_state=perturb_solution, + ) + + stream = StringIO() + + sp.report_scaling_profiles(stream=stream) + + expected = """ +============================================================================ +Scaling Profile Report +---------------------------------------------------------------------------- +Scaling Method || User Scaling || Perfect Scaling +Unscaled || 1.826E+16 | Solved 4 || +Vars Only || 4.843E+13 | Solved 4 || 2.014E+21 | Solved 4 +Harmonic || 9.974E+17 | Failed 49 || 4.443E+22 | Solved 18 +Inverse Sum || 3.001E+17 | Solved 10 || 2.399E+14 | Solved 4 +Inverse Root Sum Squares || 3.001E+17 | Solved 4 || 3.412E+14 | Solved 4 +Inverse Maximum || 3.001E+17 | Solved 4 || 4.809E+14 | Solved 4 +Inverse Minimum || 9.974E+17 | Failed 49 || 4.455E+22 | Solved 18 +Nominal L1 Norm || 2.365E+09 | Solved 4 || 2.842E+14 | Solved 4 +Nominal L2 Norm || 1.648E+09 | Solved 4 || 3.755E+14 | Solved 4 +Actual L1 Norm || 8.636E+08 | Solved 4 || 5.461E+13 | Solved 4 +Actual L2 Norm || 7.902E+08 | Solved 4 || 6.491E+13 | Solved 4 +============================================================================ +""" + + assert stream.getvalue() == expected + + +@pytest.mark.requires_idaes_solver +@pytest.mark.unit +def test_scaling_profiler_with_iscale(): + sp = ScalingProfiler( + build_model=build_model, + user_scaling=scale_vars_with_iscale, + perturb_state=perturb_solution, + ) + + stream = StringIO() + + sp.report_scaling_profiles(stream=stream) + + expected = """ +============================================================================ +Scaling Profile Report +---------------------------------------------------------------------------- +Scaling Method || User Scaling || Perfect Scaling +Unscaled || 1.826E+16 | Solved 4 || +Vars Only || 8.948E+12 | Solved 4 || 2.014E+21 | Solved 4 +Harmonic || 1.044E+17 | Solved 57 || 4.443E+22 | Solved 18 +Inverse Sum || 5.247E+17 | Failed 50 || 2.399E+14 | Solved 4 +Inverse Root Sum Squares || 5.220E+17 | Failed 55 || 3.412E+14 | Solved 4 +Inverse Maximum || 5.208E+17 | Failed 52 || 4.809E+14 | Solved 4 +Inverse Minimum || 2.103E+17 | Solved 65 || 4.455E+22 | Solved 18 +Nominal L1 Norm || 7.817E+09 | Solved 4 || 2.842E+14 | Solved 4 +Nominal L2 Norm || 1.278E+10 | Solved 4 || 3.755E+14 | Solved 4 +Actual L1 Norm || 3.950E+09 | Solved 3 || 5.461E+13 | Solved 4 +Actual L2 Norm || 4.339E+09 | Solved 3 || 6.491E+13 | Solved 4 +============================================================================ +""" + + assert stream.getvalue() == expected