Inaccurate simulation of systems with noEvent and saturations using the C runtime

[Francesco Casella]

When signal saturations are modelled using noEvent(), the C runtime has some issues, as exemplified by the Modelica.Electrical.Analog.Examples.OpAmps.Multivibrator model.

When simulating the model with OMEdit, some numerical errors are reported during the simulation, around some of the events. In fact, the simulation reaches stopTime successfully, but for some reason the final solver statistics and the end message "The simulation finished successfully." are not displayed. This should be corrected.

The testsuite reports a verification failure, see, e.g., the plot of ​opAmp.out.v, which shows how the period of the OMC simulation is significantly shorter than the one of the reference plot.

The experiment annotation sets StopTime = 1, Tolerance = 1e-6, and Interval = 0.001.

If the simulation is run with these parameters (as Jenkins does), the last commutation of opAmp.out.v takes place at time 0.975, much earlier than the value of 0.9835 of the reference plot, which was computed by Dymola using the same parameters.

If I now decrease the Interval to 0.0002, the last commutation takes place at time 0.9826, which is much closer to the reference plot. However, the precision of the location of time events shouldn't really depend so much on the reporting interval length, so there is probably something fishy going on here.

On the other hand, if I turn off Equidistant time grid off, and keep the tolerance to 1e-6, the last commutation takes place even earlier than with 0.001 communication interval, at time 0.97266. Only if I decrease the tolerance to 1e-10 do I get the last commutation to take place at time 0.9843, which is probably very close to the exact result, and is the same that I get with Dymola using the same settings.

Summing up, even with a fairly tight tolerance setting (1e-6), it seems that the actual choice of communication interval (or the absence thereof) has an excessivly high impact on the accuracy of the solution of a simple system (it is a first-order ODE) that has state events. On the other hand, Dymola (which basically uses the same ODE solver, DASSL) doesn't show such a dependency, which probably means there is something wrong with our implementation.

This is immediately apparent if you check the plots of ​c.v, which is the only state variable of the system. If you zoom in around time = 0.034, it is apparent how the OMC solution change direction at a value of -7.4, much earlier than the correct value of -7.5 which is shown by the reference solution.

It looks like the event is not triggered at the zero-crossing time, but rather at the previous step (or communication interval). This seems quite a serious issue for the simulation of systems that heavily depend on state events.

Karim, Andreas, Lennart, can you please have a look?

[Andreas Heuermann]

If im not mistaken, there are no events in this model.
When setting both Dymola and OpenModelica to a tolerance of 1e-06 for DASSL and the same output interval (e.g. 0.001) I get the (nearly) exact same results.

[Andreas Heuermann]

Never mind, I was using the Cpp runtime with the new frontend, not the C runtime.
Only for C runtime the error gets bigger over time. I'm not shure why though.

[Francesco Casella]

In fact you are right, with the current parameter settings, the opAmp model uses this equation to compute the output voltage:

v_out = smooth(0, noEvent(if V0*v_in>vps then vps else if V0*v_in<vns then vns else V0*v_in));

so there are in fact no events in this simulation. Unfortunately, the solver statistics are not shown (see issue 1.), so I didn't notice.

Anyway, it still seems that the C runtime deals with the if-expression in a somwhat wrong way, even in the presence of both smooth() and noEven() operators.

I tried to improve the description of the ticket.

[Andreas Heuermann]

I looked further into the Problem.
The generated code looks correct to me.

/*
equation index: 6
type: SIMPLE_ASSIGN
opAmp._v_out = homotopy(smooth(0, if noEvent(opAmp.simplifiedExpr > opAmp.Vps) then opAmp.Vps else if noEvent(opAmp.simplifiedExpr < opAmp.Vns) then opAmp.Vns else opAmp.simplifiedExpr), opAmp.simplifiedExpr)
*/
void Multivibrator_eqFunction_6(DATA *data, threadData_t *threadData) {
  modelica_boolean tmp0;
  modelica_boolean tmp1;
  modelica_boolean tmp2;
  modelica_real tmp3;
  tmp0 = Greater(data->localData[0]->realVars[11] /* opAmp.simplifiedExpr variable */,data->simulationInfo->realParameter[10] /* opAmp.Vps PARAM */);
  tmp2 = (modelica_boolean)tmp0;
  if(tmp2){
    tmp3 = data->simulationInfo->realParameter[10] /* opAmp.Vps PARAM */;
  } else {
    tmp1 = Less(data->localData[0]->realVars[11] /* opAmp.simplifiedExpr variable */,data->simulationInfo->realParameter[9] /* opAmp.Vns PARAM */);
    tmp3 = (tmp1?data->simulationInfo->realParameter[9] /* opAmp.Vns PARAM */:data->localData[0]->realVars[11] /* opAmp.simplifiedExpr variable */);
  }
  data->localData[0]->realVars[13] /* opAmp.v_out variable */ = homotopy(tmp3, data->localData[0]->realVars[11] /* opAmp.simplifiedExpr variable */);
}

Indeed it looks like opAmp.v_out is switched always one time step earlier then in Dymola, but I can't see why.
I would say that DASSL in Dymola and OMC are simply different implementations and because of that behave different for the combination of homotopy, smooth and noEvent.

I don't think we can solve this up to the next release and I don't knwo if it is realy an error.

[Andreas Heuermann]

Since for other solvers (CVode, Euler) Dymola and OMC the computed solutions is nearly equal we think the problem is independent of homotopy, smooth.
So our DASSL implementation seems to do something different for noEvent. Maybe our step size adaption is different from Dymola and thats the reason for the differences.

[Francesco Casella]

Replying to AnHeuermann:

Since for other solvers (CVode, Euler) Dymola and OMC the computed solutions is nearly equal we think the problem is independent of homotopy, smooth.

I agree. In fact, homotopy is totally irrelevant during simulation (which uses the actual expression by definition) and the presence of noEvent means that the expression value should be taken literally, which shouldn't be such a big deal.

So our DASSL implementation seems to do something different for noEvent. Maybe our step size adaption is different from Dymola and thats the reason for the differences.

Dymola uses a modified version of DASSL, though I don't know where the modifications are. However, the issue is not that we unfavourably compare to Dymola's result, but rather that we unfavourably compare to results computed with tighter tolerance and larger number of intervals, i.e. to the "exact" solution.

The problem with this model is that the opamp output changes very sharply with the input voltage. There are no discontinuities, but due to the very high gain of the opamp, the transition is very fast. Apparently, this creates some issues. If accuracy is an issue, we should probably keep event detection in this case. This would change the results significantly, though for good.

I think the issue here is in the influence of the communication interval on the solution. In principle, I would assume that the variable step-size DASSL algorithm chooses whatever time steps are required to keep the estimated error amplitude below the tolerance, and then it re-interpolates the result on the equidistant time grid. In fact, this turns out not to be the case, since we get very different solutions based on the choice of the communication interval.

I guess the problem is somehow in there. Do you have a clear explanation of how the equidistant time grid output is managed by DASSL and by the C runtime? I could never figure that out fully...

[Francesco Casella]

@AnHeuermann, any progress on this issue? Should we reschedule it to 2.0.0?

[Andreas Heuermann]

Yes, we should reschedule it.