Set input derivatives when simulating CS FMUs

[greenwoodms06]

@t-sommer the title of my issue was going to be exactly this so I'll chance commenting here first instead of creating a new issue first.

When using the simulate_fmu function in CS mode, I expected the set_input_derivatives=True would change the input variables from values to derivatives. However, there doesn't seem to be a difference when I switch between True and False. For testing I am using a simple Lotka Volterra model with a control variable on the y (predator) term:

  # Lotka Volterra Model FMUd via Modelica/Dymola
  # der(x) = alpha*x - beta*x*y;
  # der(y) = -gamma*y + delta*x*y + u;

    set_input_derivatives = True #False # The image below looks the same regardless. The issue discussed in #227 seems to also occurs on step changes (i.e., some overshoot thing occurs) when you switch between boolean values. That may be a separate issue though
    dtype = [('time', np.double), ('u', np.double)]
    signals = np.array([(0.0, 1)], dtype=dtype)

Here is what I believe thought would happen when I switched between True and False:

So questions:

1. Is my assumption that the set_input_derivatives flag is supposed to change the inputs variable from values to derivatives of the variables?
2. If so, do you have any idea what is occurring? If not, is there something that will help be better understand what set_input_derivatives=True is doing? I expected it to be doing fmu.setRealInputDerivatives().
3. Is there a way to mix and match derivatives and variables while still leverageing simulate_fmu? For example if you had two inputs, one you want to specify the value and the other you want to specify the derivative? I know you can't plan for all scenarios and that I can do it manually myself. Just thinking out loud! :)

[t-sommer]

Can you share the Modelica model?

[greenwoodms06]

Yup. You can get it here: https://github.com/ORNL-Modelica/ModelicaPy/tree/master/modelicapy/raven/fmpy

The simulate.py should run standalone (no RAVEN needed).

[t-sommer]

With this code

from fmpy import *
import numpy as np

filename = 'simulator.fmu'

signals = np.array([(-1, -1), (1, 1), (3, 4)], dtype=[('time', np.double), ('u', np.double)])

result = simulate_fmu(filename, input=signals, output=['u', 'x', 'y'], stop_time=2, output_interval=0.1,
                      set_input_derivatives=True,
                      fmi_call_logger=print)

plot_result(result)

I get the following FMI log

fmi2Instantiate(instanceName="simulator", fmuType=1, guid="{1461401a-4053-40d3-b99f-fff118f4703c}", resourceLocation="file:///C:/Users/xxx/AppData/Local/Temp/tmp63eb9_va/resources", callbacks=<fmpy.fmi2.fmi2CallbackFunctions object at 0x000001AAA56354C0>, visible=0, loggingOn=0) -> 0x1aaa517e720
fmi2SetupExperiment(component=0x1aaa517e720, toleranceDefined=1, tolerance=0.0001, startTime=0.0, stopTimeDefined=1, stopTime=2.0) -> OK
fmi2EnterInitializationMode(component=0x1aaa517e720) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.0]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2ExitInitializationMode(component=0x1aaa517e720) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[10.0, 5.0, 0.0]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.0]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=0.0, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[5.934365381785415, 8.256522858103269, 0.1]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.10000000000000003]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=0.1, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[2.7983573763100504, 9.334239942337462, 0.20000000000000004]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.19999999999999996]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=0.2, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[1.3238786223919625, 8.432111123465882, 0.3]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.30000000000000004]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=0.30000000000000004, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.7127780620673684, 6.9182173908263005, 0.4]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.39999999999999997]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=0.4, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.4464642803344967, 5.461853441937712, 0.49999999999999994]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.5]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=0.5, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.3198200143246329, 4.249688610020823, 0.6]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.6000000000000001]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=0.6000000000000001, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.25510318122886877, 3.2957302659570518, 0.7000000000000001]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.7000000000000002]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=0.7000000000000001, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.2213477386530397, 2.5652427540480853, 0.8000000000000002]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.8]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=0.8, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.20473990467595654, 2.015097082913305, 0.9]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[0.8999999999999999]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=0.9, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.198651136158947, 1.6058927910504204, 0.9999999999999999]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[1.0]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.0]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=1.0, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.19967530122902769, 1.3050898912718183, 1.1]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[1.1500000000000001]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.5]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=1.1, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.20588631274241542, 1.0935806761396247, 1.3000000000000003]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[1.3000000000000003]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.5]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=1.2000000000000002, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.21609916482436667, 0.9476950049598687, 1.4500000000000004]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[1.4500000000000002]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.5]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=1.3, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.22955068737062076, 0.851233552649409, 1.6000000000000003]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[1.6]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.5]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=1.4000000000000001, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.24572108252478383, 0.7923540826355314, 1.7500000000000002]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[1.75]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.5]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=1.5, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.2641905255579365, 0.761993194257294, 1.9000000000000001]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[1.9000000000000001]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.5]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=1.6, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.28459165048585894, 0.7533979118400448, 2.0500000000000003]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[2.0500000000000003]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.5]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=1.7000000000000002, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.3065658538727145, 0.7615115450533201, 2.2]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[2.2]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.5]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=1.8, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.3297489858029899, 0.7823397893659345, 2.3500000000000005]) -> OK
fmi2SetReal(component=0x1aaa517e720, vr=[352321536], nvr=1, value=[2.35]) -> OK
fmi2SetRealInputDerivatives(c=0x1aaa517e720, vr=[352321536], nvr=1, order=<fmpy.simulation.c_long_Array_1 object at 0x000001AAA56358C0>, value=[1.5]) -> OK
fmi2DoStep(component=0x1aaa517e720, currentCommunicationPoint=1.9000000000000001, communicationStepSize=0.1, noSetFMUStatePriorToCurrentPoint=1) -> OK
fmi2GetReal(component=0x1aaa517e720, vr=[335544320, 335544321, 352321536], nvr=3, value=[0.3537587716068322, 0.8129561481162643, 2.499999999999999]) -> OK
fmi2Terminate(component=0x1aaa517e720) -> OK
fmi2FreeInstance(component=0x1aaa517e720)

where the first order derivative is calculated from the input and set via fmi2SetRealInputDerivatives().

[greenwoodms06]

Hmm. ok so it looks like set_input_derivative=True calculates the slope between the provided signals data? To me this behavior seems really strange. I expected when set_input_derivative=True that instead of setting a value ('u') you use the input data to set the derivative of the value.

Below is what I expected the behavior to be: signals = np.array([(-1, -1), (1, 1), (3, 4)], dtype=[('time', np.double), ('u', np.double)])

*PseudoLog*

// Initialize
fmi2SetRealInputDerivatives(order=???,value=[-1])

// Simulation
fmi2DoStep(currentCommunicationPoint = 0.0, communicationStepSize=0.1)
.
.
.
fmi2SetRealInputDerivatives(order=???,value=[1])
fmi2DoStep(currentCommunicationPoint = 1.0, communicationStepSize=0.1)
.
.
.
fmi2SetRealInputDerivatives(order=???,value=[4])
fmi2DoStep(currentCommunicationPoint = 3.0, communicationStepSize=0.1)
.
.
.

What are your thoughts on this approach? I can't quite grasp how the current approach would be used.

P.S. I left order=??? because it doesn't seem that simulate_fmu permits one to specify the order of the derivative... Assumes 1 I suppose?

If the behavior I am showing is more preferred. I would think the set_input_derivatives variable could potentially be dropped and instead a new input variable be included (e.g., inputs_der) which allows for the same format as inputs but is tied to fmi2SetRealInputDerivatives. though the data format would need to be modified to permit one to specify the order as well.

[t-sommer]

ok so it looks like set_input_derivative=True calculates the slope between the provided signals data?

That's correct. If you know the higher order derivatives for your input signals, you have to write a custom simulation loop and set them manually. You can use https://github.com/CATIA-Systems/FMPy/blob/master/fmpy/examples/custom_input.py as a starting point or call setReal() and setRealInputDerivatives() in the step_finished callback.

[greenwoodms06]

Yes I know I can do it manually... I've done it before... Is the current treatment of specifying derivatives for simulate_fmu really the way you want it to be? I have a hard time understanding how it is useful though I imagine there was a reason behind its current behavior. Is looking to switch it from its current behavior to the one I propose of potential interest or is it a bad decision for some reason?