diff --git a/docs/explanations/scaling_toolbox/scaling_theory.rst b/docs/explanations/scaling_toolbox/scaling_theory.rst index 6c8a9a4300..52b2a00e9d 100644 --- a/docs/explanations/scaling_toolbox/scaling_theory.rst +++ b/docs/explanations/scaling_toolbox/scaling_theory.rst @@ -14,7 +14,7 @@ Types of Scaling * Constraint residual scaling refers to the magnitude of the residual of all constraints in the model. This is important as this is required to determine whether or not a model has converged, thus it is important that a residual equal to the solver tolerance does not significantly alter the solution to the model. * E.g. consider a constraint A=B. If the magnitude of A and B are 1e6 and the solver tolerance is 1e-6, this means that A and B need to be solved to 12 significant figures of precision to converge the constraint (which may be unnecessarily onerous). Similarly, if A and B were of order 1e-6, then they will be considered equal if A=1e-6 and B=0; i.e. ±100% error. - + * Jacobian scaling refers to the overall scaling and conditioning of the problem Jacobian. This is important as this determines the ability of the solver to find a path from the initial state to the final solution. It is important to ensure that the Jacobian is well-conditioned both in order to get good numerical behavior from the solver and also to ensure that floating point round-off error does not result in large changes to the variables at the solution. These aspects are not always complimentary, and there are often cases where improving one aspect of the model scaling can negatively affect another aspect (e.g., focusing too much on Jacobian condition number can often result in poor constraint residual tolerances). Additionally, each aspect affects different parts of the solver routine; Jacobian scaling is important for determining the step direction and size at each iteration, whilst constraint residual scaling is important for determining if and when a solution is found. The relative importance of each of these depends on the solver being used (e.g., for IPOPT constraint residual scaling is more important than Jacobian scaling as IPOPT does its own internal scaling of the Jacobian), however in general all of these will have an impact upon the solver behavior and users should endeavor to have good scaling for all aspects of the model. diff --git a/docs/reference_guides/model_libraries/generic/unit_models/equilibrium.rst b/docs/reference_guides/model_libraries/generic/unit_models/equilibrium.rst index d34c60f844..9f7536639c 100644 --- a/docs/reference_guides/model_libraries/generic/unit_models/equilibrium.rst +++ b/docs/reference_guides/model_libraries/generic/unit_models/equilibrium.rst @@ -51,3 +51,9 @@ EquilibriumReactorData Class .. autoclass:: EquilibriumReactorData :members: + +EquilibriumReactorScaler Class +------------------------------ + +.. autoclass:: EquilibriumReactorScaler + :members: diff --git a/idaes/core/base/property_base.py b/idaes/core/base/property_base.py index 2e0743f994..b7afb011bb 100644 --- a/idaes/core/base/property_base.py +++ b/idaes/core/base/property_base.py @@ -593,6 +593,14 @@ def include_inherent_reactions(self): # pylint: disable-next=protected-access return self.parent_component()._include_inherent_reactions() + @property + def default_initializer(self): + return self.parent_component().default_initializer + + @property + def default_scaler(self): + return self.parent_component().default_scaler + def build(self): """ General build method for StateBlockDatas. diff --git a/idaes/core/base/reaction_base.py b/idaes/core/base/reaction_base.py index 7a42c36f0f..5d43c5485c 100644 --- a/idaes/core/base/reaction_base.py +++ b/idaes/core/base/reaction_base.py @@ -267,6 +267,14 @@ def phase_list(self): def phase_component_set(self): return self.state_ref.phase_component_set + @property + def default_initializer(self): + return self.parent_component().default_initializer + + @property + def default_scaler(self): + return self.parent_component().default_scaler + def lock_attribute_creation_context(self): """Returns a context manager that does not allow attributes to be created while in the context and allows attributes to be created normally outside diff --git a/idaes/core/scaling/__init__.py b/idaes/core/scaling/__init__.py index 0d1b13ee0f..e3e3661d2b 100644 --- a/idaes/core/scaling/__init__.py +++ b/idaes/core/scaling/__init__.py @@ -11,7 +11,7 @@ # for full copyright and license information. ################################################################################# from .autoscaling import AutoScaler -from .custom_scaler_base import CustomScalerBase +from .custom_scaler_base import CustomScalerBase, ConstraintScalingScheme from .scaler_profiling import ScalingProfiler from .util import ( scaling_factors_from_json_file, diff --git a/idaes/core/scaling/custom_scaler_base.py b/idaes/core/scaling/custom_scaler_base.py index 4da3cdb20b..0d47673805 100644 --- a/idaes/core/scaling/custom_scaler_base.py +++ b/idaes/core/scaling/custom_scaler_base.py @@ -571,6 +571,7 @@ def call_submodel_scaler_method( submodel: local name of submodel to be scaled as str submodel_scalers: user provided ComponentMap of Scalers to use for submodels method: name of method to call from submodel (as string) + overwrite: whether to overwrite existing scaling factors Returns: None diff --git a/idaes/core/scaling/scaler_profiling.py b/idaes/core/scaling/scaler_profiling.py index fe8b72e187..9bef06ead3 100644 --- a/idaes/core/scaling/scaler_profiling.py +++ b/idaes/core/scaling/scaler_profiling.py @@ -19,7 +19,11 @@ from pyomo.common.tempfiles import TempfileManager from idaes.core.util.scaling import jacobian_cond -from idaes.core.scaling import AutoScaler, CustomScalerBase +from idaes.core.scaling.autoscaling import AutoScaler +from idaes.core.scaling.custom_scaler_base import ( + CustomScalerBase, + ConstraintScalingScheme, +) from idaes.core.solvers import get_solver @@ -83,23 +87,23 @@ def __init__( "Vars Only": (None, {}), "Harmonic": ( cscaler.scale_constraint_by_nominal_value, - {"scheme": "harmonic_mean"}, + {"scheme": ConstraintScalingScheme.harmonicMean}, ), "Inverse Sum": ( cscaler.scale_constraint_by_nominal_value, - {"scheme": "inverse_sum"}, + {"scheme": ConstraintScalingScheme.inverseSum}, ), "Inverse Root Sum Squares": ( cscaler.scale_constraint_by_nominal_value, - {"scheme": "inverse_root_sum_squared"}, + {"scheme": ConstraintScalingScheme.inverseRSS}, ), "Inverse Maximum": ( cscaler.scale_constraint_by_nominal_value, - {"scheme": "inverse_maximum"}, + {"scheme": ConstraintScalingScheme.inverseMaximum}, ), "Inverse Minimum": ( cscaler.scale_constraint_by_nominal_value, - {"scheme": "inverse_minimum"}, + {"scheme": ConstraintScalingScheme.inverseMinimum}, ), "Nominal L1 Norm": ( cscaler.scale_constraint_by_nominal_derivative_norm, diff --git a/idaes/core/scaling/tests/test_custom_scaling_integration.py b/idaes/core/scaling/tests/test_custom_scaling_integration.py index 548846f786..68c5e3af5c 100644 --- a/idaes/core/scaling/tests/test_custom_scaling_integration.py +++ b/idaes/core/scaling/tests/test_custom_scaling_integration.py @@ -38,8 +38,6 @@ from idaes.models.properties.activity_coeff_models.methane_combustion_ideal import ( MethaneParameterBlock as MethaneCombustionParameterBlock, ) -from idaes.core.util.testing import PhysicalParameterTestBlock, initialization_tester -from idaes.core.solvers import get_solver from idaes.core.util import to_json, from_json, StoreSpec from idaes.core.util.scaling import jacobian_cond from idaes.core.scaling import AutoScaler, CustomScalerBase, set_scaling_factor diff --git a/idaes/models/properties/examples/saponification_reactions.py b/idaes/models/properties/examples/saponification_reactions.py index a855868104..ffd2334685 100644 --- a/idaes/models/properties/examples/saponification_reactions.py +++ b/idaes/models/properties/examples/saponification_reactions.py @@ -31,11 +31,13 @@ from idaes.core.util.misc import add_object_reference from idaes.core.util.constants import Constants as const import idaes.logger as idaeslog +from idaes.core.scaling import CustomScalerBase # Some more information about this module __author__ = "Andrew Lee" +from idaes.core.util.scaling import get_scaling_factor # Set up logger _log = idaeslog.getLogger(__name__) @@ -111,12 +113,84 @@ def define_metadata(cls, obj): ) +class SaponificationReactionScaler(CustomScalerBase): + """ + Scaler for saponification reaction package. + + Variables are scaled by nominal order of magnitude, and constraints + using the inverse maximum scheme. + """ + + DEFAULT_SCALING_FACTORS = {"reaction_rate": 1e2} + + def variable_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: dict = None + ): + if model.is_property_constructed("k_rxn"): + # First check to see if k_rxn is already scaled + sf = get_scaling_factor(model.k_rxn) + + if sf is not None and not overwrite: + # k_rxn already has a scaling factor and we are not set to overwrite - move on + pass + else: + # Hopefully temperature has been scaled, so we can get the nominal value of k_rxn + # by walking the expression in the constraint. + nominals = self.get_expression_nominal_values(model.arrhenius_eqn) + + # We should get two values, k_rxn (LHS) and the Arrhenius equation (RHS) + # As of 10/3/2024, the LHS will be the 0-th element of the list, and the RHS the 1st + # However, we cannot assume this will always be the case + + # If LHS has been scaled, nominal will be 1/sf, otherwise it will be k_rxn.value + # Find the value which does NOT match this - guess that this is the 1st element + if nominals[1] != model.k_rxn.value and sf is None: + # This is the most likely case, so check it first + nominal = nominals[1] + elif sf is not None and nominals[1] != 1 / sf: + # Next, check for case where k_rxn was already scaled + nominal = nominals[1] + else: + # Otherwise we have the case where something changed in Pyomo since 10/3/2024 + nominal = nominals[0] + + self.set_variable_scaling_factor( + model.k_rxn, + 1 / nominal, + overwrite=overwrite, + ) + + if model.is_property_constructed("reaction_rate"): + for j in model.reaction_rate.values(): + self.scale_variable_by_default(j, overwrite=overwrite) + + def constraint_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: dict = None + ): + if model.is_property_constructed("arrhenius_eqn"): + self.scale_constraint_by_nominal_value( + model.arrhenius_eqn, + scheme="inverse_maximum", + overwrite=overwrite, + ) + + if model.is_property_constructed("rate_expression"): + for j in model.rate_expression.values(): + self.scale_constraint_by_nominal_value( + j, + scheme="inverse_maximum", + overwrite=overwrite, + ) + + class _ReactionBlock(ReactionBlockBase): """ This Class contains methods which should be applied to Reaction Blocks as a whole, rather than individual elements of indexed Reaction Blocks. """ + default_scaler = SaponificationReactionScaler + def initialize(blk, outlvl=idaeslog.NOTSET, **kwargs): """ Initialization routine for reaction package. diff --git a/idaes/models/properties/examples/saponification_thermo.py b/idaes/models/properties/examples/saponification_thermo.py index e21ed6094d..7086858fe7 100644 --- a/idaes/models/properties/examples/saponification_thermo.py +++ b/idaes/models/properties/examples/saponification_thermo.py @@ -44,6 +44,7 @@ from idaes.core.util.model_statistics import degrees_of_freedom from idaes.core.util.initialization import fix_state_vars, revert_state_vars import idaes.logger as idaeslog +from idaes.core.scaling import CustomScalerBase # Some more information about this module __author__ = "Andrew Lee" @@ -135,12 +136,57 @@ def define_metadata(cls, obj): ) +class SaponificationPropertiesScaler(CustomScalerBase): + """ + Scaler for saponification properties package. + + Flow and concentration are scaled by default value (if no user input provided), + pressure is scaled assuming order of magnitude of 1e5 Pa, and temperature is + scaled using the average of the bounds. Constraints using the inverse maximum + scheme. + """ + + UNIT_SCALING_FACTORS = { + # "QuantityName: (reference units, scaling factor) + "Pressure": (units.Pa, 1e-5), + } + + DEFAULT_SCALING_FACTORS = { + "flow_vol": 1e2, + "conc_mol_comp": 1e-2, + } + + def variable_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: dict = None + ): + self.scale_variable_by_default(model.flow_vol, overwrite=overwrite) + self.scale_variable_by_units(model.pressure, overwrite=overwrite) + self.scale_variable_by_bounds(model.temperature, overwrite=overwrite) + for k, v in model.conc_mol_comp.items(): + if k == "H2O": + self.set_variable_scaling_factor(v, 1e-4, overwrite=overwrite) + else: + self.scale_variable_by_default(v, overwrite=overwrite) + + def constraint_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: dict = None + ): + if model.is_property_constructed("conc_water_eqn"): + self.set_constraint_scaling_factor( + model.conc_water_eqn, + 1e-4, + overwrite=overwrite, + ) + + class _StateBlock(StateBlock): """ This Class contains methods which should be applied to Property Blocks as a whole, rather than individual elements of indexed Property Blocks. """ + default_scaler = SaponificationPropertiesScaler + def fix_initialization_states(self): """ Fixes state variables for state blocks. diff --git a/idaes/models/properties/examples/tests/test_saponification_reaction.py b/idaes/models/properties/examples/tests/test_saponification_reaction.py index 1da9acbd77..fdb879b76e 100644 --- a/idaes/models/properties/examples/tests/test_saponification_reaction.py +++ b/idaes/models/properties/examples/tests/test_saponification_reaction.py @@ -14,18 +14,20 @@ Tests for saponification property package example. Authors: Andrew Lee """ +from math import exp import pytest -from pyomo.environ import ConcreteModel, Constraint, Param, units, value, Var -from idaes.core import MaterialFlowBasis +from pyomo.environ import ConcreteModel, Constraint, Param, Suffix, units, value, Var + +from idaes.core import MaterialFlowBasis from idaes.models.properties.examples.saponification_reactions import ( SaponificationReactionParameterBlock, ReactionBlock, + SaponificationReactionScaler, ) from idaes.models.properties.examples.saponification_thermo import ( SaponificationParameterBlock, ) - from idaes.core.solvers import get_solver @@ -100,6 +102,8 @@ def test_build(self, model): assert model.rxns[1].temperature_ref is model.props[1].temperature assert model.rxns[1].dh_rxn is model.rparams.dh_rxn + assert model.rxns[1].default_scaler is SaponificationReactionScaler + @pytest.mark.unit def test_rate_constant(self, model): assert isinstance(model.rxns[1].k_rxn, Var) @@ -124,3 +128,90 @@ def test_initialize(self, model): def check_units(self, model): units.assert_units_consistent(model) + + +class TestSaponificationReactionScaler(object): + @pytest.mark.unit + def test_variable_scaling_routine(self): + model = ConcreteModel() + model.pparams = SaponificationParameterBlock() + model.rparams = SaponificationReactionParameterBlock( + property_package=model.pparams + ) + + model.props = model.pparams.build_state_block([1]) + model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) + + # Trigger build of reaction properties + model.rxns[1].reaction_rate + + scaler = model.rxns[1].default_scaler() + assert isinstance(scaler, SaponificationReactionScaler) + + scaler.variable_scaling_routine(model.rxns[1]) + + assert isinstance(model.rxns[1].scaling_factor, Suffix) + + sfx = model.rxns[1].scaling_factor + assert len(sfx) == 2 + assert sfx[model.rxns[1].k_rxn] == pytest.approx( + 1 / (3.132e6 * exp(-43000 / (8.31446262 * 298.15))), rel=1e-8 + ) + assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8) + + @pytest.mark.unit + def test_constraint_scaling_routine(self): + model = ConcreteModel() + model.pparams = SaponificationParameterBlock() + model.rparams = SaponificationReactionParameterBlock( + property_package=model.pparams + ) + + model.props = model.pparams.build_state_block([1]) + model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) + + # Trigger build of reaction properties + model.rxns[1].reaction_rate + + scaler = model.rxns[1].default_scaler() + assert isinstance(scaler, SaponificationReactionScaler) + + scaler.constraint_scaling_routine(model.rxns[1]) + + assert isinstance(model.rxns[1].scaling_factor, Suffix) + + sfx = model.rxns[1].scaling_factor + assert len(sfx) == 2 + assert sfx[model.rxns[1].arrhenius_eqn] == pytest.approx(1, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx(1e-4, rel=1e-8) + + @pytest.mark.unit + def test_scale_model(self): + model = ConcreteModel() + model.pparams = SaponificationParameterBlock() + model.rparams = SaponificationReactionParameterBlock( + property_package=model.pparams + ) + + model.props = model.pparams.build_state_block([1]) + model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) + + # Trigger build of reaction properties + model.rxns[1].reaction_rate + + scaler = model.rxns[1].default_scaler() + assert isinstance(scaler, SaponificationReactionScaler) + + scaler.scale_model(model.rxns[1]) + + assert isinstance(model.rxns[1].scaling_factor, Suffix) + + sfx = model.rxns[1].scaling_factor + k_rxn_sf = 1 / (3.132e6 * exp(-43000 / (8.31446262 * 298.15))) + assert len(sfx) == 4 + assert sfx[model.rxns[1].k_rxn] == pytest.approx(k_rxn_sf, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].arrhenius_eqn] == pytest.approx(k_rxn_sf, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx( + 1e-4 * k_rxn_sf, rel=1e-8 + ) diff --git a/idaes/models/properties/examples/tests/test_saponification_thermo.py b/idaes/models/properties/examples/tests/test_saponification_thermo.py index a5d7b8ff00..2d298df3c5 100644 --- a/idaes/models/properties/examples/tests/test_saponification_thermo.py +++ b/idaes/models/properties/examples/tests/test_saponification_thermo.py @@ -16,13 +16,14 @@ """ import pytest -from pyomo.environ import ConcreteModel, Constraint, Param, value, Var +from pyomo.environ import ConcreteModel, Constraint, Param, Suffix, value, Var from pyomo.util.check_units import assert_units_consistent from idaes.core import MaterialBalanceType, EnergyBalanceType, MaterialFlowBasis from idaes.models.properties.examples.saponification_thermo import ( SaponificationParameterBlock, SaponificationStateBlock, + SaponificationPropertiesScaler, ) from idaes.core.solvers import get_solver from idaes.core.initialization import ( @@ -102,6 +103,8 @@ def test_build(self, model): assert isinstance(model.props[1].conc_water_eqn, Constraint) assert len(model.props[1].flow_vol) == 1 + assert model.props[1].default_scaler is SaponificationPropertiesScaler + @pytest.mark.unit def test_build_defined_state(self): model = ConcreteModel() @@ -313,3 +316,89 @@ def test_initializer(): initializer.initialize(model.props) assert initializer.summary[model.props]["status"] == InitializationStatus.Ok + + +class TestSaponificationPropertiesScaler: + @pytest.mark.unit + def test_variable_scaling_routine(self): + model = ConcreteModel() + model.params = SaponificationParameterBlock() + + model.props = model.params.build_state_block([1], defined_state=False) + + scaler = model.props[1].default_scaler() + assert isinstance(scaler, SaponificationPropertiesScaler) + + scaler.variable_scaling_routine(model.props[1]) + + assert isinstance(model.props[1].scaling_factor, Suffix) + + sfx = model.props[1].scaling_factor + assert len(sfx) == 8 + assert sfx[model.props[1].flow_vol] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.props[1].pressure] == pytest.approx(1e-5, rel=1e-8) + assert sfx[model.props[1].temperature] == pytest.approx(1 / 310.65, rel=1e-8) + for k, v in model.props[1].conc_mol_comp.items(): + if k == "H2O": + assert sfx[v] == pytest.approx(1e-4, rel=1e-8) + else: + assert sfx[v] == pytest.approx(1e-2, rel=1e-8) + + @pytest.mark.unit + def test_constraint_scaling_routine(self): + model = ConcreteModel() + model.params = SaponificationParameterBlock() + + model.props = model.params.build_state_block([1], defined_state=False) + + scaler = model.props[1].default_scaler() + assert isinstance(scaler, SaponificationPropertiesScaler) + + scaler.constraint_scaling_routine(model.props[1]) + + assert isinstance(model.props[1].scaling_factor, Suffix) + + sfx = model.props[1].scaling_factor + assert len(sfx) == 1 + assert sfx[model.props[1].conc_water_eqn] == pytest.approx(1e-4, rel=1e-8) + + @pytest.mark.unit + def test_constraint_scaling_routine_defined_state(self): + model = ConcreteModel() + model.params = SaponificationParameterBlock() + + model.props = model.params.build_state_block([1], defined_state=True) + + scaler = model.props[1].default_scaler() + assert isinstance(scaler, SaponificationPropertiesScaler) + + scaler.constraint_scaling_routine(model.props[1]) + + # No constraints, so there should be no Suffix + assert not hasattr(model.props[1], "scaling_factor") + + @pytest.mark.unit + def test_scale_model(self): + model = ConcreteModel() + model.params = SaponificationParameterBlock() + + model.props = model.params.build_state_block([1], defined_state=False) + + scaler = model.props[1].default_scaler() + assert isinstance(scaler, SaponificationPropertiesScaler) + + scaler.scale_model(model.props[1]) + + assert isinstance(model.props[1].scaling_factor, Suffix) + + sfx = model.props[1].scaling_factor + assert len(sfx) == 9 + assert sfx[model.props[1].flow_vol] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.props[1].pressure] == pytest.approx(1e-5, rel=1e-8) + assert sfx[model.props[1].temperature] == pytest.approx(1 / 310.65, rel=1e-8) + for k, v in model.props[1].conc_mol_comp.items(): + if k == "H2O": + assert sfx[v] == pytest.approx(1e-4, rel=1e-8) + else: + assert sfx[v] == pytest.approx(1e-2, rel=1e-8) + assert sfx[model.props[1].conc_water_eqn] == pytest.approx(1e-4, rel=1e-8) diff --git a/idaes/models/unit_models/equilibrium_reactor.py b/idaes/models/unit_models/equilibrium_reactor.py index 4015e4ea76..a8a92c96ae 100644 --- a/idaes/models/unit_models/equilibrium_reactor.py +++ b/idaes/models/unit_models/equilibrium_reactor.py @@ -16,7 +16,7 @@ # Import Pyomo libraries from pyomo.common.config import ConfigBlock, ConfigValue, In, Bool -from pyomo.environ import Reference +from pyomo.environ import Constraint, Reference, units # Import IDAES cores from idaes.core import ( @@ -32,16 +32,175 @@ is_physical_parameter_block, is_reaction_parameter_block, ) +from idaes.core.scaling import CustomScalerBase __author__ = "Andrew Lee" +class EquilibriumReactorScaler(CustomScalerBase): + """ + Default modular scaler for Equilibrium reactors. + + This Scaler relies on modular the associated property and reaction packages, + either through user provided options (submodel_scalers argument) or by default + Scalers assigned to the packages. + + Reaction generation terms are scaled based on component flow rates, whilst + extents of reaction are unscaled. Heat duty is scaled to kW and pressure drop + to 0.1 bar. All constraints are scaled using the inverse maximum scheme. + """ + + UNIT_SCALING_FACTORS = { + # "QuantityName: (reference units, scaling factor) + "Pressure Change": (units.bar, 10), + } + + def variable_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: dict = None + ): + """ + Routine to apply scaling factors to variables in model. + + Submodel Scalers are called for the property and reaction blocks. + Reaction generation terms are scaled based on component flow rates, whilst + extents of reaction are unscaled. Heat duty is scaled to kW and pressure drop + to 0.1 bar. + + 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( + model=model, + submodel="control_volume.properties_in", + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.propagate_state_scaling( + target_state=model.control_volume.properties_out, + source_state=model.control_volume.properties_in, + overwrite=overwrite, + ) + + self.call_submodel_scaler_method( + model=model, + submodel="control_volume.properties_out", + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + model=model, + submodel="control_volume.reactions", + method="variable_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scaling control volume variables + # Reaction generation and extent are hard to know a priori + # A bad guess is worse than no guess, so leave these unscaled + # AL 10/2024: Tried scaling generation by component flow, but that was bad + + # Pressure drop - optional + if hasattr(model.control_volume, "deltaP"): + for t in model.flowsheet().time: + self.scale_variable_by_units( + model.control_volume.deltaP[t], overwrite=overwrite + ) + + # Heat transfer - optional + # Scale heat based on enthalpy flow entering reactor + if hasattr(model.control_volume, "heat"): + for t in model.flowsheet().time: + h_in = 0 + for p in model.control_volume.properties_in.phase_list: + h_in += sum( + self.get_expression_nominal_values( + model.control_volume.properties_in[ + t + ].get_enthalpy_flow_terms(p) + ) + ) + # Scale for heat is general one order of magnitude less than enthalpy flow + self.set_variable_scaling_factor( + model.control_volume.heat[t], 1 / (0.1 * h_in) + ) + + 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 shceme. + + 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( + model=model, + submodel="control_volume.properties_in", + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + model=model, + submodel="control_volume.properties_out", + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + model=model, + submodel="control_volume.reactions", + method="constraint_scaling_routine", + submodel_scalers=submodel_scalers, + overwrite=overwrite, + ) + + # Scale control volume constraints + for c in model.control_volume.component_data_objects( + Constraint, descend_into=False + ): + self.scale_constraint_by_nominal_value( + c, + scheme="inverse_maximum", + overwrite=overwrite, + ) + + # Scale unit level constraints + if hasattr(model, "rate_reaction_constraint"): + for c in model.rate_reaction_constraint.values(): + self.scale_constraint_by_nominal_value( + c, + scheme="inverse_maximum", + overwrite=overwrite, + ) + + @declare_process_block_class("EquilibriumReactor") class EquilibriumReactorData(UnitModelBlockData): """ Standard Equilibrium Reactor Unit Model Class """ + default_scaler = EquilibriumReactorScaler + CONFIG = ConfigBlock() CONFIG.declare( "dynamic", diff --git a/idaes/models/unit_models/tests/test_equilibrium_reactor.py b/idaes/models/unit_models/tests/test_equilibrium_reactor.py index 71fdb23ff8..881bf70df6 100644 --- a/idaes/models/unit_models/tests/test_equilibrium_reactor.py +++ b/idaes/models/unit_models/tests/test_equilibrium_reactor.py @@ -15,9 +15,18 @@ Authors: Andrew Lee """ +from math import exp import pytest -from pyomo.environ import check_optimal_termination, ConcreteModel, value, units +from pyomo.environ import ( + assert_optimal_termination, + ComponentMap, + ConcreteModel, + Suffix, + TransformationFactory, + units, + value, +) from idaes.core import ( FlowsheetBlock, @@ -25,7 +34,10 @@ EnergyBalanceType, MomentumBalanceType, ) -from idaes.models.unit_models.equilibrium_reactor import EquilibriumReactor +from idaes.models.unit_models.equilibrium_reactor import ( + EquilibriumReactor, + EquilibriumReactorScaler, +) from idaes.models.properties.examples.saponification_thermo import ( SaponificationParameterBlock, ) @@ -37,6 +49,10 @@ number_total_constraints, number_unused_variables, ) +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) from idaes.core.util.testing import ( PhysicalParameterTestBlock, ReactionParameterTestBlock, @@ -49,6 +65,7 @@ InitializationStatus, ) from idaes.core.util import DiagnosticsToolbox +from idaes.core.scaling import set_scaling_factor # ----------------------------------------------------------------------------- # Get default solver for testing @@ -86,6 +103,7 @@ def test_config(): assert m.fs.unit.config.reaction_package is m.fs.reactions assert m.fs.unit.default_initializer is SingleControlVolumeUnitInitializer + assert m.fs.unit.default_scaler is EquilibriumReactorScaler # ----------------------------------------------------------------------------- @@ -167,7 +185,7 @@ def test_solve(self, sapon): results = solver.solve(sapon) # Check for optimal solution - assert check_optimal_termination(results) + assert_optimal_termination(results) @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @@ -370,3 +388,383 @@ def test_block_triangularization(self, model): assert not model.fs.unit.inlet.temperature[0].fixed assert not model.fs.unit.inlet.pressure[0].fixed + + +class DummyScaler: + def variable_scaling_routine(self, model, **kwargs): + model._dummy_var_scaler = True + + def constraint_scaling_routine(self, model, **kwargs): + model._dummy_con_scaler = True + + +class TestEquilibriumReactorScaler: + @pytest.fixture + def model(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = SaponificationParameterBlock() + m.fs.reactions = SaponificationReactionParameterBlock( + property_package=m.fs.properties + ) + + m.fs.unit = EquilibriumReactor( + property_package=m.fs.properties, + reaction_package=m.fs.reactions, + has_equilibrium_reactions=False, + has_heat_transfer=True, + has_heat_of_reaction=True, + has_pressure_change=True, + ) + + m.fs.unit.inlet.flow_vol[0].set_value(1.0e-03) + m.fs.unit.inlet.conc_mol_comp[0, "H2O"].set_value(55388.0) + m.fs.unit.inlet.conc_mol_comp[0, "NaOH"].set_value(100.0) + m.fs.unit.inlet.conc_mol_comp[0, "EthylAcetate"].set_value(100.0) + m.fs.unit.inlet.conc_mol_comp[0, "SodiumAcetate"].set_value(0.0) + m.fs.unit.inlet.conc_mol_comp[0, "Ethanol"].set_value(0.0) + + m.fs.unit.inlet.temperature[0].set_value(303.15) + m.fs.unit.inlet.pressure[0].set_value(101325.0) + + m.fs.unit.heat_duty.fix(0) + m.fs.unit.deltaP.fix(0) + + return m + + @pytest.mark.component + def test_variable_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, EquilibriumReactorScaler) + + scaler.variable_scaling_routine(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.control_volume.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 8 + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].flow_vol + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].pressure + ] == pytest.approx(1e-5, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].temperature + ] == pytest.approx(1 / 310.65, rel=1e-8) + for k, v in model.fs.unit.control_volume.properties_in[0].conc_mol_comp.items(): + if k == "H2O": + assert sfx_in[v] == pytest.approx(1e-4, rel=1e-8) + else: + assert sfx_in[v] == pytest.approx(1e-2, rel=1e-8) + + # Outlet state - should be the same as the inlet + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 8 + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].flow_vol + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].pressure + ] == pytest.approx(1e-5, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].temperature + ] == pytest.approx(1 / 310.65, rel=1e-8) + for k, v in model.fs.unit.control_volume.properties_out[ + 0 + ].conc_mol_comp.items(): + if k == "H2O": + assert sfx_out[v] == pytest.approx(1e-4, rel=1e-8) + else: + assert sfx_out[v] == pytest.approx(1e-2, rel=1e-8) + + # Reaction block + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 2 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].k_rxn + ] == pytest.approx( + 1 / (3.132e6 * exp(-43000 / (8.31446262 * 310.65))), rel=1e-8 + ) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R1"] + ] == pytest.approx(1e2, rel=1e-8) + + # Check that unit model has scaling factors + sfx_cv = model.fs.unit.control_volume.scaling_factor + assert isinstance(sfx_cv, Suffix) + assert len(sfx_cv) == 2 + assert sfx_cv[model.fs.unit.control_volume.heat[0]] == pytest.approx( + 1.917448e-05, rel=1e-3 + ) + assert sfx_cv[model.fs.unit.control_volume.deltaP[0]] == pytest.approx( + 1e-4, rel=1e-3 + ) + + # No unit level variables to scale, so no suffix + assert not hasattr(model.fs.unit, "scaling_factor") + + @pytest.mark.component + def test_variable_scaling_routine_submodel_scaler(self, model): + scaler = model.fs.unit.default_scaler() + + scaler_map = ComponentMap() + scaler_map[model.fs.unit.control_volume.properties_in] = DummyScaler + scaler_map[model.fs.unit.control_volume.properties_out] = DummyScaler + scaler_map[model.fs.unit.control_volume.reactions] = DummyScaler + + scaler.variable_scaling_routine( + model.fs.unit, + submodel_scalers=scaler_map, + ) + + # Should call DummyScaler submethod for each submodel + # Should add _dummy_var_scaler = True to all submodels + assert model.fs.unit.control_volume.properties_in[0]._dummy_var_scaler + assert model.fs.unit.control_volume.properties_out[0]._dummy_var_scaler + assert model.fs.unit.control_volume.reactions[0]._dummy_var_scaler + + @pytest.mark.component + def test_constraint_scaling_routine(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, EquilibriumReactorScaler) + + scaler.constraint_scaling_routine(model.fs.unit) + + # Check that sub-models have suffixes - we will assume they are right at this point + sfx_in = model.fs.unit.control_volume.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert ( + len(sfx_in) == 0 + ) # inlet has no constraints. Not quite sure why the Suffix exists + + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 1 + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0.0].conc_water_eqn + ] == pytest.approx(1e-4, rel=1e-8) + + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 2 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].arrhenius_eqn + ] == pytest.approx(1, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R1"] + ] == pytest.approx(1e-4, rel=1e-8) + + # Check that unit model has scaling factors + sfx_cv = model.fs.unit.control_volume.scaling_factor + assert isinstance(model.fs.unit.control_volume.scaling_factor, Suffix) + assert len(sfx_cv) == 12 + assert sfx_cv[ + model.fs.unit.control_volume.enthalpy_balances[0.0] + ] == pytest.approx(8.03894082e-10, rel=1e-8) + assert sfx_cv[ + model.fs.unit.control_volume.pressure_balance[0.0] + ] == pytest.approx(9.86923267e-6, rel=1e-8) + for c in model.fs.unit.control_volume.material_balances.values(): + assert sfx_cv[c] == pytest.approx(1e-2, rel=1e-8) + for ( + c + ) in ( + model.fs.unit.control_volume.rate_reaction_stoichiometry_constraint.values() + ): + assert sfx_cv[c] == pytest.approx(1, rel=1e-8) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 1 + assert sfx_unit[ + model.fs.unit.rate_reaction_constraint[0.0, "R1"] + ] == pytest.approx(1, rel=1e-8) + + @pytest.mark.component + def test_constraint_scaling_routine_submodel_scaler(self, model): + scaler = model.fs.unit.default_scaler() + + scaler_map = ComponentMap() + scaler_map[model.fs.unit.control_volume.properties_in] = DummyScaler + scaler_map[model.fs.unit.control_volume.properties_out] = DummyScaler + scaler_map[model.fs.unit.control_volume.reactions] = DummyScaler + + scaler.constraint_scaling_routine( + model.fs.unit, + submodel_scalers=scaler_map, + ) + + # Should call DummyScaler submethod for each submodel + # Should add _dummy_con_scaler = True to all submodels + assert model.fs.unit.control_volume.properties_in[0]._dummy_con_scaler + assert model.fs.unit.control_volume.properties_out[0]._dummy_con_scaler + assert model.fs.unit.control_volume.reactions[0]._dummy_con_scaler + + @pytest.mark.component + def test_scale_model(self, model): + scaler = model.fs.unit.default_scaler() + + assert isinstance(scaler, EquilibriumReactorScaler) + + scaler.scale_model(model.fs.unit) + + # Inlet state + sfx_in = model.fs.unit.control_volume.properties_in[0].scaling_factor + assert isinstance(sfx_in, Suffix) + assert len(sfx_in) == 8 + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].flow_vol + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].pressure + ] == pytest.approx(1e-5, rel=1e-8) + assert sfx_in[ + model.fs.unit.control_volume.properties_in[0].temperature + ] == pytest.approx(1 / 310.65, rel=1e-8) + for k, v in model.fs.unit.control_volume.properties_in[0].conc_mol_comp.items(): + if k == "H2O": + assert sfx_in[v] == pytest.approx(1e-4, rel=1e-8) + else: + assert sfx_in[v] == pytest.approx(1e-2, rel=1e-8) + + # Outlet state - should be the same as the inlet + sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor + assert isinstance(sfx_out, Suffix) + assert len(sfx_out) == 9 + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].flow_vol + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].pressure + ] == pytest.approx(1e-5, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0].temperature + ] == pytest.approx(1 / 310.65, rel=1e-8) + for k, v in model.fs.unit.control_volume.properties_out[ + 0 + ].conc_mol_comp.items(): + if k == "H2O": + assert sfx_out[v] == pytest.approx(1e-4, rel=1e-8) + else: + assert sfx_out[v] == pytest.approx(1e-2, rel=1e-8) + assert sfx_out[ + model.fs.unit.control_volume.properties_out[0.0].conc_water_eqn + ] == pytest.approx(1e-4, rel=1e-8) + + # Reaction block + sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor + assert isinstance(sfx_rxn, Suffix) + assert len(sfx_rxn) == 4 + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].k_rxn + ] == pytest.approx( + 1 / (3.132e6 * exp(-43000 / (8.31446262 * 310.65))), rel=1e-8 + ) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0].reaction_rate["R1"] + ] == pytest.approx(1e2, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].arrhenius_eqn + ] == pytest.approx(5.4240896, rel=1e-8) + assert sfx_rxn[ + model.fs.unit.control_volume.reactions[0.0].rate_expression["R1"] + ] == pytest.approx(5.4240896e-4, rel=1e-8) + + # Check that unit model has scaling factors + sfx_cv = model.fs.unit.control_volume.scaling_factor + assert isinstance(sfx_cv, Suffix) + assert len(sfx_cv) == 14 + assert sfx_cv[model.fs.unit.control_volume.heat[0]] == pytest.approx( + 1.917448e-05, rel=1e-3 + ) + assert sfx_cv[model.fs.unit.control_volume.deltaP[0]] == pytest.approx( + 1e-4, rel=1e-3 + ) + assert sfx_cv[ + model.fs.unit.control_volume.enthalpy_balances[0.0] + ] == pytest.approx(7.71546823e-08, rel=1e-8) + assert sfx_cv[ + model.fs.unit.control_volume.pressure_balance[0.0] + ] == pytest.approx(1e-5, rel=1e-8) + for (_, _, j), c in model.fs.unit.control_volume.material_balances.items(): + if j == "H2O": + assert sfx_cv[c] == pytest.approx(1e-2, rel=1e-8) + else: + assert sfx_cv[c] == pytest.approx(1, rel=1e-8) + for ( + _, + _, + j, + ), c in ( + model.fs.unit.control_volume.rate_reaction_stoichiometry_constraint.items() + ): + assert sfx_cv[c] == pytest.approx(1, rel=1e-8) + + sfx_unit = model.fs.unit.scaling_factor + assert isinstance(sfx_unit, Suffix) + assert len(sfx_unit) == 1 + assert sfx_unit[ + model.fs.unit.rate_reaction_constraint[0.0, "R1"] + ] == pytest.approx(1e2, rel=1e-8) + + @pytest.mark.integration + def test_example_case(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = SaponificationParameterBlock() + m.fs.reactions = SaponificationReactionParameterBlock( + property_package=m.fs.properties + ) + + m.fs.equil = EquilibriumReactor( + property_package=m.fs.properties, + reaction_package=m.fs.reactions, + has_equilibrium_reactions=False, + has_heat_of_reaction=True, + ) + + m.fs.equil.inlet.flow_vol.fix(1.0e-03) + m.fs.equil.inlet.conc_mol_comp[0, "H2O"].fix(55388.0) + m.fs.equil.inlet.conc_mol_comp[0, "NaOH"].fix(100.0) + m.fs.equil.inlet.conc_mol_comp[0, "EthylAcetate"].fix(100.0) + m.fs.equil.inlet.conc_mol_comp[0, "SodiumAcetate"].fix(1e-8) + m.fs.equil.inlet.conc_mol_comp[0, "Ethanol"].fix(1e-8) + + m.fs.equil.inlet.temperature.fix(303.15) + m.fs.equil.inlet.pressure.fix(101325.0) + + initializer = BlockTriangularizationInitializer() + initializer.initialize(m.fs.equil) + + set_scaling_factor(m.fs.equil.control_volume.properties_in[0].flow_vol, 1e3) + + scaler = EquilibriumReactorScaler() + scaler.scale_model(m.fs.equil) + + m.fs.equil.inlet.flow_vol.fix(1) + m.fs.equil.inlet.conc_mol_comp[0, "NaOH"].fix(200.0) + m.fs.equil.inlet.conc_mol_comp[0, "EthylAcetate"].fix(100.0) + m.fs.equil.inlet.conc_mol_comp[0, "SodiumAcetate"].fix(50) + m.fs.equil.inlet.conc_mol_comp[0, "Ethanol"].fix(1e-8) + + m.fs.equil.inlet.temperature.fix(320) + + solver = get_solver( + "ipopt_v2", writer_config={"linear_presolve": True, "scale_model": True} + ) + results = solver.solve(m, tee=True) + assert_optimal_termination(results) + + # 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.987e05, rel=1e-3 + )