Set solution domain of an equation

Question

One of the advantages of Modelica is that you can define your model non-causally. For example you can write h=10 - g/2 * t^2 to define the relation between the height of a falling object (starting at 10m above the ground) and time, but you can also let modelica solve time as function of height, without having to rewrite the equation to a sqrt. This works as long as the system has only one possible solution.

The reverse function in this case however has 2 solutions, -t gives the same height as +t. When the object was thrown up in the air and reached its highest point (=10m) at t=0, the negative solution is when the object was thrown up.

The subject is mentioned on modelica university, but no solution is given there. In general it applies to any situation where multiple solutions are possible, e.g. solve alpha in y=sin(alpha) and define -pi/2<alpha<pi/2 or pi/2<alpha<3pi/2. This is relevant because the slopes in both intervals are opposite.

The problem is demonstrated in the model below

model test_range
  Real x,x_rev(start=2),y;
equation 
  der(x)=1000;
  when x>=3*Modelica.Constants.pi/2 then
    reinit(x,Modelica.Constants.pi/2);
  end when;
  y=sin(x);
  y=sin(x_rev);
  annotation (experiment(
      __Dymola_NumberOfIntervals=4000,
      Tolerance=1e-06,
      __Dymola_Algorithm="Dassl"));
end test_range;

Dymola simulation with solution interval jump: it starts in the interval around pi, but jumps to the interval around 2pi.

Simulation of same model in OpenModelica: the solution for x_rev jumps back and forth between the -pi/2--pi/2 interval and the pi/2--3pi/2 interval. (Note that the solution depends on the number of intervals, =100 in this plot)

Is it possible to set a solution domain for variables when there are multiple possible solutions?

Note that Dymola solves some equations symbolically (like inverted sin); and ensures that the new solution is the closest one to the previous solution (so not just a simple arcsin call). In general there is no such guarantee as it depends on the non-linear solver. — Hans Olsson, Aug 17 '22 at 08:05

score 0 · Answer 1 · answered Aug 08 '22 at 07:05

0

There are quite a number of cases in which multiple solutions are possible. One from the electrical domain is described here. Here, the start attribute of the variable is used to "push" the solver in to the right direction. Some more details on the attribute can be found here.

An alternative possibility is to add min and max attributes (also described on the above website). The issue with this is, that those attributes are by default not checked in every simulation step. As an example, the check needs to be activated in Dymola as shown below:

I would try to use the start attribute first. Not fully sure if there are other possibilities though...

answered Aug 08 '22 at 07:05

Markus A.

6,300
10
26

The suggestion to use min and max raises an assertion and possibly stops the simulation with an error rather than trigger the solver to look for an alternative valid solution and continue. Using 'start' is only valid at t=0. My problem is that when the solution comes near a peak during the simulation, it may land in an other interval afterwards. I'll edit the question with a demonstration model. – paul Aug 09 '22 at 11:00
If I'm understanding the issue correctly, something like `Modelica.Math.atan3` (especially the input `y0`) for a sinus function could help, correct? This would include reformulating the model a bit though... – Markus A. Aug 10 '22 at 06:32
It could help somewhat, but as a general solution to this type of problem it is not practical. It would require a specific function for each of these cases. In the case of a simple x^2 function that becomes awkward very quick because x^2 is more intuitive than some generic poly(factors and limits) function. I would like modelica to try to respect min and max when searching for a solution rather than picking the wrong solution and throw an error. If that is (currently) impossible, than the answer is 'impossible' and I should propose this t – paul Aug 11 '22 at 07:54
(sorry, was times out) ....propose this to the standard. – paul Aug 11 '22 at 08:03

score 0 · Answer 2 · answered Aug 17 '22 at 08:46

There is no general solution for this in Modelica. However, in many cases it is possible to do as follows:

Find a solution
Turn it into the unique one

So in this case:

  Real x_rev2=smallestSinInverse(x_rev);
  function smallestSinInverse
    input Real x;
    output Real z;
  algorithm 
    z:=x - floor(x/(Modelica.Constants.pi*2))*Modelica.Constants.pi*2;
    if z>3*Modelica.Constants.pi/2 then
      z:=z-Modelica.Constants.pi*2;
    elseif z>Modelica.Constants.pi/2 then
      z:=Modelica.Constants.pi-z;
    end if;
  end smallestSinInverse;

That also works for the simple t^2 case, but I agree that it is not a general solution - and it is messy to get right.

Alternatively in some cases it is possible to use Modelica.Math.Nonlinear.solveOneNonlinearEquation - but it is slower - and is less equation-based; and requires knowledge of accurate min/max:

  Real x_rev3=Modelica.Math.Nonlinear.solveOneNonlinearEquation(
    function mySin(u2=y), -Modelica.Constants.pi/2, Modelica.Constants.pi/2);
  function mySin
    input Real u;
    input Real u2;
    output Real y;
  algorithm 
    y:=Modelica.Math.sin(u)-u2;
  end mySin;

Set solution domain of an equation

2 Answers2