Generated FMU missing events in other tools

[Andreas Heuermann]

When exporting BouncingBall

model BouncingBallFMI20
  parameter Real e=0.7 "coefficient of restitution";
  parameter Real g=9.81 "gravity acceleration";
  Real h(start=1) "height of ball";
  Real v "velocity of ball";
  Real v_new;
equation
  der(v) = -g;
  der(h) = v;

  when h <= 0.0 then
    v_new = -e*pre(v);
    reinit(v, v_new);
  end when;
end BouncingBallFMI20;

as ModelExchange 2.0 FMU with omc and simulating with Dymola or FMUComplianceChecker the event triggered by when h <= 0.0 is not found by those tools.

At least OMC and OMSimulator are simulating correctly.
Maybe our fmi2GetEventIndicators function is not correct any more.

[Andreas Heuermann]

Bouncing Ball FMU

[Francesco Casella]

@AnHeuermann, can you discuss with Bernhard who can take care of this and update the owner accordingly? A blocker ticket without a responsible person has little use...

Thanks!

[Francesco Casella]

Releasing 1.14.0 which is stable and has many improvements w.r.t. 1.13.2.

This issue, previously marked as blocker for 1.14.0, is rescheduled to 1.15.0

[Andreas Heuermann]

The problem is still there but only for ModelExchange. For CoSimulation it is working just fine.
My test script to compare with FMI Compliance Checker:

// Clean up
system("rm -rf BouncingBallFMI20*");

// Load model
loadString("model BouncingBallFMI20
  parameter Real e=0.7 \"coefficient of restitution\";
  parameter Real g=9.81 \"gravity acceleration\";
  Real h(start=1) \"height of ball\";
  Real v \"velocity of ball\";
  Real v_new;
equation
  der(v) = -g;
  der(h) = v;

  when h <= 0.0 then
    v_new = -e*pre(v);
    reinit(v, v_new);
  end when;
end BouncingBallFMI20;
"); getErrorString();

// Build FMU
translateModelFMU(BouncingBallFMI20, version="2.0", fmuType="me"); getErrorString();

// Simulate with OMSimulator
system(getInstallationDirectoryPath() + "/bin/OMSimulator BouncingBallFMI20.fmu --mode=me --tolerance=1e-6 --resultFile=\"BouncingBallFMI20_OMS_res.mat\" --stopTime=0.5 --suppressPath=true --tempDir=\"BouncingBallFMI20-tmp-1\"", "BouncingBallFMI20_cs_systemCall.log"); getErrorString();
readFile("BouncingBallFMI20_cs_systemCall.log"); getErrorString();

// Simulate with FMI Compliance Checker
mkdir("BouncingBallFMI20-tmp-2"); getErrorString();
system("/path/to/FMUComplianceChecker/install/fmuCheck.linux64 -d -f -h 0.002 -l 2 -n 0 -o \"BouncingBallFMI20_fmuCheck_res.csv\" -s 0.5 -t \"BouncingBallFMI20-tmp-2\" -k me BouncingBallFMI20.fmu", "BouncingBallFMI20_cs_fmuCheck_systemCall.log"); getErrorString();
readFile("BouncingBallFMI20_cs_fmuCheck_systemCall.log"); getErrorString();

So comparing CS and ME exported FMU functions should get us to a solution.

[Andreas Heuermann]

Bouncing ball example: Red=OMSimulator result for h, Blue=FMUCheck result for h

[Andreas Heuermann]

A bit more digging into log files shows what FMUCheck is doing:

Directly after the event continuous mode will be entered with fmi2EnterContinuousTimeMode. FMUCheck will read fmi2GetReal: v = 3.117618000000008 but then sets fmi2SetContinuousStates: v = -4.473360000000011.

[VERBOSE][FMICAPI] Calling fmi2EnterEventMode
	[FMU][logEvents][FMU status:OK] fmi2EnterEventMode
[VERBOSE][FMICAPI] Calling fmi2NewDiscreteStates
	[FMU][logFmi2Call][FMU status:OK] fmi2NewDiscreteStates
	[FMU][logFmi2Call][FMU status:OK] fmi2EventUpdate: Start Event Update! Next Sample Event 0
	[FMU][logFmi2Call][FMU status:OK] fmi2EventUpdate: Need to iterate(discrete changes)!
	[FMU][logFmi2Call][FMU status:OK] fmi2EventUpdate: newDiscreteStatesNeeded true
	[FMU][logFmi2Call][FMU status:OK] fmi2EventUpdate: Checked for Sample Events! Next Sample Event 0
[VERBOSE][FMICAPI] Calling fmi2NewDiscreteStates
	[FMU][logFmi2Call][FMU status:OK] fmi2NewDiscreteStates
	[FMU][logFmi2Call][FMU status:OK] fmi2EventUpdate: Start Event Update! Next Sample Event 0
	[FMU][logFmi2Call][FMU status:OK] fmi2EventUpdate: newDiscreteStatesNeeded false
	[FMU][logFmi2Call][FMU status:OK] fmi2EventUpdate: Checked for Sample Events! Next Sample Event 0
	[FMU][logFmi2Call][FMU status:OK] fmi2GetEventIndicators: z0 = 1
[VERBOSE][FMICAPI] Calling fmi2EnterContinuousTimeMode
	[FMU][logFmi2Call][FMU status:OK] fmi2EnterContinuousTimeMode
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: h = -0.006545240000002877
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: v = 3.117618000000008
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: der(h) = 3.117618000000008
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: der(v) = -9.81
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: v_new = 3.117618000000008
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: e = 0.7
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: g = 9.81
	[FMU][logFmi2Call][FMU status:OK] fmi2GetDerivatives: der(h) = 3.117618000000008
	[FMU][logFmi2Call][FMU status:OK] fmi2GetDerivatives: der(v) = -9.81
[VERBOSE][FMUCHK] Simulation time: 0.456
	[FMU][logFmi2Call][FMU status:OK] fmi2SetTime: time=0.4560000000000003
	[FMU][logFmi2Call][FMU status:OK] fmi2SetContinuousStates: h = -0.0003100040000028562
	[FMU][logFmi2Call][FMU status:OK] fmi2SetContinuousStates: v = -4.473360000000011
	[FMU][logFmi2Call][FMU status:OK] fmi2GetEventIndicators: z0 = 1
	[FMU][logFmi2Call][FMU status:OK] fmi2CompletedIntegratorStep
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: h = -0.0003100040000028562
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: v = -4.473360000000011
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: der(h) = -4.473360000000011
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: der(v) = -9.81
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: v_new = 3.117618000000008
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: e = 0.7
	[FMU][logFmi2Call][FMU status:OK] fmi2GetReal: g = 9.81
	[FMU][logFmi2Call][FMU status:OK] fmi2GetDerivatives: der(h) = -4.473360000000011
	[FMU][logFmi2Call][FMU status:OK] fmi2GetDerivatives: der(v) = -9.81

To trigger fmi2GetContinuousStates in fmuCheck eventInfo->valuesOfContinuousStatesChanged must be true. OMSimulator is calling fmi2GetContinuousStates in any case and therefore working with this.

I would argue that FmuCheck is not correct. The FMI documentation states

If valuesOfContinuousStatesChanged = fmi2True. then at least one element of the
continuous state vector has changed its value due to the function call. The new
values of the states can be retrieved with fmi2GetContinuousStates or individually
for each state for which reinit = "true" by calling getReal.

For the first call to fmi2NewDiscreteStates valuesOfContinuousStatesChanged=fmi2True and for the second one it's fmi2False.
So the integrator needs to keep track of the values of eventInfo if it calling fmi2NewDiscreteStates multiple times.

[Andreas Heuermann]

I'm resetting eventInfo now only at fmi2EnterEventMode. In fmi2NewDiscreteStates those values will only be updated.
So if valuesOfContinuousStatesChanged becomes fmi2True it will stay that way until the next Event a.k.a fmi2EnterEventMode get's called.

​PR 959 should fix this.

[Francesco Casella]

Release 1.15.0 was scrapped, because replaceable support eventually turned out to be more easily implemented in 1.16.0. Hence, all 1.15.0 tickets are rescheduled to 1.16.0

[Andreas Heuermann]

Fixed with ​commit e374ec0.

[Andreas Heuermann]

Testing a different version of Bouncing Ball revealed that the problem is still not solved.

model BouncingBall
  parameter Real e=0.7 "coefficient of restitution";
  parameter Real g=9.81 "gravity acceleration";
  output Real h(fixed=true, start=1) "height of ball";
  output Real v(fixed=true) "velocity of ball";
  Boolean flying(fixed=true, start=true) "true, if ball is flying";
  Boolean impact;
  Real v_new(fixed=true);
  Integer foo;
equation
  impact = h <= 0.0;
  foo = if impact then 1 else 2;
  der(v) = if flying then -g else 0;
  der(h) = v;
  when {h <= 0.0 and v <= 0.0,impact} then
    v_new = if edge(impact) then -e*pre(v) else 0;
    flying = v_new > 0;
    reinit(v, v_new);
  end when;
annotation(
    experiment(StartTime = 0, StopTime = 1, Tolerance = 1e-6, Interval = 0.006));
end BouncingBall;

But this time the FMU will miss the first event for all FMU solvers I tested (OMSimulator, FMUCrossCheck, FMPy).

See attached screenshot for height h.

[Andreas Heuermann]

Height "h" of BouncingBall

[Francesco Casella]

I see two problems with this bouncing ball model.

A) The first condition is redundant

​Section 8.3.5 of the specification states:

The statements within a when-equation are activated when the scalar expression or any of the elements of the vector expression becomes true.

According to Boolean logic, given that impact = h <= 0.0, the first element of the activation vector is redundant, because I can't see any case whereby h <= 0.0 and v <= 0.0 becomes true without also impact becoming true.

B) Removing a redundant condition affects the simulation results

If I remove the first element of the when triggering condition, which is apparently redundant, and only keep when impact, the simulation (even with the standard C runtime) proceeds correcty until about time = 2.5, then the ball falls through the table. From a mathematical standpoint, this is weird.

I guess this means there is something fishy going on in the backend well ahead of FMI code generation. Do I miss something?

I'm not sure if addressing this fishiness will solve the specific issue raised by the ticket, but I would at least give it a try.

[Francesco Casella]

BTW, this model works with the C runtime, you can check if it also works with FMUs

model BouncingBall
  parameter Real e=0.7 "coefficient of restitution";
  parameter Real g=9.81 "gravity acceleration";
  parameter Real v_stick = 1e-4 "below this velocity, the ball is stuck to the table";
  
  output Real h(fixed=true, start=1) "height of ball";
  output Real v(fixed=true) "velocity of ball";
  Boolean flying(fixed=true, start=true) "true, if ball is flying";
  Boolean impact;
  Real v_new(fixed=true);
  Integer foo;
equation
  impact = h <= 0.0;
  foo = if impact then 1 else 2;
  der(v) = if flying then -g else 0;
  der(h) = v;
  when impact then
    v_new = -e*pre(v);
    flying = v_new > v_stick;
    reinit(v, if flying then v_new else 0);
  end when;
annotation(
    experiment(StartTime = 0, StopTime = 4, Tolerance = 1e-6, Interval = 0.006));
end BouncingBall;

Unfortunately, the ball eventually falls through the table if you take v_stick < 4.89e-5. I couldn't figure out why. Maybe investigating this could also shed some light on why some solvers do not behave as expected.

[Andreas Heuermann]

Replying to casella:

If I remove the first element of the when triggering condition, which is apparently redundant, and only keep when impact, the simulation (even with the standard C runtime) proceeds correcty until about time = 2.5, then the ball falls through the table. From a mathematical standpoint, this is weird.

I guess this means there is something fishy going on in the backend well ahead of FMI code generation. Do I miss something?

Mhhh... Yes (A and B) or A <=> A, but we have a when around that, so we actually have:

when{A and B, A} <==> edge(A and B) or edge(A)
                 <==> [(A and B) and not(A_pre and B_pre)] or [A and not A_pre]
case 1: [(A and B) and not(A_pre and B_pre)] is true
    ==> A=true, B=true, (not(A_pre) or not(B_pre))=true
case 2: [A and not A_pre] is true
    ==> A=true, A_pre=false, B doesn't matter, B_pre doesn't matter

So the when clause can become true if and only if one of the following cases is happening:

case 1_1: A=true, A_pre=false, B=true,  B_pre= true
case 1_2: A=true, A_pre=true,  B=true,  B_pre= false
case 1_3: A=true, A_pre=false, B=true,  B_pre= false

case 2_1: A=true, A_pre=false, B=true,  B_pre=true
case 2_2: A=true, A_pre=false, B=false, B_pre=true
case 2_3: A=true, A_pre=false, B=true,  B_pre=false
case 2_4: A=true, A_pre=false, B=false, B_pre=false

We have a few cases that are duplicate but not all case 1_x are equal to case 2_x, so we can't reduce the when condition to when{A}.

I actually needed pen and paper for this one... But I'm glad I did it. I was roughly aware that it's not so straight forward but could not pin the difference down right away.
And I hope I did it correctly :-P

Replying to casella:

BTW, this model works with the C runtime, you can check if it also works with FMUs

No, it has the same problem, that the FMU is missing the first event.

[Andreas Heuermann]

I noticed that changing the causality of h and v from independent to output is causing the changed results.
So it seems that this problem is related to #6040. Not only get the names changed but also the event handling will be changed for those alias variables.

[Francesco Casella]

I'll look at it ASAP, sorry

[Francesco Casella]

Replying to AnHeuermann:

Mhhh... Yes (A and B) or A <=> A, but we have a when around that, so we actually have:

when{A and B, A} <==> edge(A and B) or edge(A)
                 <==> [(A and B) and not(A_pre and B_pre)] or [A and not A_pre]
case 1: [(A and B) and not(A_pre and B_pre)] is true
    ==> A=true, B=true, (not(A_pre) or not(B_pre))=true
case 2: [A and not A_pre] is true
    ==> A=true, A_pre=false, B doesn't matter, B_pre doesn't matter

So the when clause can become true if and only if one of the following cases is happening:

case 1_1: A=true, A_pre=false, B=true,  B_pre= true
case 1_2: A=true, A_pre=true,  B=true,  B_pre= false
case 1_3: A=true, A_pre=false, B=true,  B_pre= false

case 2_1: A=true, A_pre=false, B=true,  B_pre=true
case 2_2: A=true, A_pre=false, B=false, B_pre=true
case 2_3: A=true, A_pre=false, B=true,  B_pre=false
case 2_4: A=true, A_pre=false, B=false, B_pre=false

We have a few cases that are duplicate but not all case 1_x are equal to case 2_x, so we can't reduce the when condition to when{A}.

Of course you are right.

I'm still wondering about the physical meaning of those conditions. The event should be triggered when the height becomes negative and the speed is already negative. Unfortunately there is no way to express that in Modelica, because you can't combine when ad if. On the other hand, some of those conditions are not feasible from a physical point of view, given the governing differential equations for the system.

These two cases are covered with both formulations

A=true, A_pre=false, B=true, B_pre= true
A=true, A_pre=false, B=true, B_pre= false

the first meaning that the ball is trying to go through the table with a negative velocity (of course it should bounce then), the second that the ball is trying to go through the table exactly when its speed is going from negative to positive, which is not possible because der(v) is either -g or zero, so the speed can only go from positive to negative, not the other way round. So, only the first case is reachable and relevant.

This case is only covered by the first formulation:

A=true, A_pre=true, B=true, B_pre= false

meaning that the ball is already firmly below the table, when the ball speed becomes negative. This condition is also unreachable; in order to get below the table, the ball has to go through it with negative speed, and this should cause it to bounce (see previous firs case), so you'll never get to this point.

These cases are only covered by the second formulation

A=true, A_pre=false, B=false, B_pre=true
A=true, A_pre=false, B=false, B_pre=false

The first corresponds to the ball trying to go through the table at the same time that the speed becomes negative. In principle this shouldn't happen, because at that point in time the speed is still zero, so h can't go from positive to negative, since der(h) = v. But I guess that's precisely what happens at the end of the Zeno transient, when the ball should rest on the table. And in fact, through it it goes. In principle, it shouldn't happen, but I guess numerical errors and tolerances play a bad role here. It is interesting that this singular case is ruled out by your original formulation, which I found quite odd. In fact, it makes sense from this point of view.

The second corresponds to the ball trying to go through the table when the speed is positive. This is for sure not possible, since der(h) = v.

As far as I am concerned, I think it's safer and easier to understand if you put in a sticking condition in the model, which avoids getting to a point of singularity of the combination between ODEs and events when you have h = 0 and very small v, which apparently makes the second formulation to fail. That's what I did with my alternative model. But I understand this has nothing to do with the topic of this ticket, which involved the first event, not the behaviour after the end of the Zeno transient.

[Andreas Heuermann]

Should be fixed when PR from #6050 is merged.

[Andreas Heuermann]

Fixed with #6050