idaes-pse
idaes-pse copied to clipboard
Remove instances of `log(x)` where x can possibly be zero.
The new Pyomo NL writer is far less permissive and will raise exceptions when ever it encounters terms with a log(0), and we have a lot of cases where this is possible in our code. One issue that has been identified in testing the new NL writer is cases of mole fractions appearing in log terms where a mole fraction has been fixed to 0 (e.g. entropy calculations which contain a x*log(x) term).
Due to this (and the fact that it is generally good practice), we should look at reformulating as many of these instances as possible to remove the logs. Additionally, we should be much more careful about fixing mole fractions to 0 (and making sure they are bounded to be strictly positive).
Do we know how this interacts with IPOPT relaxing bounds by 1e-8 by default? I'd think that would happen after the NL is written.
Edit: Thinking about this more, exceptions are probably being thrown when simplifying terms while writing the NL file. Fixed Vars don't even appear as variables in IPOPT, but as numerical constants.
Edit 2: Thinking even more, replacing log(x) with a log_x variable might prevent Pyomo from yelling at you, but the constraint x = exp(log_x) is still going to be degenerate for extremely small x. IPOPTs regularization has a good chance of saving you, but in order to resolve the problem you shouldn't have the component x in that property package to begin with if it's going to be fixed to 0 or 1e-16.
@dallan-keylogic
- Yes, this issue is happening during the NL write so this is before IPOPT sees it.
- Unfortunately that is difficult to achieve whilst still having general property packages; the simplest case being a reactor where the inlet concentration of a component is 0. I do generally encourage people to think like this and to either use multiple property package or set the concentration to a small number, but there are cases where this is not feasible (electrolytes) or where users don;t listen to me.
The reactor case is important. In something like the Gibbs Reactor, things work out just fine, because log_x gets defined through the equations using chemical potential, while x is defined through the material balances.
Where we get real issues, though, are isentropic compressors/expanders. There you end up with a degenerate subsystem because x is still near zero on both sides of the mass balance, and changes to log_x on both sides of the entropy balance cancel each other out. If this starts to trip users up, we could create special entropy terms excluding entropy of mixing in the equations of state, and using those when the phase_change=False for the PressureChanger.
@dallan-keylogic Also, the underlying issue here is actually not numeric per se, but IEE standards. Mathematically, it can be shown the 0*log(0) == 0, but IEE standard state that log(0) == NAN and that anything multiplied by NAN should also be NAN. This is what causes the NL writer to raise exceptions. That said, it probably should raise an exception because this means you've fixed a value to 0 which probably should no be.