OrdinaryDiffEq: unexpected results in DEDataArray

Specific Domains Modelling & Simulations
,
1

I am trying to solve a ODE using a DEDataVector to pass and store parameters that come from evaluating functions taking the state vector as argument. After reading the documentation I expected to obtain the value of the fields in SimData at the same times as the solution for the state vector. However, the following example shows that this is not always true:

using StaticArrays
using OrdinaryDiffEq


mutable struct SimData{T, N} <: DEDataVector{T}
    x::SVector{N, T}
    t::T
    α::T
end


foo(x) = cos(x) - 1

function fun(x, p, t)
    α = foo(x[1])
    xd1 = α * x[1] + 0.5 * x[2]
    xd2 = x[1]
    xd = @SVector [xd1, xd2]
    return SimData(xd, t, α)
end

x0_arr = @SVector [1.0, 0.0]
x0 = SimData(x0_arr, 0.0, foo(x0_arr[1]))

prob = ODEProblem{false}(fun, x0, (0.0, 1.0))

Case working as expected:

julia> sol = solve(prob, Tsit5());
julia> [sol[ii].t for ii in 1:length(sol)] - sol.t
9-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

However:

julia> sol = solve(prob, RK4());
julia> [sol[ii].t for ii in 1:length(sol)] - sol.t
13-element Vector{Float64}:
  0.0
 -0.0004995004468281904
 -0.004303836649044169
 -0.011007355182091345
 -0.017836652779592538
 -0.025106222288033136
 -0.03475184792175873
 -0.04537449406429234
 -0.05776887934747188
 -0.07134099643887737
 -0.08652208123574157
 -0.10318679719026902
 -0.04230133645599987

I also tried:

sol = solve(prob, RK4(), saveat=0:0.1:1)
sol = solve(prob, RK4(), dt=0.1)

but I couldn’t obtain the same result as with Tsit5. As it seems to me that this is related to automatic timestep determination, I finally arrived to:

julia> sol = solve(prob, RK4(), dt=0.1, adaptive=false);
julia> [sol[ii].t for ii in 1:length(sol)] - sol.t
11-element Vector{Float64}:
  0.0
 -0.05
 -0.04999999999999999
 -0.050000000000000044
 -0.04999999999999999
 -0.04999999999999999
 -0.04999999999999993
 -0.04999999999999993
 -0.04999999999999993
 -0.04999999999999993
 -0.050000000000000044

which makes sense to me (t at which f is called is the one from the previous solution point).

At this point, my question is: is there a way to obtain a result equivalent to Tsit5() example where sol[ii].t is evaluated at the same time as sol.u[ii]?

2

I don’t have a solution to the DEDataArray problem, but if you just want to log intermediate variables like α, I’ve got a new package called SimulationLogs.jl that might be worth checking out. Here’s how I’d handle this with SimulationLogs:

using StaticArrays
using OrdinaryDiffEq
using SimulationLogs

foo(x) = cos(x) - 1

function fun(x, p, t)
    @log α = foo(x[1])
    xd1 = α * x[1] + 0.5 * x[2]
    xd2 = x[1]
    return  @SVector [xd1, xd2]
end

x0_arr = @SVector [1.0, 0.0]

prob = ODEProblem{false}(fun, x0_arr, (0.0, 1.0))

sol_tsit5 = solve(prob, Tsit5());
sol_rk4 = solve(prob, RK4());

log_tsit5 = get_log(sol_tsit5)
log_rk4 = get_log(sol_rk4)

julia> log_tsit5.α - foo.(x[1] for x in sol_tsit5.u)
9-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

julia> log_rk4.α - foo.(x[1] for x in sol_rk4.u)
13-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

If you actually want to log the t value here, you can also put a

@log t t

somewhere in your fun function (although there really isn’t any point because SimulationLogs works by going back and calculating things for each time point of the solution or supplied time vector). The one caveat to all of this is that your parameters p can’t change over the course of the simulation or your results might be incorrect.

1 Like
3

Thanks for the answer, @jonniedie. I didn’t know SimulationLogs.jl. It seems quite a useful package! Yes, the caveat would be a problem for me. Maybe I didn’t choose the most representative MWE, but I need to deal with changing parameters

4

Ah, gotcha. Then DEDataArrays are probably the right choice here. FWIW, I am working on ways to handle parameters that are changed in callbacks (that’s really the only good way to handle changing parameters anyway). This would be nice to have for, say, discrete time signals in a control system simulation. It should be doable, I just haven’t had the time to get to it yet.

5

You’re relying on behavior that cannot and should not be guaranteed. Basically, you’re relying on the assumption that the final f call of the Runge-Kutta method is at t + dt. That is not the case in many methods (for example, it’s not required in explicit or implicit Runge-Kutta methods). So it’s returning the t that corresponds to that final f call, which then is the “bug”.

Why is RK4 of all algorithms doing that? Because, if you look at the documentation, this adaptive RK4 is:

  • RK4 - The canonical Runge-Kutta Order 4 method. Uses a defect control for adaptive stepping using maximum error over the whole interval.

Defect control is done via:

Reference: http://faculty.smu.edu/shampine/residuals.pdf

You can see two f calls in there. This algorithm has many nice properties of course since it can bound the interval error, and adaptive=false then returns you to the simple RK4 that most people know. But that right there is the reason for why the t is not what you think it is.

In general, this kind of thing is the reason why DEDataArray is being deprecated, because relying on hidden values implicit in the array type also relies on other assumptions about the algorithms which are not necessarily correct. Instead, for all of these use cases with observables and such, ModelingToolkit.jl is a much better solution that we are now focusing on.

https://mtk.sciml.ai/dev/

1 Like
6

Thanks Chris for such a detailed answer. Given that DEDataArrays are going to be deprecated I will look for another solution to save the problem parameters. Do you think that a ContinuousCallback might be useful for the purpose of saving those parameters at the same times where the solution is obtained?
I already had a look at ModelingToolkit.jl, but found it a little bit overwhelming, probably need to check examples more carefully.

BTW I am trying to use this for a flight mechanics simulator Simulation results by AlexS12 · Pull Request #39 · AlexS12/FlightMechanicsSimulator.jl · GitHub
You encouraged me to use OrdinaryDiffEq.jl some time ago

https://github.com/AlexS12/FlightMechanics.jl/issues/46

but I was having some validation problems and I had already messed the project, so I started from scratch migrating the original Fortran code and introducing more complex physical models step by step now that I know a little bit more of Julia coding

7

We probably need better examples in MTK.

DiscreteCallback. But MTK would be cheaper by simply reconstructing observables lazily.

ModelingToolkit is just a symbolic framework, so it’s still using DifferentialEquations.jl and OrdinaryDiffEq.jl under the hood. It just gives a higher level description of the problem in a symbolic form, and allows the symbolic compiler to generate the most efficient code for the problem. So you can say, A is this, B is this, and C is this, and C is the only differential equation", and it can automatically eliminate the A and B code, but make sol[A] and sol[B] generate their timeseries on demand. This issue with DEDataArrays here is that it’s tying implementation to the representation, i.e. it’s coupling how you’re numerically implementing the model’s discretization to the way you’re describing the model. This is inherently something we then cannot do much with, because if you say you want A numerically at t=-.67, it has to create it, even if what you really mean is “when I want something like A, make it and plot it, but do whatever is most efficient”.

1 Like
8

@AlexS12 So it turns out your sim is compatible with SimulationLogs. When I said you can’t have changing parameters, I just meant that you can’t actually change values in your p variable within your simulation function–which is something you would only want to do with callbacks anyway. But using lookup tables or function calls within the sim to calculate parameters like local gravity, aero forces/moments, etc. is totally fine.

I went ahead and forked your repo because 1. I’m a sucker for a good 6-DOF model, and 2. I wanted to make sure SimulationLogs worked with it as expected. You can see here where I sprinkled some @logs on variables of interest. I also tested the logged state variables against the simulation output to make sure there wasn’t anything funny happening.

Also, as a bonus, the logging won’t slow your simulation down because it happens lazily. It’s not until you call get_log that it runs and calculates the values.

1 Like
9

Yes! :smiley: at the moment it is compatible because the structs I was storing in DEDataArray only depend on the current state vector, so it is feasible to reconstruct them afterwards. However, that is only the simplest model. For example: it wouldn’t be valid anymore if I wanted to take aircraft consumed fuel into account: fuel flow would be calculated at each time step and and total consumed fuel would be stored in Aircraft. So Aircraft would not be a function of time and current state only. That’s my limitation with SimulationLogs.jl

Anyway the example seems so good! :tada: Thanks for the effort! I think I could use it as a preliminary implementation for saving some interesting quantities as long as it remains valid. Eventually, I could try to use MTK.

I think I will use SimulationLogs.jl also in other projects. If you like the simulator package, any other recommendation/comment is really welcomed. Feel free to open an issue if you want!

10

For this example we’d actually be good too. The only way to keep track of consumed fuel is to make it a state variable and integrate fuel flow rate (unless you know the closed-form solution of the integral, such as a constant flow rate, but even then its better to just make a state variable and handle it numerically for generality’s sake). Then instead of making a DEDataVector that holds your changing mass properties, you just calculate them on the fly as a function of the current fuel mass state variable. If you want, you can then wrap them back up into an Aircraft so you can pass things around more easily.

I know this was really just meant as an example, but the point kinda holds in general. Any value you have that depends continuously on some sort of accumulated previous state is going to need to be included in your state variables for integration. If you have something that needs to be updated at discrete time points (like a digital controller) SimulationLogs can actually now handle updates from DiscreteCallbacks. They’ll have to be in the p struct instead of in a DEDataVector, but that’s a better idea in general anyway. Also ContinousCallbacks are coming soon, I just haven’t got around to them. That way if you want to, say, switch autopilot modes at a certain altitude, that could be handled too.

It looks really nice so far. I’ll play around with it some more. And actually, it would be really nice if there was a modular set of ModelingToolkit subsystems for aerospace similar to MATLAB’s Aerospace Blockset. Maybe this could turn into that if you wanted to go down the MTK route.

1 Like
11

Definitely!!! :heart_eyes: I was thinking that while having a deeper look at MTK yesterday