Invalid evaluation of auxiliary wrapFunctionCalls variable in conditional equations/expressions

[anonymous]

I was advised by sjoelund.se on Stackoverflow to report an error i had with "acos(variable) outside domain -1.0 <= -1 <= +1.0"

The error occured with OM version 1.13, but is also present with 1.16.

I dont have enough knowledge to determine if it is something to fix/ if it is already handled the appropriate way.

I'd like to add, that since i've been doing a few trigonometry calculation models, that it can be common to run into asin/acos out of domain blockers, despite equations that can mathematically never get out of domain (only numerically, i guess).

In his words "[...]OpenModelica added a bad equation for common subexpression elimination someTmpVar = acos(cosbeta) which always runs regardless of the if-equations you use to guard this[...]", see:

​https://stackoverflow.com/questions/55064752/math-acos-out-of-domain-error-when-combined-with-backlash-element

Detailed description of the model/purpose, see stack overflow link.
For convenience repeated:
Model Diagram: ​https://i.stack.imgur.com/IjOS4.jpg
Model code is:

package test  
  model LD_Abwickel_ZF_VAL  
     Modelica.Blocks.Sources.Ramp ramp2(duration = 1.24, height = 122, offset = -2);
     Modelica.Mechanics.Rotational.Components.ElastoBacklash2 elastoBacklash(b = 1.85005, c = 1e4, d = 1e2, phi_rel0 = -0.925025);
     Modelica.Mechanics.Rotational.Components.Fixed fixed2;
     Modelica.Mechanics.Rotational.Sources.Position position1(exact = true);
     Modelica.Blocks.Math.UnitConversions.From_deg from_deg1;
     test.LD_Pendelarm_ZugfederF2 lD_Pendelarm_ZugfederF21;
  equation
     connect(position1.flange, lD_Pendelarm_ZugfederF21.flange);
     connect(elastoBacklash.flange_a, lD_Pendelarm_ZugfederF21.flange);
     connect(ramp2.y, from_deg1.u);
     connect(from_deg1.y, position1.phi_ref);
     connect(elastoBacklash.flange_b, fixed2.flange);
  end LD_Abwickel_ZF_VAL;


    model LD_Pendelarm_ZugfederF2  
      Modelica.Mechanics.Rotational.Interfaces.Flange_a flange;
      .Modelica.SIunits.Angle phi(displayUnit = "deg");
      .Modelica.SIunits.Torque M;
      parameter .Modelica.SIunits.Length a(displayUnit = "mm") = 25.49e-3;
      parameter .Modelica.SIunits.Length b(displayUnit = "mm") = 23.38e-3;
      parameter .Modelica.SIunits.Length d(displayUnit = "mm") = 43.89e-3;
      parameter .Modelica.SIunits.Length L0(displayUnit = "mm") = 59.5e-3 "Ungespannte Laenge d. Feder";
      parameter .Modelica.SIunits.Length Lk(displayUnit = "mm") = 47.19e-3 "Laenge unbelasteter Federkörper";
      parameter .Modelica.SIunits.Force F0 = 4.1 "innere Vorspannkraft Feder";
      parameter .Modelica.SIunits.TranslationalSpringConstant R(displayUnit = "N/mm") = 0.868e+3 "Federkonstante";
      parameter Integer z = 2 "Anzahl Federn";
      .Modelica.SIunits.Length c(displayUnit = "mm");
      .Modelica.SIunits.Angle delta(displayUnit = "deg", min = -2 * .Modelica.Constants.pi, max = 2 * .Modelica.Constants.pi, nominal = 0);
      .Modelica.SIunits.Angle alpha(displayUnit = "deg", min = -2 * .Modelica.Constants.pi, max = 2 * .Modelica.Constants.pi, nominal = 0.5 * .Modelica.Constants.pi);
      .Modelica.SIunits.Angle beta(displayUnit = "deg", min = -4 * .Modelica.Constants.pi, max = 4 * .Modelica.Constants.pi, nominal = .Modelica.Constants.pi);
      .Modelica.SIunits.Angle psi(displayUnit = "deg", min = -4 * .Modelica.Constants.pi, max = 4 * .Modelica.Constants.pi, nominal = 0);
      .Modelica.SIunits.Length L_gerade(displayUnit = "mm");
      .Modelica.SIunits.Length L_bogen(displayUnit = "mm") = 46.1e-3;
      .Modelica.SIunits.Length L_c(displayUnit = "mm") "rel. Federlaenge (gesamt)";
      .Modelica.SIunits.Force F_c "Summe Federkraft";
      Real cosalpha(min = -1, max = 1, nominal = 0);
      Real cosdelta(min = -1, max = 1, nominal = 0);
      Real cosbeta(min = -1, max = 1, nominal = 0);
    equation
      flange.phi = phi;
      flange.tau = M;
      if noEvent(phi <= 0) then
         c = a - b;
         cosdelta = 1;
         delta = 0;
         L_gerade = sqrt(c ^ 2 + d ^ 2);
         cosalpha = -1;
         alpha = .Modelica.Constants.pi;
      else
         c = sqrt(a ^ 2 + b ^ 2 - 2 * a * b * cos(phi));
         cosdelta = (b ^ 2 - a ^ 2 - c ^ 2) / (-2 * a * c);
         if noEvent(cosdelta <= 0.999999 and cosdelta >= (-0.999999)) then
            cos(delta) = cosdelta;
         else
            delta = if noEvent(sign(cosdelta) > 0) then 0 else .Modelica.Constants.pi;
         end if;
         L_gerade = sqrt(c ^ 2 + d ^ 2 - 2 * c * d * cos(delta + .Modelica.Constants.pi * 0.5));
         cosalpha = (a ^ 2 - b ^ 2 - c ^ 2) / (-2 * b * c);
         if noEvent(cosalpha <= 0.999999 and cosalpha >= (-0.999999)) then
            cos(alpha) = cosalpha;
         else
            alpha = if noEvent(sign(cosalpha) > 0) then 0 else .Modelica.Constants.pi;
         end if;
      end if;
      cosbeta = (d ^ 2 - c ^ 2 - L_gerade ^ 2) / (-2 * c * L_gerade);
      if noEvent(cosbeta <= 0.999999 and cosbeta >= (-0.999999)) then
         cos(beta) = cosbeta;
      else
         beta = if noEvent(sign(cosbeta) > 0) then 0 else .Modelica.Constants.pi;
      end if;
      psi = alpha - beta;
      L_c = L_gerade + L_bogen - (L0 - Lk) - Lk;
      F_c = z * (R * L_c + F0);
      M = F_c * sin(psi);
    end LD_Pendelarm_ZugfederF2;
end test;
{{{

}}}

[Francesco Casella]

The following MWE reproduces the issue in a nutshell

model TestConditionalAcos
  Real x(start = 0, fixed = true);
  Real y;
equation
  der(x) = 1;
  if noEvent(x < 1) then 
    x = cos(y);
  else 
    y = 0;
  end if;
  annotation(experiment(StopTime = 2));
end TestConditionalAcos;

@Karim, @AnHeuermann, I'm not sure how to do this, but we have to figure out how to guard the evaluation of auxiliary variables created by wrapFunctionCalls when the branch they belong to is not active.

[Karim Adbdelhak]

I can try to limit the wrapping of functions to specific branches or to exclude them entirely. Having a full on smart way to do this requires more research. One would have to check conditions for equality and such.

[Francesco Casella]

I would go for a simple solution that addresses this specific case.

[Karim Adbdelhak]

I dug deeper into this and found out that it is not wrapFunctionCalls which messes up here but simply the solving module. It creates those intermediate variables to check for certain things. There is actually so much on wrong here that i will summarize what i found and propose possible changes that we need to discuss.

Our matching module correctly detected to solve the equation

  if noEvent(x < 1) then
    x = cos(y);
  else
    y = 0;
  end if;

for variable y. Now the solver tries to split it up and comes up with the following (i shortened the names for the sake of simplicity because our compiler spits out quite big messy names, furthermore some formatting):

1: $ACOS = acos(x) [Real]
2: $PRE  = pre(y) [Real]
3: $k1   = $_round(($PRE + $ACOS) / 6.283185307179586) [Real]
4: $k2   = $_round(($PRE - $ACOS) / 6.283185307179586) [Real]
5: $x1   = $k1 * 6.283185307179586 - $ACOS [Real]
6: $x2   = $ACOS + $k2 * 6.283185307179586 [Real]
7: y     = if noEvent(x < 1.0) then if noEvent(abs($x1 - $PRE) < abs($x2 - $PRE)) then $x1 else $x2 else 0.0 [Real]

Apart from acos(x) being computed outside the if condition (which has to be fixed) i want to point out the following:

As we can see ($k1) in (3) resolves to be (pre(y) + acos(x))/(2*PI) rounded to the closest integer, and ($k2) in (4) resolves to be (pre(y) - acos(x))/(2*PI) rounded to the closest integer. (both still marked as Real though)

If we now have a look at (5) we see that (x1) would be equal to (pre(y)) if it was not for the integer rounding. If we rounded down this will be slightly less then pre(y) and vice versa.
The inverse is true for (x2) in (6).

So in the end it comes to (y) being solved in (7). First of all, if we do not meet the first condition noEvent(x<1.0) then we just have y=0.0 which is just what we expect. So far so good, but if we do meet that condition we compare the absolute differences of $x1 - pre(y) and $x2 - pre(y) which should just be the integer rounding differences. Then we return the one with less absolute difference.

This inner check is (according to one inline comment in the code) based on this transformation:

  cos(y) = x -> y = acos(x) + 2*pi*k

Maybe i am missing something, but shouldn't k represent the phase and instead of rounding to the closes integer we need something like the rest of the division?

Some intel on this would be great, because i am really not sure if this is correct as is! For our problem i will try to prevent implicit solving of such functions in if branches.

[Philip Hannebohm]

Replying to casella:
Thanks for the pleasantly small MWE!

Replying to Karim.Abdelhak:
(I swapped the quotes for the sake of argument)

This inner check is (according to one inline comment in the code) based on this transformation:

  cos(y) = x -> y = acos(x) + 2*pi*k

Maybe i am missing something, but shouldn't k represent the phase and instead of rounding to the closes integer we need something like the rest of the division?

That k looks like the winding number to keep y continuous instead of wrapping around at 2pi, in which case it must be an integer.

As we can see ($k1) in (3) resolves to be (pre(y) + acos(x))/(2*PI) rounded to the closest integer, and ($k2) in (4) resolves to be (pre(y) - acos(x))/(2*PI) rounded to the closest integer. (both still marked as Real though)

This is probably because the inverse of cosine is inherently ambiguous even in the -pi to pi region, because it only gives values between 0 and pi, so we cannot otherwise distinguish if y is on the positive or negative side. Hence the two choices.

Especially around x = 1 this seems to me like a very unstable system since the derivative of acos(x) is infinite (or 0 w.r.t. y) at that point. Basically like when a steam engine is starting with the piston fully extended, the wheel can't decide whether it should rotate clockwise or anticlockwise...

Some intel on this would be great, because i am really not sure if this is correct as is! For our problem i will try to prevent implicit solving of such functions in if branches.

We could think about better ways to do this, but the principle looks correct to me.
If only there was a way to execute this whole reasoning exclusively inside the if branch...

[Karim Adbdelhak]

Replying to phannebohm:

That k looks like the winding number to keep y continuous instead of wrapping around at 2pi, in which case it must be an integer.

Ah yes totally, thank you! I shouldn't start working at midnight :D

We could think about better ways to do this, but the principle looks correct to me.
If only there was a way to execute this whole reasoning exclusively inside the if branch...

we have to keep it implicit and solve it like that.

[Vitalij Ruge]

Maybe helpful ​page13/21

[Francesco Casella]

Replying to vitalij:

Maybe helpful ​page13/21

Thanks @vitalij!

This presentation is in fact very interesting, I think we should add a section to the openmodelica.org website where we reference all publicly available material like this.

Do you have a ready-made collection of such presentations?

[Francesco Casella]

Replying to phannebohm:

Replying to casella:
Thanks for the pleasantly small MWE!

Some people play Sudoku in their spare time, I prefer to write MWEs for OMC :)

[Andreas Heuermann]

Related ticket (with a much simpler example): #6201

[Philip Hannebohm]

Should be fixed with ​this PR.

[Francesco Casella]

Unfortunately, while the simpler test case of #6201 now runs fine after ​PR 6990, the case in comment:1 fails with

JSON value expected, got: ]]},
{"eqIndex":11,"

Simulation process failed. Exited with code 3.

No idea why this happens.

[Philip Hannebohm]

Replying to casella:

Unfortunately, while the simpler test case of #6201 now runs fine after ​PR 6990, the case in comment:1 fails with

JSON value expected, got: ]]},
{"eqIndex":11,"

Simulation process failed. Exited with code 3.

No idea why this happens.

In the terminal it runs fine. I guess there is something broken in the result file.
Vim suggests to me that line 29 in the .json file has a comma too much in the last ["true",]?

{"eqIndex":10,"section":"initial","tag":"if-equation","display":"if-equation","equation":[["noEvent(x < 1.0)",3,4,5,6,7,8,9],["true",]]},

Possibly I forgot something when making the auxiliary equation. No idea what that is exactly. I'll look into it if I have time.

[Philip Hannebohm]

Yeah so the construct i made to support the conditional evaluation of acos(x) has an else branch without equations. It may sound wrong at first, but those auxiliary variables are only used in one of the two branches and in the other case nothing needs to be computed.

Now, the json writer expects a non empty list. I changed that in ​PR 7187.
I think that was the issue. I don't know how to test it in OMEdit though.

[Francesco Casella]

I can test it tomorrow with the Windows nightly

[Francesco Casella]

Works fine, thanks @phannebohm!