modelica: compute minimum/maximum of continuous variable over time

Question

As stated above: I wish to compute the minimum (and/or maximum) of a continuous variable over time. Here is a minimal example to demonstrate:

model MinMaxTest
  Real u;
  Real u_min(start = 10);
  Real u_max(start = -10);
equation
  u = sin(time / 180 * Modelica.Constants.pi);
  u_min = min(u, u_min);
  u_max = max(u, u_max);
  annotation(experiment(StartTime = 0, StopTime = 360, Tolerance = 1e-06, Interval = 1));
end MinMaxTest;

u is the arbitrary continuous variable (for demo purposes a simple sinus wave). u_min/u_max is the minimum/maximum over time.

Obviously the expected result is u_min=-1 and u_max=1. Unfortunately the simulation crashes with a "Matrix singular!" error. Can anyone direct me how to avoid that?

---

EDIT 1

I'm using OpenModelica 1.15 (was 1.9.2)

EDIT 2

As I'm quite new to Modelica, I'm struggling to understand the differences between the following approaches:

1. u_min = if noEvent(u < u_min) then u else pre(u_min);
2. if noEvent(u < u_min) then u_min = u; else u_min = pre(u_min); end if;
3. u_min = if noEvent(u < u_min) then u else u_min;
4. u_min = if u < u_min then u else pre(u_min);
5. u_min = if u < u_min then u else u_min;
6. when u < u_min then u_min = u; end when;
7. u_min + T*der(u_min) = if u <= u_min then u else u_min;

1 and 2 are equivalent and result in the expected behavior.

3 produces the desired result but gives a "Translation Notification" about an "algebraic loop", why?

4 fails in so far, that the resulting u_min curve is identical to u?! why?

5 combines 3 and 4.

6 fails to compile with Sorry - Support for Discrete Equation Systems is not yet implemented

7 I'm unclear what the idea behind this is, but it works if T is of the suggested size.

If I'm understanding the Modelica documentation correctly then 1-5 have in common that exactly one equation is active at all times. noEvent suppresses event generation at the specified zero crossing. I had the impression that this is mostly an efficiency improvement. Why does leaving it out cause 4 to fail? pre refers to the previous value of the variable, so I guess that makes sense if we want to keep a variable constant, but why does 7 work without it? My understanding of when was, that its equation is only active at that precise event, and otherwise keeps the previous value, which is why I tried using it in 6. It seems to work if I compare against constant values (which is of no use for this particular problem).

EDIT3

8. u_min = smooth(0, if u < u_min then u else pre(u_min));

Interestingly, this works also.

Answer 1

For your initial question, what seems to work correctly for me in OpenModelica is this:

u_min = min(u, pre(u_min));
u_max = max(u, pre(u_max));

For me that compiles, simulates, and gives the expected results, but also does say "Matrix singular!". On the other hand, if I change the initial declaration for u_max to this:

Real u_max(start = 0);

Then, the "Matrix singular!" goes away.
I don't know why, but that does seems to do the job, and I would suggest is more straightforward then the other options you have listed.

Answer 2

I tested your model with Dymola 2016 and it works, however you can try to use an alternative approach. In Modelica you have to think in terms of equations and not in terms of assignments.

u_min = min(u, u_min);

Is what you would do if the code were to be executed as a sequence of instructions. Under the hood the Modelica tool converts this equation into a nonlinear system that is solved as the simulation proceed.

These are the statistics I get when simulating your model

Statistics

Original Model
Number of components: 1
Variables: 3
Unknowns: 3 (3 scalars)
Equations: 3
Nontrivial: 3

Translated Model
Time-varying variables: 3 scalars
Number of mixed real/discrete systems of equations: 0
Sizes of linear systems of equations: { }
Sizes after manipulation of the linear systems: { }
Sizes of nonlinear systems of equations: {1, 1}
Sizes after manipulation of the nonlinear systems: {1, 1}
Number of numerical Jacobians: 0

As you can see there are two nonlinear systems, one for u_min and one for u_max.

An alternative solution to you problem is the following

model Test
  Real x;
  Real y;
  Real u_min;
  Real u_max;
  parameter Real T = 1e-4;
equation 
  x = sin(time) + 0.1*time;
  y = sin(time) - 0.1*time;
  u_min + T*der(u_min) = if y <= u_min then y else u_min;
  u_max + T*der(u_max) = if x >= u_max then x else u_max;

end Test;

In this case u_min and u_max are two state variables and they follow the variables x and y, depending on their values. For example, when x is lower than u_max then u_max gets "stuck" to the maximum value reached up to that point in time.

Sorry but I can't post an image of the model running since this is my first reply.

Answer 3

The main issue here is that you get an equation which is singular, since you try to solve the equation u_min = min(u,u_min). Where u_min depends on u and u_min and every value of u_min that is smaller than u would fit in that equation, also a tool might try to use a non-linear solver for that. An other solution for that could be perhaps the delay operator:

  u_min = min(u, delay(u_min,0));
  u_max = max(u, delay(u_max,0));

Some notes on the different approaches:

1. u_min = if noEvent(u < u_min) then u else pre(u_min);
2. if noEvent(u < u_min) then u_min = u; else u_min = pre(u_min); end if;

These both are semantically identical, so the result is should be the same. Also the usage of the pre operator solves the issue, since here u_min depends on u and pre(u_min), so there is no need for a non-linear solver.

3. u_min = if noEvent(u < u_min) then u else u_min;

Like above where min() is used here the solution of u_min depends on u and u_min, what leads to a non-linear solution.

4. u_min = if u < u_min then u else pre(u_min);

The semantic of the noEvent() operator results into literally usage the if-expression, in this case here an event u < u_min is triggered and the expression u_min = u is used all the time.

5. u_min = if u < u_min then u else u_min;

Yes, it combines the problems of 3 and 4.

6. when u < u_min the u_min = u; end when;

Here again the solution of u_min depends on u_min and u.

7. u_min + T*der(u_min) = if u <= u_min then u else u_min;

Here u_min is a state and so the calculation of u_min is done by the integrator and this equation is now solved for der(u_min), which then effects u_min.