Troubleshooting stability of PFR model with enthalpy balances enabled and reactions proceeding to completion

[jasonmbray-p66]

System specs:
python 3.11.5
pyomo 6.7.1
idaes-pse 2.2.0

I am modeling a complex reaction network (>100 reactions) using the PFR reactor model and starting to have problems with the solver not finding feasible solutions and/or taking a very long time to solve. I initially had it working isothermally (energy_balance_type = none) in what seemed to be a stable manner and reasonable initialization and solve times of a few seconds. Then I have recently turned on the enthalpy balance (energy_balance_type = enthalpyTotal, default for FTPx state definition) to incorporate heats of reaction and now I am having serious problems getting this to initialize or solve.

After a lot of troubleshooting, I've narrowed things down somewhat to where I suspect the infeasibility issues are related to some species being reacted to extinction. Some of the reactions happen very quickly, and as a result the reaction rates are attempting to drive the mole fraction of those reactants below the 1e-20 lower bound for mole_frac_comp variables. I think this is at least part of the problem because trying some workarounds like removing the reaction from the model or changing the lower bound on the mole_frac_comp variable have allowed the model to solve in some cases.

Question 1: Is there a recommended or preferred way in IDAES to handle reactions that run to completion, that is, one of the reactants is completely consumed?

Solve time and reliability were not as much of a problem until I introduced the energy balance to the reactor. For example, it was not having any issues solving, even with the very small concentrations of reactant, for the isothermal case with energy_balance_type = none. What is more, even if I deactivate the enthalpy_balances constraint and fix the reactor temperatures to be exactly the same as the isothermal case, there is still a 10-fold slow-down in solve times, if the reactor solves at all.

Question 2: What happens when the enthalpy balance equations are added to make the model so much slower and less stable in solving what should otherwise be the exact same problem?

[adowling2]

Hello,

I recommend using the diagnostics toolbox to confirm the numerics of your model are reasonable. You may find there is a modeling mistake or the equations are poorly scaled.

https://idaes-examples.readthedocs.io/en/2.4.0/docs/diagnostics/diagnostics_toolbox_doc.html

Cheers,
Alex

[dallan-keylogic]

Question 1: Is there a recommended or preferred way in IDAES to handle reactions that run to completion, that is, one of the reactants is completely consumed?

If you have reactions that are completely consumed, you might want to use log-form variables in your kinetic expressions. They've dramatically improved the convergence of the GibbsReactor when a single component is consumed to near-exhaustion.

Also, consider if you have highly-reactive intermediates you can eliminate from your model entirely. This paper is very helpful in that regard. However, if you're in a breakthrough situation in which the component consumed to exhaustion is a function of feed conditions, the QSSA will not help.

Solve time and reliability were not as much of a problem until I introduced the energy balance to the reactor. For example, it was not having any issues solving, even with the very small concentrations of reactant, for the isothermal case with energy_balance_type = none. What is more, even if I deactivate the enthalpy_balances constraint and fix the reactor temperatures to be exactly the same as the isothermal case, there is still a 10-fold slow-down in solve times, if the reactor solves at all.

Question 2: What happens when the enthalpy balance equations are added to make the model so much slower and less stable in solving what should otherwise be the exact same problem?

Which equation of state are you using for your system?

I imagine if you have ~100 reactions, you have a similar number of components. Performing enthalpy balances on such a system would involve creating a bunch of Expressions for the enthalpy in terms of temperature, depending on what correlation you're using. When you used energy_balance_type = enthalpyTotal, the ControlVolume1D generates Vars for the enthalpy and the spatial derivative of the enthalpy as well as discretization equations for the derivative. Even when the enthalpy balance itself is deactivated, these additional constraints still exist and can still cause problems.

Enthalpy is the poster child for a badly scaled quantity because it tends to be extremely large in SI units. If variable and equation scaling are not performed correctly, IPOPT will spend all its efforts reducing the residual in the enthalpy balance at the expense of everything else. Unfortunately, the new scaling toolbox has not been rolled out for the ControlVolume1D yet. You'll either have to use the old scaling tools or manually scale all the balance equations you're using. If you use the new toolbox, remember to perform a Pyomo scaling transformation before solving the model. The diagnostics toolbox that @adowling2 suggested can help you identify badly-scaled variables and constraints.

For more about scaling, check out this paper.

[jasonmbray-p66]

Thanks to both of you for your quick responses and insights. It sounds like short of some undiscovered error, scaling is the most likely culprit. I'll provide a little more detail along those lines and see if you have any further suggestions.

I recommend using the diagnostics toolbox to confirm the numerics of your model are reasonable. You may find there is a modeling mistake or the equations are poorly scaled.

I have been using the diagnostic tools, so I know there are no structural issues in the model. The numerical issues list is lengthy, but almost all related to the VLE aspect of the system. I spent a lot of time looking at scaling in the isothermal case a few months ago when I was having my first round of solver troubles. In the end, I found a different error that was causing the problems, so it's hard to tell what, if any, impact the scaling had on getting the isothermal case to work, although I was able to knock a few variables off the list at that time. The condition number did not improve at all, really, and is still much too large... but again, this was already the case for my working isothermal reactor model.

I would estimate 99% of the variables on the large Jacobian list involve mole_frac_phase_comp of species that can be in both vapor and liquid, and these correspond to StateBlock equilibrium_constraint constraints with large Jacobian values throughout the model, not just in the reactor. I think it comes down to the fact that pressures, including vapor pressures, are extremely large in the base units, but since vapor pressures are expressions, not variables, I haven't found a good way to scale them. And many of the mole_frac_phase_comp's are also quite small, as most species strongly prefer one phase or the other, with H2 dominating the vapor phase and making fractions of everything else pretty small and H2 in liquid also being very small as you might expect. All mole_frac_phase_comp's are already scaled throughout the model, though I could explore picking out the specific ones in the list to see if I can modify scaling just for those phase-component pairs in question. I'm fairly new to Pyomo and model scaling, so I'm certainly open to suggestions here as it can only help.

Comparing the isothermal reactor numerical issues list with that of the reactor with enthalpy balances, I have the same long list of VLE-related variables and constraints, but with the control volume's enthalpy_flow_linking_constraint added near the top of the list. This constraint seems to also correspond to mole_frac_phase_comp of literally every species in the system showing up in the large Jacobian variables list. This constraint is being auto-scaled by whatever the 1d control volume's calculate_scaling_factors method is doing. Right now it looks like it's maybe only getting scaled by 0.01, so maybe not enough to move the needle. Would you say that built-in scaling method is inadequate and I need to give this constraint special treatment?

If you have reactions that are completely consumed, you might want to use log-form variables in your kinetic expressions. They've dramatically improved the convergence of the GibbsReactor when a single component is consumed to near-exhaustion.

I already explored this a little bit - not re-writing reaction rate expressions in terms of log variables, but just hoping I might be able to relax the 1e-20 lower bound on mole_frac_comp's and instead use the log form of the composition variable to control the true lower limits. Unfortunately, I couldn't even get the feed streams to initialize, so was never able to test how this might help in the reactor. The model seems extremely temperamental, maybe due to the VLE scaling issues above, so I got stuck at this point.

Also, consider if you have highly-reactive intermediates you can eliminate from your model entirely. This paper is very helpful in that regard. However, if you're in a breakthrough situation in which the component consumed to exhaustion is a function of feed conditions, the QSSA will not help.

I think I mentioned I tried removing a few obviously misbehaving reactions from the model with some success. It is challenging to do this in a general way because while some reactions might be able to be eliminated completely, there will always be some reactions that are on the border where it does depend on the feed or reactor conditions. I don't think the issue is related to transient intermediates, rather the main reactive species present in the inlet being consumed completely.

Unfortunately, the new scaling toolbox has not been rolled out for the ControlVolume1D yet. You'll either have to use the old scaling tools or manually scale all the balance equations you're using.

I guess I'm using the old scaling tools as you say.

If you use the new toolbox, remember to perform a Pyomo scaling transformation before solving the model

I did attempt transforming the model post-scaling following the generic Pyomo documentation before learning how to ask ipopt to use the scaling factors internally. It completely messes up initialization in IDAES because there are physical assumptions built into the initialization routines that are violated when certain variables become scaled. And since I can't initialize the model without at least some scaling, using a transformed model does not work. Hopefully you all are taking initialization into account as you re-work the scaling tools.

[dallan-keylogic]

I would estimate 99% of the variables on the large Jacobian list involve mole_frac_phase_comp of species that can be in both vapor and liquid, and these correspond to StateBlock equilibrium_constraint constraints with large Jacobian values throughout the model, not just in the reactor.

Here the devil is in the implementation details. What equation of state and equilibrium form are you using? We've found that equating log_fugacity is better conditioned than equating fugacity directly. If you're using the cubic equation of state, that naturally leads to use of log mole fraction variables, but the ideal equation of state takes the log function of the normal mole fraction variables. We really should update the ideal EOS to use log variables as well, but it never came up in any of our applications (generally if we need to do VLE we're using something more complex than an ideal EOS). There's even more technical debt involving handling of Henry's Law components here.

Comparing the isothermal reactor numerical issues list with that of the reactor with enthalpy balances, I have the same long list of VLE-related variables and constraints, but with the control volume's enthalpy_flow_linking_constraint added near the top of the list. This constraint seems to also correspond to mole_frac_phase_comp of literally every species in the system showing up in the large Jacobian variables list. This constraint is being auto-scaled by whatever the 1d control volume's calculate_scaling_factors method is doing. Right now it looks like it's maybe only getting scaled by 0.01, so maybe not enough to move the needle. Would you say that built-in scaling method is inadequate and I need to give this constraint special treatment?

The problem is probably not with the built-in scaling method, but rather that it's looking for a scaling factor that it assumed was going to be set by the user but wasn't. Here we see it looking for a scaling factor on _enthalpy_flow[t, x, p]. Backing up, that scaling factor is taken from properties[t, x].get_enthalpy_flow_terms(p). Figuring out where that got a scaling factor took a bit of legwork, but it's here in the modular property base file. There, it's getting scaling factors from the phase molar flow rate and phase molar enthalpy. The molar flow rate's default scaling factor is set by default based on what the user specified as the initial molar flow rate. (If this isn't a good scaling factor, the user should specify the default scaling factor elsewhere). There is no default scaling factor for phase molar enthalpy, but get_scaling_factor is specified to warn the user that there is no scaling factor set for enth_mol_phase, so you're probably seeing warnings that no scaling factor has been set.

So there is a well-thought-out procedure for scaling in place, but it's opaque to all but the most experienced users. The new scaling tools are supposed to help with this problem, but I'm not sure that they'll actually end up doing all that much. The response most users have to getting a ton of logged warnings about unset scaling factors is "How do I hide those messages?" not "I should set scaling factors." If you have any suggestions about how to make the process more visible to users, I'd love to hear them.

For more detail about how to use default scaling factors in the modular property framework, look at this example.

I did attempt transforming the model post-scaling following the generic Pyomo documentation before learning how to ask ipopt to use the scaling factors internally. It completely messes up initialization in IDAES because there are physical assumptions built into the initialization routines that are violated when certain variables become scaled. And since I can't initialize the model without at least some scaling, using a transformed model does not work. Hopefully you all are taking initialization into account as you re-work the scaling tools.

Yes, we are aware of the problem with initialization of scaling-transformed models. The Pyomo function calculate_variable_from_constraint and the block triangularization initializer work equally well on the original and transformed model. However, for more complicated initialization procedures like the Rachford-Rice method we use for initializing VLE, handling scaling is much more difficult. The solution we eventually want to apply is to have Pyomo perform the scaling transformation while writing the .nl file to pass to IPOPT (or your solver of choice). While that writer option currently exists, it also has bugs in it causing the scaling factors to be applied incorrectly. We don't have an example simple enough to troubleshoot, though.

Note that even if you have IPOPT use user-set scaling factors, you still can run into trouble because it only uses them for computing and inverting the Jacobian. That means things like the constraint residual and computation of how close variables are to bounds don't benefit from scaling (and how close variables are to bounds is very important to the interior point method implemented in IPOPT).

[jasonmbray-p66]

Thanks a lot for unpacking the details of how the scaling works here. I was starting to figure out some things on my own before I asked this question, but you've explained several issues where I had hit walls before, like how the ideal EOS doesn't support log variables.

Fortunately, after spending time looking closely at each constraint on the extreme jacobian list, I was able to get my model to solve with all the enthalpy constraints active! More below about what I learned from this exercise.

The molar flow rate's default scaling factor is set by default based on what the user specified as the initial molar flow rate. (If this isn't a good scaling factor, the user should specify the default scaling factor elsewhere).

There is a bug I reported a few months ago with how default scaling factors get applied, and the response I got was something like "we aren't going to fix that right now because the whole scaling approach is being overhauled soon". I can appreciate that is a valid reason, so I had to make a hot fix in my own code to be able to define my own default scaling factors and have them take effect.

So there is a well-thought-out procedure for scaling in place, but it's opaque to all but the most experienced users.

I completely agree! It is a challenging topic and not easy to fix without adding cleaner examples to the documentation, and making it part of the discussion from the very beginning, just like initialization. It's just another step in the model building process.

The response most users have to getting a ton of logged warnings about unset scaling factors is "How do I hide those messages?" not "I should set scaling factors." If you have any suggestions about how to make the process more visible to users, I'd love to hear them.

My experience is that when I first heard about scaling, I viewed it as something that I could ignore as long as my model is solving ok (I think ipopt has some default scaling behavior that is good enough for some things). Only if I'm having trouble getting feasible solutions should I need to think about scaling, and the next step is to try using the defaults only and see what happens. And only after that doesn't work, do I really need to try and understand it and get it right. Obviously, that's not the right mindset.

Your new Scaling Theory doc page is actually really good, but it's hard to find, and I never found and read it until recently after I had already struggled through figuring things out on my own. Maybe it didn't exist until recently, I don't know. It could be helpful to add a more realistic example from a state block with a few variables and a constraint with realistic numbers instead of only talking in vague generalities. Maybe even demonstrating what I describe below of breaking a constraint down piece by piece. However, I don't see anything in the latest scaling doc pages about where/how to set different default scaling factors or how to manually set them when not using the defaults. Another snippet that could be added, even though it generally seems too simple to even include in documentation, is just to explain that scaling factors are multiplied by the variable or constraint they are associated with. That would have been really helpful for me starting out when trying to decide what scaling factor to choose for a given variable - I had to read the Pyomo source code just to figure out whether variables are multiplied or divided by their scaling factors (i.e. should the scaling factor be set equal to the magnitude of my variable value or 1 / magnitude?).

To figure out the scaling for my model, I ended up having to look at each constraint on the extreme jacobian list one by one (or groups of identical indexed constraints) and print out the scaled values of every variable in the constraint as well as the scaled left hand and right hand sides of the equality. I wrote a couple functions to help me, but it was very tedious. This let me see where my variables needed to be scaled differently and then determine what kind of scaling the constraint needed additionally. I found that there was not one set of default scaling factors I could use for all state blocks, because some of my streams had extremely (10+ orders of magnitude) different component concentrations for which I had to set scaling manually for the special cases.

Something else I had trouble with relating to scaling is getting the DiagnosticsToolbox to run the numerical checks. The jacobian calcs depend on the PyNumero extension in Pyomo, but this was not installed by default. The first time I used the diagnostics I got to trying to run the checks and got an error about missing the PyNumero package, but I struck out over and over trying to find a Python package named PyNumero to add to my conda environment. Eventually I found the Pyomo documentation explaining how to install it. But this dependency should not be a secret. That should be written all over the IDAES installation pages and especially the DiagnosticsToolbox doc pages that explain the diagnostics workflow, with links to the PyNumero installation instructions.

The biggest breakthrough, though, was actually discovering that the DiagnosticsToolbox has a major bug in all of the display_*_with_extreme_jacobians functions. These functions are (incorrectly) written to always calculate the Jacobian based on the unscaled model, so no matter what scaling factors I introduced in the model, every time I reviewed the diagnostics I got the exact same list of problem variables and constraints. It was incredibly demoralizing to try to understand if I was doing scaling correctly and see absolutely no change in the condition number or extreme values lists. This bug literally cost me months of lost productivity, and once I discovered and fixed the issue with my own patch, it was a matter of a few days for me to have my model scaling problems mostly fixed and feeling like I actually knew what was going on.

I've prepared a branch with some proposed changes to the DiagnosticsToolbox class and will submit a PR to hopefully incorporate these changes in a future release. I hope you'll take a look.