How to implement a real time DAE solver

Specific Domains Modelling & Simulations
,
1

I have the following solver:

differential_vars =  ones(Bool, 36)
prob = DAEProblem(residual!, yd0, y0, tspan, differential_vars=differential_vars)
solver = IDA(linear_solver=:Dense)

@time sol = solve(prob, solver, saveat=0.025, abstol=0.000001, reltol=0.001)

This works fine. But now I want to change the simulation parameters (control inputs) in a fixed interval. Pseudocode:

dt=0.05
t_start = 0.0
t_end = 0.0 + dt
while true
    tspan = (t_start, t_end)
    solve(prob, solver, ...)
    # execute controller
    # update control parameters
    t_start += dt
    t_end += dt
end

How can I achieve that? Is there a way to continue solving a differential equation without completely restarting it?

I am using the package Sundials.

2

Do you really wan to stop integration in each dt step? There is a callback specifically for this: PeriodicCallback. Or, if you want to really stop solver, you can use integrator interface and step into each time step.

3

Ok, I tried:

differential_vars =  ones(Bool, 36)
solver = IDA(linear_solver=:Dense)
dt = 0.5
t_start = 0.0
t_end   = 5*dt
tspan = (t_start, t_end) 

prob = DAEProblem(residual!, yd0, y0, tspan, differential_vars=differential_vars)
integrator = init(prob, solver, abstol=0.000001, reltol=0.001)

for i in 1:5
    step!(integrator, dt, stop_at_tdt=true)
    for (u,t) in tuples(integrator)
        @show u[18], t
    end
end

But the parameter stop_at_tdt is not accepted:

julia> include("src/RTSim.jl")
[728.4833579505422, 734.950422647022, 741.5051811137938, 748.1406855648668, 754.8497626818295, 761.6991795961667]
ERROR: LoadError: MethodError: no method matching step!(::Sundials.IDAIntegrator{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Sundials.Handle{Sundials.IDAMem}, DAESolution{Float64, 2, Vector{MVector{36, Float64}}, Nothing, Nothing, Nothing, Vector{Float64}, DAEProblem{MVector{36, Float64}, MVector{36, Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DAEFunction{true, typeof(residual!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, IDA{:Dense, Nothing, Nothing}, SciMLBase.HermiteInterpolation{Vector{Float64}, Vector{MVector{36, Float64}}, Vector{MVector{36, Float64}}}, DiffEqBase.DEStats}, IDA{:Dense, Nothing, Nothing}, Sundials.var"#104#108"{DAEProblem{MVector{36, Float64}, MVector{36, Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DAEFunction{true, typeof(residual!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, Tuple{Int64}, Tuple{Int64}}, Sundials.FunJac{Sundials.var"#104#108"{DAEProblem{MVector{36, Float64}, MVector{36, Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DAEFunction{true, typeof(residual!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, Tuple{Int64}, Tuple{Int64}}, Nothing, Nothing, SciMLBase.NullParameters, Nothing, Nothing, MVector{36, Float64}, Vector{Float64}, Nothing, Nothing}, Nothing, Sundials.DEOptions{DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, CallbackSet{Tuple{}, Tuple{}}, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)}, Vector{Float64}, Tuple{Int64}, Tuple{Int64}, MVector{36, Float64}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Nothing}, ::Float64; stop_at_tdt=true)
Closest candidates are:
  step!(::SciMLBase.DEIntegrator, ::Any) at /home/ufechner/.julia/packages/SciMLBase/Z1NtH/src/integrator_interface.jl:576 got unsupported keyword argument "stop_at_tdt"
  step!(::SciMLBase.DEIntegrator, ::Any, ::Any) at /home/ufechner/.julia/packages/SciMLBase/Z1NtH/src/integrator_interface.jl:576 got unsupported keyword argument "stop_at_tdt"
  step!(::Sundials.AbstractSundialsIntegrator) at /home/ufechner/.julia/packages/Sundials/vDmSJ/src/common_interface/integrator_types.jl:236 got unsupported keyword argument "stop_at_tdt"
  ...
Stacktrace:
 [1] top-level scope

You can find the complete code here: KiteViewer/RTSim.jl at sim · ufechner7/KiteViewer · GitHub

Any idea?

4

step!(integrator, dt, true), it’s documented as an optional positional argument IIRC.

5

Ok, this fixed the syntax error. But the result is not what I would expect:

julia> include("src/RTSim.jl")
[728.4833579505422, 734.950422647022, 741.5051811137938, 748.1406855648668, 754.8497626818295, 761.6991795961667]
(u[18], t) = (368.74700170881215, 0.6359999999999999)
(u[18], t) = (368.7470017992227, 0.7719999999999998)
(u[18], t) = (368.7470020404195, 1.0439999999999996)
(u[18], t) = (368.74700225092147, 1.3159999999999994)
(u[18], t) = (368.74700290175184, 1.859999999999999)
(u[18], t) = (368.7470033622797, 2.4039999999999986)
(u[18], t) = (368.74700352032835, 2.5)

I want to know the results at 0.5, 1, 1.5, 2.0 and 2.5 seconds simulation time. Is there a way to fix that?

6

Open an issue with the code. That seems like a Sundials thing. I can dig in, but not until tonight.

7

I am trying now:

for i in 1:5
    step!(integrator, dt, true)
    ts = range(i*dt, stop=i*dt, length=1)
    for (u,t) in TimeChoiceIterator(integrator,ts)
        @show u[18],t
    end
end

The result looks promising, but I still need to check if it is correct:

julia> include("src/RTSim.jl")
[728.4833579505422, 734.950422647022, 741.5051811137938, 748.1406855648668, 754.8497626818295, 761.6991795961667]
(u[18], t) = (368.7470014438838, 0.05)
(u[18], t) = (368.7470014812494, 0.1)
(u[18], t) = (368.7470015049025, 0.15000000000000002)
(u[18], t) = (368.7470015230225, 0.2)
(u[18], t) = (368.7470015381994, 0.25)

But if fails for larger loops:

julia> include("src/RTSim.jl")
[-78.38285932922957, -78.38285932922957, -78.38285932922957, -78.38285932936324, -78.38285932856124, -78.38285932989791]
(u[18], t) = (370.05608255364035, 0.1)
(u[18], t) = (370.1557855621726, 0.2)
(u[18], t) = (370.2440497926481, 0.30000000000000004)
(u[18], t) = (370.3201838088508, 0.4)
(u[18], t) = (370.428361513166, 0.5)
(u[18], t) = (370.428361513166, 0.6000000000000001)
(u[18], t) = (370.428361513166, 0.7000000000000001)
(u[18], t) = (370.428361513166, 0.8)

[IDAS ERROR]  IDAGetDky
  Illegal value for t. t = 0.9 is not between tcur - hu = 0.90002 and tcur = 0.900892.

┌ Warning: IDAGetDky failed with error code = 
│   flag = -26
└ @ Sundials ~/.julia/packages/Sundials/vDmSJ/src/simple.jl:20
(u[18], t) = (370.428361513166, 0.9)

[IDAS ERROR]  IDAGetDky
  Illegal value for t. t = 1 is not between tcur - hu = 1.00046 and tcur = 1.00089.

┌ Warning: IDAGetDky failed with error code = 
│   flag = -26
└ @ Sundials ~/.julia/packages/Sundials/vDmSJ/src/simple.jl:20
(u[18], t) = (370.428361513166, 1.0)

julia>
8

You need to be looping over step!(integrator, dt, true). It’s just:

for i in 1:5
    step!(integrator, dt, true)
    @show integrator.u, integrator.t
end
9

I opened an issue:

10

Well, it seams to work now:

for i in 1:round(t_end/dt)
    step!(integrator, dt, true)
    u = integrator.u
    t = integrator.t
    @show u[18], t
end
1 Like
11

Your problem is that iterating the integrator makes it step.

for i in take(integrator,n) end

makes it step n times for example. So notice

for i in 1:5
    step!(integrator, dt, true)
    for (u,t) in TimeChoiceIterator(integrator,ts)
        @show u[18],t
    end
end

your inner loop code is making it step a bunch, which I assume isn’t what you wanted (but correctly steps it to the end)

12

These iterators that not only iterate over a result set, but actually calculate the result step-by-step still confuse me. But thanks for explaining!

Just don’t expect that everybody who uses DifferentialEquations is a Julia expert.

13

Yes, that’s why we have the step! function-based setup instead, since not everyone might be accustomed to using iterator interfaces (though iterators can be quite nice for handling things generically.

1 Like