FMI 3 co-simulation timestepping offset by one timestep

steve-biggs-fox commented 5 months ago

I have a simple Python co-simulation FMU that I built with pythonFMU3 for testing, i.e. can I build and run a basic FMU before moving on to more advanced things.

The simple FMU simulates exponential decay, y = exp(-t).

The problem is that when I run the FMU using fmpy, the current time actually uses the time of the previous timestep, thus there is an offset of one timestep in the results. Full details and steps to reproduce are given below...

The code that defines the FMU is as follows:

exp_decay_fmu.py:

# Script to create a basic FMU where the output decays to zero exponentially (y = exp(-t))

import numpy as np
from pythonfmu3 import Fmi3Causality, Fmi3Slave, Fmi3Variability, Float64, Fmi3Initial

class ExpDecay(Fmi3Slave):

    def __init__(self, **kwargs):
        super().__init__(**kwargs)

        self.description = (
            "A basic FMU where the output decays to zero exponentially (y = exp(-t))"
        )

        self.time = 0.0
        self.analytical = 1.0
        self.forward_euler = 1.0
        self.register_variable(
            Float64(
                "time",
                causality=Fmi3Causality.independent,
                variability=Fmi3Variability.continuous,
            )
        )
        self.register_variable(
            Float64(
                "analytical",
                causality=Fmi3Causality.output,
                start=1.0,
                variability=Fmi3Variability.continuous,
                initial=Fmi3Initial.exact,
                description="y = exp(-t) calculated analytically",
            ),
            nested=False,
        )
        self.register_variable(
            Float64(
                "forward_euler",
                causality=Fmi3Causality.output,
                start=1.0,
                variability=Fmi3Variability.continuous,
                initial=Fmi3Initial.exact,
                description="y = exp(-t) calculated via forwards Euler algorithm",
            ),
            nested=False,
        )

        self.default_experiment()

    def do_step(self, current_time, step_size):
        self.analytical = np.exp(-current_time)
        self.forward_euler -= self.forward_euler * step_size
        return True

Notice that the FMU has two output variables, analytical and forward_euler. These compute y = exp(-t) analytically and via a forward Euler algorithm respectively, i.e. calculating the result directly using the current time, and incrementing the value based on the previous value and the time derivative at the previous timestep.

Next, build the FMU using pythonFMU3 from the command line:

pythonfmu3 build -f exp_decay_fmu.py

This generates ExpDecay.fmu, which is an FMI Version 3.0 Co-Simulation FMU.

Then, run the following commands via a Python script or interactive Python session to import fmpy and simulate the FMU:

import fmpy
result = fmpy.simulate_fmu(exp_decay_fmu, stop_time=10, output_interval=0.5)
print(result)

This produces the following output:

[( 0. , 1.00000000e+00, 1.00000000e+00)
 ( 0.5, 1.00000000e+00, 5.00000000e-01)
 ( 1. , 6.06530660e-01, 2.50000000e-01)
 ( 1.5, 3.67879441e-01, 1.25000000e-01)
 ( 2. , 2.23130160e-01, 6.25000000e-02)
 ( 2.5, 1.35335283e-01, 3.12500000e-02)
 ( 3. , 8.20849986e-02, 1.56250000e-02)
 ( 3.5, 4.97870684e-02, 7.81250000e-03)
 ( 4. , 3.01973834e-02, 3.90625000e-03)
 ( 4.5, 1.83156389e-02, 1.95312500e-03)
 ( 5. , 1.11089965e-02, 9.76562500e-04)
 ( 5.5, 6.73794700e-03, 4.88281250e-04)
 ( 6. , 4.08677144e-03, 2.44140625e-04)
 ( 6.5, 2.47875218e-03, 1.22070312e-04)
 ( 7. , 1.50343919e-03, 6.10351562e-05)
 ( 7.5, 9.11881966e-04, 3.05175781e-05)
 ( 8. , 5.53084370e-04, 1.52587891e-05)
 ( 8.5, 3.35462628e-04, 7.62939453e-06)
 ( 9. , 2.03468369e-04, 3.81469727e-06)
 ( 9.5, 1.23409804e-04, 1.90734863e-06)
 (10. , 7.48518299e-05, 9.53674316e-07)]

First column is time, t. Notice that the middle column (analytical result) is offset by one timestep, i.e. result at t = 0.5 is actually exp(0), at t = 1 is actually exp(0.5), etc.

On inspecting the code in fmpy.simulation.simulateCS (line 1171 onward), it appears that for FMI 3, fmpy is doing something fancy with the time; this seems to me the most likely source of the error. Unfortunately, this code is quite difficult to understand (maybe more comments here would help), so I cannot work out how to fix it myself.

Please investigate and fix this bug. Many thanks :)

from line 1231 onward:

    # simulation loop
    while True:

        recorder.sample(time, force=True)

        ....

        input.apply(time)

        if is_fmi1:

            ....

        elif is_fmi2:

            ...

        else:

            t_input_event = input.nextEvent(time)

            t_next = (n_steps + 1) * output_interval

            input_event = t_next > t_input_event

            step_size = t_input_event - time if input_event else t_next - time

            event_encountered, terminate_simulation, early_return, last_successful_time = fmu.doStep(currentCommunicationPoint=time, communicationStepSize=step_size)

            ...

            if early_return and last_successful_time < t_next:
                time = last_successful_time
            else:
                time = t_next
                n_steps += 1

            if use_event_mode and (input_event or event_encountered):

                recorder.sample(last_successful_time, force=True)

                fmu.enterEventMode()

                input.apply(last_successful_time, after_event=True)

                update_discrete_states = True

                while update_discrete_states and not terminate_simulation:
                    update_discrete_states, terminate_simulation, _, _, _, _ = fmu.updateDiscreteStates()

                ...

                fmu.enterStepMode()

        if step_finished is not None and not step_finished(time, recorder):
            recorder.sample(time, force=True)
            break

t-sommer commented 5 months ago

Can you run your FMU with simulate_fmu(..., fmi_call_logger=print) and show the problem based on the FMI call sequence? You can also compare it to the call sequence of fmusim by running fmusim --log-fmi-calls ExpDecay.fmu.

steve-biggs-fox commented 5 months ago

Output from simulate_fmu(..., fmi_call_logger=print) as requested. Discussion below.

fmi3InstantiateCoSimulation(instanceName="ExpDecay", instantiationToken="d63f337f-fd8c-11ee-8e9c-8cec4bc18684", resourcePath="C:\Users\steve\AppData\Local\Temp\tmpvdmj__bp\resources\", visible=False, loggingOn=False, eventModeUsed=False, earlyReturnAllowed=False, requiredIntermediateVariables=NULL, nRequiredIntermediateVariables=0, instanceEnvironment=0x0, logMessage=<CFunctionType object at 0x00000194592CDB10>, intermediateUpdate=<CFunctionType object at 0x00000194592CDBE0>) -> 0x194595fcd70
fmi3EnterInitializationMode(instance=0x194595fcd70, toleranceDefined=False, tolerance=0.0, startTime=0.0, stopTimeDefined=True, stopTime=10.0) -> OK
fmi3ExitInitializationMode(instance=0x194595fcd70) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[1.0, 1.0], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=0.0, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[1.0, 0.5], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=0.5, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.6065306597126334, 0.25], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=1.0, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.36787944117144233, 0.125], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=1.5, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.22313016014842982, 0.0625], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=2.0, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.1353352832366127, 0.03125], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=2.5, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.0820849986238988, 0.015625], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=3.0, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.049787068367863944, 0.0078125], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=3.5, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.0301973834223185, 0.00390625], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=4.0, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.01831563888873418, 0.001953125], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=4.5, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.011108996538242306, 0.0009765625], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=5.0, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.006737946999085467, 0.00048828125], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=5.5, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.004086771438464067, 0.000244140625], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=6.0, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.0024787521766663585, 0.0001220703125], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=6.5, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.0015034391929775724, 6.103515625e-05], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=7.0, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.0009118819655545162, 3.0517578125e-05], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=7.5, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.0005530843701478336, 1.52587890625e-05], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=8.0, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.00033546262790251185, 7.62939453125e-06], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=8.5, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.00020346836901064417, 3.814697265625e-06], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=9.0, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[0.00012340980408667956, 1.9073486328125e-06], nValues=2) -> OK
fmi3DoStep(instance=0x194595fcd70, currentCommunicationPoint=9.5, communicationStepSize=0.5, noSetFMUStatePriorToCurrentPoint=True, eventHandlingNeeded=False, terminateSimulation=False, earlyReturn=False, lastSuccessfulTime=0.0) -> OK
fmi3GetFloat64(instance=0x194595fcd70, valueReferences=[1, 2], nValueReferences=2, values=[7.48518298877006e-05, 9.5367431640625e-07], nValues=2) -> OK
fmi3Terminate(instance=0x194595fcd70) -> OK
fmi3FreeInstance(instance=0x194595fcd70)

So, it seems this is what happens:

The simulation starts at t = 0.0, as expected (fmi3EnterInitializationMode(..., startTime=0.0, ...)).
The outputs also initialize as expected to 1.0 (analytical) and 1.0 (forward Euler) (fmi3GetFloat64(..., values=[1.0, 1.0], ...)).
Then, I think this is where it goes wrong. On the first call to fmi3DoStep(...), we have currentCommunicationPoint=0.0 and communicationStepSize=0.5. The currentCommunicationPoint=0.0 is the problem, I think. This step is unnecessary as we already have the values at t=0.0, i.e. the initialization values. So the first step should have currentCommunicationPoint=startTime+communicationStepSize. As it is, it will calculate y_i = exp(t_i) = exp(0.0) = 1.0 (analytical) and y_i = y_i-1 + (dy/dt)_i-1 dt = 1.0 - 1.0 0.5 = 0.5 (forward Euler). So we end up with an analytical result for t = 0.0 and a forward Euler result for t = 0.5. It appears these then get stored in the results alongside t = 0.5, so it looks like it's the analytical results (i.e. currentCommunicationPoint) that is offset.

t-sommer commented 5 months ago

The communicationPoint for the first call to fmi3DoStep() is startTime, so the time of the FMU is advanced from startTime to startTime + communicationStepSize. How to deal with this information is completely up to the FMU.

For reference you can simulate your FMU with fmusim --log-fmi-calls ExpDecay.fmu from the Reference FMUs which will give you the same call sequence.

steve-biggs-fox commented 5 months ago

Fair enough. I will just have to deal with it.

FWIW, I still think it is wrong. The results for currentCommunicationPoint=0.0 are stored next to t=0.5. Really, it should be called previousCommunicationPoint.

steve-biggs-fox commented 5 months ago

Hmmm... No, actually, you are right. I guess it's how the standard says that time gets advanced. The standard effectively says by running doStep, you are trying to advance time from currentCommunicationPoint by communicationStepSize. If successful, then you have made it to currentCommunicationPoint + communicationStepSize, so that is where the result should be stored. Right? I guess this is just potentially non-obvious if you're new to FMUs. Thanks for your help :)

CATIA-Systems / FMPy

FMI 3 co-simulation timestepping offset by one timestep #664