Acos in implicit function, failing to solve non-linear system

[AnHeuermann]

For the following model the runtime needs to solve x implicit:

model impAcos
  Real x(start=1);
equation
  time = acos(x);
end impAcos;

but will fail to do so. The equation should be well defined (​print from WolframAlpha), but the solution at time=0 is directly on the edge of the domain where acos is well defined.

The error message:

record SimulationResult
    resultFile = "",
    simulationOptions = "startTime = 0.0, stopTime = 1.0, numberOfIntervals = 500, tolerance = 1e-06, method = 'dassl', fileNamePrefix = 'impAcos', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags = ''",
    messages = "Simulation execution failed for model: impAcos
assert            | warning | The following assertion has been violated during initialization at time 0.000000
assert            | debug   | Model error: Argument of acos(x) outside the domain -1.0 <= 1 <= 1.0
assert            | warning | The following assertion has been violated during initialization at time 0.000000
assert            | debug   | Model error: Argument of acos(x) outside the domain -1.0 <= 1 <= 1.0
assert            | warning | The following assertion has been violated during initialization at time 0.000000
assert            | debug   | Model error: Argument of acos(x) outside the domain -1.0 <= 1 <= 1.0
assert            | warning | The following assertion has been violated during initialization at time 0.000000
assert            | debug   | Model error: Argument of acos(x) outside the domain -1.0 <= 1 <= 1.0
stdout            | warning | While solving non-linear system an assertion failed during initialization.
|                 | |       | | The non-linear solver tries to solve the problem that could take some time.
|                 | |       | | It could help to provide better start-values for the iteration variables.
|                 | |       | | For more information simulate with -lv LOG_NLS_V
assert            | warning | The following assertion has been violated during initialization at time 0.000000
assert            | debug   | Model error: Argument of acos(x) outside the domain -1.0 <= 1 <= 1.0
assert            | warning | The following assertion has been violated during initialization at time 0.000000
assert            | debug   | Model error: Argument of acos(x) outside the domain -1.0 <= 1.01 <= 1.0

[...]

assert            | debug   | Model error: Argument of acos(x) outside the domain -1.0 <= 1 <= 1.0
assert            | warning | The following assertion has been violated during initialization at time 0.000000
assert            | debug   | Model error: Argument of acos(x) outside the domain -1.0 <= 1 <= 1.0
assert            | warning | The following assertion has been violated during initialization at time 0.000000
assert            | debug   | Model error: Argument of acos(x) outside the domain -1.0 <= 1 <= 1.0
assert            | debug   | Solving non-linear system 2 failed at time=0.
|                 | |       | For more information please use -lv LOG_NLS.
assert            | info    | simulation terminated by an assertion at initialization

[AnHeuermann]

Related ticket #6171

[AnHeuermann]

The default solver hybrid can't find the solution because it starts at the solution 1 and want's to check values slightly bigger then 1 which is not a good idea.

Changing the non-linear solver to kinsol gets us stuck in an endless-loop.
Using newton results in a segmentation fault. (That should be fixed as well)
Solver homotopy will give up a bit more early.

So what I need to do is create a new strategy for hybrid and change the iteration direction after some failed tries. Or maybe add a flag to pass to the non-linear solvers to start in a different direction.

[phannebohm]

Sorry, what am I missing?
The equation should be solved symbolically for x to be x = cos(time), right?

[AnHeuermann]

Nope, it should get solved numerically.

[casella]

Replying to AnHeuermann:

Nope, it should get solved numerically.

Why?

In any case, the derivative of acos(x) at the solution is undefined (it approaches infinity), so I doubt any derivative-based numerical solver will ever make it.

I would consider this model as ill conditioned, is there any reason why we want it to succeed?

[AnHeuermann]

I would consider this model as ill conditioned, is there any reason why we want it to succeed?

Yeah, that was my thought as well. If we are very close to invalid range of a function there isn't much we can do about it.

But while I was investigating the segmentation fault for newton I couldn't reproduce the issue today. And then I noticed, that newton produces a correct solution for this example. I think I switched from Windows to Ubuntu while testing. On my Windwos install it still get's a segmentation fault.
The reason for this will probably be, that the direction newton starts with a new guess is opposite to the other solvers.

[casella]

I'm not surprised that the segfault is OS-dependent. Very hard to determine how a Newton solver will behave when trying to compute a solution at a singularity point...

That said, I agree we should get a proper error message, not a hard crash. Regarding this issue, it looks like there is some uncatched invalid operation which wreaks havoc somewhere.

d/dx(acos(x)) = -1/sqrt(1 - x^2)

so when x approaches 1 it goes bad for two reasons: division by zero and negative square root. How exactly is the derivative of acos(x) defined by the backend, and how is it computed in the generated code?

[AnHeuermann]

That is exactly the definition we are using in Differentiate.mo

    // der(arccos(x)) = -der(x)/sqrt(1-x^2)
    case ("acos",_)
      equation
        tp = Expression.typeof(exp);
        (exp_1, funcs) = differentiateExp(exp, inDiffwrtCref, inInputData,inDiffType,inFuncs, maxIter);
        exp_2 = Expression.makePureBuiltinCall("sqrt", {DAE.BINARY(DAE.RCONST(1.0),DAE.SUB(tp),DAE.BINARY(exp,DAE.MUL(tp),exp))}, tp);
      then (DAE.UNARY(DAE.UMINUS(tp),DAE.BINARY(exp_1,DAE.DIV(tp),exp_2)), funcs);

In the C code we call a function called acos(). Not sure where it is defined though, could be from cmath library.

I will try to fix the segmentation fault and leave the rest as it is.

[casella]

Replying to AnHeuermann:

That is exactly the definition we are using in Differentiate.mo

    // der(arccos(x)) = -der(x)/sqrt(1-x^2)
    case ("acos",_)
      equation
        tp = Expression.typeof(exp);
        (exp_1, funcs) = differentiateExp(exp, inDiffwrtCref, inInputData,inDiffType,inFuncs, maxIter);
        exp_2 = Expression.makePureBuiltinCall("sqrt", {DAE.BINARY(DAE.RCONST(1.0),DAE.SUB(tp),DAE.BINARY(exp,DAE.MUL(tp),exp))}, tp);
      then (DAE.UNARY(DAE.UMINUS(tp),DAE.BINARY(exp_1,DAE.DIV(tp),exp_2)), funcs);

The question is, how are those "sqrt" and DAE.DIV operators eventually mapped into C code, and how does this C code handle the case of sqrt(neg) and div(0)?

In the C code we call a function called acos(). Not sure where it is defined though, could be from cmath library.

I guess so. acos(x) is a builtin function in Modelica, and what it means is well known at least from Leonhard Euler's times :)

I will try to fix the segmentation fault and leave the rest as it is.

Sounds good.

[phannebohm]

Replying to casella:

Replying to AnHeuermann:

Nope, it should get solved numerically.

Why?

In any case, the derivative of acos(x) at the solution is undefined (it approaches infinity), so I doubt any derivative-based numerical solver will ever make it.

I would consider this model as ill conditioned, is there any reason why we want it to succeed?

I still don't get why we can't solve this directly. The equation

cos(time) = x

has no problem at time = 0 and there is no ambiguity about the value of x. We can of course always have an assert(time <= pi, "not sure what the exact problem would be").

In general there are more complicated equations like

time = acos(x) + x

for example, which we can not solve directly. Still, in that case the alternative

cos(time - x) = x

would make it less ill conditioned, so why don't we try that? In a NLS or equation that can't be solved explicitly (which is a NLS of size 1, I guess), let's avoid functions with branch points in favor of their inverses.

[phannebohm]

Replying to phannebohm:

In a NLS or equation that can't be solved explicitly (which is a NLS of size 1, I guess), let's avoid functions with branch points in favor of their inverses.

Okay, there are cases where it is not obvious how to do that, for example

acos(x) = asin(x)

Which side shall we prefer in that situation?

EDIT: Oh well, in this case the equation can be simplified to x = sqrt(1 - x^2) and ultimately x = 1/sqrt(2). Bad example, but you get my point.

[casella]

Replying to phannebohm:

I still don't get why we can't solve this directly. The equation

cos(time) = x

has no problem at time = 0

This is true, but the actual problem that you posted in the ticket

time = acos(x)

indeed has a big problem at time = 0. It is easy to argue that based on Hadamard's definition of "​well-posed problem". If you perturb the data (i.e. the left-hand-side) slightly

time - eps = acos(x)

all hell breaks loose, no matter how small eps is. For a problem to be well-posed it is not enough to prove that you have one exact solution on paper.

The reason why this happens is that the derivative of the acos(x) function is singular in x = 0, which also means that the inverse function theorem does not hold there.

I am well aware that the following identity holds

0 = acos(1)

but that doesn't mean that our inverse problem is well-posed. It isn't. For well-posed problems, you can add eps anywhere and nothing bad happens provided that eps is small enough. Not here. Incidentally, eps may represent numerical errors due to finite-precision floating point representation, which are ubiquitous in numerical system simulation. It doesn't make sense to solve ill-posed problems in a simulation software.

The situation is similar to fluid problems, where you'd like to write

m_flow = k*sqrt(dp)

Unfortunately, when dp = 0 this thing is ill-posed and numerically fragile, exactly for the same reason, i.e., the derivative approaches infinity at dp = 0, which is also the boundary of function domain. In fact, people invented the "thermo-fluid square root" specifically to avoid all the bad consequences of that, see ​Modelica.Fluid.Utilities.regRoot.

[casella]

Incidentally, for the same reasons I would consider even this problem to be ill-posed:

model M
  Real x, y;
equation
  y = time;
  y = x^2;
end M;

There are two exact solutions of this problem, x = +/-sqrt(time). But, again, if you subtract eps from time in the first equation, or if you subtract eps from y in the second equation, boom!

Also in this case, you can see there is something fishy here because the Jacobian of the implicit equation y = x^2 is zero when time = 0. Even if you don't solve it numerically, this singularity is a revealing sign of ill-posedness, because it tells you that the function is not locally invertible when time = 0.

[phannebohm]

Replying to casella:

Replying to phannebohm:

I still don't get why we can't solve this directly. The equation

cos(time) = x

has no problem at time = 0

This is true, but the actual problem that you posted in the ticket

time = acos(x)

indeed has a big problem at time = 0.

How can that be? The functional relationship is exactly the same!

It is easy to argue that based on Hadamard's definition of "​well-posed problem". If you perturb the data (i.e. the left-hand-side) slightly

time - eps = acos(x)

all hell breaks loose, no matter how small eps is. For a problem to be well-posed it is not enough to prove that you have one exact solution on paper.

This is only because of the way the equation is handled. If you look at other branches of acos it can easily give negative values for whichever eps you like. The x will just pretty much stay constant.

Replying to casella:

Incidentally, for the same reasons I would consider even this problem to be ill-posed:

model M
  Real x, y;
equation
  y = time;
  y = x^2;
end M;

There are two exact solutions of this problem, x = +/-sqrt(time). But, again, if you subtract eps from time in the first equation, or if you subtract eps from y in the second equation, boom!

Also in this case, you can see there is something fishy here because the Jacobian of the implicit equation y = x^2 is zero when time = 0. Even if you don't solve it numerically, this singularity is a revealing sign of ill-posedness, because it tells you that the function is not locally invertible when time = 0.

I agree 100% with you on that one. First, the solution is not unique and on top of that you start at the singular point, and second, perturbing y gets outside the range of x^2. In other words it's outside the domain of the sqrt function. On top of that the Jacobian is zero so all your argumentation above holds here.

The difference to the other problem

model M
  Real x, y;
equation
  y = time;
  y = acos(x);
end M;

is you don't solve for the ambiguous part of the equation, which in this case is y, but for x and that gives the unique solution x = cos(time).

Here the directionality is the other way around! acos plays the role of ^2, restricting the range of x not y, and cos, which is perfectly fine for any argument, plays the role of sqrt. So your second point does not apply here either, unless you insist that the range of acos is to be taken between 0 and pi. But then the sqrt would need to be non-negative for the same reason.

In the original formulation the Jacobian of the implicit equation y = acos(x) is actually infinite, not zero, as you pointed out earlier. The explicit one has zero Jacobian, but who cares since we can calculate cos(time) directly, no inverse function theorem required.

Maybe I'm not making myself clear. But I believe it's not a complicated issue.

Also I really hope my head isn't too far up in the theoretical clouds right now. Sorry If I'm wrong, but you haven't convinced me yet :D

[sjoelund.se]

If you write 0.0 = acos(x), there should be only 1 answer (x=1). If you write x = cos(0.0) there is only 1 answer as well (x=1).

Now the reason we fail to solve it at time=0.0 is simply that the non-linear solver will never guess that x=1... It prints it as 1 above, but it's actually 1+eps that is guessed. If it actually did try x=1 from the start it should get a 0.0000 residual and stop without trying to improve the answer :) We have had lots of issues for these kinds of NLS problems where there is a min/max/bounds/etc. Ideally all of them should handle this, but it's hard to give a good hint to the solver that this is outside the domain.

Is it ever a problem to invert acos? I think not. So we should just do it even if it hides this problem with NLS solvers.

Inverting cos could be problematic because if you have 0.0 = cos(x) then x has infinite solutions and it could be that the nonlinear system only accepts one of them (2.0 = cos(x) is not a problem because acos(2.0) gives an error). Also see comment:14

[sjoelund.se]

Is it ever a problem to invert acos? I think not.

Actually... thinking about it more... acos(...) only returns 0 .. pi so we should add the assertion in that case as well.

assert(time >= 0 and time <= pi, "...");
x = cos(time);

Even though we could obviously easily solve it for all values of time if inverted.

[casella]

I agree that the case of comment:13 is different, that may have brought some confusion.

But I still claim that inverting acos(x) exactly at the boundary of it domain of definition (i.e. x = 1, y = 0), where the function is not locally invertible because it has no finite derivative, is not a well-posed problem.

If you look at the ​graph of acos(x), it is obvious why that is the case: if you perturb the zero value that time has at the beginning of the simulation and make it slighly negative, say -1e-16, all of a sudden there is no solution at all.

I understand that in this specific case the initial value of time is exactly zero and is not -1e-16. But I still think that we can't expect a tool that ultimately uses double-precision numbers to solve problems that are potentially sensitive to errors that are comparable to machine-epsilon.

For the very same reason, writing if x == 0.0 or if x <> 0.0 is not allowed in Modelica. Even though, technically, you can set x = 0 with an equation and that of course will have no approximation. However, double-precision algebra is simply not good enough to deal with such conditions, in general.

In some cases, you can carry out some parts of the computation symbolically. That's fine. But sooner or later, eventually you'll have to resort to numerical algebra. Hence, I think we shouldn't care much about problems that can only be solved symbolically.

Bottom line, if OMC can't solve this problem, I'm not worried at all. I think we have much more important issues to solve than this.

My 2 cts :)

[phannebohm]

Replying to casella:

If you look at the ​graph of acos(x), it is obvious why that is the case: if you perturb the zero value that time has at the beginning of the simulation and make it slightly negative, say -1e-16, all of a sudden there is no solution at all.

I was aware of the graph, but I thought we are interested in the relation and not the restricted function. Thank you! I get it now. If the modeler wants to have the full relation they can always write it in terms of cos, so writing acos implies the restriction to the real interval [0, pi].
So then let's do as comment:16 says.

In some cases, you can carry out some parts of the computation symbolically. That's fine. But sooner or later, eventually you'll have to resort to numerical algebra. Hence, I think we shouldn't care much about problems that can only be solved symbolically.

In particular I was wondering when to apply a transformation, so as to make the resulting expression numerically easier. But you're right, this is rather esoteric.

In general I believe our symbolic approach has the benefit of suffering from less numerical issues by exploiting symbolic analyses. Sometimes I overthink that :D

[casella]

Replying to phannebohm:

I was aware of the graph, but I thought we are interested in the relation and not the restricted function. Thank you! I get it now. If the modeler wants to have the full relation they can always write it in terms of cos, so writing acos implies the restriction to the real interval [0, pi].

Absolutely. We don't have multi-valued functions in Modelica.

So then let's do as comment:16 says.

Sounds good.

In general I believe our symbolic approach has the benefit of suffering from less numerical issues by exploiting symbolic analyses. Sometimes I overthink that :D

In many cases, yes, but you have to be very careful about the assumptions, particularly the implicit ones :)

[casella]

After ​PR 6990 was merged, the test case of this ticket runs fine.