DifferentialEquations.jl saving p vector solution

I’m currently using Functional Mockup Units (fmus) in my work and have been combining it with DifferentialEquations.jl in simulations. An FMU consists of an input, state and output, where the state comprises all variables that must be integrated. I use p as a buffer to write fmu io signals to and use u to integrate the state variables of the fmu.

A DifferentialEquations.jl simulation does not keep track of p automatically, so I use a SavingCallback from the DiffEqCallback library to do that. The recorded p vector is not always correct though. Why is this, and is there a better way to do this?

MWE (without fmu)

using DifferentialEquations
using DiffEqCallbacks
using Plots

start = 1.0
stop = 100.0
save_period = 1.0
control_period = 10.0

# 1. store the fmu output in p[2] (here output = state, but that is not always so)
# 2. use the input p[1] (and the fmu state u) to calculate derivatives
f(u, p, t) = (p[2] = u; 1e-2.*p[1])

# controller providing fmu input
control_cb = PeriodicCallback(control_period) do integrator
    integrator.p[1] = (0.5 - integrator.p[2])

# save output of the fmu for later use
# note that the output of an fmu is not always a state variable
plog = SavedValues(Float64, Vector{Float64})
save_cb = SavingCallback((_, _, int) -> copy(int.p), plog, saveat=start:save_period:stop)

u0 = 0.0
p0 = [0.0, 0.0] # [input, output], ideally p[2] is always equal to u

sol = solve(ODEProblem(f, u0, (start, stop), p0), callback=CallbackSet(control_cb, save_cb))

plot(sol) # nice continuous solution
scatter!(plog.t, last.(plog.saveval)) # should be the same as sol, but isn't
vline!(sol.t) # end of each integration step

Thanks for the advice

I have a solution coming in like 2 weeks or so, so I’m just not going to give workarounds because I want people to use the real solution when it’s out :sweat_smile: