# How to implement a real time DAE solver

**URL:** <https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719>\
**Category:** Modelling & Simulations\
**Tags:** question, diffeq\
**Created:** [May 24, 2021, 12:35pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719 "2021-05-24T12:35:52Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 24, 2021, 12:35pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/1 "2021-05-24T12:35:52Z")

</div>

I have the following solver:

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [May 24, 2021, 12:59pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/2 "2021-05-24T12:59:54Z")

</div>

Do you really wan to stop integration in each `dt` step? There is a callback specifically for this: [`PeriodicCallback`](https://diffeq.sciml.ai/latest/features/callback_library/#PeriodicCallback). Or, if you want to really stop solver, you can use [integrator interface](https://diffeq.sciml.ai/latest/basics/integrator/) and step into each time step.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 24, 2021, 1:54pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/3 "2021-05-24T13:54:42Z")

</div>

Ok, I tried:

```julia
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
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](https://github.com/ufechner7/KiteViewer/blob/sim/src/RTSim.jl)

Any idea?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 24, 2021, 1:58pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/4 "2021-05-24T13:58:57Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 24, 2021, 2:01pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/5 "2021-05-24T14:01:51Z")

</div>

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 24, 2021, 2:12pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/6 "2021-05-24T14:12:21Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 24, 2021, 2:15pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/7 "2021-05-24T14:15:34Z")

</div>

I am trying now:

```julia
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
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
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> 

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 24, 2021, 2:33pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/8 "2021-05-24T14:33:50Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 24, 2021, 2:47pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/9 "2021-05-24T14:47:59Z")

</div>

> [@ChrisRackauckas](#):
>
> Open an issue with the code.

I opened an issue:

> <https://github.com/SciML/Sundials.jl/issues/307>
>
> The following code gives an unexpected output:
> 
> \`\`\`Julia
> 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, true)
> for (u,t) in tuples(integrator)
> @show u\[18\], t
> end
> end
> \`\`\`
> Output:
> \`\`\`
> (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 would expect the results at 0.5, 1, 1.5, 2.0 and 2.5 seconds simulation time.
> 
> Full example: https://github.com/ufechner7/KiteViewer/blob/sim/src/RTSim\_bug.jl

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 24, 2021, 2:49pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/10 "2021-05-24T14:49:32Z")

</div>

Well, it seams to work now:

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 24, 2021, 3:17pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/11 "2021-05-24T15:17:46Z")

</div>

Your problem is that iterating the integrator makes it step.

```julia
for i in take(integrator,n) end

```

makes it step `n` times for example. So notice

```julia
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)

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 24, 2021, 3:58pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/12 "2021-05-24T15:58:21Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 24, 2021, 4:13pm UTC](https://discourse.julialang.org/t/how-to-implement-a-real-time-dae-solver/61719/13 "2021-05-24T16:13:27Z")

</div>

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.