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?
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?
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
end
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.
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â.
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
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?
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? - #3 by asprionj) 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}
x::Array{T,1}
ctrl_x::T # controller state
ctrl_u::T # controller output
end
# 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]
end
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
end
end
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
ctrl_fun(integrator)
step!(integrator,Ts,true)
end
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:
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?
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â.
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:
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:
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.
SoâŚ
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.
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 ) 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.
