Now that we have already introduced both time events and state events, let’s examine some important complications associated with state events. Surprisingly, these complications can be introduced by even simple models.
Consider the following almost trivial model:
model Decay1
Real x;
initial equation
x = 1;
equation
der(x) = -sqrt(x);
end Decay1;
If we attempt to simulate this model for 5 seconds, we find that the simulation terminates after about 2 seconds with the following trajectory:
Again, numerical issues creep in. Even though mathematically it
should not be possible for the value of x
to drop below zero,
using numerical integration techniques it is possible for small
amounts of error to creep in and drive x
below zero. When that
happens, the sqrt(x)
expression generates a floating point
exception and the simulation terminates.
To prevent this, we might introduce an if
expression to guard
against evaluating the square root of a negative number, like this:
model Decay2
Real x;
initial equation
x = 1;
equation
der(x) = if x>=0 then -sqrt(x) else 0;
end Decay2;
Simulating this model we get the following trajectory [1]:
Again, the simulation fails. But why? It fails for the same reason, a numerical exception that results from taking the square root of a negative number.
Most people are quite puzzled when they see an error message about a
floating point exception like this (or, for example division by zero)
after they have introduced a guard expression as we have done. They
naturally assume that there is no way that sqrt(x)
can be
evaluated if x
is less than zero. But this assumption is incorrect.
Given the if
expression:
der(x) = if x>=0 then sqrt(x) else 0;
it is entirely possible that sqrt
will be called with a negative
argument. The reason is related to the fact that this is a state
event. Remember, the time at which a time event will occur is known
in advance. But this is not the case for a state event. In order to
determine when the event will occur, we have to search the solution
trajectory to see when the condition (e.g., x>=0
becomes
false).
The important thing to understand is that until the event occurs,
the behavior doesn’t change. In other words, the two sides of this
if
expression represent two types of behavior, der(x)=sqrt(x)
and der(x)=0
. Since x
is initially greater than zero, the
initial behavior is der(x)=sqrt(x)
. The solver will continue
using this equation until it has determined the time of the event
represented by x>=0
. In order to determine the time of that
event, it must go past the point where the value of the conditional
expression changes. This means that while attempting to determine
precisely when the condition x>=0
changes from true to false, it
will continue to use the equation der(x)=sqrt(x)
even though x
is negative.
Most users initially assume that each time der(x)
is evaluated,
the if
expression is evaluated (specifically the conditional
expression in the if
expression). Hopefully the previous
paragraph has made it clear that this is not the case.
This time spent trying and retrying integration steps can be saved
thanks to the fact that Modelica can extract a so-called “zero
crossing” function from the if
expression. This function is
called a zero crossing function because it is normally constructed to
have a root at the point where the event will occur. For example, if
we had the following if
expression:
y = if a>b then 1 else 0;
The zero crossing function would be . This function is
chosen because it changes from positive to negative precisely at the
point where a>b
.
Recall our previous equation:
der(x) = if x>=0 then sqrt(x) else 0;
In this case, the zero crossing function is simply since the event occurs when itself crosses zero.
The Modelica compiler collects all the zero crossing functions in the model for the integrator to use. During integration, the integrator checks to see if any of the zero crossing functions have changed sign. If they have, it uses the solution it computed during that step to interpolate the zero crossing function to find the location, in time, of the root of the zero crossing function and this is the point in time where the event occurs. This process is much more efficient because the root finding algorithms have more information to help them identify to location of the root (information like the derivative of the zero crossing function) and evaluation is very cheap because it doesn’t involve taking additional integration steps, only evaluating the interpolation functions from the triggering integration step.
But after all this, it still isn’t clear how to avoid the problems we
saw in the Decay1
and Decay2
models. The answer is a special
operator called noEvent
. The noEvent
operator suppresses this
special event handling. Instead, it does what most users expected
would happen in the first place, which is to evaluate the conditional
expression for every value of x
. We can see the noEvent
operator in action in the following model:
model Decay3
Real x;
initial equation
x = 1;
equation
der(x) = if noEvent(x>=0) then -sqrt(x) else 0;
end Decay3;
and the results can be seen here:
Now the simulation completes without any problem. This is because the
use of noEvent
ensures that sqrt(x)
is never called with a
negative value of x
.
It might seems strange that we have to explicitly include the
noEvent
operator in order to get what we consider the most
intuitive behavior. Why not make the default behavior the most
intuitive one? The answer is performance. Using conditional
expressions to generate events improves the performance of the
simulations by giving the solver clues about when to expect abrupt
changes in behavior. Most of the time, this approach doesn’t cause
any problem. The examples we have presented in this chapter were
designed to highlight this issue, but they are not representative of
most cases. For this reason, noEvent
is not the default, but must
be used explicitly. It should be noted that the noEvent
operator
should only be used when there is a smooth transition in behavior,
otherwise it can create performance issues.
There is a common effect known as “chattering” that you will run into sooner or later with Modelica. Consider the following model:
model WithChatter
Real x;
initial equation
x = 2;
equation
der(x) = if x>=1 then -1 else 1;
end WithChatter;
Effectively, the behavior of this model is that for any initial value
of x
, it will progress linearly toward 1. Mathematically
speaking, once the value of x
gets to 1, it should just stay
there. This is because any deviation away from 1, either greater than
or less than, will immediately cause it to go back to 1.
But we will not be solving these equations in a strictly mathematical
way. We’ll be using floating point representations and using
numerical integrators. As such, we have limited precision and
integration error to content with. The net effect will be that the
trajectory of x
will not remain exactly 1 but will deviate
slightly above and below. Each time this happens, it will generate an
event.
Simulating this model gives us the following results:
This kind of model can introduce an effect known as “chattering”.
Chattering is simply the degradation in simulation performance due to
a large number of events occurring that artificially shorten the time
steps taken by the solver. The impact on simulation performance is
clear if we look at the CPU time taken during the simulation. It
starts to rise dramatically once x
is close to 1. This is because
behind the scenes the events are causing lots of very small time steps
which dramatically increases the number of computations being
performed. The important thing about the WithChatter
example is
that it has a seemingly obvious mathematical solutions but still
suffers from degraded simulation performance because of the high
frequency of events.
This is another case where the noEvent
operator can help us out.
We can suppress the events being generated by the conditional
expression by using the noEvent
operator as follows:
model WithoutChatter
Real x;
initial equation
x = 2;
equation
der(x) = if noEvent(x>=1) then -1 else 1;
end WithoutChatter;
In doing so, we will get approximately the same solution, but with better simulation performance:
Note how there is no perceptible change in the slope of the CPU usage
before and after x
arrives at 1. Contrast this with the
significant difference seen in the case of WithChatter
.
In reality, equations like this are uncommon. In this case, we’ve used an extreme case in an attempt to clearly show the impact of chattering. The behavior being described here is not particularly realistic or physical. In this case, we’ve exaggerated the effect to clearly demonstrate the impact on simulation performance.
More typical examples of chattering in the real world will feature a
conditional expression that sits at some stable point (Decay2
is a
good example). In such cases, chattering occurs because the system
tends to naturally settle at or near the point where the conditional
expression occurs. But because of precision and numerical
considerations, the event associated with the conditional expression
is frequently triggered. The effect is exacerbated in cases where
there are many components in the system all sitting at or near such
equilibrium points.
Hopefully the discussion so far has made it clear why it is necessary to suppress events in some cases. But one might reasonably ask, why not skip events and just evaluate conditional expressions all the time? So let’s take some time to explore this question and explain why, on the whole, associating events with conditional expressions is very good idea [2].
Without event detection, the integrator will simply step right over events. When this happens, the integrator will miss important changes in behavior and this will have a significant impact on the accuracy of the simulation. This is because the accuracy of most integration routines is based on assumptions about the continuity of the underlying function and its derivatives. If those assumptions are violated, we need to let the integration routines know so they can account these changes in behavior.
This is where events come in. They force the integration to stop at the point where a behavior change occurs and then restart again after the behavior change has occurred. The result is greater accuracy but at the cost of slower simulations. Let’s look at a concrete example. Consider the following simple Modelica model:
model WithEvents "Integrate with events"
parameter Real freq = 1.0;
Real x(start=0);
Real y = time;
equation
der(x) = if sin(2*Modelica.Constants.pi*freq*time)>0 then 2.0 else 0.0;
end WithEvents;
Looking at this system, we can see that half the time the derivative
of x
will be 2
and the other half of the time the derivative
of x
will be 0
. So over each of these cycles, the average
derivative of x
should be 1
. This means at the end of each
cycle, x
and y
should be equal.
If we simulate the WithEvents
model, we get the following results:
Note how, at the end of each cycle, the trajectories of x
and
y
meet. This is a visual indication of the accuracy of the
underlying integration. Even if we increase the frequency of the
underlying cycle, we see that this property holds true:
However, now let us consider the case where we use exactly the same
integration parameters but suppress events by using the noEvents
operator as follows:
model WithNoEvents "Integrate without events"
parameter Real freq = 1.0;
Real x(start=0);
Real y = time;
equation
der(x) = noEvent(if sin(2*Modelica.Constants.pi*freq*time)>0 then 2.0 else 0.0);
end WithNoEvents;
In this case, the integrator is blind to the changes in behavior. It
does its best to integrate accurately but without explicit knowledge of
where the behavior changes occur, it will blindly continue using the
wrong value of the derivative and extrapolate well beyond the change
in behavior. If we simulate the WithNoEvents
model, using the
same integrator settings, we can see how significantly different our
results will be:
Note how quickly the integrator introduces some pretty significant error.
The integration settings used in these examples were chosen to
demonstrate the impact that the noEvent
operator can have on
accuracy. However, the settings were admittedly chosen to accentuate
these differences. Using more typical settings, the differences in
the results probably would not have been so dramatic. Furthermore,
the impact of using noEvent
are impossible to predict or quantify
since they will vary significantly from one solver to another. But
the underlying point is clear, using the noEvent
operator can have
a significant impact on accuracy simulation results.
[1] | This model will not always fail. The failure depends on how much integration error is introduced and this, in turn, depends on the numerical tolerances used. |
[2] | A special thanks to Lionel Belmon for challenging my original discussion and identifying several unsubstantiated assumptions on my part. As a result, this explanation is much better and includes results to support the conclusions drawn. |