From d18ade632a3001c4de74e42d574862147736d86d Mon Sep 17 00:00:00 2001 From: MarcusHolly Date: Fri, 20 Dec 2024 11:39:55 -0500 Subject: [PATCH] Add scalers and testing for thickener & dewaterer --- watertap/unit_models/dewatering.py | 108 ++++- .../unit_models/tests/test_dewatering_unit.py | 366 ++++++++++++++- .../unit_models/tests/test_thickener_unit.py | 415 +++++++++++++++++- watertap/unit_models/thickener.py | 106 +++++ 4 files changed, 991 insertions(+), 4 deletions(-) diff --git a/watertap/unit_models/dewatering.py b/watertap/unit_models/dewatering.py index a5c14b00b2..ccade111b8 100644 --- a/watertap/unit_models/dewatering.py +++ b/watertap/unit_models/dewatering.py @@ -34,16 +34,18 @@ import idaes.logger as idaeslog from pyomo.environ import ( + Constraint, Param, - units as pyunits, Var, NonNegativeReals, + units as pyunits, ) from pyomo.common.config import ConfigValue, In from idaes.core.util.exceptions import ( ConfigurationError, ) +from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme from watertap.costing.unit_models.dewatering import cost_dewatering __author__ = "Alejandro Garciadiego, Adam Atia" @@ -65,12 +67,116 @@ class ActivatedSludgeModelType(Enum): modified_ASM2D = auto() +class DewatererScaler(CustomScalerBase): + """ + Default modular scaler for the dewatering 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.mixed_state, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.propagate_state_scaling( + target_state=model.underflow_state, + source_state=model.mixed_state, + overwrite=overwrite, + ) + self.propagate_state_scaling( + target_state=model.overflow_state, + source_state=model.mixed_state, + overwrite=overwrite, + ) + + self.call_submodel_scaler_method( + submodel=model.underflow_state, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.overflow_state, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level variables + self.scale_variable_by_default(model.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.mixed_state, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.underflow_state, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.overflow_state, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level constraints + # Constraint( + for c in model.component_data_objects(Constraint, descend_into=False): + self.scale_constraint_by_nominal_value( + c, + scheme=ConstraintScalingScheme.inverseMaximum, + overwrite=overwrite, + ) + + @declare_process_block_class("DewateringUnit") class DewateringData(SeparatorData): """ Dewatering unit block for BSM2 """ + default_scaler = DewatererScaler + CONFIG = SeparatorData.CONFIG() CONFIG.outlet_list = ["underflow", "overflow"] CONFIG.split_basis = SplittingType.componentFlow diff --git a/watertap/unit_models/tests/test_dewatering_unit.py b/watertap/unit_models/tests/test_dewatering_unit.py index 50a157ae62..ccbd61880b 100644 --- a/watertap/unit_models/tests/test_dewatering_unit.py +++ b/watertap/unit_models/tests/test_dewatering_unit.py @@ -18,6 +18,9 @@ ConcreteModel, value, assert_optimal_termination, + Suffix, + TransformationFactory, + Var, units as pyunits, ) @@ -27,6 +30,12 @@ MomentumBalanceType, ) +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) +from idaes.core.scaling.scaling_base import ScalerBase + from idaes.models.unit_models.separator import SplittingType from pyomo.environ import ( @@ -49,9 +58,14 @@ ConfigurationError, ) -from watertap.unit_models.dewatering import DewateringUnit, ActivatedSludgeModelType +from watertap.unit_models.dewatering import ( + DewateringUnit, + ActivatedSludgeModelType, + DewatererScaler, +) from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( ASM1ParameterBlock, + ASM1PropertiesScaler, ) from watertap.property_models.unit_specific.activated_sludge.asm2d_properties import ( @@ -821,3 +835,353 @@ def test_du_costing_config_err(): "dewatering_type": "foo", }, ) + + +class TestThickenerScaler: + @pytest.fixture + def model(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = ASM1ParameterBlock() + + m.fs.unit = DewateringUnit(property_package=m.fs.props) + + m.fs.unit.inlet.flow_vol.fix(178.4674 * 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(130.867 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(258.5789 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix( + 17216.2434 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(2611.4843 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(626.0652 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix( + 1442.7882 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.54323 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(100.8668 * units.mg / units.liter) + m.fs.unit.inlet.alkalinity.fix(97.8459 * units.mol / units.m**3) + + m.fs.unit.hydraulic_retention_time.fix() + + return m + + @pytest.mark.component + def test_variable_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, DewatererScaler) + + scaler.variable_scaling_routine(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.mixed_state[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.mixed_state[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_in[model.fs.unit.mixed_state[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_in[model.fs.unit.mixed_state[0].temperature] == pytest.approx( + 1e-1, rel=1e-8 + ) + + # Outlet state - should be the same as the inlet + sfx_underflow = model.fs.unit.underflow_state[0].scaling_factor + assert isinstance(sfx_underflow, Suffix) + assert len(sfx_underflow) == 3 + assert sfx_underflow[ + model.fs.unit.underflow_state[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_underflow[ + model.fs.unit.underflow_state[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_underflow[ + model.fs.unit.underflow_state[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + sfx_overflow = model.fs.unit.overflow_state[0].scaling_factor + assert isinstance(sfx_overflow, Suffix) + assert len(sfx_overflow) == 3 + assert sfx_overflow[model.fs.unit.overflow_state[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_overflow[model.fs.unit.overflow_state[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_overflow[ + model.fs.unit.overflow_state[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Check that unit model has scaling factors + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 1 + assert sfx_unit[model.fs.unit.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, DewatererScaler) + + scaler.constraint_scaling_routine(model.fs.unit) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 62 + assert sfx_unit[model.fs.unit.eq_electricity_consumption[0]] == pytest.approx( + 1, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_hydraulic_retention[0]] == pytest.approx( + 1.14755273e-6, rel=1e-8 + ) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_I"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_S"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_P"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_BH"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_BA"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_ND"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "H2O"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_I"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_S"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_O"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_NO"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_NH"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_ND"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_ALK"] + ] == pytest.approx(1, rel=1e-8) + + @pytest.mark.component + def test_scale_model(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, DewatererScaler) + + scaler.scale_model(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.mixed_state[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.mixed_state[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_in[model.fs.unit.mixed_state[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_in[model.fs.unit.mixed_state[0].temperature] == pytest.approx( + 1e-1, rel=1e-8 + ) + + # Outlet state - should be the same as the inlet + sfx_underflow = model.fs.unit.underflow_state[0].scaling_factor + assert isinstance(sfx_underflow, Suffix) + assert len(sfx_underflow) == 3 + assert sfx_underflow[ + model.fs.unit.underflow_state[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_underflow[ + model.fs.unit.underflow_state[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_underflow[ + model.fs.unit.underflow_state[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + sfx_overflow = model.fs.unit.underflow_state[0].scaling_factor + assert isinstance(sfx_overflow, Suffix) + assert len(sfx_overflow) == 3 + assert sfx_overflow[model.fs.unit.underflow_state[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_overflow[model.fs.unit.underflow_state[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_overflow[ + model.fs.unit.underflow_state[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Check that unit model has scaling factors + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 63 + assert sfx_unit[model.fs.unit.volume[0]] == pytest.approx(1e-3, rel=1e-3) + assert sfx_unit[model.fs.unit.eq_electricity_consumption[0]] == pytest.approx( + 1, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_hydraulic_retention[0]] == pytest.approx( + 2.06559491e-6, rel=1e-8 + ) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_I"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_S"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_P"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_BH"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_BA"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_ND"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "H2O"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_I"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_S"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_O"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_NO"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_NH"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_ND"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_ALK"] + ] == pytest.approx(1, 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 = ASM1ParameterBlock() + + m.fs.unit = DewateringUnit(property_package=m.fs.props) + + m.fs.unit.inlet.flow_vol.fix(178.4674 * 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(130.867 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(258.5789 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix( + 17216.2434 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(2611.4843 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(626.0652 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix( + 1442.7882 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.54323 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(100.8668 * units.mg / units.liter) + m.fs.unit.inlet.alkalinity.fix(97.8459 * units.mol / units.m**3) + + m.fs.unit.hydraulic_retention_time.fix() + + 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( + 3.3758830009e5, rel=1e-3 + ) + + @pytest.mark.integration + def test_example_case_scaler(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = ASM1ParameterBlock() + + m.fs.unit = DewateringUnit(property_package=m.fs.props) + + m.fs.unit.inlet.flow_vol.fix(178.4674 * 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(130.867 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(258.5789 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix( + 17216.2434 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(2611.4843 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(626.0652 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(1e-6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix( + 1442.7882 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.54323 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(100.8668 * units.mg / units.liter) + m.fs.unit.inlet.alkalinity.fix(97.8459 * units.mol / units.m**3) + + m.fs.unit.hydraulic_retention_time.fix() + + sb = ScalerBase() + + scaler = DewatererScaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.mixed_state: ASM1PropertiesScaler, + m.fs.unit.underflow_state: ASM1PropertiesScaler, + m.fs.unit.overflow_state: ASM1PropertiesScaler, + }, + ) + + # 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( + 2.100583344e4, rel=1e-3 + ) diff --git a/watertap/unit_models/tests/test_thickener_unit.py b/watertap/unit_models/tests/test_thickener_unit.py index d64cb3ae00..7ef024751d 100644 --- a/watertap/unit_models/tests/test_thickener_unit.py +++ b/watertap/unit_models/tests/test_thickener_unit.py @@ -19,6 +19,9 @@ from pyomo.environ import ( ConcreteModel, units, + Suffix, + TransformationFactory, + Var, ) from idaes.core import FlowsheetBlock @@ -28,14 +31,23 @@ from watertap.core.solvers import get_solver import idaes.core.util.scaling as iscale - +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) +from idaes.core.scaling.scaling_base import ScalerBase from idaes.core.util.exceptions import ( ConfigurationError, ) -from watertap.unit_models.thickener import Thickener, ActivatedSludgeModelType +from watertap.unit_models.thickener import ( + Thickener, + ActivatedSludgeModelType, + ThickenerScaler, +) from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( ASM1ParameterBlock, + ASM1PropertiesScaler, ) from watertap.property_models.unit_specific.activated_sludge.asm2d_properties import ( @@ -429,3 +441,402 @@ def test_list_error(): m.fs.unit = Thickener( property_package=m.fs.props, outlet_list=["outlet1", "outlet2"] ) + + +class TestThickenerScaler: + @pytest.fixture + def model(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = ASM1ParameterBlock() + + m.fs.unit = Thickener( + property_package=m.fs.props, + outlet_list=["underflow", "overflow"], + split_basis=SplittingType.componentFlow, + ) + + m.fs.unit.inlet.flow_vol.fix(300 * 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(28.0643 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(0.67336 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(3036.2175 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(63.2392 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix( + 4442.8377 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(332.5958 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(1922.8108 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(1.3748 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(9.1948 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(0.15845 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.55943 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(4.7411 * units.mg / units.liter) + m.fs.unit.inlet.alkalinity.fix(4.5646 * units.mol / units.m**3) + + m.fs.unit.hydraulic_retention_time.fix() + m.fs.unit.diameter.fix() + + return m + + @pytest.mark.component + def test_variable_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ThickenerScaler) + + scaler.variable_scaling_routine(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.mixed_state[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.mixed_state[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_in[model.fs.unit.mixed_state[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_in[model.fs.unit.mixed_state[0].temperature] == pytest.approx( + 1e-1, rel=1e-8 + ) + + # Outlet state - should be the same as the inlet + sfx_underflow = model.fs.unit.underflow_state[0].scaling_factor + assert isinstance(sfx_underflow, Suffix) + assert len(sfx_underflow) == 3 + assert sfx_underflow[ + model.fs.unit.underflow_state[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_underflow[ + model.fs.unit.underflow_state[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_underflow[ + model.fs.unit.underflow_state[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + sfx_overflow = model.fs.unit.overflow_state[0].scaling_factor + assert isinstance(sfx_overflow, Suffix) + assert len(sfx_overflow) == 3 + assert sfx_overflow[model.fs.unit.overflow_state[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_overflow[model.fs.unit.overflow_state[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_overflow[ + model.fs.unit.overflow_state[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Check that unit model has scaling factors + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 1 + assert sfx_unit[model.fs.unit.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, ThickenerScaler) + + scaler.constraint_scaling_routine(model.fs.unit) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 63 + assert sfx_unit[model.fs.unit.eq_electricity_consumption[0]] == pytest.approx( + 1, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_hydraulic_retention[0]] == pytest.approx( + 5.5555555556e-4, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_volume[0]] == pytest.approx( + 5.5555555556e-4, rel=1e-8 + ) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_I"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_S"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_P"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_BH"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_BA"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_ND"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "H2O"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_I"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_S"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_O"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_NO"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_NH"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_ND"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_ALK"] + ] == pytest.approx(1, rel=1e-8) + + @pytest.mark.component + def test_scale_model(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ThickenerScaler) + + scaler.scale_model(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.mixed_state[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.mixed_state[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_in[model.fs.unit.mixed_state[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_in[model.fs.unit.mixed_state[0].temperature] == pytest.approx( + 1e-1, rel=1e-8 + ) + + # Outlet state - should be the same as the inlet + sfx_underflow = model.fs.unit.underflow_state[0].scaling_factor + assert isinstance(sfx_underflow, Suffix) + assert len(sfx_underflow) == 3 + assert sfx_underflow[ + model.fs.unit.underflow_state[0].flow_vol + ] == pytest.approx(1e1, rel=1e-8) + assert sfx_underflow[ + model.fs.unit.underflow_state[0].pressure + ] == pytest.approx(1e-6, rel=1e-8) + assert sfx_underflow[ + model.fs.unit.underflow_state[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + sfx_overflow = model.fs.unit.underflow_state[0].scaling_factor + assert isinstance(sfx_overflow, Suffix) + assert len(sfx_overflow) == 3 + assert sfx_overflow[model.fs.unit.underflow_state[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_overflow[model.fs.unit.underflow_state[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_overflow[ + model.fs.unit.underflow_state[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + # Check that unit model has scaling factors + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 64 + assert sfx_unit[model.fs.unit.volume[0]] == pytest.approx(1e-3, rel=1e-3) + assert sfx_unit[model.fs.unit.eq_electricity_consumption[0]] == pytest.approx( + 1, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_hydraulic_retention[0]] == pytest.approx( + 0.001, rel=1e-8 + ) + assert sfx_unit[model.fs.unit.eq_volume[0]] == pytest.approx(0.001, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_I"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_S"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_P"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_BH"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_BA"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.overflow_particulate_fraction[0, "X_ND"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "H2O"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_I"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_S"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_O"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_NO"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_NH"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_ND"] + ] == pytest.approx(1, rel=1e-8) + assert sfx_unit[ + model.fs.unit.non_particulate_components[0, "S_ALK"] + ] == pytest.approx(1, 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 = ASM1ParameterBlock() + + m.fs.unit = Thickener( + property_package=m.fs.props, + outlet_list=["underflow", "overflow"], + split_basis=SplittingType.componentFlow, + ) + + m.fs.unit.inlet.flow_vol.fix(300 * 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(28.0643 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(0.67336 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(3036.2175 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(63.2392 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix( + 4442.8377 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(332.5958 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(1922.8108 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(1.3748 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(9.1948 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(0.15845 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.55943 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(4.7411 * units.mg / units.liter) + m.fs.unit.inlet.alkalinity.fix(4.5646 * units.mol / units.m**3) + + m.fs.unit.hydraulic_retention_time.fix() + m.fs.unit.diameter.fix() + + # Set scaling factors for badly scaled variables + iscale.set_scaling_factor(m.fs.unit.underflow_state[0.0].flow_vol, 1e4) + iscale.set_scaling_factor(m.fs.unit.underflow_state[0.0].pressure, 1e-6) + iscale.set_scaling_factor( + m.fs.unit.underflow_state[0.0].conc_mass_comp["S_S"], 1e4 + ) + iscale.set_scaling_factor( + m.fs.unit.underflow_state[0.0].conc_mass_comp["S_NH"], 1e4 + ) + iscale.set_scaling_factor( + m.fs.unit.underflow_state[0.0].conc_mass_comp["S_ND"], 1e4 + ) + iscale.set_scaling_factor(m.fs.unit.overflow_state[0.0].pressure, 1e-6) + iscale.set_scaling_factor( + m.fs.unit.overflow_state[0.0].conc_mass_comp["S_S"], 1e4 + ) + iscale.set_scaling_factor( + m.fs.unit.overflow_state[0.0].conc_mass_comp["S_NH"], 1e4 + ) + iscale.set_scaling_factor( + m.fs.unit.overflow_state[0.0].conc_mass_comp["S_ND"], 1e4 + ) + iscale.set_scaling_factor( + m.fs.unit.overflow_state[0.0].conc_mass_comp["X_ND"], 1e4 + ) + + 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.63071968e8, rel=1e-3 + ) + + @pytest.mark.integration + def test_example_case_scaler(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = ASM1ParameterBlock() + + m.fs.unit = Thickener( + property_package=m.fs.props, + outlet_list=["underflow", "overflow"], + split_basis=SplittingType.componentFlow, + ) + + m.fs.unit.inlet.flow_vol.fix(300 * 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(28.0643 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(0.67336 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(3036.2175 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(63.2392 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix( + 4442.8377 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(332.5958 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(1922.8108 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(1.3748 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(9.1948 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(0.15845 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(0.55943 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(4.7411 * units.mg / units.liter) + m.fs.unit.inlet.alkalinity.fix(4.5646 * units.mol / units.m**3) + + m.fs.unit.hydraulic_retention_time.fix() + m.fs.unit.diameter.fix() + + sb = ScalerBase() + + # Apply scaling to poorly scaled variables + for var in m.fs.component_data_objects(Var, descend_into=True): + if "flow_vol" in var.name: + sb.set_variable_scaling_factor(var, 1) + if "conc_mass_comp" in var.name: + sb.set_variable_scaling_factor(var, 1e1) + + scaler = ThickenerScaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.mixed_state: ASM1PropertiesScaler, + m.fs.unit.underflow_state: ASM1PropertiesScaler, + m.fs.unit.overflow_state: ASM1PropertiesScaler, + }, + ) + + # 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( + 4.8463788512e2, rel=1e-3 + ) diff --git a/watertap/unit_models/thickener.py b/watertap/unit_models/thickener.py index 6d3f818f37..ed9790770f 100644 --- a/watertap/unit_models/thickener.py +++ b/watertap/unit_models/thickener.py @@ -34,6 +34,7 @@ import idaes.logger as idaeslog from pyomo.environ import ( + Constraint, Param, Var, NonNegativeReals, @@ -44,6 +45,7 @@ from idaes.core.util.exceptions import ( ConfigurationError, ) +from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme from watertap.costing.unit_models.thickener import cost_thickener __author__ = "Alejandro Garciadiego, Adam Atia" @@ -65,12 +67,116 @@ class ActivatedSludgeModelType(Enum): modified_ASM2D = auto() +class ThickenerScaler(CustomScalerBase): + """ + Default modular scaler for the thickener 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.mixed_state, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.propagate_state_scaling( + target_state=model.underflow_state, + source_state=model.mixed_state, + overwrite=overwrite, + ) + self.propagate_state_scaling( + target_state=model.overflow_state, + source_state=model.mixed_state, + overwrite=overwrite, + ) + + self.call_submodel_scaler_method( + submodel=model.underflow_state, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.overflow_state, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level variables + self.scale_variable_by_default(model.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.mixed_state, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.underflow_state, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.overflow_state, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level constraints + # Constraint( + for c in model.component_data_objects(Constraint, descend_into=False): + self.scale_constraint_by_nominal_value( + c, + scheme=ConstraintScalingScheme.inverseMaximum, + overwrite=overwrite, + ) + + @declare_process_block_class("Thickener") class ThickenerData(SeparatorData): """ Thickener unit model for BSM2 """ + default_scaler = ThickenerScaler + CONFIG = SeparatorData.CONFIG() CONFIG.outlet_list = ["underflow", "overflow"] CONFIG.split_basis = SplittingType.componentFlow