diff --git a/idaes/config.py b/idaes/config.py index b698966954..9ce1ebb2d3 100644 --- a/idaes/config.py +++ b/idaes/config.py @@ -385,7 +385,7 @@ def _new_idaes_config_block(): "scale_model", pyomo.common.config.ConfigValue( domain=Bool, - default=False, # TODO: Change to true once transition complete + default=True, description="Whether to apply model scaling in writer", ), ) diff --git a/idaes/core/scaling/tests/test_scaling_profiler.py b/idaes/core/scaling/tests/test_scaling_profiler.py index bb6f770d11..7f1ceb9be5 100644 --- a/idaes/core/scaling/tests/test_scaling_profiler.py +++ b/idaes/core/scaling/tests/test_scaling_profiler.py @@ -18,6 +18,7 @@ from io import StringIO import os +import re import pytest import numpy as np @@ -30,6 +31,7 @@ MethaneParameterBlock as MethaneCombustionParameterBlock, ) from idaes.core.util import from_json, StoreSpec +from idaes.core.util.scaling import constraint_scaling_transform from idaes.core.scaling import set_scaling_factor from idaes.core.scaling.scaler_profiling import ScalingProfiler @@ -207,6 +209,62 @@ def test_solve_perturbed_state(self): assert res["iterations"] == 0 +def demo_model_cst(): + m = ConcreteModel() + + m.v1 = Var(initialize=2) + m.v1.fix() + m.v2 = Var(initialize=4) + m.v3 = Var(initialize=-6) + + m.c1 = Constraint(expr=m.v2 == m.v1**2) + m.c2 = Constraint(expr=0 == m.v1 + m.v2 + m.v3) + + constraint_scaling_transform(m.c2, 0.5) + + return m + + +def demo_scaling_cst(model): + set_scaling_factor(model.v2, 0.5) + constraint_scaling_transform(model.c1, 0.5) + + +@pytest.mark.unit +def test_constraint_scaling_transform(): + # In this case the user has put old-style scaling in the build_model method + sp = ScalingProfiler( + build_model=demo_model_cst, + user_scaling=demo_scaling, + perturb_state=demo_pertubration, + ) + with pytest.raises( + RuntimeError, + match=re.escape( + "Attempted to set constraint scaling factor for transformed constraint c2. " + "Please use only one of set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ), + ): + _ = sp.profile_scaling_methods() + + # Here the user includes old-style scaling in the scaling method + sp = ScalingProfiler( + build_model=demo_model, + user_scaling=demo_scaling_cst, + perturb_state=demo_pertubration, + ) + with pytest.raises( + RuntimeError, + match=re.escape( + "Attempted to set constraint scaling factor for transformed constraint c1. " + "Please use only one of set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ), + ): + _ = sp.profile_scaling_methods() + + # Case study using Gibbs reactor model # Get solution json from scaling tests FILENAME = "gibbs_solution.json" diff --git a/idaes/core/scaling/tests/test_util.py b/idaes/core/scaling/tests/test_util.py index f11cb1d804..4dff80cb7c 100644 --- a/idaes/core/scaling/tests/test_util.py +++ b/idaes/core/scaling/tests/test_util.py @@ -45,6 +45,7 @@ unscaled_variables_generator, unscaled_constraints_generator, ) +from idaes.core.util.scaling import constraint_scaling_transform import idaes.logger as idaeslog currdir = this_file_dir() @@ -1742,6 +1743,52 @@ def test_set_scaling_factor_overwrite_false(self, caplog): ) assert m.scaling_factor[m.v] == 10 + @pytest.mark.unit + def test_set_scaling_factor_transformed_constraint(self): + m = ConcreteModel() + m.x = Var(initialize=500) + m.c1 = Constraint(expr=m.x <= 1e3) + m.c2 = Constraint(expr=m.x == 1e3) + m.c3 = Constraint(expr=m.x >= 1e3) + + def indexed_constraint_rule(b, idx): + return b.x == 1e3 + + m.c4 = Constraint([1, 2, 3], rule=indexed_constraint_rule) + + def match(cname): + return re.escape( + f"Attempted to set constraint scaling factor for transformed constraint {cname}. " + "Please use only one of set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ) + + constraint_scaling_transform(m.c1, 1e-3) + constraint_scaling_transform(m.c2, 1e-3) + constraint_scaling_transform(m.c3, 1e-3) + for idx in m.c4: + constraint_scaling_transform(m.c4[idx], 1e-3) + + with pytest.raises(RuntimeError, match=match("c1")): + set_scaling_factor(m.c1, 1e-3) + with pytest.raises(RuntimeError, match=match("c2")): + set_scaling_factor(m.c2, 1e-3) + with pytest.raises(RuntimeError, match=match("c3")): + set_scaling_factor(m.c3, 1e-3) + with pytest.raises( + RuntimeError, + match=re.escape( + "Attempted to set constraint scaling factor for indexed constraint c4 " + "with transformed ConstraintData children. Please use only one of " + "set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ), + ): + set_scaling_factor(m.c4, 1e-3) + for idx in m.c4: + with pytest.raises(RuntimeError, match=match(f"c4[{idx}]")): + set_scaling_factor(m.c4[idx], 1e-3) + @pytest.fixture def model_expr(self): m = ConcreteModel() diff --git a/idaes/core/scaling/util.py b/idaes/core/scaling/util.py index 8c84d31c8e..bfb9d46476 100644 --- a/idaes/core/scaling/util.py +++ b/idaes/core/scaling/util.py @@ -45,11 +45,13 @@ from pyomo.core.base.constraint import ConstraintData from pyomo.core.base.expression import ExpressionData from pyomo.core.base.param import ParamData +from pyomo.core.base.constraint import ConstraintData from pyomo.core import expr as EXPR from pyomo.common.numeric_types import native_types from pyomo.core.base.units_container import _PyomoUnit from idaes.core.util.exceptions import BurntToast +from idaes.core.util.scaling import get_constraint_transform_applied_scaling_factor import idaes.logger as idaeslog @@ -509,14 +511,31 @@ def set_scaling_factor(component, scaling_factor: float, overwrite: bool = False elif math.isnan(scaling_factor): raise ValueError(f"scaling factor for {component.name} is NaN.") + if isinstance(component, ConstraintData): + xform_factor = get_constraint_transform_applied_scaling_factor(component) + if xform_factor is not None: + raise RuntimeError( + f"Attempted to set constraint scaling factor for transformed constraint {component.name}. " + "Please use only one of set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ) + elif isinstance(component, Constraint): + for condata in component.values(): + xform_factor = get_constraint_transform_applied_scaling_factor(condata) + if xform_factor is not None: + raise RuntimeError( + "Attempted to set constraint scaling factor for indexed constraint " + f"{component.name} with transformed ConstraintData children. Please " + "use only one of set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ) + if component.is_indexed(): - # What if a scaling factor already exists for the indexed component? - # for idx in component: - # set_scaling_factor(component[idx], scaling_factor=scaling_factor, overwrite=overwrite) raise TypeError( f"Component {component.name} is indexed. Set scaling factors for individual indices instead." ) + # Get suffix and assign scaling factor sfx = get_component_scaling_suffix(component) if not overwrite and component in sfx: diff --git a/idaes/core/solvers/tests/test_solvers.py b/idaes/core/solvers/tests/test_solvers.py index 9dc21c00a7..4f135c1bfa 100644 --- a/idaes/core/solvers/tests/test_solvers.py +++ b/idaes/core/solvers/tests/test_solvers.py @@ -299,7 +299,7 @@ def test_get_solver_ipopt_v2(): assert solver.options.max_iter == 200 assert solver.config.writer_config.linear_presolve - assert not solver.config.writer_config.scale_model + assert solver.config.writer_config.scale_model @pytest.mark.skipif(not pyo.SolverFactory("ipopt").available(False), reason="no Ipopt") @@ -308,7 +308,7 @@ def test_get_solver_ipopt_v2_w_options(): solver = get_solver( "ipopt_v2", options={"tol": 1e-5, "foo": "bar"}, - writer_config={"linear_presolve": False}, + writer_config={"linear_presolve": False, "scale_model": False}, ) assert isinstance(solver, LegacySolverWrapper) diff --git a/idaes/core/util/scaling.py b/idaes/core/util/scaling.py index 8d0aa3aef8..e0cdbd31da 100644 --- a/idaes/core/util/scaling.py +++ b/idaes/core/util/scaling.py @@ -41,7 +41,7 @@ from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP from pyomo.contrib.pynumero.asl import AmplInterface from pyomo.common.modeling import unique_component_name -from pyomo.core.base.constraint import ConstraintData +from pyomo.core.base.constraint import ConstraintData, Constraint from pyomo.common.collections import ComponentMap, ComponentSet from pyomo.dae import DerivativeVar from pyomo.dae.flatten import slice_component_along_sets @@ -49,11 +49,14 @@ from pyomo.core import expr as EXPR from pyomo.common.numeric_types import native_types from pyomo.core.base.units_container import _PyomoUnit +from pyomo.core.base.suffix import SuffixFinder import idaes.logger as idaeslog _log = idaeslog.getLogger(__name__) +sfFinder = SuffixFinder("scaling_factor") + def __none_left_mult(x, y): """PRIVATE FUNCTION, If x is None return None, else return x * y""" @@ -230,6 +233,26 @@ def set_scaling_factor(c, v, data_objects=True, overwrite=True): # doesn't exist. This handles the case where you get a constant 0 and # need its scale factor to scale the mass balance. return 1 + + if isinstance(c, ConstraintData): + xform_factor = get_constraint_transform_applied_scaling_factor(c) + if xform_factor is not None: + raise RuntimeError( + f"Attempted to set constraint scaling factor for transformed constraint {c.name}. " + "Please use only one of set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ) + elif isinstance(c, Constraint): + for condata in c.values(): + xform_factor = get_constraint_transform_applied_scaling_factor(condata) + if xform_factor is not None: + raise RuntimeError( + f"Attempted to set constraint scaling factor for indexed constraint {c.name}" + "with transformed ConstraintData children. Please use only one of " + "set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ) + try: suf = c.parent_block().scaling_factor @@ -460,6 +483,14 @@ def constraint_scaling_transform(c, s, overwrite=True): s = pyo.value(s) if not isinstance(c, ConstraintData): raise TypeError(f"{c} is not a constraint or is an indexed constraint") + scaling_factor = sfFinder.find(component_data=c) + if scaling_factor is not None: + raise RuntimeError( + f"Attempted to transform constraint {c.name} with existing scaling factor. " + "Please use only one of set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ) + st = get_constraint_transform_applied_scaling_factor(c, default=None) if not overwrite and st is not None: diff --git a/idaes/core/util/tests/test_scaling.py b/idaes/core/util/tests/test_scaling.py index f07ad129cd..8799d1a40c 100644 --- a/idaes/core/util/tests/test_scaling.py +++ b/idaes/core/util/tests/test_scaling.py @@ -464,6 +464,53 @@ def c2(b, i): assert sc.get_scaling_factor(m.z[i]) is None +@pytest.mark.unit +def test_set_scaling_factor_transformed_constraint(): + m = pyo.ConcreteModel() + m.x = pyo.Var(initialize=500) + m.c1 = pyo.Constraint(expr=m.x <= 1e3) + m.c2 = pyo.Constraint(expr=m.x == 1e3) + m.c3 = pyo.Constraint(expr=m.x >= 1e3) + + def indexed_constraint_rule(b, idx): + return b.x == 1e3 + + m.c4 = pyo.Constraint([1, 2, 3], rule=indexed_constraint_rule) + + def match(cname): + return re.escape( + f"Attempted to set constraint scaling factor for transformed constraint {cname}. " + "Please use only one of set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ) + + sc.constraint_scaling_transform(m.c1, 1e-3) + sc.constraint_scaling_transform(m.c2, 1e-3) + sc.constraint_scaling_transform(m.c3, 1e-3) + for idx in m.c4: + sc.constraint_scaling_transform(m.c4[idx], 1e-3) + + with pytest.raises(RuntimeError, match=match("c1")): + sc.set_scaling_factor(m.c1, 1e-3) + with pytest.raises(RuntimeError, match=match("c2")): + sc.set_scaling_factor(m.c2, 1e-3) + with pytest.raises(RuntimeError, match=match("c3")): + sc.set_scaling_factor(m.c3, 1e-3) + with pytest.raises( + RuntimeError, + match=re.escape( + "Attempted to set constraint scaling factor for indexed constraint c4" + "with transformed ConstraintData children. Please use only one of " + "set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ), + ): + sc.set_scaling_factor(m.c4, 1e-3) + for idx in m.c4: + with pytest.raises(RuntimeError, match=match(f"c4[{idx}]")): + sc.set_scaling_factor(m.c4[idx], 1e-3) + + @pytest.mark.unit def test_set_and_get_scaling_factor(): m = pyo.ConcreteModel() @@ -698,6 +745,42 @@ def test_remove_units(self, model): assert model.constraint_transformed_scaling_factor[model.c2] == 1e-3 +@pytest.mark.unit +def test_constraint_scaling_transform_existing_scaling_factor(): + m = pyo.ConcreteModel() + m.x = pyo.Var(initialize=500) + m.c1 = pyo.Constraint(expr=m.x <= 1e3) + m.c2 = pyo.Constraint(expr=m.x == 1e3) + m.c3 = pyo.Constraint(expr=m.x >= 1e3) + + def indexed_constraint_rule(b, idx): + return b.x == 1e3 + + m.c4 = pyo.Constraint([1, 2, 3], rule=indexed_constraint_rule) + + def match(cname): + return re.escape( + f"Attempted to transform constraint {cname} with existing scaling factor. " + "Please use only one of set_scaling_factor and constraint_scaling_transform " + "per constraint to avoid double scaling." + ) + + sc.set_scaling_factor(m.c1, 1e-3) + sc.set_scaling_factor(m.c2, 1e-3) + sc.set_scaling_factor(m.c3, 1e-3) + sc.set_scaling_factor(m.c4, 1e-3) + + with pytest.raises(RuntimeError, match=match("c1")): + sc.constraint_scaling_transform(m.c1, 1e-3) + with pytest.raises(RuntimeError, match=match("c2")): + sc.constraint_scaling_transform(m.c2, 1e-3) + with pytest.raises(RuntimeError, match=match("c3")): + sc.constraint_scaling_transform(m.c3, 1e-3) + for idx in range(1, 4): + with pytest.raises(RuntimeError, match=match(f"c4[{idx}]")): + sc.constraint_scaling_transform(m.c4[idx], 1e-3) + + class TestScaleSingleConstraint: @pytest.fixture(scope="class") def model(self): diff --git a/idaes/models/flowsheets/demo_flowsheet.py b/idaes/models/flowsheets/demo_flowsheet.py index 7e461bcbe9..debfa758d2 100644 --- a/idaes/models/flowsheets/demo_flowsheet.py +++ b/idaes/models/flowsheets/demo_flowsheet.py @@ -141,19 +141,6 @@ def set_scaling(m): ) iscale.set_scaling_factor(m.fs.M01.mixed_state[0].enthalpy_flow_terms["Liq"], 1) iscale.set_scaling_factor(m.fs.M01.mixed_state[0].enthalpy_flow_terms["Vap"], 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_total, 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_sum_mol_frac, 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_mol_frac_out, 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_comp["benzene"], 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_comp["toluene"], 1) - iscale.set_scaling_factor( - m.fs.M01.mixed_state[0].eq_phase_equilibrium["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.M01.mixed_state[0].eq_phase_equilibrium["toluene"], 1 - ) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_P_vap["benzene"], 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_P_vap["toluene"], 1) iscale.set_scaling_factor(m.fs.M01.mixed_state[0].pressure, 1e-5) iscale.set_scaling_factor(m.fs.M01.mixed_state[0].pressure_sat_comp, 1e-5) @@ -214,42 +201,6 @@ def set_scaling(m): m.fs.H02.control_volume.properties_in[0].enth_mol_phase_comp["Vap", "toluene"], 1e-4, ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].eq_comp["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].eq_comp["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].eq_comp["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].eq_comp["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].eq_phase_equilibrium["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].eq_phase_equilibrium["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].eq_phase_equilibrium["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].eq_phase_equilibrium["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].eq_P_vap["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].eq_P_vap["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0.0].eq_P_vap["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0.0].eq_P_vap["toluene"], 1 - ) iscale.set_scaling_factor( m.fs.H02.control_volume.properties_out[0].enth_mol_phase_comp["Liq", "benzene"], 1e-4, @@ -266,17 +217,6 @@ def set_scaling(m): m.fs.H02.control_volume.properties_out[0].enth_mol_phase_comp["Vap", "toluene"], 1e-4, ) - iscale.set_scaling_factor(m.fs.H02.control_volume.properties_in[0.0].eq_total, 1) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].eq_sum_mol_frac, 1 - ) - iscale.set_scaling_factor(m.fs.H02.control_volume.properties_out[0.0].eq_total, 1) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0.0].eq_sum_mol_frac, 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0.0].eq_mol_frac_out, 1 - ) iscale.set_scaling_factor( m.fs.H02.control_volume.properties_in[0.0].material_flow_terms[ @@ -336,18 +276,6 @@ def set_scaling(m): m.fs.F03.control_volume.properties_in[0].mole_frac_comp["toluene"], 1e1 ) - iscale.set_scaling_factor(m.fs.F03.control_volume.properties_in[0.0].eq_total, 1) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].eq_sum_mol_frac, 1 - ) - iscale.set_scaling_factor(m.fs.F03.control_volume.properties_out[0.0].eq_total, 1) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0.0].eq_sum_mol_frac, 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0.0].eq_mol_frac_out, 1 - ) - iscale.set_scaling_factor( m.fs.F03.control_volume.properties_in[0].enth_mol_phase_comp["Liq", "benzene"], 1e-4, @@ -380,43 +308,6 @@ def set_scaling(m): m.fs.F03.control_volume.properties_out[0].enth_mol_phase_comp["Vap", "toluene"], 1e-4, ) - - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].eq_comp["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].eq_comp["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].eq_comp["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].eq_comp["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].eq_phase_equilibrium["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].eq_phase_equilibrium["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].eq_phase_equilibrium["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].eq_phase_equilibrium["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].eq_P_vap["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].eq_P_vap["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0.0].eq_P_vap["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0.0].eq_P_vap["toluene"], 1 - ) iscale.set_scaling_factor( m.fs.F03.control_volume.properties_in[0.0].material_flow_terms[ "Liq", "benzene"