Generic property package: estimating bubble and dew temperatures

[alma-walmsley]

Summary

I was having problems with the "Bubble, dew, and critical point initialization" step in the generic property package initialization. I ended up finding out that the initial guess for temperature_bubble was higher than that of temperature_dew, which I suspect was the key cause of infeasibility. I then had a look through the estimate_Tbub and estimate_Tdew methods in idaes.models.properties.modular_properties.base.utility. These use Newton's method to converge on a guess for temperature_bubble and temperature_dew respectively. I eventually realized that I was hitting the iteration limit of 30, which meant I got a bad guess which would cause major issues.

I was able to resolve this by increasing the iteration limit via idaes.models.properties.modular_properties.base.utility.MAX_ITER = 1000. This ended up giving me a much better guess since I was actually converging rather than stopping early. Note, I don't have much chemical engineering knowledge, so I am not sure if this is a sign of a modelling error.

I decided to open this issue since it took me a while to debug.

Potential actions:

a) Log a warning:

• At a minimum, I think a warning should be logged if MAX_ITER is reached in the estimate_Tbub and estimate_Tdew methods.

b) Increase default max iterations:

• Assume that failing to converge in MAX_ITER iterations will cause major issues, so increase the default number of max iterations. The argument against this is that a poorly-defined model (ie. one that never manages to converge) would take extra time to end up reaching the same result. Also, I can imagine if MAX_ITER were set to something like 100 instead of 30, I never would have realized there was an issue (and then we get "flaky" initialization when a more complex model is used, which is not fun). However, some combination of a) and b) may be useful.

c) Switch to a solver instead of using Newton's method

• Despite the added overhead, I was able to get increased performance (time-wise) by defining the model in pyomo and solving externally rather than solving directly in python. ipopt could solve in ~"3-5 iterations" compared to what was ~60 iterations when using Newton's method in python. I imagine it is more scalable too (for complex models).

For example,

def _custom_estimate_Tbub(blk, T_units, raoult_comps, henry_comps, liquid_phase):
    """
    Function to estimate bubble point temperature
    Args:
        blk: StateBlock to use
        T_units: units of temperature
        raoult_comps: list of components that follow Raoult's Law
        henry_comps: list of components that follow Henry's Law
        liquid_phase: name of liquid phase
    Returns:
        Estimated bubble point temperature as a float.
    """
    # Use lowest component temperature_crit as starting point
    # Starting high and moving down generally works better,
    # as it under-predicts next step due to exponential form of
    # Psat.
    # Subtract 1 to avoid potential singularities at Tcrit
    Tbub_initial = (
        min(blk.params.get_component(j).temperature_crit.value for j in raoult_comps)
        - 1
    )
    blk.tbub_initialize = Block()
    m = blk.tbub_initialize
    m.Tbub = Var(initialize=Tbub_initial)
    m.f = Expression(expr=(
        sum(
            get_method(blk, "pressure_sat_comp", j)(
                blk, blk.params.get_component(j), m.Tbub * T_units
            )
            * blk.mole_frac_comp[j]
            for j in raoult_comps
        )
        + sum(
            blk.mole_frac_comp[j]
            * blk.params.get_component(j)
            .config.henry_component[liquid_phase]["method"]
            .return_expression(blk, liquid_phase, j, m.Tbub * T_units)
            for j in henry_comps
        )
        - blk.pressure
    ))
    m.df = Expression(expr=(
        sum(
            get_method(blk, "pressure_sat_comp", j)(
                blk, blk.params.get_component(j), m.Tbub * T_units, dT=True
            )
            * blk.mole_frac_comp[j]
            for j in raoult_comps
        )
        + sum(
            blk.mole_frac_comp[j]
            * blk.params.get_component(j)
            .config.henry_component[liquid_phase]["method"]
            .dT_expression(blk, liquid_phase, j, m.Tbub * T_units)
            for j in henry_comps
        )
    ))
    m.tbub_constraint = Constraint(rule=(
        m.Tbub == m.Tbub - m.f / m.df
    ))
    solver = get_solver("ipopt")
    solver.options["tol"] = TOL  # default currently 1e-1
    solver.options["max_iter"] = MAX_ITER  # default currently 30
    solver.solve(m)

    # probably need to check for infeasibility here
    # ...

    del blk.tbub_initialize  # cleanup
    return m.Tbub.value

Feedback would be useful 🙂

[andrewlee94]

@alma-walmsley Hitting the maximum iteration limit is not necessarily an issue as long as it gets "close enough" - the idea is that you just need to get close enough that the main solve can take over and get the final solution. Thus, the iteration limit was kept relatively low to avoid doing too much computation in the initialization sequence (in fact it was written on the basis that it might hit the maximum iterations regularly).

As such, I think instead of raising the warning if the iteration limit is hit that it would be better to check that Tbub an Tdew make sense and raise a warning/exception if you get Tbub > Tdew with a suggestion to change the iteration limit. Bumping the iteration limit a bit might not be bad either.

However, if option c) is reliable, then that sounds like a good way to go.

On a side note however, I have some ideas for an alternative formulation that might manage to avoid some of these issues by removing the need for the explicit calculation of the bubble and dew points (using a pseudo bubble and dew point instead).

[alma-walmsley]

Cool as, thank you for your comments. These make sense to me. A warning for Tbub > Tdew is a much better approach than just seeing if max iterations is reached. The question then is what to do going forward:

• It would be quick and easy to add a warning for Tbub > Tdew if max iterations is reached.
• If you want, it may be worth looking into that alternative formulation. This is SmoothVLE v1 though, and I'm not sure if v2 is the preferred option now...?
• I am not sure how reliable solving in ipopt might be, though I imagine ipopt can easily handle a problem like this...?

[andrewlee94]

The new complementarity formulation is only applicable for systems that use a cubic EoS for both phases, so it has some limitations. Thus, we do need to keep the older version around (at least for now).

As for solving the problem with IPOPT, if I recall the initialization is solved for ideal conditions so the problem should be fairly straight forward.

[dallan-keylogic]

Do you have any more information about the specific case in which this method fails?

IPOPT is probably performing better because it's actually using a line search to determine how big a step is necessary instead of the extremely heuristic step limiter

        # Limit temperature step to avoid excessive overshoot
        if f / df < -50:
            Tbub1 = Tbub0 + 50
        elif f / df > 50:
            Tbub1 = Tbub0 - 50
        else:
            Tbub1 = Tbub0 - f / df

The problem with IPOPT, though, is that calls to it are expensive (due to writing the .nl file), so we want to limit the number of times we call it. For example, one thing I do in a similar situation is concatenate constraints from an indexed StateBlock along the index and then make one single call. The individual problems are independent, so the only way they interact is through determining line search step length.

Are you using units of Kelvin?

[alma-walmsley]

@dallan-keylogic Here is an example:

from pyomo.environ import ConcreteModel
from idaes.core import FlowsheetBlock
import idaes.models.properties.modular_properties.base.utility as utility
from property_packages.build_package import build_package

m = ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False)
m.fs.properties = build_package("peng-robinson", ["benzene", "toluene"], ["Liq","Vap"])
m.fs.state = m.fs.properties.build_state_block([0], defined_state=True)

m.fs.state[0].flow_mol.fix(1)
m.fs.state[0].pressure.fix(500000)
m.fs.state[0].temperature.fix(300)
m.fs.state[0].mole_frac_comp["benzene"].fix(0.5)
m.fs.state[0].mole_frac_comp["toluene"].fix(0.5)

utility.MAX_ITER = 30  # default is 30

m.fs.state.initialize()

It uses the build_package function from the Ahuora property packages repo, which essentially pulls from a compound database to dynamically generate the configuration for a generic property package.

Again, I am not sure if this is indicative of a modelling error, but I have had no issues since increasing MAX_ITER to 1000