DiffEqs: Hybrid (continuous+discrete) system / periodic callback

I want to model the combination of a physical (and thus continuous-time) system with a (discrete-time) controller. The controller algorithm is executed, say, every second. The output of the controller (that thus changes every second) is the input to the physical system.

I tried the straightforward way: Model the physical system as an ODE, then run a loop that evaluates the controller, changes parameters of that model (that represent the controller output), then run the simulation for one second. The problem is that the initialisation phase of the solve command takes about 5ms, whereas the simulation of the physical system over 30 minutes simulated time takes around 60ms in total. So if I call solve 1800 times (every simulated second), this already takes 9s.

What would be a good way to implement this? Basically, it would have to be something like tstops every second, but then pause the solver and let the discrete function (the controller algorithm) execute and change parameters and/or some state variables with zero autonomous dynamics. Is a DiscreteCallback the way to go?

1 Like

Yup, tstops + DiscreteCallback is a good way to do it. Or just use the integrator interface.

1 Like

Just found out about PeriodicCallback – seems to be exactly what I’m looking for…?
Edit: yes, that’s exactly it. Can even modify the parameters; awesome.


Hrm, one more thing though: CVODE_BDF works exceptionally well for my problem – it’s like 10x faster than any other solver. But it does not implement the callback interface, right?

It does implement the callback interface.

(If you’d like to help donate this as a test problem in the future to help us keep improving benchmarks, we would appreciate it!)


Here’s a “minimal-as-possible” example:

using DifferentialEquations
using Plots

f(x, p, t) = p .- x # the system ("unit first-order lag")
ctrl_f(e, x) = 0.85(e + 1.2x) # the controller (proportional-integral)
ctrl_x = [0.0] # array to record (discrete) state of controller

function ctrl_cb(int) # the callback function
    e = 1 - int.u[1] # control error
    push!(ctrl_x, ctrl_x[end] + e) # integrate the control error
    int.p = ctrl_f(e, ctrl_x[end]) # update controller output

pcb  = PeriodicCallback(ctrl_cb, 1.0)
prob = ODEProblem(f, [0.0], (0.0, 8.0), 0.0, callback=pcb)
sol  = solve(prob)

plot(sol, label="system")

ctrl_t = 0:(length(ctrl_x)-2)
plot!(0:(length(ctrl_x)-2), ctrl_f.(1 .- sol(ctrl_t, idxs=1).u, ctrl_x[2:end]),
    color=:black, linetype=:steppost, label="controller")


When using CVODE_BDF() as solver, the callback is never invoked (and thus still ctrl_x == [0.0]), and the system stays flat at 0.0…


I cannot donate the original problem, sorry – it comprises quite a lot of critical know-how from my job. I can try to come up with a generic example that has similar “numerical features” as the original one, but that will take some time. Basically, the majority of state variables of the problem stem from 1D approximations of flows, both by “sequence of volumes” and spectral methods. Right now autodiff fails on the built-in stiff solvers, but that’s another point I have to look at.

BTW, today I played around with pretty much every solver. TRBDF2 for low accuracy, and KenCarp4 or Cash4 for higher accuracy work well – all within a factor 2 from CVODE_BDF.

1 Like

Sundials doesn’t initialize its callbacks. That should get changed.

Callback initialization was added to Sundials master @asprionj. I’m just going to finish up its progress bar support and add a release.

Awesome! :+1:

Tested the use-case here with CVODE_BDF() from the latest Sundials Release 3.5.0, it works. Thanks again Chris! Will apply it to the original problem tomorrow and see how performance is compared to the native solvers. (Larger, computationally intensive model, and much longer simulation.)

I’m curious, did you consider Modia, and if so, what made you decide against using it for this application? It sounds like your application is exactly the kind of stuff modia aims to handle, but I have not tried it myself yet, and I would be interested in knowing if you looked elsewhere due to lack of features or documentation etc.?

Yes, I did consider it, but did not invest that much time. Frankly, I newer could warm to Modelica, don’t know the exact reason. For this specific task/project, just to have one more language / environment besides the “main language for the everyday technical computing” was something I did not want, and it is quite different from a “normal programming language”, so not that quick to learn.

Ok, Modia is written in Julia, but it still adheres to Modelica’s philosophy (which, again, is perfectly fine in itself) and provides a single way of doing things (e.g. always using the IDA solver, or PyPlot for plotting). Also, it was (and still is) in a somewhat “experimental” state and does not provide any systematic documentation (AFAIK). All of this is no problem for a new, mostly cummunity-driven project trying out new ideas… but it’s not something you’d use to solve a real-world problem as fast as possible and need to be flexible to adapt to whatever may get in your way. So I decided to just use DiffEqs as base simulation framework and build a minimal “framework” around it, providing exactly (and only) what I needed (see this post).

In that post I also mention that ModelingToolkit.jl seems to be a very promising approach to “unify” such modeling frameworks by providing a “core engine”.

1 Like

Did this now. CVODE_BDF is about 5x slower than Cash4 and KenCarp4 on the original hybrid problem for medium accuracies (tolerance 1e-4 to 1e-6). Seems that either

  • enforcing a stop every second (thus limiting the max. step size) and/or
  • introducing a discontinuity (the step change of the control signal) each second, thus exciting the system’s stiffness,

does no good to CVODE_BDF. It is, however, rather insensitive to increasing the accuracy, i.e. lowering the tolerance, whereas the other two take around 10x longer when going from 1e-6 to 1e-8, and thus again are around 2x slower than CVODE_BDF. So the limiting of the step size seems to be the main show stopper for CVODE_BDF.

Since 1e-4 is already sufficient for my application / model / use-case, I’ll probably stick to said native solvers for hybrid-system simulation.

Yes, that is definitely the case. BDF methods are multistep methods, so it needs X previous steps for the X order method. If you keep having discontinuities, then it keeps going back to implicit Euler and shrinking the time steps to accommodate for the low order. This is demonstrated as an example in one of our papers as to why no single integrator is perfect for all cases.

Makes sense… nice paper, could just skim over it, Fig. 2 exactly illustrates the problem at hand. I’ll give a shot at the Rodas* solvers as well :wink:

BTW, would there be a way to enforce some fixed (low) number of steps per interval between discrete events? (Once it is experimentally tested how many usually suffice.) How large could the performance benefits be?

Tstops. It shouldn’t make it perform better, just worse

Just wanted to try the MWE also with the “clean” solution of subtyping DEDataVector and letting DiffEq do everything (saving the time trajectories, interpolation, etc.). Besides the problems I described here (What does "user_cache: method has not been implemented for the integrator" mean?) and the necessity to loop over the cache, there’s another issue with the interpolation of the “discrete” user fields. Here’s the code for the MWE:

mutable struct CtrlSimTypeScalar{T} <: DEDataVector{T}
    ctrl_x::T # controller state
    ctrl_u::T # controller output

# NOTE: non-in-place formulation does not work; should mention this in docs
function f(dx,x,p,t)
    dx[1] = x.ctrl_u - x[1]

ctrl_f(e, x) = 0.85(e + 1.2x)

x0 = CtrlSimTypeScalar([0.0], 0.0, 0.0)
prob = ODEProblem(f, x0, (0.0, 8.0))

function ctrl_fun(int)
    e = 1 - int.u[1] # control error

    # pre-calculate values to avoid overhead in iteration over cache
    x_new = int.u.ctrl_x + e
    u_new = ctrl_f(e, x_new)

    # iterate over cache...
    for c in full_cache(int)
        c.ctrl_x = x_new
        c.ctrl_u = u_new

integrator = init(prob, Tsit5()) # only Tsit5 seems to work

# take 8 steps with a sampling time of 1s
Ts = 1.0
for i in 1:8

Then, for illustrating what I mean, do some plotting:

# extract solution, plot...
sol = integrator.sol
plot(sol, label="system")

# actual full steps, with steppre line-style ("correct")
plot!(sol.t, [u.ctrl_u for u in sol.u], linetype=:steppre, color=:black, marker=true, label="u[.].ctrl_u")

# interpolated control signal ("wrong")
t_ctrl = 0:0.01:sol.t[end]
plot!(t_ctrl, [sol(t).ctrl_u for t in t_ctrl], color=:red, label="u(t).ctrl_u")


It seems that:

  1. The cache does not contain the current stop, because changing the cache only affects future points. Is this always the case? Is there a closer description of the role and inner workings of the cache?

  2. The interpolation method used by DiffEq for the user fields is “last value”, which in this context is wrong. Can this behaviour be changed by some kwarg or similar? If not, would adding this be a simple addition or would we have to change the internals a lot?

Why bother? Isn’t this just cosmetics? Nope: When calculating intermediate quantities from the overall solution (that is, continuous state and the discrete signals), the correct value for the discrete signals at any given time has to be used. And I’m always against implementing a “core functionality” such as interpolating the solution again “on the outside”.

1 Like

Let me just start by saying, I’m not convinced this array-based approach is a good idea at all, and have been against it for about the past year now. I think we need to do something drastically differently and just remove this from the docs, and hopefully this post will explain why.

First of all,

What you did is not the recommended approach, precisely because of issues with discontinuities like this. What’s recommended is that you use a DiscreteCallback that updates the value. That will save two values at the discontinuity, one pre and one post discontinuity. The reason is because that makes sure the two exact (non-unique) values are known at that point, which ensures that you can always know the left-continuous and right-continuous behavior (assuming you keep the solution array ordered). So if you followed the documentation, you would not have this issue.

However, given that you seem very interested in this topic and possibly want to try and reduce memory, let me dig into why trying to assume that a point in time with a discontinuity has a unique value that can be saved, and also why DEDataArray is somewhat problematic.

The cache is just literally the cache. When you do f(du,u,p,t), what’s du? The cache variable. You’re just updating the values in the future cache variable. Of course, the user both interacts with the cache variables and doesn’t. The user never “explicitly” wants to interact with the cache, but we do hand it to the user as the du. But then when we do the broadcasts in the next step, those values need to be correct or the next u will not have the correct discrete values. That’s why they need to be updated.

It’s basically random due to the broadcast. Let me explain a bit. The interpolation is just calculated here:

Inside of that broadcast expression

@inbounds @.. out = y₀ + dt*(k[1]*b1Θ + k[2]*b2Θ + k[3]*b3Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ)

you have DEDataArrays of both the forward and backward values. The array values are easy: you perform broadcast. But… what do you do with the discrete values?

Let’s unpack that. Notice that there is no in-place broadcast overload, which means that the values of the DEDataArray are updated in-place but the discrete values are unchanged. The values are thus whatever was already put into the DEDataArray. Which comes from:

So technically this issue can be fixed by generating the in-place part with the right-hand side, but… that’s somewhat odd? I guess we can do that if all tests pass.

That description of broadcast gives a full scope of the issue with handling these objects. When you broadcast in-place, which is what all of the internal parts of the algorithm do, you can control what the value of the discrete variables are because they don’t change in broadcast, so you just need to set them correctly before hitting the algorithm which is what the cache loop is doing. However, for out-of-place broadcast, it’s going to just grab it from some array in the broadcast expression (at least broadcast always lowers the same way). However, if you dig in you’ll see:

y₀.ctrl_u = 0.4938506200402129
(k[1]).ctrl_u = 1.2099094768961067

at that point, so this issue of why you have to update the cache directly comes into play here.

So we could have it similar the right-most point, but that means that DEDataArray will work in DiffEq the way you’d like, which we can throw a test on and it should be okay-ish.

That seems a little snarky. I tried tens of explicit RKs and they seem to work. What don’t work are methods that use linear algebra because DEDataArray doesn’t have BLAS overloads. In theory doing linear algebra is just telling BLAS how to lower it, like https://github.com/SciML/LabelledArrays.jl/blob/master/src/larray.jl#L88 . If you’re willing to test it out, try adding that one line for DEDataArray and see if that fixes linear algebra. If it does, send a PR to DiffEqBase.


In conclusion, if you donate some test cases we can probably get this fixed up for you with a really tricky fix as I described there. But again, trying to save one unique point at a discontinuity is just a problematic idea, and DEDataArray will always need a bit of manual intervention because what to do with the discrete variables when doing the array operations is ill-defined. And if you test out linear algebra a bit we can get that overload in there too. But I generally don’t think that DEDataArray is the best way to do explicit algebraic variables, but don’t have a much better idea right now.

1 Like

Thanks a lot for this very detailed explanation. Couldn’t 100% follow it on the first read-through. Didn’t want to be snarky with the “only Tsit5”, it’s just the only one out of the ones I tried that worked. And yes I suspected something like that all the others (I tried) are implicit ones and thus involve linear algebra since my problem seems to be quite stiff…

The thing is, I’m using Julia with DiffEq and some of the other awesome packages mainly at work, to “get problems solved”. Currently, I have to literally squeeze some spare minutes between much necessary work to work on that cool stuff and get/keep it running. Bad thing is that I cannot afford the time to work myself thoroughly into every topic / package / approach involved… which I hate because I’m the kind of guy always wanting to understand at least as much as to select the appropriate way of doing things. So I try to squeeze out some other minutes from my currently very sparse spare time to play around a bit (can’t call this anything else).

I do have a solution for my “in production” models / simulations, and I don’t want you to put in any extra work (on “unfunded mandates”, as you used to call it once :wink:) to fix anything I asked / pointed to here. It’s just wondering what (and what not) and how stuff works.

The approach I currently use “in production” is to just abuse parameters as algebraic variables. First, I did use discrete callbacks to change them (as “actions” of the controller). But in the controller (the discrete part of the “closed-loop model”), I have to keep track of “controller events” (e.g. some actions take effect only after some variable time, but no actual controller code is executed at that point in time) to know what exactly to do at each tstop. That is, I had to keep my list of irregularly occurring events (which also store data) consistent with the stopping-time list of the integrator (expanded on the go using add_tstop!). It did work, but it was very cumbersome to maintain and work on. Also, opening this approach up to be able to also use the continuous part (DEProblem) in another “environment” (e.g. without the controller) would have been quite tedious. So, soon I decided to use the integrator interface, have a little “outer loop” keeping track of the events, and just stepping the integrator from event to event. It works like a charm, is conceptionally simple, and very flexible.

As parameters, I use a named tuple and basically merge to update some of them with almost no overhead. I keep track of all those changes by an external, custom data structure only storing the actual events on each control signal. I also split my physical model automatically in a dynamic and an algebraic part (which of course may overlap) by sorting a digraph and cutting off the irrelevant parts. The dynamic part is executed during simulation, while the algebraic one (to calculate interesting intermediate and pure “output” values) is calculated by interpolating the solution. And now, the problem is that I would have to match my externally recorded “discrete data” with the DE solution (the “continuous data”), to reconstruct the parameter tuple at that time instance, in order to have all the correct signals to calculate said intermediate / output variables.

This approach works for a closed-loop system of “discrete controller” + “continuous physics” since there is no interaction in between controller events. So it’s easy to just run the integrator until the time of the next event, make it stop there, change something, then run the integrator until the next event. I also eliminated any delays by approximating transport effects (by spatial discretisation, barycentric interpolation, Bessel filters). Having something available which is able to handle this all (discrete controller, delays, hierarchic / component-based model, etc.) would be absolutely awesome, thus my interest in the stuff above, ModelingToolkit, etc.

1 Like