Asserts in conditional equations are not handled correctly

[Francesco Casella]

Please consider the following MWE

model M3
  Real x;
equation
  if time > 0.5 then
    x = time;
    assert(x > 0, "error");
  else
    x = 0;
  end if;
end M3;

At time = 0.501 the runtime fails with

The following assertion has been violated at time 0.502000
if time > 0.5 then x > 0.0 else true

Obviously there is something wrong with the way asserts are handled, apparently they still use some "old value".

[Andreas Heuermann]

I think the question is at which super-dense time do asserts need to be true?
So does it need to be true during all of event-iteration or only at the end of it?

But the problem here is:
At time=0.502 we evaluate (a bit simplified)

[...]
/* x = if time > 0.5 then time else 0.0 */
RELATIONHYSTERESIS(tmp0, data->localData[0]->timeValue, 0.5, 0, Greater);
data->localData[0]->realVars[0] /* x variable */ = (tmp0?data->localData[0]->timeValue:0.0);
[...]
/* assert(if time > 0.5 then x > 0.0 else true, "error");*/
tmp1 = Greater(data->localData[0]->timeValue,0.5);
if (tmp1) {
  tmp2 = Greater(data->localData[0]->realVars[0] /* x variable */,0.0);
}
if (!tmp2) {
/* Print error */
}

We are still in continuous-time-mode, so the value of x is not allowed to change and we get the pre-value x=0.0.
So the assert checks the value of x and issues an error.

The problem is, that we don't save the actual value of x during continuous simulation mode, but only use the pre value of it.
I'm not sure if can just use that value for asserts. I think this would cause more issues but can't come up with an example for that yet.

But to implement it we can check if the variable in the assert is discrete and if so use a new value that has its actual value and not the one from data->localData[0]->realVars.

Shouldn't be hard to implement.

[Francesco Casella]

Replying to AnHeuermann:

I think the question is at which super-dense time do asserts need to be true?
So does it need to be true during all of event-iteration or only at the end of it?

I would say only at the end, for two reasons. The first is that we never check assertions during iterations, otherwise we may have too many failures, but only on (provisionally) accepted steps, when a solution has finally be found.

The second is that, in practice, we have models with if-equations where one branch corresponds to the "normal" operating state and the other one corresponds to some "parked" or "inactive" state, usually with dummy equations for variables that do not have any meaning. If I need to make some assertions on the "normal" operating state, they should never be called in the "parked" state, not even during the transitions, because they are very likely to fail.

[Francesco Casella]

I'm not sure I'm getting the details of your comment:3. My intuitive understanding (based on Appendix C of the Specification) of event handling when you have if-equations is the following

• During continuous integration, you stick to the currently active branch of the if-equation to compute the next step
• After computing the step, you check the zero-crossing functions.
• If one of them changes, you use the solver to compute the zero-crossing instant, while keeping the same branch active
• Once at the crossing, you process the event, possibly causing event iteration in super-dense time
• When the event iteration reaches a fixed point, you switch the branch of the if-equation according to the condition(s) that actually changed
• You restart integration, keeping the new branch active

Is this consistent with your view?

If so, I find it weird that the assert of the now inactive branch is still computed, once the simulation has restarted on the new branch. This is not supposed to be the case. How you can implement that, that's another story, and I have no strong opinions about it.

[Francesco Casella]

Replying to AnHeuermann:

We are still in continuous-time-mode, so the value of x is not allowed to change and we get the pre-value x=0.0.

Specifically, I don't get this. The assert equation depends on x. So, it should be evaluated after x has been computed. Why is the pre value used? Aren't we considering asserts when reordering equations during causalization?

[Francesco Casella]

BTW, x is not a discrete variable, since it changes continuously during some parts of the simulation.

Discrete variables only change at events and remain constant during continuous-time integration, that's not the case of x in this example.

[Francesco Casella]

Maybe related to #6201

[Philip Hannebohm]

Replying to casella:

Maybe related to #6201

I think that one has another focus as it deals more with design questions an numeric stability.

This reminded me more of #6171 where the focus is on the problem of branching (as in computer science :) ) and evaluating.
The solution might also be similar, see ​this PR.

[Philip Hannebohm]

Replying to phannebohm:

The solution might also be similar, see ​this PR.

Scrap that. On second thought I don't really know what's going on here. Sorry :)

[Andreas Heuermann]

My previous explanation is bad. Let me try again:

We call the solver from 0.5 to 0.502 (continuous time mode):
1.   Evaluate x = if time > 0.5 then time else 0.0
       x := 0 (not in event mode, so using pre value of relation time > 0.5, a.k.a false
2.   Check assert(if time > 0.5 then x > 0.0 else true, "error")
       tmp1 = Greater(data->localData[0]->timeValue,0.5);  <--- this is true
       if tmp1 then
         tmp2 = Greater(data->localData[0]->realVars[0] /* x variable */,0.0);   <--- this is true as well
         if not tmp2 then
           Throw assert

We are using RELATIONHYSTERESIS to update x only when we are in event mode. When we have an assert depending on a relation we should do the same.

Doing that in the code seems to work, so changing that in the code generation could solve this problem. I hope it's easy to find in the templates.

[Francesco Casella]

Replying to AnHeuermann:

My previous explanation is bad. Let me try again:

We are using RELATIONHYSTERESIS to update x only when we are in event mode. When we have an assert depending on a relation we should do the same.

Doing that in the code seems to work, so changing that in the code generation could solve this problem. I hope it's easy to find in the templates.

Agreed. Mabe this would fix other issues that were revealed by the verification process in the library testuite, there are some glitches just after events take place.

[Francesco Casella]

After some discussion with @AnHeuermann, we thought that maybe introducing some auxiliary boolean variables to compute the conditions would solve the problem, because the BLT causalization process would handle them appropriately.

I tried to figure out how to do that, but it turned out not to be easy. If I do it like this:

model M4
  Real x;
  Boolean condition;
equation
  if time > 0.5 then
    x = time;
    //assert(condition, "error");
  else
    x = 0;
  end if;
  condition = x > 0;
  annotation(experiment(StopTime = 1, Interval = 0.1),
             __OpenModelica_simulationFlags(lv = "LOG_EVENTS_V"));
end M4;

I still get the same problem. Unfortunately, I can't move the condition = x > 0 equation inside the first branch of the if-equation, because that would violate the rule that all branches of it have the same number of equations.

I guess that's where the crux of the problem is.

If I temporarily comment out the assert statement from M4, I get this simulation log

Integrator attempt to handle a problem with a called assert.
The following assertion has been violated at time 0.600000
if time > 0.5 then condition else true
error
check for discrete changes at time=0.5
Integrator attempt to handle a problem with a called assert.
The following assertion has been violated at time 0.550000
if time > 0.5 then condition else true
error
model terminate | Simulation terminated by an assert at time: 0.55


model terminate | Simulation terminated by an assert at time: 0.55


The initialization finished successfully without homotopy method.
state event at time=0.500000000151
[1] time > 0.5
check for discrete changes at time=0.500000000151
status of relations at time=0.500000000151
status of zero crossings at time=0.500000000151
check for discrete changes at time=0.500000000151
terminal event at stop time 1

Everything works smoothly in this case; the event is correctly located, within the tolerance set for the zero-crossing function solutions. We should handle time events more efficiently, see #2152, but that is another story.

If I now uncomment the assert, I get the following simulation log

Integrator attempt to handle a problem with a called assert.
The following assertion has been violated at time 0.550000
if time > 0.5 then condition else true
error
model terminate | Simulation terminated by an assert at time: 0.55
The initialization finished successfully without homotopy method.
The following assertion has been violated at time 0.600000
if time > 0.5 then condition else true
error
check for discrete changes at time=0.5
Integrator attempt to handle a problem with a called assert.
The following assertion has been violated at time 0.550000
if time > 0.5 then condition else true
error
model terminate | Simulation terminated by an assert at time: 0.55

It looks like the event-location mechanism was never triggered, nor the event processed, because the assert killed it first.

I tried to figure out what happens, by running the simulation in my mind. This model gets run by fixed time step Euler, because it has no differential equations, so we get steps at time 0.0, 0.1, 0.2, 0.3, 0.4, and 0.5, all being computed with the active else branch, which is selected during initialization. So far so good.

Then, a step at time = 0.6 is attemped. According to ​Appendix B of the specification, the following actions should be performed

1. The DAE (B.1a) is solved by a numerical integration method. In this phase the conditions c of theif- and when-clauses, as well as the discrete-time variables z and m are kept constant.
2. During integration, all relations from (B.1d) are monitored. If one of the relations changes its value an event is triggered, i.e., the exact time instant of the change is determined and the integration is halted.
3. At an event instant, (B.1) is a mixed set of algebraic equations which is solved for the Real, Boolean and Integer unknowns
4. After an event is processed, the integration is restarted with 1.

Note, that both the values of the conditions c as well as the values of z and m (all discrete-time Real, Boolean and Integer variables) are only changed at an event instant and that these variables remain constant during continuous integration.

Apparently, our runtime is doing something different.

According to point 1. above, as the step at time t = 0.6 is attempted, we should still consider the else branch of the if-equation, so x would still be zero, condition should remain false, and the assert should not even be called, because it is on the other branch.

The change of sign of zcf = time - 0.5 would then be detected, and the iterations to find the event instant performed, all by remaining on the else branch of the if-equation.

The event would then be processed, and only at that point we would switch to the other branch, but then, thanks to the causalization, x = 0.5 should be computed first, then condition = x > 0 = true, then assert(condition, ...) would not fail.

I'm not too much into the details of the current implementation, but is seems to me that it is not following the specification, which would be a problem. Do I miss something?

[Francesco Casella]

Ticket retargeted after milestone closed