Bad strategy to solve rational equations

[Francesco Casella]

Please consider the package contained in the attached .mo file. It contains a model with an implicit rational equation to compute mu given p and T

T =   r_mu * (p *1e-5) / (mu * 1e7 - c_mu)
    - a_mu / (mu * 1e7 * (mu * 1e7 + b_mu)
    + b_mu * (mu * 1e7 - b_mu));

The equation is equivalent to a polynomial equation, which is obtained by multiplying both sides of the equation by the product of all the denominators. The polynomial equation has the same roots, but its residual is more well-behaved, as it doesn't have singular points (vertical asymptotes). So, it might be convenient to turn the original equation into a cubic one for improved convergence robustness.

However, OMC does it by fully expanding the polynomials, thus obtaining this very complex expression that I got from the info.json file

-2e+016 * p * r_mu * b_mu * 10000000.0 * mu ^ 4.0
+ c_mu * p * r_mu * b_mu * -2000000000.0 * c_mu * mu ^ 2.0
+ 10000000.0 * p * r_mu * -1e+016 * (10000000.0 * mu - c_mu) * mu ^ 4.0
+ 100000000000000.0 * T * b_mu * (10000000.0 * mu - c_mu) ^ 3.0 * mu ^ 2.0
+ T * 1e+021 * (10000000.0 * mu - c_mu) ^ 3.0 * mu ^ 2.0 * mu
+ 100000000000000.0 * b_mu * T * 10000000.0 * ((10000000.0 * mu - c_mu) * mu) ^ 2.0 * mu
+ (10000000.0 * mu - c_mu) * ((-c_mu) * (a_mu + (-b_mu) * b_mu * T) - ((-a_mu) - (-b_mu) * b_mu * T) * 10000000.0 * mu) * (10000000.0 * mu - c_mu) ^ 2.0
+ b_mu * T * c_mu * -100000000000000.0 * ((10000000.0 * mu - c_mu) * mu) ^ 2.0
+ p * r_mu * b_mu ^ 2.0 * 1e-005 * (10000000.0 * mu - c_mu) ^ 3.0
+ c_mu ^ 2.0 * -200.0 * p * r_mu * b_mu * (10000000.0 * mu - c_mu) * mu 
+ c_mu ^ 2.0 * -1000000000.0 * p * r_mu * (10000000.0 * mu - c_mu) * mu ^ 2.0
+ (-c_mu) * -10000000.0 * b_mu * T * c_mu * 10000000.0 * (10000000.0 * mu - c_mu) * mu ^ 2.0
+ -100000000000000.0 * c_mu * 10000000.0 * b_mu * T * (10000000.0 * mu - c_mu) * mu ^ 2.0 * mu
+ 1e+016 * r_mu * p * c_mu * mu ^ 3.0 * (10000000.0 * mu - c_mu)
+ (-c_mu) * 10000000.0 * T * b_mu * (10000000.0 * mu - c_mu) * mu * (10000000.0 * mu - c_mu) ^ 2.0 
+ (-c_mu) * 100000000000000.0 * T * (10000000.0 * mu - c_mu) ^ 3.0 * mu ^ 2.0
+ (-c_mu) * 10000000.0 * b_mu * T * -10000000.0 * c_mu * (10000000.0 * mu - c_mu) * mu ^ 2.0
+ c_mu * -10000000.0 * b_mu * T * c_mu ^ 2.0 * (10000000.0 * mu - c_mu) * mu
+ 2000000000.0 * c_mu * p * r_mu * b_mu * (10000000.0 * mu - c_mu) * mu ^ 2.0
+ 1e+016 * c_mu * p * r_mu * (10000000.0 * mu - c_mu) * mu ^ 3.0 
+ mu ^ 3.0 * c_mu * 4.000000000000001e+016 * p * r_mu * b_mu

This expression is about 10 times more time consuming to compute. Furthermore it is numerically ill-conditioned, since it involves the summation of very large terms that may be both positive and negative, which is well known as a numerical no-no. Polynomials should be computed using Horner's rule, not by expanding them to individual power terms.

As a consequence, the nonlinear solver behaves erratically, while the original equation has a rather well-behaved shape in a large neighbourhood of the solution. The attached test.mos file provides two test cases demonstrating the problem.

In the first case, the default solver is used. The log file shows that the solver fumbles around for a while, then eventually accepts the start value as the solution (!!), probably because of really bad scaling of the residual.

This is really unacceptable, because no error is issued, though the solution is completely bogus. It is OK if we cannot solve some problems, but claiming they're solved when they aren't is really dangerous.

Please also note that all the involved variables have nominal values for scaling.

If the kinsol solver is used, the simulation fails at initialization, which is at least a consistent behaviour.

Dymola can solve this equation without problems. Apparently, it does not transform it to an equivalent polynomial form at all.

Q1: What is the rationale behind the symbolic transformation carried out by the back-end?

Q2: Is it possible to turn it off and see what happens?

Q3: Why does the default solver accept such a bad solution for the equation?

[Francesco Casella]

Test package

[Francesco Casella]

Test script

[Vitalij Ruge]

Q2: using modelica function with Inline=false
e.g.

  DynamicViscosity mu(start = 10e-5) "Dynamic viscosity";
function foo
    input Real a_mu;
    input Real b_mu;
    input Real c_mu;
    input Real r_mu;
    input DynamicViscosity mu;
    input Pressure p;
    output Temperature T;
   algorithm
      T := r_mu * (p *1e-5) / (mu * 1e7 - c_mu) - a_mu / (mu * 1e7 * (mu * 1e7 + b_mu) + b_mu * (mu * 1e7 - b_mu));
      annotation(Inline=false);
  end foo;
equation
//Coefficients
  T = (B - A) * time + A;
  T1=T-273.14;

  r_mu = rc_mu / (1 + Q1_mu * ((p / p_crit * T / T_crit) ^ 0.5 - 1)) ^ 2;
  c_mu = b_mu * (exp(Q2_mu * (T / T_crit) ^ 0.5 - 1) + Q3_mu * ((p / p_crit) ^ 0.5 - 1) ^ 2);
  T = foo(a_mu, b_mu ,c_mu , r_mu, mu, p);

[Francesco Casella]

Replying to vitalij:

Q2: using modelica function with Inline=false

This is a workaround which I cannot use because the model will then be used in JModelica/Optimica and that does not allow using functions. I really need an equation-based fluid model by design.

In any case, a good OOM tool should not force people to do this kind of things. If symbolic transformations are applied, they should make the solution more robust and more effficient, not less.

[Lennart Ochel]

As far as I can tell from a first quick look, this transformation is performed by the optimization module solveSimpleEquations. You can turn it off an see what happens.

[Francesco Casella]

I tried that. If I turn off solveSimpleEquations at initialization, the transformation is not performed and the initial equation is solved correctly. At that point, even if the module is not turned off for simulation, it seems that the solver latches onto the correct solution ofthe polynomial version of the equation and works fine through the simulation.

Bottom line 1: the optimization of this type of equations carried out by solveSimpleEquations doesn't seem to be a good idea, and should be either fixed or removed. Can you please point me out the rationale behind it?

Bottom line 2: there is something wrong in the acceptance criterion of the default nonlinear solver.

Let me stress once more that having a simulation tool producing a completely wrong result without any warning is the worst possible thing that can happen, so I'd take this very seriously. Not because it's my problem, but because it may happen (or have already happened) to other people that happen to write rational equations in their models.

[Francesco Casella]

Moved to milestone 1.13.0 as milestone 1.12.0 was closed

[Francesco Casella]

Rescheduled to 1.14.0

[Francesco Casella]

I guess this should be addressed, but it's not a blocker for 1.14.0

[Francesco Casella]

@Karim, I now need to handle models with Peng-Robinson's equation (see above) and I'm running again into this issue.

Of course I can disable solveSimpleEquations, but that's rather drastic and it is preventing other useful optimizations.

Could you check the part of the code that turns a rational equation into a polinomial one and just disable that?

[Karim Adbdelhak]

I have a pending PR that fixes the issue for this. Four models from the testsuite are failing, i will add more information later. ​PR6877

[Karim Adbdelhak]

The PR has been pushed and it is resolved. One module had to be changed to not be verified because the results changed (to arguably better results). Copy/Paste from PR:

has to change /simulation/libraries/msl32/Modelica.Fluid.Examples.Explanatory.MomentumBalanceFittings.mos from verified to simple simulation mathematically the initialization has multiple solutions and our results display the second possible solution. Taking the description into account the new solution is the intended one and the old one is not. We cannot change the reference file because it is provided by modelica association and we have to comply. Updating the start value rightBoundary1.ports[1].m_flowleftBoundary1.ports[1].m_flow to something like 40 will produce the inteded solution in dymola. Once that and the result file have been updated by the modelica association we can change it back.

​Relevant MSL PR