Opened 8 years ago
Last modified 4 years ago
#4293 assigned defect
Bad strategy to solve rational equations — at Initial Version
| Reported by: | Francesco Casella | Owned by: | Lennart Ochel |
|---|---|---|---|
| Priority: | blocker | Milestone: | 2.0.0 |
| Component: | Backend | Version: | |
| Keywords: | Cc: | Andreas Heuermann, Vitalij Ruge |
Description
Please consider the models contained in the attached .mo file. They contain 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.
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?
Test package