idaes-pse icon indicating copy to clipboard operation
idaes-pse copied to clipboard

Remove instances of `log(x)` where x can possibly be zero.

Open andrewlee94 opened this issue 3 years ago • 4 comments
trafficstars

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).

andrewlee94 avatar Oct 05 '22 15:10 andrewlee94

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 avatar Oct 06 '22 19:10 dallan-keylogic

@dallan-keylogic

  1. Yes, this issue is happening during the NL write so this is before IPOPT sees it.
  2. 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.

andrewlee94 avatar Oct 06 '22 19:10 andrewlee94

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 avatar Oct 06 '22 19:10 dallan-keylogic

@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.

andrewlee94 avatar Oct 06 '22 20:10 andrewlee94