###############################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES).
#
# Copyright (c) 2018-2026 by the software owners: The Regents of the
# University of California, through Lawrence Berkeley National Laboratory,
# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
# University, West Virginia University Research Corporation, et al.
# All rights reserved.  Please see the files COPYRIGHT.md and LICENSE.md
# for full copyright and license information.
###############################################################################

Advanced Scaling Techniques

Author: Doug Allan
Maintainer: Doug Allan
Updated: 2026-04-30

Introduction

This notebook is intended to give the user further insight into how the scaling toolbox can be used to improve robustness. It is a followup to a previous notebook about how to create scaler objects in IDAES, and it assumes some familiarity with the material contained therein. A large portion of the material is an adaptation of Allan, Ostace, and Polly (2025). By the end of this study, the user should understand:

• How a default scaler object can be specified for a property package

• How to set default scaling factors for property package variables

• How to use the DiagnosticsToolbox to troubleshoot saling issues

• The strengths and weaknesses of the NominalValueExpressionWalker used in CustomScalerBase.get_sum_terms_nominal_values and CustomScalerBase.get_expression_nominal_value for use in constraint scaling

• How to use the SVDToolbox to troubleshoot scaling issues

Step 1: Set Up Test Case

We will use the SolventReboiler unit model as a case study. It is a component in a stripping column used for solvent regeneration. It is modeled as a single equilibrium stage with an external heat duty. The inlet stream is entirely liquid, while it has both boilup and bottoms outlet streams. The liquid and vapor phases use different configurations of the modular property framework.

First, we set up and attempt to initialize the unit model.

import logging
import traceback

from pyomo.environ import (
    ComponentMap,
    ConcreteModel,
    Constraint,
    exp,
    Param,
    Reference,
    units as pyunits,
    value,
)

from idaes.core import (
    ControlVolume0DBlock,
    FlowsheetBlock,
    MaterialFlowBasis,
    MomentumBalanceType,
)
import idaes.logger as idaeslog
from idaes.core.scaling.util import get_scaling_factor
from idaes.core.util.exceptions import ConfigurationError, InitializationError

from idaes.models.properties.modular_properties.base.generic_property import (
    GenericParameterBlock,
)

from idaes.models_extra.column_models.solvent_reboiler import SolventReboiler
from idaes.models_extra.column_models.properties.MEA_solvent import (
    configuration as aqueous_mea,
)
from idaes.models_extra.column_models.properties.MEA_vapor import flue_gas

logging.getLogger("pyomo.repn.plugins.nl_writer").setLevel(logging.ERROR)


def create_model():
    m = ConcreteModel()
    m.fs = FlowsheetBlock(dynamic=False)

    m.fs.liquid_properties = GenericParameterBlock(**aqueous_mea)
    m.fs.vapor_properties = GenericParameterBlock(**flue_gas)

    m.fs.unit = SolventReboiler(
        liquid_property_package=m.fs.liquid_properties,
        vapor_property_package=m.fs.vapor_properties,
    )

    m.fs.unit.inlet.flow_mol[0].fix(83.89)
    m.fs.unit.inlet.temperature[0].fix(392.5)
    m.fs.unit.inlet.pressure[0].fix(183700)
    m.fs.unit.inlet.mole_frac_comp[0, "CO2"].fix(0.0326)
    m.fs.unit.inlet.mole_frac_comp[0, "H2O"].fix(0.8589)
    m.fs.unit.inlet.mole_frac_comp[0, "MEA"].fix(0.1085)

    m.fs.unit.heat_duty.fix(430.61e3)
    return m


m = create_model()
reboiler_init = m.fs.unit.default_initializer(always_estimate_states=True)
try:
    reboiler_init.initialize(m.fs.unit, output_level=idaeslog.DEBUG)
except InitializationError:
    traceback.print_exc()

Here we see that initialization failed due to the unit model solve terminating after reaching the maximum number of iterations. However, when solution of a square problem reaches 200-300 iterations without converging, increasing the number of iterations usually does not help. In this case, it terminates with a Restoration Failed status after about 500 iterations. When confronted with a failing solve, we should first head to the DiagnosticsToolbox.

Step 2: Run Model Diagnostics

While initialization fails in every version of Python supported presently (3.10, 3.11, 3.12 3.13), the model termination point(s) that Python 3.10 and 3.11 come to are different than the model termination point that 3.12 and 3.13 come to. In order to ensure that the diagnostic output is the same in every supported version, we load the model state from a json file.

from idaes.core.util import from_json, StoreSpec

from_json(m, fname="advanced_scaling_failed_solution.json", wts=StoreSpec.value())

First, we will use the DiagnosticsToolbox to check to see if there are structural issues with the model. (More information about the DiagnosticsToolbox can be found in its documentation and tutorial notebook

from idaes.core.util import DiagnosticsToolbox

dt = DiagnosticsToolbox(m)
dt.report_structural_issues()

The density relationship is being flagged for potential evaluation errors. The interval arithmetic used to identify potential evaluation errors can produce relatively loose bounds, especially for the sort of large polynomial expressions favored in engineering literature. Independent testing of the density relationship over a wide range of values has demonstrated that no evaluation errors occur, so we can ignore these warnings.

The next step is determining which numerical issues exist.

dt.report_numerical_issues()

In order to help interpret this output, we need to review the solution of systems of nonlinear equations.

When we try to solve a square problem (i.e., a problem with zero degrees of freedom) using IPOPT, we are asking it to solve an equation of the form $\( \mathbf{f}(\mathbf{x}) = \mathbf{0} \)\( in which \)\mathbf{f}(\cdot)\( is a vector-valued function of length \)n\(, and \)\mathbf{x}\( is a vector also of length \)n\(. Due to the finite accuracy of floating point arithmetic, the most we can ask for is \)||\mathbf{f}(\mathbf{x})|| < \varepsilon\(, for some small value of \)\varepsilon\(. The IDAES default solver configuration sets \)\varepsilon = 10^{-6}\(. For some constraint \)j\(, the value of \)f_j(x_k)\( is the called its *residual*, i.e., the remaining error in constraint satisfaction. IPOPT displays the infinity norm of the constraint residual, \)\( \lvert\lvert \mathbf{f}(\mathbf{x}_k)\rvert\rvert_{\infty} := \max_j \lvert f_j(\mathbf{x}_k)\rvert \)\( i.e., the constraint residual with largest magnitude, for each step \)k$ in the primal infeasibility (inf_pr) column.

Problems of this form are typically solved with some variation on Newton’s method. We denote the Jacobian matrix for function \(\mathbf{f}(\mathbf{x})\) as $\( \mathbf{J}(\mathbf{x}) := \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \dots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial x_1} & \dots & \frac{\partial f_n}{\partial x_n} \end{bmatrix} \)\( Note that each row of \)\mathbf{J}(\cdot)$ is associated with a constraint and every column is associated with a variable.

In the classical form of Newton’s method, we start with an initial guess \(\mathbf{x}_0\), then iterate using the relationship

\begin{align} \mathbf{J}(\mathbf{x}_k)\cdot\mathbf{\delta x}_k &= -\mathbf{f}(\mathbf{x}k) \ \mathbf{x}{k+1} &= \mathbf{x}_k + \mathbf{\delta x}_k \end{align}

IPOPT incorporates additional features, such as a line search, in order to enhance its robustness to bad initial guesses, but it is designed to take full Newton steps in the neighborhood of a solution (for the full details, refer to Wächter and Biegler (2006)).

The condition number \(\kappa_2(\cdot)\) (for the Euclidean norm) is a worst-case error bound in the solution of a system of linear equations. In this case, we have $\( \lvert\lvert\mathbf{\delta x}_k\rvert\rvert_2 \leq \kappa(\mathbf{J}(\mathbf{x}_k)) \cdot \lvert\lvert \mathbf{f}(\mathbf{x}_k)\rvert\rvert_2 \)\( If we terminate under the condition \)||\mathbf{f}(\mathbf{x}_f)||_2 < \varepsilon\(, then \)\( \lvert\lvert\mathbf{\delta x}_f\rvert\rvert_2 \leq \kappa(\mathbf{J}(\mathbf{x}_f)) \cdot \varepsilon \)\( With a condition number on the order of \)10^{13}\( and the default tolerance \)\varepsilon = 10^{-6}\(, then the difference between the computed solution \)x_f\( and the true solution \)x\( could be on the order of \)10^7$.

In practice, however, we frequently end up with much better results than those predicted by this worst-case bound. The condition number of the unscaled Jacobian does not taken into account that some variables are much bigger than others. For example, the liquid phase mole fraction of \(\text{CO}_2\) is on the order of \(10^{-5}\), whereas the heat duty is on the order of \(10^6\). Unless we scale the model, we have no way of determining which changes to the solution are significant and which are not. Additionally, we frequently find that most constraints have residuals on the order of \(10^{-10}\) or less, with only a few problematic constraints having larger residuals. That is the case here.

import idaes.core.util.model_statistics as mstat

print(f"Number of constraints: {mstat.number_activated_constraints(m)}")
dt.config.constraint_residual_tolerance = 1e-10
dt.display_constraints_with_large_residuals()

Even with a non-optimal termination, we have only three constraints with residuals greater than \(10^{-5}\) for a model with 77 active constraints. Two of those constraints are on the order of \(10^{-5}\), but the third is on the order of \(10^4\).

Nevertheless, a condition number on the order of \(10^{13}\) is significantly higher than we’d like, because large condition numbers can degrade the convergence of Newton’s method or even prevent it from converging entirely. Ideally, we want flowsheets to have condition numbers less than \(10^8\). Connecting unit models together with things like recycle streams can amplify the overall flowsheet condition number to be multiple orders of magnitude greater than that of its individual components. Consequently, we want unit models to have condition numbers less than \(10^4\).

Model scaling can both reduce the Jacobian’s condition number and help ensure that constraints are satisfied to an appropriate precision. Both of those effects can help the model initialization to converge. Before beginning to write a scaler object, however, we should look at the structure of the unit model.

def build(self):
    """Build the model.

    Args:
        None
    Returns:
        None
    """
    # Call UnitModel.build to setup dynamics
    super().build()

    # Check phase lists match assumptions
    if self.config.vapor_property_package.phase_list != ["Vap"]:
        raise ConfigurationError(
            f"{self.name} SolventReboiler model requires that the vapor "
            f"phase property package have a single phase named 'Vap'"
        )
    if self.config.liquid_property_package.phase_list != ["Liq"]:
        raise ConfigurationError(
            f"{self.name} SolventReboiler model requires that the liquid "
            f"phase property package have a single phase named 'Liq'"
        )

    # Check for at least one common component in component lists
    if not any(
        j in self.config.vapor_property_package.component_list
        for j in self.config.liquid_property_package.component_list
    ):
        raise ConfigurationError(
            f"{self.name} SolventReboiler model requires that the liquid "
            f"and vapor phase property packages have at least one "
            f"common component."
        )

    # ---------------------------------------------------------------------
    # Add Control Volume for the Liquid Phase
    self.liquid_phase = ControlVolume0DBlock(
        dynamic=self.config.dynamic,
        has_holdup=self.config.has_holdup,
        property_package=self.config.liquid_property_package,
        property_package_args=self.config.liquid_property_package_args,
    )

    self.liquid_phase.add_state_blocks(has_phase_equilibrium=True)

    # Separate liquid and vapor phases means that phase equilibrium will
    # be handled at the unit model level, thus has_phase_equilibrium is
    # False, but has_mass_transfer is True.
    self.liquid_phase.add_material_balances(
        balance_type=self.config.material_balance_type,
        has_mass_transfer=True,
        has_phase_equilibrium=False,
    )

    # Need to include enthalpy transfer term for the mass transfer
    self.liquid_phase.add_energy_balances(
        balance_type=self.config.energy_balance_type,
        has_heat_transfer=True,
        has_enthalpy_transfer=True,
    )

    self.liquid_phase.add_momentum_balances(
        balance_type=self.config.momentum_balance_type,
        has_pressure_change=self.config.has_pressure_change,
    )

    # ---------------------------------------------------------------------
    # Add single state block for vapor phase
    tmp_dict = dict(**self.config.vapor_property_package_args)
    tmp_dict["has_phase_equilibrium"] = False
    tmp_dict["defined_state"] = False
    self.vapor_phase = self.config.vapor_property_package.build_state_block(
        self.flowsheet().time, doc="Vapor phase properties", **tmp_dict
    )

    # ---------------------------------------------------------------------
    # Check flow basis is compatible
    t_init = self.flowsheet().time.first()
    if (
        self.vapor_phase[t_init].get_material_flow_basis()
        != self.liquid_phase.properties_out[t_init].get_material_flow_basis()
    ):
        raise ConfigurationError(
            f"{self.name} vapor and liquid property packages must use the "
            f"same material flow basis."
        )

    # ---------------------------------------------------------------------
    # Add Ports for the reboiler
    self.add_inlet_port(name="inlet", block=self.liquid_phase, doc="Liquid feed")
    self.add_outlet_port(name="bottoms", block=self.liquid_phase, doc="Bottoms stream")
    self.add_outlet_port(
        name="vapor_reboil",
        block=self.vapor_phase,
        doc="Vapor stream from reboiler",
    )

    # ---------------------------------------------------------------------
    # Add unit level constraints
    # First, need the union and intersection of component lists
    all_comps = (
        self.vapor_phase.component_list
        | self.liquid_phase.properties_out.component_list
    )
    common_comps = (
        self.vapor_phase.component_list
        & self.liquid_phase.properties_out.component_list
    )

    # Get units for unit conversion
    vunits = self.config.vapor_property_package.get_metadata().get_derived_units
    lunits = self.config.liquid_property_package.get_metadata().get_derived_units
    flow_basis = self.vapor_phase[t_init].get_material_flow_basis()
    if flow_basis == MaterialFlowBasis.molar:
        fb = "flow_mole"
    elif flow_basis == MaterialFlowBasis.mass:
        fb = "flow_mass"
    else:
        raise ConfigurationError(
            f"{self.name} SolventReboiler only supports mass or molar "
            f"basis for MaterialFlowBasis."
        )

    if any(j not in common_comps for j in self.vapor_phase.component_list):
        # We have non-condensable components present, need zero-flow param
        self.zero_flow_param = Param(
            mutable=True, default=1e-8, units=vunits("flow_mole")
        )

    # Material balances
    def rule_material_balance(blk, t, j):
        if j in common_comps:
            # Component is in equilibrium
            # Mass transfer equals vapor flowrate
            return -blk.liquid_phase.mass_transfer_term[t, "Liq", j] == pyunits.convert(
                blk.vapor_phase[t].get_material_flow_terms("Vap", j),
                to_units=lunits(fb),
            )
        elif j in self.vapor_phase.component_list:
            # Non-condensable component
            # No mass transfer term
            # Set vapor flowrate to an arbitrary small value
            return (
                blk.vapor_phase[t].get_material_flow_terms("Vap", j)
                == blk.zero_flow_param
            )
        else:
            # Non-vaporisable component
            # Mass transfer term is zero, no vapor flowrate
            return blk.liquid_phase.mass_transfer_term[t, "Liq", j] == 0 * lunits(fb)

    self.unit_material_balance = Constraint(
        self.flowsheet().time,
        all_comps,
        rule=rule_material_balance,
        doc="Unit level material balances",
    )

    # Phase equilibrium constraints
    # For all common components, equate fugacity in vapor and liquid
    def rule_phase_equilibrium(blk, t, j):
        return blk.liquid_phase.properties_out[t].fug_phase_comp[
            "Liq", j
        ] == pyunits.convert(
            blk.vapor_phase[t].fug_phase_comp["Vap", j], to_units=lunits("pressure")
        )

    self.unit_phase_equilibrium = Constraint(
        self.flowsheet().time,
        common_comps,
        rule=rule_phase_equilibrium,
        doc="Unit level phase equilibrium constraints",
    )

    # Temperature equality constraint
    def rule_temperature_balance(blk, t):
        return blk.liquid_phase.properties_out[t].temperature == pyunits.convert(
            blk.vapor_phase[t].temperature, to_units=lunits("temperature")
        )

    self.unit_temperature_equality = Constraint(
        self.flowsheet().time,
        rule=rule_temperature_balance,
        doc="Unit level temperature equality",
    )

    # Unit level energy balance
    # Energy leaving in vapor phase must be equal and opposite to enthalpy
    # transfer from liquid phase
    def rule_energy_balance(blk, t):
        return -blk.liquid_phase.enthalpy_transfer[t] == pyunits.convert(
            blk.vapor_phase[t].get_enthalpy_flow_terms("Vap"),
            to_units=lunits("energy") / lunits("time"),
        )

    self.unit_enthalpy_balance = Constraint(
        self.flowsheet().time,
        rule=rule_energy_balance,
        doc="Unit level enthalpy_balance",
    )

    # Pressure balance constraint
    def rule_pressure_balance(blk, t):
        return blk.liquid_phase.properties_out[t].pressure == pyunits.convert(
            blk.vapor_phase[t].pressure, to_units=lunits("pressure")
        )

    self.unit_pressure_balance = Constraint(
        self.flowsheet().time,
        rule=rule_pressure_balance,
        doc="Unit level pressure balance",
    )

    # Set references to balance terms at unit level
    self.heat_duty = Reference(self.liquid_phase.heat[:])

    if (
        self.config.has_pressure_change is True
        and self.config.momentum_balance_type != MomentumBalanceType.none
    ):
        self.deltaP = Reference(self.liquid_phase.deltaP[:])

As is typical with IDAES models, two submodels are created: the ControlVolume0D named liquid_phase and the StateBlock named vapor_phase. Then several unit model level constraints are written.

Step 3: Creating a New Scaler Class

To create a new scaling routine for the solvent reboiler, we start by creating a new Scaler class which inherits from the CustomScalerBase class in idaes.core.scaling.

from idaes.core.scaling import CustomScalerBase


class SolventReboilerScaler(CustomScalerBase):
    def variable_scaling_routine(
        self, model, overwrite: bool = False, submodel_scalers: dict = None
    ):
        # Empty method for now
        pass

    def constraint_scaling_routine(
        self, model, overwrite: bool = False, submodel_scalers: dict = None
    ):
        # Empty method for now
        pass

As you know from the first scaling tutorial, the variable_scaling_routine and constraint_scaling_routine methods are used to implement subroutines for scaling the variables and constraints in the model respectively. Separately, there is a scale_model method that will call each of these in sequence in order to scale an entire model by applying the following steps:

1. apply variable scaling routine,

2. apply first stage scaling fill-in,

3. apply constraint scaling routine,

4. apply second stage scaling fill-in.

The second and fourth steps are intended to allow users to provide methods to fill in missing scaling information that was not provided by the first and second steps. However, these methods are prone to perform poorly on all but the simplest models (see Step 5 for more information).

Step 4: Apply Scaling to Sub-Models

First, let’s look at how to scale the control volume and state block sub-models. Because the property package for the vapor_phase state block is specified by the user, we do not know what variables and constraints it may contain, so we cannot (and should not) scale it directly. The property package creator (hopefully) created a scaler object that we can use. The zero-dimensional control volume liquid_phase contains two state blocks, material, energy, and pressure balance constraints, and variables corresponding to additional interaction terms. Because we are choosing which terms to create in the control volume, we could, in principle, scale it at the unit model level. However, because control volumes typically should be scaled in the same way, a submodel scaler has been provided for it as well. Thus, what we want to do here is to call the variable and constraint scaling routines from the Scaler associated with each sub-model, which we can do using the call_submodel_scaler_method method.

class SolventReboilerScaler(CustomScalerBase):
    def variable_scaling_routine(
        self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None
    ):
        # Call scaling methods for sub-models
        self.call_submodel_scaler_method(
            submodel=model.liquid_phase,
            method="variable_scaling_routine",
            submodel_scalers=submodel_scalers,
            overwrite=overwrite,
        )

        self.call_submodel_scaler_method(
            submodel=model.vapor_phase,
            method="variable_scaling_routine",
            submodel_scalers=submodel_scalers,
            overwrite=overwrite,
        )

    def constraint_scaling_routine(
        self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None
    ):
        # Call scaling methods for sub-models
        self.call_submodel_scaler_method(
            submodel=model.liquid_phase,
            method="constraint_scaling_routine",
            submodel_scalers=submodel_scalers,
            overwrite=overwrite,
        )

        self.call_submodel_scaler_method(
            submodel=model.vapor_phase,
            method="constraint_scaling_routine",
            submodel_scalers=submodel_scalers,
            overwrite=overwrite,
        )

In order for the submodel scaling to function properly, default scaling factors need to be set for some variables. The property package scaler cannot know if the system is supposed to be bench scale, pilot scale, or production scale unless the user sets a default for flow rate. These defaults can be set on a unit model by unit model basis by using the submodel_scalers argument, or they can be set at the flowsheet level on the parameter block. We will do the latter.

First, we need to create objects using the modular property submodel scaler class.

liquid_properties_scaler = m.fs.liquid_properties.default_state_scaler_class()
vapor_properties_scaler = m.fs.vapor_properties.default_state_scaler_class()

Now, let’s see what default scaling factors these scaler objects have.

liquid_properties_scaler.default_scaling_factors

We can see that we are required to set a default scaling factor for flow_mol_phase. Furthermore, it is recommended to set a default scaling factor for enth_mol_phase, visc_d_phase, therm_cond_phase, dens_mol_phase, andvol_mol_phase. For these latter quantities, there is a method to estimate reasonable scaling factors, but it might not be suitible for all cases.

With an inlet flow of about 80 mol/s, a scaling factor of 1/80 is appropriate. We can ignore visc_d_phase and therm_cond_phase, because this unit model does not require dynamic viscosity or thermal conductivity. Molar enthalpy is a tougher quantity to scale.

By convention, the molar enthalpy of each element in its pure form at \(25^\circ\text{C}\) and \(1\:\text{atm}\) is set equal to zero. In some property packages, enthalpy of formation is taken into account for compounds, while in some property packages it is not. However, the scaling factor for enthalpy should not depend on this convention—only enthalpy differences are meaningful. So we need to determine how to determine what constitutes a “significant” enthalpy difference for this property package. The way the modular property scaler object estimates this is through the temperature scaling factor, which determines what a “significant” temperature difference is, and the observation that a wide range of substances have a specific heat capacity around \(2\:\frac{\text{J}}{\text{g K}}\).

This default method is good enough for our purposes now. We can revisit it later, if necessary.

liquid_properties_scaler.default_scaling_factors["flow_mol_phase"] = 1 / 80
vapor_properties_scaler.default_scaling_factors["flow_mol_phase"] = 1 / 10

m.fs.liquid_properties.default_state_scaler_object = liquid_properties_scaler
m.fs.vapor_properties.default_state_scaler_object = vapor_properties_scaler

Now we can see how much submodel scaling improves the conditioning of the unit model.

reboiler_scaler = SolventReboilerScaler()
reboiler_scaler.scale_model(m.fs.unit)
dt = DiagnosticsToolbox(m)
dt.report_numerical_issues()

We see in this case that partly scaling the model improved the condition number of the Jacobian by only two and a half orders of magnitude. This sort of marginal improvement is common for partial scaling, and sometimes partial scaling makes the model’s condition number worse. When variables are scaled without the constraints they appear in also being scaled, it generates new extreme Jacobian entries. If we look at the number of extreme Jacobian rows and columns, however, we see the progress we’ve made. Before applying submodel scaling, we had 27 variables with extreme column norms and 16 constraints with extreme row norms. Now we have only 10 variables and 4 constraints.

We are concerned with Jacobian rows and columns with extreme values because they allow us to estimate the Jacobian condition number. Determining the (Euclidean) condition number requires the computation of the singular value decomposition (SVD) of \(\mathbf{J}(\mathbf{x})\). However, that process is computationally expensive for large models, and it does not identify individual variables and constraints as problematic (instead, it identifies combinations of variables and constraints as problematic). However, we have a lower bound to the condition number in terms of the norms of rows and columns of the Jacobian matrix: $\( \max\left(\frac{||\mathbf{r}_\text{max}||_2}{||\mathbf{r}_\text{min}||_2}, \frac{||\mathbf{c}_\text{max}||_2}{||\mathbf{c}_\text{min}||_2}\right) \leq \kappa(\mathbf{J}(\mathbf{x})) \)$ Consequently, Jacobian rows or columns with extremely large or extremely small norms indicate that the Jacobian is ill-conditioned. Each row is associated with a constraint and each column is associated with a variable. Thus, the extreme Jacobian rows and columns give us a list of variables and constraints to investigate for scaling.

Just because a variable or constraint has a Jacobian column or row with an extreme norm does not mean it is badly-scaled. Frequently, either a constraint containing the variable or a variable present in a constraint is the actual problem.

NOTE

Jacobian rows and columns with extreme norms reveal that the Jacobian matrix is ill-conditioned, but the converse is not true. A Jacobian matrix can have no rows and columns with extreme norms but still be singular or ill-conditioned.

Let’s see which variables and constraints are problematic.

dt.display_variables_with_extreme_jacobians()
dt.display_constraints_with_extreme_jacobians()

These four unit model level constraints are causing our remaining problems.

Step 5: Unit Model Level Scaling

We should take a look at the constraints before determining their scaling factors. However, Pyomo’s pprint method often prints out multiple terminal pages’ worth of thermodynamic expressions, which are extremely difficult to parse. The output is so verbose because it descends into named Expressions and recursively subsititutes them into the constraint expression. Since IDAES makes heavy use of named Expressions to reduce the number of Constraints in a model, the output becomes unmanangable for all but the most simple constraints. The print_compact_form utility function usually results in much more readable output.

from idaes.core.util.misc import print_compact_form

print_compact_form(m.fs.unit.unit_enthalpy_balance[0.0])

This is a relatively straightforward constraint.

So long as we have scaling factors set for the terms in the enthalpy balance, we can scale the constraint by the same factor. Because liquid_phase.enthalpy_transfer[0.0] is a Var created by the ControlVolume0D, it already has a scaling factor assigned by the ControlVolume0DScaler.

get_scaling_factor(m.fs.unit.liquid_phase.enthalpy_transfer[0.0])

However, in the modular property framework, get_enthalpy_flow_terms(p) returns a named Expression. Expressions do not need to be scaled on their own. However, scaling hints can be assigned for them to assist in determining appropriate scaling factors for other variables and constraints. When the get_scaling_factor function is called on an Expression, it will return a scaling hint if one is assigned. (This behavior is useful, because what is a Var in one property package can be implemented as an Expression in another.)

In this case, however, no scaling hint is assigned to vapor_phase[0.0]._enthalpy_flow_term["Vap"].

print(get_scaling_factor(m.fs.unit.vapor_phase[0.0]._enthalpy_flow_term["Vap"]))

To get an estimate of the size of that expression, then, we can look at its constituent terms:

print_compact_form(m.fs.unit.vapor_phase[0.0]._enthalpy_flow_term["Vap"])

The term vapor_phase[0.0].flow_mol_phase["Vap"] is a Var and has a scaling factor applied.

get_scaling_factor(m.fs.unit.vapor_phase[0.0].flow_mol_phase["Vap"])

However, vapor_phase[0.0].enth_mol_phase["Vap"] is a named Expression. Nevertheless, it has a scaling hint.

get_scaling_factor(m.fs.unit.vapor_phase[0.0].enth_mol_phase["Vap"])

Because the scaling factor of the molar flow rate is \(0.1\), we can estimate that the molar flow rate will be about \(10\). Similarly, we can estimate the size of a significant molar enthalpy difference is around the size of 1/5.46e-5\( \approx 18,000\). Therefore the general size of the enthalpy flow terms should be around \(180,000\).

get_scaling_factor(m.fs.unit.liquid_phase.properties_out[0.0].flow_mol_phase["Liq"])

The get_sum_terms_nominal_values method in CustomScalerBase which uses an expression walker to go through an expression to return a list of the expected magnitude (or “nominal value”) of all additive terms in the expression based on the scaling factors for the variables involved. For example, suppose we have an expression $\( a + b\cdot c - d\cdot (e - f) \)\( Let \)\sigma(v)\( denote the scaling factor of the variable \)v\(. Then `get_sum_terms_nominal_values` substitutes the inverse of the variable scaling factor into each term in the sum and returns it as a list: \)\( \left[\frac{1}{\sigma(a)},\; \frac{1}{\sigma(b) \cdot \sigma(c)}, \frac{-\left(\frac{1}{\sigma(e)} - \frac{1}{\sigma(f)}\right)}{\sigma(d)} \right] \)$

This process can produce reliable estimates of expression magnitude for certain operations:

• Addition of positive numbers

• Multiplication and division

• Raising a variable to a fixed power

For these operations, only an order-of-magnitude estimate for the variable is sufficient to get an order-of-magnitude estimate of the resulting expression. However, it is far less reliable for other operations:

• Subtraction (or addition of negative numbers to positive numbers)

• Functions with a variable exponent

• Trigometric functions

• Logarithms

For these operations, an expression’s final value is highly dependent on the exact values of its included variables. For example, suppose \(\sigma(e) = \sigma(f)\). In that case, the term \(d\cdot (e - f)\) would have an estimated magnitude of \(0\), which is useless for scaling purposes. It is in cases like this that scaling hints for named expressions. If we denote \(G = e - f\) and assign \(\sigma(G) = \sigma(e)\), we have a far better estimate of the term \(d\cdot (e - f)\) in the quantity \(1/(\sigma(d)\cdot\sigma(G))\).

Logarithms are a special case. Since they convert multiplication to addition, application of a scaling factor simply shifts the zero point of the resulting expression. $\( \log(h \cdot \sigma(h)) = \log(h) + \log(\sigma(h)) \)$

The output of a logarithmic expression should be considered well-scaled by default, unless its argument is an absolutely enormous number (in which case floating point arithmetic is probably inappropriate anyway).

In the case of scaling unit.unit_enthalpy_balance[0.0], we only have two additive terms, one consisting of a single variable, the other consisting of the product of a variable with an expression with a scaling hint. In this case, get_sum_terms_nominal_values should produce a good estimate of the order-of-magnitude of each term.

reboiler_scaler.get_sum_terms_nominal_values(m.fs.unit.unit_enthalpy_balance[0.0])

The inverse of either of these numbers would be reasonable scaling factors. The scale_constraint_by_nominal value function automatically scales a constraint by using one of these numbers. Which number we use is selected using the ConstraintScalingScheme Enum. If we take the inverse of the larger number, we are saying that the constraint needs to be satisfied with the same precision as the largest term in the sum. This corresponds to ConstraintScalingScheme.inverseMaximum. If we take the inverse of the smaller number, we are saying that the constraint needs to be satisfied with the same precision as the smallest term in the sum. This corresponds to ConstraintScalingScheme.inverseMinimum.

In general, the inverseMinimum scheme is the safest to use, because it demands the most precision. Let’s use the inverseMinimum method for this constraint in order to more precisely determine the vapor phase’s temperature.

The next two constraints can be dealt with together:

print_compact_form(m.fs.unit.unit_phase_equilibrium[0.0, "CO2"])
print_compact_form(m.fs.unit.unit_phase_equilibrium[0.0, "H2O"])

These are simple vapor-liquid equilibrium constraints for the volatile components. Let’s see what get_sum_terms_nominal_values shows for them:

print(
    reboiler_scaler.get_sum_terms_nominal_values(
        m.fs.unit.unit_phase_equilibrium[0.0, "CO2"]
    )
)
print(
    reboiler_scaler.get_sum_terms_nominal_values(
        m.fs.unit.unit_phase_equilibrium[0.0, "H2O"]
    )
)

Here, we see that, while the vapor phase fugacities are associated with a nominal value of 10000, the liquid phase fugacities have wildly different nominal values. We do expect the vapor phase to be enriched in \(\text{CO}_2\), but not by a factor of \(2\cdot 10^4\). Let’s unpack these values. We are using the ideal gas law for the vapor phase, so the fugacity is simply the product of the mole fraction and vapor pressure:

print_compact_form(m.fs.unit.vapor_phase[0.0].fug_phase_comp["Vap", "CO2"])

The liquid phase, by contrast, contains a complicated Henry’s Law type relationship for \(\text{CO}_2\) and a slightly less complicated Antoine’s Law relationship for \(\text{H}_2\text{O}\):

print_compact_form(
    m.fs.unit.liquid_phase.properties_out[0.0].fug_phase_comp["Liq", "CO2"]
)

Parsing the entire pile of thermodynamic gibberish is unnecessary to plan our next move. We can see that the liquid phase fugacity expression is rife with subtraction and exponentiation, both of which are dangerous for the nominal value expression walker. Even worse, we are scaling based on a failed solution, so the values may be off anyway. It is far better to base the scaling factor on the vapor phase fugacity. The get_expression_nominal_value function allows the user to get the nominal value of a single expression.

reboiler_scaler.get_expression_nominal_value(
    m.fs.unit.vapor_phase[0.0].fug_phase_comp["Vap", "CO2"]
)

While in this case this value is positive, the nominal value can be negative, so you should use the absolute value when getting a scaling factor.

Finally, we have the last constraint: mechanical equilibrium between phases.

print_compact_form(m.fs.unit.unit_pressure_balance[0.0])

The inverseMinimum scheme will work here as well.

In addition to these problematic unit model level constraints, there are some additional constraints that need scaling factors. Although they might not be problematic under these process conditions, they may be for a scaled-up flowsheet. If we suspect they are well-scaled already, we should use the report_scaling_factors function to show what they are:

from idaes.core.scaling.util import report_scaling_factors

report_scaling_factors(m.fs.unit, descend_into=False)

We have yet to scale the unit_temperature_equality and unit_material_balanceconstraints. The temperature equality constraint is simple.

print_compact_form(m.fs.unit.unit_temperature_equality[0.0])

inverseMinimum works fine here.

The material balance constraints have different forms depending on whether a component occurs in the liquid phase, the vapor phase, or both:

print_compact_form(m.fs.unit.unit_material_balance[0.0, "CO2"])
print_compact_form(m.fs.unit.unit_material_balance[0.0, "MEA"])
print_compact_form(m.fs.unit.unit_material_balance[0.0, "N2"])

The first constraint equates the vapor flow rate to the mass transfer from the liquid phase. Either inverse maximum or inverse minimum would work. The next two are a bit trickier. For \(\text{MEA}\), which is considered nonvolatile in this property package, there is no mass transfer to the vapor phase. A common misconception is that variables should always be scaled to be order one. That is incorrect. We only need to calculate mass_transfer_term[0.0,"Liq","MEA"] to the same precision as the other (nonzero) terms of the mass balance (and the ControlVolume0D submodel scaler does just that). Similarly, flow_mol_phase_comp["Vap","N2"] is set to a tiny nonzero value (because flows that are exactly zero often cause problems for property calculations). We do not need more significant figures for either of those variables than the rest of the terms in their corresponding material balances. Therefore we want to use the inverseMaximum because inverseMinimum would demand unnecessary precision.

NOTE

If we had relied on autoscaling methods to fill in scaling factors for these variables and constraints, we would have ended up with enormous scaling factors for flow_mol_phase_comp["Vap","N2"] and flow_mol_phase_comp["Vap","O2"] because they are not-quite equal to zero at the model solution.

Similarly, m.fs.unit.unit_phase_equilibrium[0.0,"CO2"] would have been overscaled due to the enormous nominal value given by the expression walker for the liquid phase fugacity

Finally, we can complete the scaler object.

from idaes.core.scaling.custom_scaler_base import ConstraintScalingScheme


class SolventReboilerScaler(CustomScalerBase):
    def variable_scaling_routine(
        self, model, overwrite: bool = False, submodel_scalers: dict = None
    ):
        # Call scaling methods for sub-models
        self.call_submodel_scaler_method(
            submodel=model.liquid_phase,
            method="variable_scaling_routine",
            submodel_scalers=submodel_scalers,
            overwrite=overwrite,
        )

        self.call_submodel_scaler_method(
            submodel=model.vapor_phase,
            method="variable_scaling_routine",
            submodel_scalers=submodel_scalers,
            overwrite=overwrite,
        )

    def constraint_scaling_routine(
        self, model, overwrite: bool = False, submodel_scalers: dict = None
    ):
        # Call scaling methods for sub-models
        self.call_submodel_scaler_method(
            submodel=model.liquid_phase,
            method="constraint_scaling_routine",
            submodel_scalers=submodel_scalers,
            overwrite=overwrite,
        )

        self.call_submodel_scaler_method(
            submodel=model.vapor_phase,
            method="constraint_scaling_routine",
            submodel_scalers=submodel_scalers,
            overwrite=overwrite,
        )
        # Note: scaling factors cannot be applied directly to indexed objects.
        # They should be applied to VarData/ConstraintData/ExpressionData children instead.
        for condata in model.unit_material_balance.values():
            self.scale_constraint_by_nominal_value(
                condata,
                scheme=ConstraintScalingScheme.inverseMaximum,
                overwrite=overwrite,
            )
        for condata in model.unit_temperature_equality.values():
            self.scale_constraint_by_nominal_value(
                condata,
                scheme=ConstraintScalingScheme.inverseMinimum,
                overwrite=overwrite,
            )
        for condata in model.unit_enthalpy_balance.values():
            self.scale_constraint_by_nominal_value(
                condata,
                scheme=ConstraintScalingScheme.inverseMinimum,
                overwrite=overwrite,
            )
        for condata in model.unit_pressure_balance.values():
            self.scale_constraint_by_nominal_value(
                condata,
                scheme=ConstraintScalingScheme.inverseMinimum,
                overwrite=overwrite,
            )
        for (t, j), condata in m.fs.unit.unit_phase_equilibrium.items():
            nom = abs(
                self.get_expression_nominal_value(
                    m.fs.unit.vapor_phase[t].fug_phase_comp["Vap", j]
                )
            )
            self.set_constraint_scaling_factor(condata, 1 / nom, overwrite=overwrite)

Now we can see how this works on the model. We will create a new instance of the model to demonstrate that the improved scaling is sufficient to allow the model to initialize without a custom-tailored initial guess.

m = create_model()

# Since this is a new instance of the model, we need to set the default
# scaler objects again.
m.fs.liquid_properties.default_state_scaler_object = liquid_properties_scaler
m.fs.vapor_properties.default_state_scaler_object = vapor_properties_scaler

scaler_obj = SolventReboilerScaler()
scaler_obj.scale_model(m.fs.unit)

reboiler_init = m.fs.unit.default_initializer(always_estimate_states=True)
reboiler_init.initialize(m.fs.unit)

We see that scaling improved the robustness of the model enough so that initialization could succeed. Now let’s look at the numerical issues.

dt = DiagnosticsToolbox(m)
dt.report_numerical_issues()

We have reduced the Jacobian condition number by another four orders of magnitude, as well as eliminating all extreme Jacobian rows and columns. A condition number of \(5\cdot 10^6\) is still higher than we’d like for a single unit model, especially for a simple unit model like this one, because connecting units together as part of a flowsheet inevitably increases the condition number.

Step 6: Singular Value Decomposition

With no more extreme Jacobian rows or columns, we have two avenues to proceed:

1. Check the extreme Jacobian entries

2. Check the singular value decomposition (SVD) of the Jacobian

Extreme Jacobian entries only suggest that the variable and constraint might be contributing to ill-conditioning—they might also be false positives. However, extreme singular values are guaranteed to be contributing to the Jacobian’s ill-conditioning. Therefore, we will look at the SVD next.

The singular value decomposition factorizes a matrix into three parts: $\( \mathbf{J} = \mathbf{U \Sigma V}^T \)\( Both \)\mathbf{U}\( and \)\mathbf{V}\( are orthogonal matrices, i.e., \)\mathbf{U}^T \mathbf{U} = \mathbf{I}\( and \)\mathbf{V}^T\mathbf{V} = \mathbf{I}\(. If \)\mathbf{J}\( is square, then \)\mathbf{\Sigma}\( is a diagonal matrix of the form: \)\( \mathbf{\Sigma} = \begin{bmatrix} \sigma_1 & 0 & \dots & 0 \\ 0 & \sigma_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \dots & \sigma_n \end{bmatrix} \)\( in which \)\sigma_j\( are the singular values. The singular values are arranged in descending order, i.e., \)\sigma_{j+1} \leq \sigma_{j}\(. If any \)\sigma_j=0$, then the matrix is singular.

The condition number is given by: $\( \kappa(\mathbf{J}) = \sigma_1 / \sigma_n \)$ Thus, extremely large and extremely small singular values cause matrix ill-conditioning. However, extremely large singular values typically show up as rows and columns with extreme norms. Therefore we will look at the smallest singular values.

Because \(\mathbf{U}\) and \(\mathbf{V}\) are orthogonal and \(\mathbf{\Sigma}\) is diagonal, computing a Newton step is easy once the SVD is calculated. $\( \mathbf{\delta x}_k = -\mathbf{V \Sigma}^{-1} \mathbf{U}^T \cdot \mathbf{f}(\mathbf{x}_k) \)\( Since we have to refactorize \)\mathbf{J}\( at each iteration, the SVD is not an efficient method to calculate a Newton step. However, this process has a geometric interpretation that is useful for model diagnostics and scaling. First, the constraint residual is decomposed into orthogonal components by \)\mathbf{U}\(: \)\( \mathbf{s} = -\mathbf{U}^T \cdot \mathbf{f}(\mathbf{x}_k) \)\( Next, these orthogonal components are scaled by the inverse of the associated singular value: \)\( \mathbf{\tilde{s}} = \mathbf{\Sigma}^{-1} \cdot \mathbf{s} \)\( Finally, these independent components are translated into the variable space by \)\mathbf{V}\(: \)\( \mathbf{\delta x}_k = \mathbf{V} \cdot \mathbf{\tilde{s}} \)\( So each orthogonal component of the constraint residual space is associated with an orthogonal component of the variable space. In particular, each singular value \)\sigma_j\( is associated with a left singular vector \)\mathbf{u}_j\( and right singular vector \)\mathbf{v}_j\(, the \)j\text{th}\( columns of \)\mathbf{U}\( and \)\mathbf{V}\(, respectively. Any constraint error in the direction of \)\mathbf{u}_j\( *must* be satisfied by a corresponding change in the direction of \)\mathbf{v}_j\(. When \)\sigma_j\( is extremely small, tiny changes in the direction of \)\mathbf{u}_j\( translate into enormous changes in \)\mathbf{v}_j\(. This dysfunctional relationship can be a sign of bad scaling, but also can be a sign of some other model issue like numerical singularity. Typically, singular values in the range \)10^{-12}\( to \)10^{-4}\( are signs of either bad scaling or a highly sensitive model (for example, a flowsheet with a large recycle stream). Values below \)10^{-12}$ are indicative of numerical singularity.

NOTE

The singular vectors \(\mathbf{u}_j\) and \(\mathbf{v}_j\) are typically dense. In order to link them with specific variables and constraints, a thresholding scheme is applied. The size_cutoff_in_singular_vector config option for the SVDToolbox controls the magnitude of an entry in \(\mathbf{u}_j\) and \(\mathbf{v}_j\) that is sufficient to link a variable or constraint to the \(\sigma_j\). It is 0.1 by default, but typically has a value between 0.05 and 0.3. Decrease it to show fewer variables and constraints, increase it to show more variables and constraints.

So let’s create an instance of the SVDToolbox and display the variables and constraints associated with the 5 smallest singular values. (For a full description of all the methods available, look at the SVD toolbox documentation.)

st = dt.prepare_svd_toolbox()
# The SVD toolbox displays 10 singular values by default, but we are going to
# view the only the first five to reduce the amount of output.
st.display_underdetermined_variables_and_constraints(singular_values=[1, 2, 3, 4, 5])

The smallest singular value is 40 times smaller than the second smallest singular value. This usually means that we can significantly improve the matrix conditioning by scaling only a small number of variables and constraints. There is only one constraint strongly associated with the smallest singular value, so let’s start with that.

print_compact_form(
    m.fs.unit.liquid_phase.properties_out[0.0].log_conc_mol_phase_comp_true_eq[
        "Liq", "CO2"
    ]
)

This constraint converts between the \(\text{CO}_2\) concentration and its logarithm. Using these sorts of log-form variables can increase the number of significant digits when dealing with reactions involving components in small quantities. (They can also introduce degeneracy when dealing with trace quantities that do not appear in reactions.)

The SVD Toolbox also has two helper functions: display_variables_in_constraint, which displays the Jacobian row corresponding to a constraint, and display_constraints_including_variable, which displays the Jacobian column corresponding to a variable.

st.display_variables_in_constraint(
    m.fs.unit.liquid_phase.properties_out[0.0].log_conc_mol_phase_comp_true_eq[
        "Liq", "CO2"
    ]
)

NOTE

At this point of the process, the author of this notebook found a shortcoming in the ModularPropertiesScaler. No scaling hint was set for the molar density expression, so the expression walker descended into it and came back with an inappropriate nominal value. Revising the scaler to estimate the phase density resulted in a condition number improvement of 1-2 orders of magnitude.

Remember:

1. Scaling is an iterative process.

2. Submodel scalers may need to be revised to fix scaling issues.

3. Using the expression walker in an expression tree that involves subtraction or exponential functions is dangerous.

From our knowledge of the physical properties, mole_frac_phase_comp_true["Liq","CO2"] and log_conc_mol_phase_comp_true["Liq","CO2"] should have the largest effects on this equation. The other variables affect the concentration indirectly through the phase density. But while mole_frac_phase_comp_true["Liq","CO2"] does have an appropriately-sized impact, log_conc_mol_phase_comp_true["Liq","CO2"] is several orders-of-magnitude too small. A naive approach might be to decrease the scaling factor on log_conc_mol_phase_comp_true["Liq","CO2"] until it has an appropriately-sized effect in the equation. Remember, though, that log variables are well-scaled by default. Something else must be wrong.

Let’s take a look at the values on both sides of the equation.

print(
    f"Left hand side value: {value(exp(m.fs.unit.liquid_phase.properties_out[0.0].log_conc_mol_phase_comp_true['Liq','CO2'])):.4e}"
)
print(
    f"CO2 true concentration value: {value(m.fs.unit.liquid_phase.properties_out[0.0].conc_mol_phase_comp_true['Liq','CO2']):.4e}"
)

Since the left hand side and right hand side both have values of order one, we expect the equation to have a similar scaling factor.

print(
    get_scaling_factor(
        (
            m.fs.unit.liquid_phase.properties_out[0.0].log_conc_mol_phase_comp_true_eq[
                "Liq", "CO2"
            ]
        )
    )
)

However, it actually has an extremely small scaling factor. To what can we attribute this?

It is a result of not setting good default values for the true species mole fractions in the liquid phase. Although the apparent \(\text{CO}_2\) mole fraction is about 0.03, almost all of that is in the form of a carbamate salt formed through reaction with \(\text{MEA}\). The true mole fraction is significantly lower.

print(
    f"True CO2 mole fraction: {m.fs.unit.liquid_phase.properties_out[0].mole_frac_phase_comp_true['Liq','CO2'].value:4e}"
)
print(
    f"True MEACOO- mole fraction: {m.fs.unit.liquid_phase.properties_out[0].mole_frac_phase_comp_true['Liq','MEACOO_-'].value:4e}"
)

Nevertheless, the true value of the \(\text{CO}_2\) mole fraction is important because it is used in the calculation of fugacity. It should therefore be assigned a larger scaling factor. Therefore, let’s revisit the default scaling factors for the liquid properties:

default_scaling_factors = (
    m.fs.liquid_properties.default_state_scaler_object.default_scaling_factors
)
# Dictionaries are mutable, so changing its value here also changes it on the scaler object
# First scale the apparent mole fractions
default_scaling_factors["mole_frac_phase_comp[Liq, H2O]"] = 1
default_scaling_factors["mole_frac_phase_comp[Liq, CO2]"] = 10
default_scaling_factors["mole_frac_phase_comp[Liq, MEA]"] = 10
# Next scale the true mole fractions
default_scaling_factors["mole_frac_phase_comp_true[Liq, H2O]"] = 1
default_scaling_factors["mole_frac_phase_comp_true[Liq, CO2]"] = 1e4
default_scaling_factors["mole_frac_phase_comp_true[Liq, MEA]"] = 10
default_scaling_factors["mole_frac_phase_comp_true[Liq, MEA_+]"] = 10
default_scaling_factors["mole_frac_phase_comp_true[Liq, MEACOO_-]"] = 10
default_scaling_factors["mole_frac_phase_comp_true[Liq, HCO3_-]"] = 100

With these new scaling factors, let’s rescale and resolve the mode, then check the new condition number.

scaler_obj = SolventReboilerScaler(overwrite=True)
scaler_obj.scale_model(m.fs.unit)
reboiler_init.initialize(m.fs.unit)
dt = DiagnosticsToolbox(m)
dt.report_numerical_issues()

The condition number decreased by a factor of 18, which is less than what we might have hoped for. Let’s check out the SVD toolbox again.

st = dt.prepare_svd_toolbox()
st.display_underdetermined_variables_and_constraints([1, 2, 3, 4, 5])

We can see an order-of-magnitude difference between the first four smallest singular values and the fifth-smallest singular value. Inspecting the variables and constraints involved, they all involve the inherent reaction calculations for ionic speciation. Improving the scaling of those variables and constraints might improve the condition number by another one to two orders of magnitude. However, we are reaching the point of diminishing returns for scaling. We can leave the scaler here unless a solver failure forces us to revisit it again.

Step 7: Reviewing Model Scaling

Before we conclude entirely, let’s take a look at all the variable and constraint scaling factors used.

from idaes.core.scaling.util import report_scaling_factors

report_scaling_factors(m.fs.unit, descend_into=True)

Most variable values have scaled values within a couple orders of magnitude of one, which is what we like. The outliers are inherent reaction extents, which we already know are a bit problematic, and variables corresponding to the components \(\text{N}_2\) and \(\text{O}_2\). Because those components are present in the vapor property package but not the vapor stream, they should not be scaled to order one, and their current scaling factors are fine.

We have reduced the condition number to \(3 \cdot 10^5\), which is not quite as low as we would have liked, but is a vast improvement over the value \(3 \cdot 10^{13}\), which we started out with.

References

[1] Allan, D. A., Ostace, A. G., & Polley, T. (2025). Jacobian-based model diagnostics and application to equation oriented modeling of a carbon capture system. Computers & Chemical Engineering, 109312.

[2] Wächter, A., & Biegler, L. T. (2006). On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical programming, 106(1), 25-57.