diff --git a/watertap/unit_models/translators/tests/test_translator_adm1_asm1.py b/watertap/unit_models/translators/tests/test_translator_adm1_asm1.py index 50c3311f2f..70738d3e6b 100644 --- a/watertap/unit_models/translators/tests/test_translator_adm1_asm1.py +++ b/watertap/unit_models/translators/tests/test_translator_adm1_asm1.py @@ -24,9 +24,12 @@ ConcreteModel, value, assert_optimal_termination, + Suffix, + TransformationFactory, ) from idaes.core import FlowsheetBlock +import idaes.core.util.scaling as iscale from pyomo.environ import units @@ -39,21 +42,29 @@ ) from idaes.core.util.testing import initialization_tester +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) -from watertap.unit_models.translators.translator_adm1_asm1 import Translator_ADM1_ASM1 +from watertap.unit_models.translators.translator_adm1_asm1 import ( + Translator_ADM1_ASM1, + ADM1ASM1Scaler, +) from watertap.property_models.unit_specific.anaerobic_digestion.adm1_properties import ( ADM1ParameterBlock, + ADM1PropertiesScaler, ) from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( ASM1ParameterBlock, + ASM1PropertiesScaler, ) from watertap.property_models.unit_specific.anaerobic_digestion.adm1_reactions import ( ADM1ReactionParameterBlock, ) - from pyomo.util.check_units import assert_units_consistent # ----------------------------------------------------------------------------- @@ -288,3 +299,306 @@ def test_conservation(self, admasm): ) <= 1e-6 ) + + +class TestADM1ASM1Scaler: + @pytest.fixture + def model(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM1 = ASM1ParameterBlock() + m.fs.props_ADM1 = ADM1ParameterBlock() + m.fs.ADM1_rxn_props = ADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + + m.fs.unit = Translator_ADM1_ASM1( + inlet_property_package=m.fs.props_ADM1, + outlet_property_package=m.fs.props_ASM1, + reaction_package=m.fs.ADM1_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + ) + + m.fs.unit.inlet.flow_vol.fix(170 * 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_su"].fix(12.394 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_aa"].fix(5.54 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_fa"].fix(107.41 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_va"].fix(12.33 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_bu"].fix(14.00 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_pro"].fix(17.584 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ac"].fix(89.315 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_h2"].fix(2.55e-4 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ch4"].fix(55.49 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_IC"].fix( + 95.149 * units.mmol / units.liter * 12 * units.mg / units.mmol + ) + m.fs.unit.inlet.conc_mass_comp[0, "S_IN"].fix( + 94.468 * units.mmol / units.liter * 14 * units.mg / units.mmol + ) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(130.87 * units.mg / units.liter) + + m.fs.unit.inlet.conc_mass_comp[0, "X_c"].fix(107.92 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ch"].fix(20.517 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_pr"].fix(84.22 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_li"].fix(43.629 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_su"].fix(312.22 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_aa"].fix(931.67 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_fa"].fix(338.39 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_c4"].fix(335.77 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_pro"].fix(101.12 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ac"].fix(677.24 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_h2"].fix(284.84 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(17216 * units.mg / units.liter) + + m.fs.unit.inlet.cations[0].fix(1e-8 * units.mmol / units.liter) + m.fs.unit.inlet.anions[0].fix(5.21 * units.mmol / units.liter) + + return m + + @pytest.mark.component + def test_variable_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ADM1ASM1Scaler) + + scaler.variable_scaling_routine(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.properties_in[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].temperature] == pytest.approx( + 1e-1, rel=1e-8 + ) + + # Outlet state - should be the same as the inlet + sfx_underflow = model.fs.unit.properties_out[0].scaling_factor + assert isinstance(sfx_underflow, Suffix) + assert len(sfx_underflow) == 3 + assert sfx_underflow[model.fs.unit.properties_out[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_underflow[model.fs.unit.properties_out[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_underflow[ + model.fs.unit.properties_out[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + @pytest.mark.component + def test_constraint_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ADM1ASM1Scaler) + + scaler.constraint_scaling_routine(model.fs.unit) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 16 + + @pytest.mark.component + def test_scale_model(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ADM1ASM1Scaler) + + scaler.scale_model(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.properties_in[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].temperature] == pytest.approx( + 1e-1, rel=1e-8 + ) + + # Outlet state - should be the same as the inlet + sfx_underflow = model.fs.unit.properties_out[0].scaling_factor + assert isinstance(sfx_underflow, Suffix) + assert len(sfx_underflow) == 3 + assert sfx_underflow[model.fs.unit.properties_out[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_underflow[model.fs.unit.properties_out[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_underflow[ + model.fs.unit.properties_out[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + sfx_overflow = model.fs.unit.properties_out[0].scaling_factor + assert isinstance(sfx_overflow, Suffix) + assert len(sfx_overflow) == 3 + assert sfx_overflow[model.fs.unit.properties_out[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_overflow[model.fs.unit.properties_out[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_overflow[ + model.fs.unit.properties_out[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) == 16 + + # 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.props_ADM1 = ADM1ParameterBlock() + m.fs.ADM1_rxn_props = ADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + + m.fs.unit = Translator_ADM1_ASM1( + inlet_property_package=m.fs.props_ADM1, + outlet_property_package=m.fs.props_ASM1, + reaction_package=m.fs.ADM1_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + ) + + m.fs.unit.inlet.flow_vol.fix(170 * 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_su"].fix(12.394 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_aa"].fix(5.54 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_fa"].fix(107.41 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_va"].fix(12.33 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_bu"].fix(14.00 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_pro"].fix(17.584 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ac"].fix(89.315 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_h2"].fix(2.55e-4 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ch4"].fix(55.49 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_IC"].fix( + 95.149 * units.mmol / units.liter * 12 * units.mg / units.mmol + ) + m.fs.unit.inlet.conc_mass_comp[0, "S_IN"].fix( + 94.468 * units.mmol / units.liter * 14 * units.mg / units.mmol + ) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(130.87 * units.mg / units.liter) + + m.fs.unit.inlet.conc_mass_comp[0, "X_c"].fix(107.92 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ch"].fix(20.517 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_pr"].fix(84.22 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_li"].fix(43.629 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_su"].fix(312.22 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_aa"].fix(931.67 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_fa"].fix(338.39 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_c4"].fix(335.77 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_pro"].fix(101.12 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ac"].fix(677.24 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_h2"].fix(284.84 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(17216 * units.mg / units.liter) + + m.fs.unit.inlet.cations[0].fix(1e-8 * units.mmol / units.liter) + m.fs.unit.inlet.anions[0].fix(5.21 * units.mmol / units.liter) + + # Check condition number to confirm scaling + iscale.calculate_scaling_factors(m) + + 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.872983349e5, 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.props_ADM1 = ADM1ParameterBlock() + m.fs.ADM1_rxn_props = ADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + + m.fs.unit = Translator_ADM1_ASM1( + inlet_property_package=m.fs.props_ADM1, + outlet_property_package=m.fs.props_ASM1, + reaction_package=m.fs.ADM1_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + ) + + m.fs.unit.inlet.flow_vol.fix(170 * 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_su"].fix(12.394 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_aa"].fix(5.54 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_fa"].fix(107.41 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_va"].fix(12.33 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_bu"].fix(14.00 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_pro"].fix(17.584 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ac"].fix(89.315 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_h2"].fix(2.55e-4 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ch4"].fix(55.49 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_IC"].fix( + 95.149 * units.mmol / units.liter * 12 * units.mg / units.mmol + ) + m.fs.unit.inlet.conc_mass_comp[0, "S_IN"].fix( + 94.468 * units.mmol / units.liter * 14 * units.mg / units.mmol + ) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(130.87 * units.mg / units.liter) + + m.fs.unit.inlet.conc_mass_comp[0, "X_c"].fix(107.92 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ch"].fix(20.517 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_pr"].fix(84.22 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_li"].fix(43.629 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_su"].fix(312.22 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_aa"].fix(931.67 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_fa"].fix(338.39 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_c4"].fix(335.77 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_pro"].fix(101.12 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ac"].fix(677.24 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_h2"].fix(284.84 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(17216 * units.mg / units.liter) + + m.fs.unit.inlet.cations[0].fix(1e-8 * units.mmol / units.liter) + m.fs.unit.inlet.anions[0].fix(5.21 * units.mmol / units.liter) + + scaler = ADM1ASM1Scaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.properties_in: ADM1PropertiesScaler, + m.fs.unit.properties_out: 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( + 5.4739117e3, rel=1e-3 + ) diff --git a/watertap/unit_models/translators/tests/test_translator_adm1_asm2d.py b/watertap/unit_models/translators/tests/test_translator_adm1_asm2d.py index e4b532e45f..203585086f 100644 --- a/watertap/unit_models/translators/tests/test_translator_adm1_asm2d.py +++ b/watertap/unit_models/translators/tests/test_translator_adm1_asm2d.py @@ -22,6 +22,8 @@ ConcreteModel, value, assert_optimal_termination, + Suffix, + TransformationFactory, ) from idaes.core import ( @@ -39,18 +41,26 @@ number_total_constraints, number_unused_variables, ) +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) +import idaes.core.util.scaling as iscale from idaes.core.util.testing import initialization_tester from watertap.unit_models.translators.translator_adm1_asm2d import ( Translator_ADM1_ASM2D, + ADM1ASM2dScaler, ) from watertap.property_models.unit_specific.anaerobic_digestion.modified_adm1_properties import ( ModifiedADM1ParameterBlock, + ModifiedADM1PropertiesScaler, ) from watertap.property_models.unit_specific.activated_sludge.modified_asm2d_properties import ( ModifiedASM2dParameterBlock, + ModifiedASM2dPropertiesScaler, ) from watertap.property_models.unit_specific.anaerobic_digestion.modified_adm1_reactions import ( @@ -324,3 +334,303 @@ def test_conservation(self, asmadm): ) <= 1e-6 ) + + +class TestADM1ASM2dScaler: + @pytest.fixture + def model(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM2D = ModifiedASM2dParameterBlock() + m.fs.ASM2d_rxn_props = ModifiedASM2dReactionParameterBlock( + property_package=m.fs.props_ASM2D + ) + m.fs.props_ADM1 = ModifiedADM1ParameterBlock() + m.fs.ADM1_rxn_props = ModifiedADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + + m.fs.unit = Translator_ADM1_ASM2D( + inlet_property_package=m.fs.props_ADM1, + outlet_property_package=m.fs.props_ASM2D, + inlet_reaction_package=m.fs.ADM1_rxn_props, + outlet_reaction_package=m.fs.ASM2d_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + ) + + m.fs.unit.inlet.flow_vol.fix(170 * 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_su"].fix(0.034597) + m.fs.unit.inlet.conc_mass_comp[0, "S_aa"].fix(0.015037) + m.fs.unit.inlet.conc_mass_comp[0, "S_fa"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_va"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_bu"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_pro"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_ac"].fix(0.025072) + m.fs.unit.inlet.conc_mass_comp[0, "S_h2"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_ch4"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_IC"].fix(0.34628) + m.fs.unit.inlet.conc_mass_comp[0, "S_IN"].fix(0.60014) + m.fs.unit.inlet.conc_mass_comp[0, "S_IP"].fix(0.22677) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(0.026599) + + m.fs.unit.inlet.conc_mass_comp[0, "X_ch"].fix(7.3687) + m.fs.unit.inlet.conc_mass_comp[0, "X_pr"].fix(7.7308) + m.fs.unit.inlet.conc_mass_comp[0, "X_li"].fix(10.3288) + m.fs.unit.inlet.conc_mass_comp[0, "X_su"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_aa"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_fa"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_c4"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_pro"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_ac"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_h2"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(12.7727) + m.fs.unit.inlet.conc_mass_comp[0, "X_PHA"].fix(0.0022493) + m.fs.unit.inlet.conc_mass_comp[0, "X_PP"].fix(1.04110) + m.fs.unit.inlet.conc_mass_comp[0, "X_PAO"].fix(3.4655) + m.fs.unit.inlet.conc_mass_comp[0, "S_K"].fix(0.02268) + m.fs.unit.inlet.conc_mass_comp[0, "S_Mg"].fix(0.02893) + + m.fs.unit.inlet.cations[0].fix(0.04) + m.fs.unit.inlet.anions[0].fix(0.02) + + return m + + @pytest.mark.component + def test_variable_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ADM1ASM2dScaler) + + scaler.variable_scaling_routine(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.properties_in[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_in[model.fs.unit.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.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 3 + assert sfx_out[model.fs.unit.properties_out[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_out[model.fs.unit.properties_out[0].pressure] == pytest.approx( + 1e-5, rel=1e-8 + ) + assert sfx_out[model.fs.unit.properties_out[0].temperature] == pytest.approx( + 1e-2, rel=1e-8 + ) + + @pytest.mark.component + def test_constraint_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ADM1ASM2dScaler) + + scaler.constraint_scaling_routine(model.fs.unit) + + sfx_out = model.fs.unit.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 1 + + @pytest.mark.component + def test_scale_model(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ADM1ASM2dScaler) + + scaler.scale_model(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.properties_in[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_in[model.fs.unit.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.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 3 + assert sfx_out[model.fs.unit.properties_out[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_out[model.fs.unit.properties_out[0].pressure] == pytest.approx( + 1e-5, rel=1e-8 + ) + assert sfx_out[model.fs.unit.properties_out[0].temperature] == pytest.approx( + 1e-2, 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_ASM2D = ModifiedASM2dParameterBlock() + m.fs.ASM2d_rxn_props = ModifiedASM2dReactionParameterBlock( + property_package=m.fs.props_ASM2D + ) + m.fs.props_ADM1 = ModifiedADM1ParameterBlock() + m.fs.ADM1_rxn_props = ModifiedADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + + m.fs.unit = Translator_ADM1_ASM2D( + inlet_property_package=m.fs.props_ADM1, + outlet_property_package=m.fs.props_ASM2D, + inlet_reaction_package=m.fs.ADM1_rxn_props, + outlet_reaction_package=m.fs.ASM2d_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + ) + + m.fs.unit.inlet.flow_vol.fix(170 * 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_su"].fix(0.034597) + m.fs.unit.inlet.conc_mass_comp[0, "S_aa"].fix(0.015037) + m.fs.unit.inlet.conc_mass_comp[0, "S_fa"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_va"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_bu"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_pro"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_ac"].fix(0.025072) + m.fs.unit.inlet.conc_mass_comp[0, "S_h2"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_ch4"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_IC"].fix(0.34628) + m.fs.unit.inlet.conc_mass_comp[0, "S_IN"].fix(0.60014) + m.fs.unit.inlet.conc_mass_comp[0, "S_IP"].fix(0.22677) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(0.026599) + + m.fs.unit.inlet.conc_mass_comp[0, "X_ch"].fix(7.3687) + m.fs.unit.inlet.conc_mass_comp[0, "X_pr"].fix(7.7308) + m.fs.unit.inlet.conc_mass_comp[0, "X_li"].fix(10.3288) + m.fs.unit.inlet.conc_mass_comp[0, "X_su"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_aa"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_fa"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_c4"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_pro"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_ac"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_h2"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(12.7727) + m.fs.unit.inlet.conc_mass_comp[0, "X_PHA"].fix(0.0022493) + m.fs.unit.inlet.conc_mass_comp[0, "X_PP"].fix(1.04110) + m.fs.unit.inlet.conc_mass_comp[0, "X_PAO"].fix(3.4655) + m.fs.unit.inlet.conc_mass_comp[0, "S_K"].fix(0.02268) + m.fs.unit.inlet.conc_mass_comp[0, "S_Mg"].fix(0.02893) + + m.fs.unit.inlet.cations[0].fix(0.04) + m.fs.unit.inlet.anions[0].fix(0.02) + + 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( + 4.47213596e5, rel=1e-3 + ) + + @pytest.mark.integration + def test_example_case_scaler(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM2D = ModifiedASM2dParameterBlock() + m.fs.ASM2d_rxn_props = ModifiedASM2dReactionParameterBlock( + property_package=m.fs.props_ASM2D + ) + m.fs.props_ADM1 = ModifiedADM1ParameterBlock() + m.fs.ADM1_rxn_props = ModifiedADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + + m.fs.unit = Translator_ADM1_ASM2D( + inlet_property_package=m.fs.props_ADM1, + outlet_property_package=m.fs.props_ASM2D, + inlet_reaction_package=m.fs.ADM1_rxn_props, + outlet_reaction_package=m.fs.ASM2d_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + ) + + m.fs.unit.inlet.flow_vol.fix(170 * 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_su"].fix(0.034597) + m.fs.unit.inlet.conc_mass_comp[0, "S_aa"].fix(0.015037) + m.fs.unit.inlet.conc_mass_comp[0, "S_fa"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_va"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_bu"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_pro"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_ac"].fix(0.025072) + m.fs.unit.inlet.conc_mass_comp[0, "S_h2"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_ch4"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "S_IC"].fix(0.34628) + m.fs.unit.inlet.conc_mass_comp[0, "S_IN"].fix(0.60014) + m.fs.unit.inlet.conc_mass_comp[0, "S_IP"].fix(0.22677) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(0.026599) + + m.fs.unit.inlet.conc_mass_comp[0, "X_ch"].fix(7.3687) + m.fs.unit.inlet.conc_mass_comp[0, "X_pr"].fix(7.7308) + m.fs.unit.inlet.conc_mass_comp[0, "X_li"].fix(10.3288) + m.fs.unit.inlet.conc_mass_comp[0, "X_su"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_aa"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_fa"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_c4"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_pro"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_ac"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_h2"].fix(0) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(12.7727) + m.fs.unit.inlet.conc_mass_comp[0, "X_PHA"].fix(0.0022493) + m.fs.unit.inlet.conc_mass_comp[0, "X_PP"].fix(1.04110) + m.fs.unit.inlet.conc_mass_comp[0, "X_PAO"].fix(3.4655) + m.fs.unit.inlet.conc_mass_comp[0, "S_K"].fix(0.02268) + m.fs.unit.inlet.conc_mass_comp[0, "S_Mg"].fix(0.02893) + + m.fs.unit.inlet.cations[0].fix(0.04) + m.fs.unit.inlet.anions[0].fix(0.02) + + scaler = ADM1ASM2dScaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.properties_in: ModifiedADM1PropertiesScaler, + m.fs.unit.properties_out: ModifiedASM2dPropertiesScaler, + }, + ) + + # 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( + 7.180968e3, rel=1e-3 + ) diff --git a/watertap/unit_models/translators/tests/test_translator_asm1_adm1.py b/watertap/unit_models/translators/tests/test_translator_asm1_adm1.py index ff4eaafaa0..64c97e74ba 100644 --- a/watertap/unit_models/translators/tests/test_translator_asm1_adm1.py +++ b/watertap/unit_models/translators/tests/test_translator_asm1_adm1.py @@ -20,7 +20,14 @@ """ import pytest -from pyomo.environ import ConcreteModel, value, assert_optimal_termination, Param +from pyomo.environ import ( + ConcreteModel, + value, + assert_optimal_termination, + Param, + Suffix, + TransformationFactory, +) from idaes.core import FlowsheetBlock import idaes.core.util.scaling as iscale @@ -36,14 +43,24 @@ import idaes.logger as idaeslog from idaes.core.util.testing import initialization_tester +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) +from idaes.core.scaling.scaling_base import ScalerBase -from watertap.unit_models.translators.translator_asm1_adm1 import Translator_ASM1_ADM1 +from watertap.unit_models.translators.translator_asm1_adm1 import ( + Translator_ASM1_ADM1, + ASM1ADM1Scaler, +) from watertap.property_models.unit_specific.anaerobic_digestion.adm1_properties import ( ADM1ParameterBlock, + ADM1PropertiesScaler, ) from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( ASM1ParameterBlock, + ASM1PropertiesScaler, ) from watertap.property_models.unit_specific.anaerobic_digestion.adm1_reactions import ( @@ -436,3 +453,291 @@ def test_intermidiates(self, asmadm): assert pytest.approx(0.040699, rel=1e-3) == value(asmadm.fs.unit.ReqOrgNx[0]) assert pytest.approx(67.5251, rel=1e-3) == value(asmadm.fs.unit.ReqCODXc[0]) assert pytest.approx(0.04779, rel=1e-3) == value(asmadm.fs.unit.ReqCODs[0]) + + +class TestASM1ADM1Scaler: + @pytest.fixture + def model(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM1 = ASM1ParameterBlock() + m.fs.props_ADM1 = ADM1ParameterBlock() + m.fs.ADM1_rxn_props = ADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + + m.fs.unit = Translator_ASM1_ADM1( + inlet_property_package=m.fs.props_ASM1, + outlet_property_package=m.fs.props_ADM1, + reaction_package=m.fs.ADM1_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + ) + + 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(28.0665 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(48.9526 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix( + 10361.7101 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix( + 20375.0176 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix( + 10210.0698 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(553.2808 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(3204.6601 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0.25225 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(1.6871 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(28.9098 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(4.6834 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(906.0933 * units.mg / units.liter) + m.fs.unit.inlet.alkalinity.fix(7.1549 * units.mol / units.m**3) + + return m + + @pytest.mark.component + def test_variable_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ASM1ADM1Scaler) + + scaler.variable_scaling_routine(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.properties_in[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].temperature] == pytest.approx( + 1e-1, rel=1e-8 + ) + + # Outlet state - should be the same as the inlet + sfx_underflow = model.fs.unit.properties_out[0].scaling_factor + assert isinstance(sfx_underflow, Suffix) + assert len(sfx_underflow) == 3 + assert sfx_underflow[model.fs.unit.properties_out[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_underflow[model.fs.unit.properties_out[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_underflow[ + model.fs.unit.properties_out[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + @pytest.mark.component + def test_constraint_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ASM1ADM1Scaler) + + scaler.constraint_scaling_routine(model.fs.unit) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 34 + + @pytest.mark.component + def test_scale_model(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ASM1ADM1Scaler) + + scaler.scale_model(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.properties_in[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].temperature] == pytest.approx( + 1e-1, rel=1e-8 + ) + + # Outlet state - should be the same as the inlet + sfx_underflow = model.fs.unit.properties_out[0].scaling_factor + assert isinstance(sfx_underflow, Suffix) + assert len(sfx_underflow) == 3 + assert sfx_underflow[model.fs.unit.properties_out[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_underflow[model.fs.unit.properties_out[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_underflow[ + model.fs.unit.properties_out[0].temperature + ] == pytest.approx(1e-1, rel=1e-8) + + sfx_overflow = model.fs.unit.properties_out[0].scaling_factor + assert isinstance(sfx_overflow, Suffix) + assert len(sfx_overflow) == 3 + assert sfx_overflow[model.fs.unit.properties_out[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_overflow[model.fs.unit.properties_out[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_overflow[ + model.fs.unit.properties_out[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) == 34 + + # TODO: Remove test once iscale is deprecated + @pytest.mark.integration + def test_example_case_iscale(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM1 = ASM1ParameterBlock() + m.fs.props_ADM1 = ADM1ParameterBlock() + m.fs.ADM1_rxn_props = ADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + + m.fs.unit = Translator_ASM1_ADM1( + inlet_property_package=m.fs.props_ASM1, + outlet_property_package=m.fs.props_ADM1, + reaction_package=m.fs.ADM1_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + ) + + 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(28.0665 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(48.9526 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix( + 10361.7101 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix( + 20375.0176 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix( + 10210.0698 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(553.2808 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(3204.6601 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0.25225 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(1.6871 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(28.9098 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(4.6834 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(906.0933 * units.mg / units.liter) + m.fs.unit.inlet.alkalinity.fix(7.1549 * units.mol / units.m**3) + + # Check condition number to confirm scaling + iscale.calculate_scaling_factors(m) + + iscale.set_scaling_factor(m.fs.unit.properties_out[0].flow_vol, 1) + iscale.set_scaling_factor(m.fs.unit.properties_out[0].conc_mass_comp["S_IC"], 1) + iscale.set_scaling_factor(m.fs.unit.properties_out[0].conc_mass_comp["S_IN"], 1) + iscale.set_scaling_factor(m.fs.unit.properties_out[0].conc_mass_comp["X_aa"], 1) + iscale.set_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["X_I"], 1e-1 + ) + + 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.00056509896e11, 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.props_ADM1 = ADM1ParameterBlock() + m.fs.ADM1_rxn_props = ADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + + m.fs.unit = Translator_ASM1_ADM1( + inlet_property_package=m.fs.props_ASM1, + outlet_property_package=m.fs.props_ADM1, + reaction_package=m.fs.ADM1_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + ) + + 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(28.0665 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(48.9526 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix( + 10361.7101 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix( + 20375.0176 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix( + 10210.0698 * units.mg / units.liter + ) + m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(553.2808 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(3204.6601 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(0.25225 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(1.6871 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(28.9098 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(4.6834 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(906.0933 * units.mg / units.liter) + m.fs.unit.inlet.alkalinity.fix(7.1549 * units.mol / units.m**3) + + scaler = ASM1ADM1Scaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.properties_in: ASM1PropertiesScaler, + m.fs.unit.properties_out: ADM1PropertiesScaler, + }, + ) + + sb = ScalerBase() + sb.set_variable_scaling_factor(m.fs.unit.properties_out[0].flow_vol, 1) + sb.set_variable_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["S_aa"], 1e3 + ) + sb.set_variable_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["S_h2"], 1e7 + ) + + # 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.2518e4, rel=1e-3 + ) diff --git a/watertap/unit_models/translators/tests/test_translator_asm2d_adm1.py b/watertap/unit_models/translators/tests/test_translator_asm2d_adm1.py index f525eb46e2..cceb3bff7b 100644 --- a/watertap/unit_models/translators/tests/test_translator_asm2d_adm1.py +++ b/watertap/unit_models/translators/tests/test_translator_asm2d_adm1.py @@ -25,6 +25,8 @@ value, assert_optimal_termination, Param, + Suffix, + TransformationFactory, ) from idaes.core import FlowsheetBlock @@ -38,18 +40,29 @@ number_total_constraints, number_unused_variables, ) +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) +import idaes.core.util.scaling as iscale +from idaes.core.scaling.scaling_base import ScalerBase import idaes.logger as idaeslog from idaes.core.util.testing import initialization_tester -from watertap.unit_models.translators.translator_asm2d_adm1 import Translator_ASM2d_ADM1 +from watertap.unit_models.translators.translator_asm2d_adm1 import ( + Translator_ASM2d_ADM1, + ASM2dADM1Scaler, +) from watertap.property_models.unit_specific.anaerobic_digestion.modified_adm1_properties import ( ModifiedADM1ParameterBlock, + ModifiedADM1PropertiesScaler, ) from watertap.property_models.unit_specific.activated_sludge.modified_asm2d_properties import ( ModifiedASM2dParameterBlock, + ModifiedASM2dPropertiesScaler, ) from watertap.property_models.unit_specific.anaerobic_digestion.modified_adm1_reactions import ( @@ -156,8 +169,6 @@ def asmadm(self): m.fs.unit.inlet.conc_mass_comp[0, "S_K"].fix(0.01979 * units.kg / units.m**3) m.fs.unit.inlet.conc_mass_comp[0, "S_Mg"].fix(0.18987 * units.kg / units.m**3) - # constraint_scaling_transform(m.fs.unit.SIC_output[0], 1e-3) - return m @pytest.mark.build @@ -634,3 +645,297 @@ def test_conservation(self, asmadm): ) <= 1e-6 ) + + +class TestADM1ASM2dScaler: + @pytest.fixture + def model(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM2d = ModifiedASM2dParameterBlock() + m.fs.props_ADM1 = ModifiedADM1ParameterBlock() + m.fs.ADM1_rxn_props = ModifiedADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + m.fs.ASM2d_rxn_props = ModifiedASM2dReactionParameterBlock( + property_package=m.fs.props_ASM2d, + ) + + m.fs.unit = Translator_ASM2d_ADM1( + inlet_property_package=m.fs.props_ASM2d, + outlet_property_package=m.fs.props_ADM1, + inlet_reaction_package=m.fs.ASM2d_rxn_props, + outlet_reaction_package=m.fs.ADM1_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + bio_P=True, + ) + + m.fs.unit.inlet.flow_vol.fix(18446 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + eps = 1e-9 * units.kg / units.m**3 + + m.fs.unit.inlet.conc_mass_comp[0, "S_O2"].fix(eps) + m.fs.unit.inlet.conc_mass_comp[0, "S_F"].fix(0.02644 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_A"].fix(0.01766 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(0.02723 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH4"].fix(0.01858 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_N2"].fix(0.00507 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO3"].fix(0.00002 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_PO4"].fix(0.00469 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_IC"].fix(0.07899 * units.kg / units.m**3) + + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(10.96441 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(19.08476 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_H"].fix(9.47939 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_PAO"].fix(3.8622 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_PP"].fix(0.45087 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_PHA"].fix(0.02464 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_AUT"].fix(0.33379 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_K"].fix(0.01979 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_Mg"].fix(0.18987 * units.kg / units.m**3) + + return m + + @pytest.mark.component + def test_variable_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ASM2dADM1Scaler) + + scaler.variable_scaling_routine(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.properties_in[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].pressure] == pytest.approx( + 1e-5, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].temperature] == pytest.approx( + 1e-2, rel=1e-8 + ) + + # Outlet state - should be the same as the inlet + sfx_out = model.fs.unit.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 3 + assert sfx_out[model.fs.unit.properties_out[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_out[model.fs.unit.properties_out[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_out[model.fs.unit.properties_out[0].temperature] == pytest.approx( + 1e-1, rel=1e-8 + ) + + @pytest.mark.component + def test_constraint_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ASM2dADM1Scaler) + + scaler.constraint_scaling_routine(model.fs.unit) + + sfx_out = model.fs.unit.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 0 + + @pytest.mark.component + def test_scale_model(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, ASM2dADM1Scaler) + + scaler.scale_model(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 3 + assert sfx_in[model.fs.unit.properties_in[0].flow_vol] == pytest.approx( + 1e1, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].pressure] == pytest.approx( + 1e-5, rel=1e-8 + ) + assert sfx_in[model.fs.unit.properties_in[0].temperature] == pytest.approx( + 1e-2, rel=1e-8 + ) + + # Outlet state - should be the same as the inlet + sfx_out = model.fs.unit.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 3 + assert sfx_out[model.fs.unit.properties_out[0].flow_vol] == pytest.approx( + 1e5, rel=1e-8 + ) + assert sfx_out[model.fs.unit.properties_out[0].pressure] == pytest.approx( + 1e-6, rel=1e-8 + ) + assert sfx_out[model.fs.unit.properties_out[0].temperature] == pytest.approx( + 1e-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_ASM2d = ModifiedASM2dParameterBlock() + m.fs.props_ADM1 = ModifiedADM1ParameterBlock() + m.fs.ADM1_rxn_props = ModifiedADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + m.fs.ASM2d_rxn_props = ModifiedASM2dReactionParameterBlock( + property_package=m.fs.props_ASM2d, + ) + + m.fs.unit = Translator_ASM2d_ADM1( + inlet_property_package=m.fs.props_ASM2d, + outlet_property_package=m.fs.props_ADM1, + inlet_reaction_package=m.fs.ASM2d_rxn_props, + outlet_reaction_package=m.fs.ADM1_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + bio_P=True, + ) + + m.fs.unit.inlet.flow_vol.fix(18446 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + eps = 1e-9 * units.kg / units.m**3 + + m.fs.unit.inlet.conc_mass_comp[0, "S_O2"].fix(eps) + m.fs.unit.inlet.conc_mass_comp[0, "S_F"].fix(0.02644 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_A"].fix(0.01766 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(0.02723 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH4"].fix(0.01858 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_N2"].fix(0.00507 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO3"].fix(0.00002 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_PO4"].fix(0.00469 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_IC"].fix(0.07899 * units.kg / units.m**3) + + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(10.96441 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(19.08476 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_H"].fix(9.47939 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_PAO"].fix(3.8622 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_PP"].fix(0.45087 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_PHA"].fix(0.02464 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_AUT"].fix(0.33379 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_K"].fix(0.01979 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_Mg"].fix(0.18987 * units.kg / units.m**3) + + iscale.set_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["S_aa"], 1e3 + ) + + 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( + 5.745705986257e3, rel=1e-3 + ) + + @pytest.mark.integration + def test_example_case_scaler(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props_ASM2d = ModifiedASM2dParameterBlock() + m.fs.props_ADM1 = ModifiedADM1ParameterBlock() + m.fs.ADM1_rxn_props = ModifiedADM1ReactionParameterBlock( + property_package=m.fs.props_ADM1 + ) + m.fs.ASM2d_rxn_props = ModifiedASM2dReactionParameterBlock( + property_package=m.fs.props_ASM2d, + ) + + m.fs.unit = Translator_ASM2d_ADM1( + inlet_property_package=m.fs.props_ASM2d, + outlet_property_package=m.fs.props_ADM1, + inlet_reaction_package=m.fs.ASM2d_rxn_props, + outlet_reaction_package=m.fs.ADM1_rxn_props, + has_phase_equilibrium=False, + outlet_state_defined=True, + bio_P=True, + ) + + m.fs.unit.inlet.flow_vol.fix(18446 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + eps = 1e-9 * units.kg / units.m**3 + + m.fs.unit.inlet.conc_mass_comp[0, "S_O2"].fix(eps) + m.fs.unit.inlet.conc_mass_comp[0, "S_F"].fix(0.02644 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_A"].fix(0.01766 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(0.02723 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH4"].fix(0.01858 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_N2"].fix(0.00507 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO3"].fix(0.00002 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_PO4"].fix(0.00469 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_IC"].fix(0.07899 * units.kg / units.m**3) + + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(10.96441 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(19.08476 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_H"].fix(9.47939 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_PAO"].fix(3.8622 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_PP"].fix(0.45087 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_PHA"].fix(0.02464 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "X_AUT"].fix(0.33379 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_K"].fix(0.01979 * units.kg / units.m**3) + m.fs.unit.inlet.conc_mass_comp[0, "S_Mg"].fix(0.18987 * units.kg / units.m**3) + + scaler = ASM2dADM1Scaler() + scaler.scale_model( + m.fs.unit, + submodel_scalers={ + m.fs.unit.properties_in: ModifiedASM2dPropertiesScaler, + m.fs.unit.properties_out: ModifiedADM1PropertiesScaler, + }, + ) + + sb = ScalerBase() + sb.set_variable_scaling_factor( + m.fs.unit.properties_out[0].flow_vol, 1, overwrite=True + ) + sb.set_variable_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["S_aa"], 1e3 + ) + sb.set_variable_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["S_ch4"], 1e9 + ) + sb.set_variable_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["X_su"], 1e9 + ) + sb.set_variable_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["X_fa"], 1e9 + ) + sb.set_variable_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["X_c4"], 1e9 + ) + sb.set_variable_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["X_pro"], 1e9 + ) + sb.set_variable_scaling_factor( + m.fs.unit.properties_out[0].conc_mass_comp["X_ac"], 1e9 + ) + + # 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.640234319533e3, rel=1e-3 + ) diff --git a/watertap/unit_models/translators/translator_adm1_asm1.py b/watertap/unit_models/translators/translator_adm1_asm1.py index 712b22bd1f..709794af06 100644 --- a/watertap/unit_models/translators/translator_adm1_asm1.py +++ b/watertap/unit_models/translators/translator_adm1_asm1.py @@ -35,10 +35,12 @@ from watertap.core.solvers import get_solver import idaes.logger as idaeslog import idaes.core.util.scaling as iscale +from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme from idaes.core.util.exceptions import InitializationError from pyomo.environ import ( + Constraint, Param, units as pyunits, check_optimal_termination, @@ -52,12 +54,87 @@ _log = idaeslog.getLogger(__name__) +class ADM1ASM1Scaler(CustomScalerBase): + """ + Default modular scaler for the ADM1-ASM1 translator block. + 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. + """ + + 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.properties_in, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + self.call_submodel_scaler_method( + submodel=model.properties_out, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + 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.properties_in, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.properties_out, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level constraints + 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("Translator_ADM1_ASM1") class TranslatorDataADM1ASM1(TranslatorData): """ Translator block representing the ADM1/ASM1 interface """ + default_scaler = ADM1ASM1Scaler + CONFIG = TranslatorData.CONFIG() CONFIG.declare( "reaction_package", diff --git a/watertap/unit_models/translators/translator_adm1_asm2d.py b/watertap/unit_models/translators/translator_adm1_asm2d.py index 410921792a..54815872f3 100644 --- a/watertap/unit_models/translators/translator_adm1_asm2d.py +++ b/watertap/unit_models/translators/translator_adm1_asm2d.py @@ -37,9 +37,13 @@ import idaes.logger as idaeslog import idaes.core.util.scaling as iscale +from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme + + from idaes.core.util.exceptions import InitializationError from pyomo.environ import ( + Constraint, units as pyunits, check_optimal_termination, Set, @@ -53,12 +57,86 @@ _log = idaeslog.getLogger(__name__) +class ADM1ASM2dScaler(CustomScalerBase): + """ + Default modular scaler for ADM1-ASM2d translator block. + 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. + """ + + 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.properties_in, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + self.call_submodel_scaler_method( + submodel=model.properties_out, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + 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.properties_in, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.properties_out, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level constraints + 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("Translator_ADM1_ASM2D") class TranslatorDataADM1ASM2D(TranslatorData): """ Translator block representing the ADM1/ASM2D interface """ + default_scaler = ADM1ASM2dScaler + CONFIG = TranslatorData.CONFIG() CONFIG.declare( diff --git a/watertap/unit_models/translators/translator_asm1_adm1.py b/watertap/unit_models/translators/translator_asm1_adm1.py index a3e0f084ba..6f620dd269 100644 --- a/watertap/unit_models/translators/translator_asm1_adm1.py +++ b/watertap/unit_models/translators/translator_asm1_adm1.py @@ -36,8 +36,10 @@ from watertap.core.solvers import get_solver import idaes.logger as idaeslog import idaes.core.util.scaling as iscale +from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme from pyomo.environ import ( + Constraint, Param, NonNegativeReals, Var, @@ -56,12 +58,87 @@ _log = idaeslog.getLogger(__name__) +class ASM1ADM1Scaler(CustomScalerBase): + """ + Default modular scaler for the ASM1-ADM1 translator block. + 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. + """ + + 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.properties_in, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + self.call_submodel_scaler_method( + submodel=model.properties_out, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + 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.properties_in, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.properties_out, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level constraints + 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("Translator_ASM1_ADM1") class TranslatorDataASM1ADM1(TranslatorData): """ Translator block representing the ASM1/ADM1 interface """ + default_scaler = ASM1ADM1Scaler + CONFIG = TranslatorData.CONFIG() CONFIG.declare( "reaction_package", diff --git a/watertap/unit_models/translators/translator_asm2d_adm1.py b/watertap/unit_models/translators/translator_asm2d_adm1.py index 7a114ae72f..921cf3b805 100644 --- a/watertap/unit_models/translators/translator_asm2d_adm1.py +++ b/watertap/unit_models/translators/translator_asm2d_adm1.py @@ -35,8 +35,10 @@ from idaes.core.util.model_statistics import degrees_of_freedom from watertap.core.solvers import get_solver import idaes.logger as idaeslog +from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme from pyomo.environ import ( + Constraint, Param, units as pyunits, check_optimal_termination, @@ -54,12 +56,86 @@ _log = idaeslog.getLogger(__name__) +class ASM2dADM1Scaler(CustomScalerBase): + """ + Default modular scaler for ASM2d-ADM1 translator block. + 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. + """ + + 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.properties_in, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + self.call_submodel_scaler_method( + submodel=model.properties_out, + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + 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.properties_in, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.properties_out, + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale unit level constraints + 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("Translator_ASM2d_ADM1") class TranslatorDataASM2dADM1(TranslatorData): """ Translator block representing the ASM2d/ADM1 interface """ + default_scaler = ASM2dADM1Scaler + CONFIG = TranslatorData.CONFIG() # TODO: Change the default to False