InputEvent and TimeEvent at the same moment may lead to inconsistencies

[joerg-tt]

If an FMU (FMI 2.0 ME) issues a time event at a certain point in time and at that same instant an input event is to occur, the simulation loop in simulateME() may handle this in an unexpected way. We have observed two different inconsistencies that resulted in input events to be triggered one step too early (example 1) and in another case where an input event was triggered too late (example 2).

The relevant code range is

FMPy/fmpy/simulation.py

Lines 1043 to 1055 in 9e9493c

	# get the next input event
	t_input_event = input.nextEvent(time)
	# check for input event
	input_event = t_input_event <= t_next
	if input_event:
	t_next = t_input_event
	time_event = nextEventTimeDefined and nextEventTime <= t_next
	if time_event and not fixed_step:
	t_next = nextEventTime

Ex.1: input event triggered too early

Assuming t_input_event <= t_next results in True which in turn sets t_next = t_input_event. If, in the same iteration, a time event occurs this would result in a reassigned t_next = nextEventTime.

Problem: In case nextEventTime is smaller than t_input_event, the condition t_input_event <= t_next would now be Falsey, although, the flag input_event is still set to True which in turn leads to the input event being processed in

FMPy/fmpy/simulation.py

Lines 1112 to 1113 in 9e9493c

	if input_event:
	input.apply(time=time, after_event=True)

.

Proposed solution: if a time event occurs, check again if the input event is still valid (L1054):

        if time_event and not fixed_step:
            t_next = nextEventTime
            input_event = t_input_event <= t_next

Ex.2: input event not triggered when it's due

This issue is related to floating point tolerances and it occurs if t_next is slightly bigger than t_input_event (resulting in input_event=False), while nextEventTime is exactly t_next (input_event=True).

Example:
t_input_event = 1.060000
nextTimeEvent = 1.0599999999999388
t_next = 1.0599999999999388

In this scenario, an input event is not triggered, although it should be. Considering the fact that 1.06 == 1.0599999999999388 when compensating for floating point imperfections.

Proposed solution: also consider an eps when checking for input events

        input_event = t_input_event <= t_next + eps

Conclusion

Both scenarios happened when simulating a real ME-FMU. Our current workaround is to adjust the time axis of the input arrays to be slightly ahead of time by an epsilon. This at least fixes the issue with input events being skipped.

[t-sommer]

Can you share an FMU and input file to reproduce the problem? And can you also reproduce the the behavior with fmusim (https://github.com/modelica/reference-fmus)?

[joerg-tt]

The reference FMUs did not contain a suitable FMU.

But we can reproduce the issue using this simple FMU (built with OMEdit 1.22.3) which produces a time event at 1.0599999999999388:

model timeevent
  input Integer inp;  // Discrete input variable
  output Integer out;
  discrete Real xd2(start = 0);
  parameter Real eventTime = 1.0599999999999388;

equation
  when sample(eventTime, 1) then
    xd2 = pre (xd2) + 1;
  end when;
  out = inp;

end timeevent;

Example script to simulate the FMU while also having an input event at 1.06 (the value steps from 1 to 106):

from fmpy import simulate_fmu
import numpy as np

if __name__ == '__main__':

    dtype = [('time', float), ('inp', float)]
    input = np.array([(0, 0), (1, 1), (1.06, 106)], dtype=dtype)

    result = simulate_fmu(
        "timeevent.fmu",
        stop_time=1.07,
        step_size=0.001,
        fmi_type='ModelExchange',
        input=input,
        output=['out'],
        output_interval=0.001,
    )

    np.set_printoptions(precision=100, threshold=10000)
    print(result[1050:])

The output is as follows when running with fmpy 0.3.21:

...
 (1.0569999999999997,   1) (1.0579999999999998,   1)
 (1.0589999999999997,   1) (1.059999999999939 ,   1)
 (1.059999999999939 ,   1) (1.06              ,   1)
 (1.06              , 106) (1.061             , 106)
 (1.0619999999999998, 106) (1.0629999999999997, 106)
 (1.0639999999999998, 106) (1.0649999999999997, 106)
...

I have also introduced the proposed changes from above (input_event = t_input_event <= t_next + eps) which results in the following output:

...
 (1.0569999999999997,   1) (1.0579999999999998,   1)
 (1.0589999999999997,   1) (1.059999999999939 ,   1)
 (1.059999999999939 , 106) (1.06              , 106)
 (1.06              , 106) (1.061             , 106)
 (1.0619999999999998, 106) (1.0629999999999997, 106)
 (1.0639999999999998, 106) (1.0649999999999997, 106)
...

In the original version, at time t=1.059999999999939 and even t=1.06 the value still remains 1 while in the modified version, the value 106 is already written at the second t=1.059999999999939.

[t-sommer]

In the original version, at time [...] t=1.06 the value still remains 1 [...].

It changes from 1 to 106 (reformatted for better readability):

(1.059999999999939 ,   1)
(1.06              ,   1)
(1.06              , 106)
(1.061             , 106)

[joerg-tt]

This is correct. But my point was that it should already change to 106 at the second instant of 1.059999999999939 because this is basically the same value as 1.06, only off by a floating point precision loss.

[joerg-tt]

Maybe the formatting of the output is confusing (it must be read left-to-right, top-to-bottom). Here is a reformatted version (top-down):

# vanilla fmpy 0.3.21

 (1.0569999999999997,   1)
 (1.0579999999999998,   1)
 (1.0589999999999997,   1)
 (1.059999999999939 ,   1)
 (1.059999999999939 ,   1)
 (1.06              ,   1)
 (1.06              , 106)  # <- change happens too late at the third instance of 1.06=1.05999999..
 (1.061             , 106)
 (1.0619999999999998, 106)
 (1.0629999999999997, 106)
 (1.0639999999999998, 106)
 (1.0649999999999997, 106)

# modified version

 (1.0569999999999997,   1)
 (1.0579999999999998,   1)
 (1.0589999999999997,   1)
 (1.059999999999939 ,   1)
 (1.059999999999939 , 106)  # <- change happens at the second instance of 1.06=1.05999999.. (aka after_event=True)
 (1.06              , 106)
 (1.06              , 106)
 (1.061             , 106)
 (1.0619999999999998, 106)
 (1.0629999999999997, 106)
 (1.0639999999999998, 106)
 (1.0649999999999997, 106)

[t-sommer]

@joerg-tt, can you give modelica/Reference-FMUs#594 at try?

[joerg-tt]

Thanks for the suggestion. So far I haven't been using fmusim and I would need to compile it myself in order to try how it behaves there. I have given fmusim-gui 0.0.34 a try to see if it behaves similar to FMpy, but I haven't found an option to export the raw recorded data as CSV. So unfortunately I cannot even tell if the problem exists in fmusim to begin with.

Our focus has solely been FMpy and that's why we proposed a solution in this repository.

[t-sommer]

I would need to compile it myself

Just download the pre-built binaries from GitHub Actions. If these changes work for you I will apply them to FMPy as well.

[joerg-tt]

Thanks for the hint where to find the .exe

I have compared the behavior of fmusim.exe 0.0.34 with the patched variant from https://github.com/modelica/Reference-FMUs/actions/runs/11053616452

TLDR: the result looks good and the patched version does not contain the extra timestamps 1.059999999999939 in the recording.

Details

Command line:

fmusim.exe --interface-type me --stop-time 1.2 --output-interval 0.001 --input-file stim.csv --output-file rec.csv timeevent.fmu

with stim.csv:

time,inp
0.0,0
1.06,1

fmusim.exe 0.0.34

rec.csv

...
1.058,0
1.059,0
1.059999999999939,0
1.059999999999939,0
1.06,0
1.06,1
1.061,1
1.062,1
...

updated fmusim.exe

rec.csv

...
1.058,0
1.059,0
1.06,0
1.06,1
1.061,1
1.062,1
...

[t-sommer]

@joerg-tt, can you give #701 a try?

[joerg-tt]

Sorry for replying so late. I had been on a longer vacation.

I have tried the same script that I used in #692 (comment) with fmpy v0.3.22 that includes the fix and the result looks promising:

[(1.054,   1)
 (1.055,   1)
 (1.056,   1)
 (1.057,   1)
 (1.058,   1)
 (1.059,   1)
 (1.06 ,   1)
 (1.06 , 106)
 (1.061, 106)
 (1.062, 106)]

I have also had a glance at your code changes and I think it looks good. Thank you!