Importing FMUs

ModelingToolkit is able to import FMUs following the FMI Standard versions 2 and 3. This integration is done through FMI.jl and requires importing it to enable the relevant functionality. Currently Model Exchange (ME) and CoSimulation (CS) FMUs are supported. Events, non-floating-point variables and array variables are not supported. Additionally, calculating the time derivatives of FMU states/outputs is not supported.

Experimental

This functionality is currently experimental and subject to change without a breaking release of ModelingToolkit.jl.

FMUs of full models

Here, we will demonstrate the usage of an FMU of an entire model (as opposed to a single component). First, the required libraries must be imported and the FMU loaded using FMI.jl.

using ModelingToolkit, FMI, FMIZoo, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D

# This is a spring-pendulum FMU from FMIZoo.jl. It is a v2 FMU
# and we are importing it in ModelExchange format.
fmu = loadFMU("SpringPendulum1D", "Dymola", "2022x"; type = :ME)

Model name:	SpringPendulum1D
Type:		0

Following are the variables in the FMU (both states and parameters):

fmu.modelDescription.modelVariables

25-element Vector{FMICore.fmi2ScalarVariable}:
 Name: 'mass_s0' (reference: 16777216)
 Name: 'mass_v0' (reference: 16777217)
 Name: 'fixed.s0' (reference: 16777218)
 Name: 'fixed.flange.s' (reference: 234881026)
 Name: 'fixed.flange.f' (reference: 637534208)
 Name: 'spring.flange_a.s' (reference: 234881026)
 Name: 'spring.flange_a.f' (reference: 637534212)
 Name: 'spring.flange_b.s' (reference: 637534209)
 Name: 'spring.flange_b.f' (reference: 637534208)
 Name: 'spring.s_rel' (reference: 637534210)
 ⋮
 Name: 'der(mass.s)' (reference: 587202560)
 Name: 'mass.L' (reference: 16777222)
 Name: 'mass.flange_a.s' (reference: 637534209)
 Name: 'mass.flange_a.f' (reference: 637534212)
 Name: 'mass.flange_b.s' (reference: 905969669)
 Name: 'mass.flange_b.f' (reference: 100663302)
 Name: 'mass.v' (reference: 33554433)
 Name: 'der(mass.v)' (reference: 587202561)
 Name: 'mass.a' (reference: 587202561)

Next, FMIComponent is used to import the FMU as an MTK component. We provide the FMI major version as a Val to the constructor, along with the loaded FMU and the type as keyword arguments.

@named model = ModelingToolkit.FMIComponent(Val(2); fmu, type = :ME)

\[ \begin{align} \mathtt{mass\_\_vˍt}\left( t \right) &= \mathtt{mass\_\_a}\left( t \right) \\ \mathtt{mass\_\_sˍt}\left( t \right) &= \mathtt{wrapper}\_{1}\left( \left[ \begin{array}{c} \mathtt{mass\_\_s}\left( t \right) \\ \mathtt{mass\_\_v}\left( t \right) \\ \end{array} \right], \left[ \begin{array}{c} \end{array} \right], \left[ \begin{array}{c} \mathtt{mass\_s0} \\ \mathtt{fixed\_\_s0} \\ \mathtt{mass\_\_m} \\ \mathtt{spring\_\_s\_rel0} \\ \mathtt{mass\_\_L} \\ \mathtt{spring\_\_c} \\ \mathtt{mass\_v0} \\ \end{array} \right], t \right) \\ \mathtt{mass\_\_a}\left( t \right) &= \mathtt{wrapper}\_{2}\left( \left[ \begin{array}{c} \mathtt{mass\_\_s}\left( t \right) \\ \mathtt{mass\_\_v}\left( t \right) \\ \end{array} \right], \left[ \begin{array}{c} \end{array} \right], \left[ \begin{array}{c} \mathtt{mass\_s0} \\ \mathtt{fixed\_\_s0} \\ \mathtt{mass\_\_m} \\ \mathtt{spring\_\_s\_rel0} \\ \mathtt{mass\_\_L} \\ \mathtt{spring\_\_c} \\ \mathtt{mass\_v0} \\ \end{array} \right], t \right) \\ \frac{\mathrm{d} \mathtt{mass\_\_s}\left( t \right)}{\mathrm{d}t} &= \mathtt{mass\_\_sˍt}\left( t \right) \\ \frac{\mathrm{d} \mathtt{mass\_\_v}\left( t \right)}{\mathrm{d}t} &= \mathtt{mass\_\_a}\left( t \right) \end{align} \]

Note how hierarchical names in the FMU (e.g. mass.m or spring.f) are turned into flattened names, with __ being the namespace separator (mass__m and spring__f).

Note

Eventually we plan to reconstruct a hierarchical system structure mirroring the one indicated by the variables in the FMU. This would allow accessing the above mentioned variables as model.mass.m and model.spring.f instead of model.mass__m and model.spring__f respectively.

Derivative variables such as der(mass.v) use the dummy derivative notation, and are hence transformed into a form similar to mass__vˍt. However, they can still be referred to as D(model.mass__v).

equations(model)

Since the FMI spec allows multiple names to alias the same quantity, ModelingToolkit.jl creates equations to alias them. For example, it can be seen above that der(mass.v) and mass.a have the same reference, and hence refer to the same quantity. Correspondingly, there is an equation mass__vˍt(t) ~ mass__a(t) in the system.

Note

Any variables and/or parameters that are not part of the FMU should be ignored, as ModelingToolkit creates them to manage the FMU. Unexpected usage of these variables/parameters can lead to errors.

defaults(model)

Dict{Any, Any} with 8 entries:
  mass__L        => 0.0
  wrapper⋆       => FMI2InstanceWrapper(Model name:	SpringPendulum1D…
  spring__c      => 10.0
  mass__m        => 1.0
  mass_s0        => 0.5
  mass_v0        => 0.0
  spring__s_rel0 => 1.0
  fixed__s0      => 0.1

All parameters in the FMU are given a default equal to their start value, if present. Unknowns are not assigned defaults even if a start value is present, as this would conflict with ModelingToolkit's own initialization semantics.

We can simulate this model like any other ModelingToolkit system.

julia> sys = structural_simplify(model)┌ Warning: Internal error: Variable mass__s(t) was marked as being in 0 ~ (wrapper(Any[mass__s(t), mass__v(t)], Float64[], Any[mass_s0, fixed__s0, mass__m, spring__s_rel0, mass__L, spring__c, mass_v0], t))[1] - mass__sˍt(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable mass__v(t) was marked as being in 0 ~ (wrapper(Any[mass__s(t), mass__v(t)], Float64[], Any[mass_s0, fixed__s0, mass__m, spring__s_rel0, mass__L, spring__c, mass_v0], t))[1] - mass__sˍt(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable mass__s(t) was marked as being in 0 ~ (wrapper(Any[mass__s(t), mass__v(t)], Float64[], Any[mass_s0, fixed__s0, mass__m, spring__s_rel0, mass__L, spring__c, mass_v0], t))[2] - mass__a(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable mass__v(t) was marked as being in 0 ~ (wrapper(Any[mass__s(t), mass__v(t)], Float64[], Any[mass_s0, fixed__s0, mass__m, spring__s_rel0, mass__L, spring__c, mass_v0], t))[2] - mass__a(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
Model model:
Equations (2):
  2 standard: see equations(model)
Unknowns (2): see unknowns(model)
  mass__s(t)
  mass__v(t)
Parameters (8): see parameters(model)
  mass__L [defaults to 0.0]
  spring__c [defaults to 10.0]
  mass__m [defaults to 1.0]
  mass_s0 [defaults to 0.5]
  ⋮
Observed (4): see observed(model)
julia> prob = ODEProblem(sys, [sys.mass__s => 0.5, sys.mass__v => 0.0], (0.0, 5.0))┌ Warning: Internal error: Variable mass__s(t) was marked as being in 0 ~ (wrapper(SymbolicUtils.BasicSymbolic{Float64}[mass__s(t), mass__v(t)], Float64[], SymbolicUtils.BasicSymbolic{Float64}[mass_s0, fixed__s0, mass__m, spring__s_rel0, mass__L, spring__c, mass_v0], t))[1] - mass__sˍt(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable mass__v(t) was marked as being in 0 ~ (wrapper(SymbolicUtils.BasicSymbolic{Float64}[mass__s(t), mass__v(t)], Float64[], SymbolicUtils.BasicSymbolic{Float64}[mass_s0, fixed__s0, mass__m, spring__s_rel0, mass__L, spring__c, mass_v0], t))[1] - mass__sˍt(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable mass__s(t) was marked as being in 0 ~ (wrapper(SymbolicUtils.BasicSymbolic{Float64}[mass__s(t), mass__v(t)], Float64[], SymbolicUtils.BasicSymbolic{Float64}[mass_s0, fixed__s0, mass__m, spring__s_rel0, mass__L, spring__c, mass_v0], t))[2] - mass__a(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable mass__v(t) was marked as being in 0 ~ (wrapper(SymbolicUtils.BasicSymbolic{Float64}[mass__s(t), mass__v(t)], Float64[], SymbolicUtils.BasicSymbolic{Float64}[mass_s0, fixed__s0, mass__m, spring__s_rel0, mass__L, spring__c, mass_v0], t))[2] - mass__a(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 5.0)
u0: 2-element Vector{Float64}:
 0.5
 0.0
julia> sol = solve(prob, Tsit5())retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 53-element Vector{Float64}:
 0.0
 0.000166333998669328
 0.000166333998669328
 0.0018296739853626079
 0.0018296739853626079
 0.018252579723965046
 0.018252579723965046
 0.0628049056311173
 0.0628049056311173
 0.13413221122185767
 ⋮
 3.8965775480896196
 4.1865811396986174
 4.1865811396986174
 4.450238792425324
 4.450238792425324
 4.739873150383787
 4.739873150383787
 5.0
 5.0
u: 53-element Vector{Vector{Float64}}:
 [0.5, 0.0]
 [0.5000000830009954, 0.0009980039459963424]
 [0.5000000830009954, 0.0009980039459963424]
 [0.5000100430926603, 0.010977982660156045]
 [0.5000100430926603, 0.010977982660156045]
 [0.5009991925476527, 0.10945467878676397]
 [0.5009991925476527, 0.10945467878676397]
 [0.5117945229973937, 0.37435700342546413]
 [0.5117945229973937, 0.37435700342546413]
 [0.553169961473471, 0.78087710666022]
 ⋮
 [0.5175451840772557, -0.45800694663626157]
 [0.6308152690522649, 1.1838110149591914]
 [0.6308152690522649, 1.1838110149591914]
 [1.0618874727123693, 1.8943128948707495]
 [1.0618874727123693, 1.8943128948707495]
 [1.5519539450537894, 1.2493052977445305]
 [1.5519539450537894, 1.2493052977445305]
 [1.6970341025378801, -0.19764541877156047]
 [1.6970341025378801, -0.19764541877156047]

We can interpolate the solution object to obtain values at arbitrary time points in the solved interval, just like a normal solution.

julia> sol(0.0:0.1:1.0; idxs = sys.mass_a)ERROR: ArgumentError: System model: variable mass_a does not exist

FMUs following version 3 of the specification can be simulated with almost the same process. This time, we will create a model from a CoSimulation FMU.

fmu = loadFMU("SpringPendulum1D", "Dymola", "2023x", "3.0"; type = :CS)
@named inner = ModelingToolkit.FMIComponent(
    Val(3); fmu, communication_step_size = 0.001, type = :CS)

\[ \begin{align} \mathtt{\_\_mtk\_internal\_u}\left( t \right) &= \left[ \begin{array}{c} \mathtt{mass\_\_s}\left( t \right) \\ \mathtt{mass\_\_v}\left( t \right) \\ \end{array} \right] \end{align} \]

This FMU has fewer equations, partly due to missing aliasing variables and partly due to being a CS FMU. CoSimulation FMUs are bundled with an integrator. As such, they do not function like ME FMUs. Instead, a callback steps the FMU at periodic intervals in time and obtains the updated state. This state is held constant until the next time the callback triggers. The periodic interval must be specified through the communication_step_size keyword argument. A smaller step size typically leads to less error but is more computationally expensive.

This model alone does not have any differential variables, and calling structural_simplify will lead to an ODESystem with no unknowns.

structural_simplify(inner)

\[ \begin{align} \end{align} \]

Simulating this model will cause the OrdinaryDiffEq integrator to immediately finish, and will not trigger the callback. Thus, we wrap this system in a trivial system with a differential variable.

@variables x(t) = 1.0
@mtkbuild sys = ODESystem([D(x) ~ x], t; systems = [inner])

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} &= x\left( t \right) \end{align} \]

We can now simulate sys.

prob = ODEProblem(sys, [sys.inner.mass__s => 0.5, sys.inner.mass__v => 0.0], (0.0, 5.0))
sol = solve(prob, Tsit5())

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 10001-element Vector{Float64}:
 0.0
 0.0
 0.001
 0.001
 0.002
 0.002
 0.003
 0.003
 0.004
 0.004
 ⋮
 4.996
 4.996
 4.997
 4.997
 4.998
 4.998
 4.999
 4.999
 5.0
u: 10001-element Vector{Vector{Float64}}:
 [1.0]
 [1.0]
 [1.0010005001667084]
 [1.0010005001667084]
 [1.0020020013340003]
 [1.0020020013340003]
 [1.003004504503377]
 [1.003004504503377]
 [1.004008010677342]
 [1.004008010677342]
 ⋮
 [147.8206921899474]
 [147.8206921899474]
 [147.96858681712632]
 [147.96858681712632]
 [148.1166294129045]
 [148.1166294129045]
 [148.26482012532432]
 [148.26482012532432]
 [148.41315910257674]

The variables of the FMU are discrete, and their timeseries can be obtained at intervals of communication_step_size.

sol[sys.inner.mass__s]

5002-element Vector{Float64}:
 0.5
 0.5
 0.5000032557675739
 0.5000126065628254
 0.5000277374082364
 0.5000489816151542
 0.5000761776744463
 0.5001092988967224
 0.5001483452819825
 0.5001934523600614
 ⋮
 1.6980584968463255
 1.697907233814026
 1.6977499916178909
 1.6975867718588842
 1.6974175761995924
 1.6972424063642613
 1.697061264138831
 1.696874151370973
 1.696874151370973

FMUs of components

FMUs can also be imported as individual components. For this example, we will use custom FMUs used in the test suite of ModelingToolkit.jl.

fmu = loadFMU(
    joinpath(@__DIR__, "..", "..", "..", "test", "fmi", "fmus", "SimpleAdder.fmu");
    type = :ME)
fmu.modelDescription.modelVariables

7-element Vector{FMICore.fmi2ScalarVariable}:
 Name: 'c' (reference: 0)
 Name: 'der(c)' (reference: 1)
 Name: 'a' (reference: 2)
 Name: 'b' (reference: 3)
 Name: 'out' (reference: 4)
 Name: 'out2' (reference: 5)
 Name: 'value' (reference: 6)

This FMU is equivalent to the following model:

@mtkmodel SimpleAdder begin
    @variables begin
        a(t)
        b(t)
        c(t)
        out(t)
        out2(t)
    end
    @parameters begin
        value = 1.0
    end
    @equations begin
        out ~ a + b + value
        D(c) ~ out
        out2 ~ 2c
    end
end

a and b are inputs, c is a state, and out and out2 are outputs of the component.

julia> @named adder = ModelingToolkit.FMIComponent(Val(2); fmu, type = :ME);
julia> isinput(adder.a)true
julia> isinput(adder.b)true
julia> isoutput(adder.out)true
julia> isoutput(adder.out2)true

ModelingToolkit recognizes input and output variables of the component, and attaches the appropriate metadata. We can now use this component as a subcomponent of a larger system.

julia> @variables a(t) b(t) c(t) [guess = 1.0];
julia> @mtkbuild sys = ODESystem(
           [adder.a ~ a, adder.b ~ b, D(a) ~ t,
               D(b) ~ adder.out + adder.c, c^2 ~ adder.out + adder.value],
           t;
           systems = [adder])┌ Warning: Internal error: Variable adder₊c(t) was marked as being in 0 ~ (adder₊wrapper(SymbolicUtils.BasicSymbolic{Float64}[adder₊c(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊a(t), adder₊b(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊value], t))[1] - adder₊cˍt(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable adder₊a(t) was marked as being in 0 ~ (adder₊wrapper(SymbolicUtils.BasicSymbolic{Float64}[adder₊c(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊a(t), adder₊b(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊value], t))[1] - adder₊cˍt(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable adder₊b(t) was marked as being in 0 ~ (adder₊wrapper(SymbolicUtils.BasicSymbolic{Float64}[adder₊c(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊a(t), adder₊b(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊value], t))[1] - adder₊cˍt(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable adder₊c(t) was marked as being in 0 ~ (adder₊wrapper(SymbolicUtils.BasicSymbolic{Float64}[adder₊c(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊a(t), adder₊b(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊value], t))[2] - adder₊out2(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable adder₊a(t) was marked as being in 0 ~ (adder₊wrapper(SymbolicUtils.BasicSymbolic{Float64}[adder₊c(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊a(t), adder₊b(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊value], t))[2] - adder₊out2(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable adder₊b(t) was marked as being in 0 ~ (adder₊wrapper(SymbolicUtils.BasicSymbolic{Float64}[adder₊c(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊a(t), adder₊b(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊value], t))[2] - adder₊out2(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable adder₊c(t) was marked as being in 0 ~ (adder₊wrapper(SymbolicUtils.BasicSymbolic{Float64}[adder₊c(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊a(t), adder₊b(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊value], t))[3] - adder₊out(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable adder₊a(t) was marked as being in 0 ~ (adder₊wrapper(SymbolicUtils.BasicSymbolic{Float64}[adder₊c(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊a(t), adder₊b(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊value], t))[3] - adder₊out(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
┌ Warning: Internal error: Variable adder₊b(t) was marked as being in 0 ~ (adder₊wrapper(SymbolicUtils.BasicSymbolic{Float64}[adder₊c(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊a(t), adder₊b(t)], SymbolicUtils.BasicSymbolic{Float64}[adder₊value], t))[3] - adder₊out(t), but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/work/ModelingToolkit.jl/ModelingToolkit.jl/src/structural_transformation/utils.jl:262
Model sys:
Equations (4):
  4 standard: see equations(sys)
Unknowns (4): see unknowns(sys)
  a(t)
  adder₊c(t)
  b(t)
  c(t)
Parameters (2): see parameters(sys)
  adder₊value [defaults to 1.0]
  adder₊wrapper⋆ [defaults to FMI2InstanceWrapper(Model name:	SimpleAdder
Type:		0, UInt32[0x00000001], UInt32[0x00000000], UInt32[0x00000005, 0x00000004], UInt32[0x00000006], UInt32[0x00000002, 0x00000003], 1.0e-6, nothing)]
Observed (6): see observed(sys)
julia> equations(sys)4-element Vector{Equation}:
 Differential(t)(a(t)) ~ t
 Differential(t)(adder₊c(t)) ~ adder₊cˍt(t)
 Differential(t)(b(t)) ~ adder₊c(t) + adder₊out(t)
 0 ~ adder₊value + adder₊out(t) - (c(t)^2)

Note how the output adder.out is used in an algebraic equation of the system. We have also given sys.c a guess, expecting it to be solved for by initialization. ModelingToolkit is able to use FMUs in initialization to solve for initial states. As mentioned earlier, we cannot differentiate through an FMU. Thus, automatic differentiation has to be disabled for the solver.

prob = ODEProblem(sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0],
    (0.0, 1.0), [sys.adder.value => 2.0])
solve(prob, Rodas5P(autodiff = false))

retcode: Success
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 23-element Vector{Float64}:
 0.0
 1.0e-6
 1.0e-6
 1.1e-5
 1.1e-5
 0.00011099999999999999
 0.00011099999999999999
 0.0011109999999999998
 0.0011109999999999998
 0.011110999999999996
 ⋮
 0.18654232672188578
 0.35269441444276967
 0.35269441444276967
 0.599660138158417
 0.599660138158417
 0.7974985670543848
 0.7974985670543848
 1.0
 1.0
u: 23-element Vector{Vector{Float64}}:
 [1.0, 2.0, 1.0, 2.4494897532242415]
 [1.0000000000005, 2.000004000003, 1.000006000005, 2.449490967528866]
 [1.0000000000005, 2.000004000003, 1.000006000005, 2.449490967528866]
 [1.0000000000605, 2.0000440003630024, 1.0000660006050037, 2.4495032150755596]
 [1.0000000000605, 2.0000440003630024, 1.0000660006050037, 2.4495032150755596]
 [1.0000000061605, 2.000444036965507, 1.000666061608875, 2.4496256995241894]
 [1.0000000061605, 2.000444036965507, 1.000666061608875, 2.4496256995241894]
 [1.0000006171605003, 2.004447705478185, 1.0066721754922148, 2.4508514423874646]
 [1.0000006171605003, 2.004447705478185, 1.0066721754922148, 2.4508514423874646]
 [1.0000617271605003, 2.0448168885834086, 1.0672871759358749, 2.4631989166722073]
 ⋮
 [1.0173990198294076, 2.863377337118342, 2.3131380863021573, 2.7074963892884942]
 [1.0621966749895644, 3.8767560783153834, 3.8827510803308853, 2.9908106651945037]
 [1.0621966749895644, 3.8767560783153834, 3.8827510803308853, 2.9908106651945037]
 [1.1797961406480866, 5.985794247052296, 7.192273566119376, 3.51741531886956]
 [1.1797961406480866, 5.985794247052296, 7.192273566119376, 3.51741531886956]
 [1.3180019822269007, 8.409295569885535, 11.026790556661403, 4.042885501777785]
 [1.3180019822269007, 8.409295569885535, 11.026790556661403, 4.042885501777785]
 [1.500000000000003, 11.855536157589174, 16.505550563651997, 4.691041238518567]
 [1.500000000000003, 11.855536157589174, 16.505550563651997, 4.691041238518567]

CoSimulation FMUs follow a nearly identical process. Since CoSimulation FMUs operate using callbacks, after triggering the callbacks and altering the discrete state the algebraic equations may no longer be satisfied. To resolve for the values of algebraic variables, we use the reinitializealg keyword of FMIComponent. This is a DAE initialization algorithm to use at the end of every callback. Since CoSimulation FMUs are not directly involved in the RHS of the system - instead operating through callbacks - we can use a solver with automatic differentiation.

fmu = loadFMU(
    joinpath(@__DIR__, "..", "..", "..", "test", "fmi", "fmus", "SimpleAdder.fmu");
    type = :CS)
@named adder = ModelingToolkit.FMIComponent(
    Val(2); fmu, type = :CS, communication_step_size = 1e-3,
    reinitializealg = BrownFullBasicInit())
@mtkbuild sys = ODESystem(
    [adder.a ~ a, adder.b ~ b, D(a) ~ t,
        D(b) ~ adder.out + adder.c, c^2 ~ adder.out + adder.value],
    t;
    systems = [adder])
prob = ODEProblem(sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0],
    (0.0, 1.0), [sys.adder.value => 2.0])
solve(prob, Rodas5P())

retcode: Success
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 2004-element Vector{Float64}:
 0.0
 0.0
 1.0e-6
 1.1e-5
 0.00011099999999999999
 0.001
 0.001
 0.002
 0.002
 0.003
 ⋮
 0.996
 0.996
 0.997
 0.997
 0.998
 0.998
 0.999
 0.999
 1.0
u: 2004-element Vector{Vector{Float64}}:
 [1.0, 1.0, 2.2360679786694857]
 [1.0, 1.0, 2.2360679786694857]
 [1.0000000000005, 1.000005, 2.23606797749979]
 [1.0000000000605, 1.0000550000000001, 2.23606797749979]
 [1.0000000061605, 1.000555, 2.23606797749979]
 [1.0000005, 1.005, 2.23606797749979]
 [1.0000005, 1.005, 2.45051025298814]
 [1.000002, 1.0110090055004999, 2.45051025298814]
 [1.000002, 1.0110090055004999, 2.4517363246283344]
 [1.0000045, 1.0170280325175058, 2.4517363246280173]
 ⋮
 [1.4960079999999851, 16.35754369983386, 4.671293088935029]
 [1.4960079999999851, 16.35754369983386, 4.674778251409642]
 [1.497004499999985, 16.389173203064352, 4.674778251407637]
 [1.497004499999985, 16.389173203064352, 4.6782665275808215]
 [1.498001999999985, 16.420855218475776, 4.678266527578814]
 [1.498001999999985, 16.420855218475776, 4.681757919680873]
 [1.499000499999985, 16.45258983225982, 4.681757919678864]
 [1.499000499999985, 16.45258983225982, 4.685252429942133]
 [1.4999999999999851, 16.484377130747994, 4.68525242994012]