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!
Tested the usecase 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 cummunitydriven project trying out new ideas… but it’s not something you’d use to solve a realworld 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 1e4
to 1e6
). 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 1e6
to 1e8
, 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 1e4
is already sufficient for my application / model / usecase, I’ll probably stick to said native solvers for hybridsystem 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?
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?  #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: noninplace 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
# precalculate 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 linestyle ("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 arraybased 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 (nonunique) values are known at that point, which ensures that you can always know the leftcontinuous and rightcontinuous 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?
https://github.com/SciML/DiffEqBase.jl/blob/master/src/data_array.jl#L148L149
Let’s unpack that. Notice that there is no inplace broadcast overload, which means that the values of the DEDataArray are updated inplace but the discrete values are unchanged. The values are thus whatever was already put into the DEDataArray. Which comes from:
https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/dense/generic_dense.jl#L367L378
So technically this issue can be fixed by generating the inplace part with the righthand 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 inplace, 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 outofplace 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 rightmost 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 okayish.
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 illdefined. 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 readthrough. 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 “closedloop 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 stoppingtime 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 closedloop 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 / componentbased model, etc.) would be absolutely awesome, thus my interest in the stuff above, ModelingToolkit, etc.
Yeah no worries, I thought I’d lay it all out there because you seem fairly interested in DEDataArrays, and so it’s probably good for me to write down somewhere why I think they may not be a long term solution. But I don’t have a better long term solution, and maybe a simple change is all that’s needed to fix 99.99% of cases, so I’ll probably play around with those ideas I proposed tonight. I really wanted this text somewhere because after I do this PR, I assume that I’ll need to have a post to point to in order to explain to future contributors why it’s so crucial that interpolation uses similar
on the right and not the left end point, even though it overwrites all of the values. Writing generic code is fun
You’re too strict on your typing. Always be careful with types if you want AD to work. Here, AD only occurs on continuous variables, so the types of continuous and discrete variables should be separated.
mutable struct CtrlSimTypeScalar{T,T2} <: DEDataVector{T}
x::Vector{T}
ctrl_x::T2 # controller state
ctrl_u::T2 # controller output
end
with this, your test works now, and it also works on outofplace and implicit methods:
function f(x,p,t)
x.ctrl_u . x
end
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{false}(f, x0, (0.0, 8.0))
function ctrl_fun(int)
e = 1  int.u[1] # control error
# precalculate 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...
if DiffEqBase.isinplace(int.sol.prob)
for c in full_cache(int)
c.ctrl_x = x_new
c.ctrl_u = u_new
end
end
end
integrator = init(prob, Tsit5())
# take 8 steps with a sampling time of 1s
Ts = 1.0
for i in 1:8
ctrl_fun(integrator)
DiffEqBase.step!(integrator,Ts,true)
end
sol = integrator.sol
@test [u.ctrl_u for u in sol.u[2:end]] == [sol(t).ctrl_u for t in sol.t[2:end]]
prob = ODEProblem{true}(f, x0, (0.0, 8.0))
integrator = init(prob, Rosenbrock23())
# take 8 steps with a sampling time of 1s
Ts = 1.0
for i in 1:8
ctrl_fun(integrator)
DiffEqBase.step!(integrator,Ts,true)
end
sol = integrator.sol
@test [u.ctrl_u for u in sol.u[2:end]] == [sol(t).ctrl_u for t in sol.t[2:end]]
with PRs:
https://github.com/SciML/DiffEqBase.jl/pull/481
https://github.com/SciML/OrdinaryDiffEq.jl/pull/1110
Awesome(ly quick), thanks!
The outofplace example doesn’t actually do anything  everything stays at 0.0 and ctrl_fun
has no effect on the integrator.
which example?
The outofplace test in the post to which I replied:
mutable struct CtrlSimTypeScalar{T,T2} <: DEDataVector{T}
x::Vector{T}
ctrl_x::T2 # controller state
ctrl_u::T2 # controller output
end
function f(x,p,t)
x.ctrl_u . x
end
ctrl_f(e, x) = 0.85(e + 1.2x)
x0 = CtrlSimTypeScalar([0.0], 0.0, 0.0)
prob = ODEProblem{false}(f, x0, (0.0, 8.0))
function ctrl_fun(int)
e = 1  int.u[1] # control error
# precalculate 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...
if DiffEqBase.isinplace(int.sol.prob)
for c in full_cache(int)
c.ctrl_x = x_new
c.ctrl_u = u_new
end
end
end
integrator = init(prob, Tsit5())
# take 8 steps with a sampling time of 1s
Ts = 1.0
for i in 1:8
ctrl_fun(integrator)
DiffEqBase.step!(integrator,Ts,true)
end
sol = integrator.sol
# nothing happened
getfield.(sol.u,:x) == [[0.0] for i in 1:length(sol)]
I know it’s a test for the interpolation of the discrete variables, but what would make this a fully functional example?
Adding the else
statement to ctrl_fun almost makes it work
function ctrl_fun(int)
e = 1  int.u[1] # control error
# precalculate 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...
if DiffEqBase.isinplace(int.sol.prob)
for c in full_cache(int)
c.ctrl_x = x_new
c.ctrl_u = u_new
end
else
int.u.ctrl_x = x_new
int.u.ctrl_u = u_new
u_modified!(int,true)
end
end
but the effect is delayed until 1 second.