Co-simulation FMU generated by openmodelica not working as designed

[cyrosd]

Description

My co-simulation FMU is not working as intended while the Model exchange part of the FMU works as intended.

The model is a simple PID with 2nd order transfer block activated by a 60 seconds period 50% pulse between 0 and 1 as an input. (Model in the additional files section)

What happens is that when the FMU is loaded in co-simulation (It is supposed to contain both Co-Simulation and Model Exchange which has been confirmed by the official FMU checker ) the output derivates exponentially, alternating between positives and negatives instead of staying in the [-1;2] range it usually stays within.

I have tested the model on 3 different FMI implementations with similar results and came to the conclusion the problem was with OpenModelica. An output csv generated by Daccosim demonstrates what happens in the files section.

Steps to reproduce

Open the file PIDtest_noinput.mo with OMEdit. (Available in the files section) Generate a FMU with File > Export > FMU

Load the FMU in a FMI-compatible software and run a simulation for 100 seconds

Expected behavior

Output generated natively by OpenModelica with a 100 seconds simulation: Output generated by Simulink with Model Exchange FMU with a 100 seconds simulation

Screenshots

Output generated by Simulink with Co-Simulation FMU with a 100 seconds simulation

Additional files

This archive contains:

• PIDtest_noinput.mo : the original model
• PIDtest_noinput.flu : the generated FMU containing both Model Exchange and Co-simulation tools
• PIDtest_CS_output.csv : the csv of the output I got running the FMU through Daccosim (coherent with Simulink graph output) PIDtest.zip

Version and OS

OpenModelica v1.21.0 (64-bit) OMSimulator v2.1.1.post188-gaf996ad-mingw

OS Windows 10 Pro V. 22H2 64 bits

[Arinomo]

I try your fmu with python 3.8.8 and fmpy module as well as exporting the PIDtest_noinput.mo using OpenModelica v1.22.0-dev-39-g23fdc84c46 (64-bit) resulting the expected results

import fmpy
import matplotlib.pyplot as plt
import pandas as pd

result = fmpy.simulate_fmu(PathtoFMU, stop_time=100, output_interval=0.001)
result_df = pd.DataFrame(result)
fig, ax = plt.subplots()
ax.plot(result_df['time'], result_df['pulse1'], result_df['time'], result_df['y'])

[cyrosd]

Are you sure you're using co-simulation? using fmpy.gui, I get this result

[Arinomo]

import pyfmi
import matplotlib.pyplot as plt
fmu = pyfmi.FMUModelCS2('X:/PIDtest/PIDtest_noinput.fmu')
opts = fmu.simulate_options()
opts['ncp'] = int(100/0.001)
result = fmu.simulate(final_time=100, options=opts)
ax.plot(result['time'],result['pulse1'], result3'time'],result['y'])

import matplotlib.pyplot as plt
import pyfmi
fmu = pyfmi.FMUModelCS2('X:/PIDtest/PIDtest_noinput.fmu')
fig, ax = plt.subplots()
result = fmu.simulate(final_time=100)
ax.plot(result['time'],result['pulse1'], result['time'],result['y'])

I dont understand how ncp change the behavior of model. IMHO: since euler is a fixed step size, one should take into account, how big/small the step should be.

Edit: yields the same problem with above example using pyfmi.

import fmpy, pandas
import matplotlib.pyplot as plt
fmu = fmpy.fmi2('X:/PIDtest/PIDtest_noinput.fmu')
res = fmpy.simulate_fmu('X:/PIDtest/PIDtest_noinput.fmu', solver='Euler', stop_time=100)
res2= fmpy.simulate_fmu('X:/PIDtest/PIDtest_noinput.fmu', solver='Euler', stop_time=100, output_interval=0.001)

Edit2: using following Model in OM, where the model was "forced" to used Euler also yield an erroneous result until the interval cahnged to 10ms:

Model with Euler

```.mo model PIDtest Modelica.Blocks.Continuous.PID pid(Td = 0, Ti = 10, k = 10) annotation( Placement(visible = true, transformation(origin = {-2, 4}, extent = {{-10, -10}, {10, 10}}, rotation = 0))); Modelica.Blocks.Continuous.SecondOrder secondOrder(D = 0.7, k = 2, w = 1) annotation( Placement(visible = true, transformation(origin = {28, 4}, extent = {{-10, -10}, {10, 10}}, rotation = 0))); Modelica.Blocks.Math.Feedback feedback annotation( Placement(visible = true, transformation(origin = {-34, 4}, extent = {{-10, -10}, {10, 10}}, rotation = 0))); Modelica.Blocks.Interfaces.RealOutput y annotation( Placement(visible = true, transformation(origin = {70, 4}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {94, 34}, extent = {{-10, -10}, {10, 10}}, rotation = 0))); Modelica.Blocks.Sources.Pulse pulse(period = 60) annotation( Placement(visible = true, transformation(origin = {-92, 4}, extent = {{-10, -10}, {10, 10}}, rotation = 0))); Modelica.Blocks.Interfaces.RealOutput pulse1 annotation( Placement(visible = true, transformation(origin = {66, 46}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {66, 46}, extent = {{-10, -10}, {10, 10}}, rotation = 0))); equation connect(pid.u, feedback.y) annotation( Line(points = {{-14, 4}, {-25, 4}}, color = {0, 0, 127})); connect(secondOrder.u, pid.y) annotation( Line(points = {{16, 4}, {9, 4}}, color = {0, 0, 127})); connect(feedback.u2, secondOrder.y) annotation( Line(points = {{-34, -4}, {-34, -16}, {44, -16}, {44, 4}, {39, 4}}, color = {0, 0, 127})); connect(y, secondOrder.y) annotation( Line(points = {{70, 4}, {39, 4}}, color = {0, 0, 127})); connect(feedback.u1, pulse.y) annotation( Line(points = {{-42, 4}, {-81, 4}}, color = {0, 0, 127})); connect(pulse1, pulse.y) annotation( Line(points = {{66, 46}, {-60.5, 46}, {-60.5, 4}, {-81, 4}}, color = {0, 0, 127})); annotation( uses(Modelica(version = "4.0.0"), Modelica_DeviceDrivers(version = "2.1.1")), Diagram, experiment(StartTime = 0, StopTime = 100, Tolerance = 1e-6, Interval = 0.2), __OpenModelica_commandLineOptions = "--matchingAlgorithm=PFPlusExt --indexReductionMethod=dynamicStateSelection -d=initialization,NLSanalyticJacobian -d=initialization -d=force-fmi-attributes -d=infoXmlOperations ", __OpenModelica_simulationFlags(lv = "stdout,assert,LOG_STATS", s = "euler", variableFilter = ".*")); end PIDtest; ```

[casella]

@cyrosd by default the co-simulation model uses explicit Euler for integration. If the model has a time constant that is smaller than half the sampling time, it will become numerically unstable.

There are two ways out of this: either using a small enough time step, or using an implicit solver inside the CS FMU.