# How to reinit! an integrator in a non-allocating and consistent way?
**URL:** https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771
**Category:** Modelling & Simulations
**Tags:** modelingtoolkit, ordinarydiffeq
**Created:** [April 6, 2025, 7:35pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771 "2025-04-06T19:35:33Z")
**Posts on this page:** 11
**Page:** 1
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [April 6, 2025, 7:35pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/1 "2025-04-06T19:35:33Z")
I am trying to create a control function from a ModelingToolkit.jl model:
f(x, u, p) → x\_next
I want to do that by re-initializing the [integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#SciMLBase.reinit!) with the current state x and input u, and then stepping the integrator:
```julia
function f!(xnext, x, u, _, p) # inplace version of f(x, u, p)
(integrator, set_x, set_u, _) = p
set_x(integrator, x)
set_u(integrator, u)
step!(integrator)
xnext .= integrator.u # integrator.u is the integrator state, called x in the function
nothing
end
```
How can I reinitialize the integrator in a non-allocating way? I am using `setu` and `setp` now, but there seems to be some caches in the integrator that are not being reset, as `xnext` is not consistent when calling the control function multiple times with the same inputs:
```julia
xnext = [1.000003514961648, 0.9999867038583831]
xnext = [1.0000016858995648, 0.9999936227203033]
xnext = [1.0000168585118323, 0.9999362274452563]
xnext = [1.0001685367381081, 0.999362298719702]
xnext = [1.0016805307279277, 0.9936254586572106]
```
And using `reinit!(integrator)` works, but that uses over 1000 allocations, and performance is critical for this function.
Here is the MWE I created:
```julia
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: D_nounits as D, t_nounits as t, setu, setp, getu
Ts = 0.1
@mtkmodel Pendulum begin
@parameters begin
g = 9.8
L = 0.4
K = 1.2
m = 0.3
τ = 0.0 # input
end
@variables begin
θ(t) = 0.0 # state
ω(t) = 0.0 # state
y(t) # output
end
@equations begin
D(θ) ~ ω
D(ω) ~ -g/L*sin(θ) - K/m*ω + τ/m/L^2
y ~ θ * 180 / π
end
end
@mtkbuild mtk_model = Pendulum()
prob = ODEProblem(mtk_model, nothing, (0.0, Ts))
integrator = OrdinaryDiffEq.init(prob, Tsit5(); dt=Ts, abstol=1e-8, reltol=1e-8, save_on=false, save_everystep=false)
set_x = setu(mtk_model, unknowns(mtk_model))
set_u = setp(mtk_model, [mtk_model.τ])
get_h = getu(mtk_model, [mtk_model.y])
p = (integrator, set_x, set_u, get_h)
function f!(xnext, x, u, _, p)
(integrator, set_x, set_u, _) = p
set_x(integrator, x)
set_u(integrator, u)
step!(integrator)
xnext .= integrator.u # integrator.u is the integrator state, called x in the function
nothing
end
for i in 1:10
local xnext = zeros(2)
f!(xnext, ones(2), 1.0, nothing, p)
@show xnext
end
for i in 1:10
local xnext = zeros(2)
f!(xnext, zeros(2), 1.0, nothing, p)
@show xnext
end
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 7, 2025, 3:11pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/2 "2025-04-07T15:11:28Z")
If you add `reset_dt = true` you get what you want?
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [April 7, 2025, 3:22pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/3 "2025-04-07T15:22:59Z")
```julia
julia> @time reinit!(integrator; reset_dt=true)
0.000716 seconds (1.05 k allocations: 30.953 KiB)
```
That still results in 1.05k allocations. `reinit!(integrator)` works, but has a lot of allocations. While `setu` is inconsistent, but isn’t allocating.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 7, 2025, 3:30pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/4 "2025-04-07T15:30:40Z")
Open an issue to look into those allocations. My guess is that it has to do with the initial mapping, so it’s related specifically to the interaction of MTK and reinit!. I presume non-MTK models are non-allocating on this?
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [April 7, 2025, 3:48pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/5 "2025-04-07T15:48:58Z")
This works:
```julia
julia> @time reinit!(integrator; reinit_dae=false)
0.000011 seconds (7 allocations: 240 bytes)
```
But `reinit_dae` is not documented as far as I know. I don’t really know what it does, and if it could be problematic to set it to false.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 7, 2025, 3:56pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/6 "2025-04-07T15:56:53Z")
It is problematic to set it to false because it won’t run the parameter inits and such. But yeah, we can make an issue to make the MTK dae reinit leaner. That would be an MTK issue.
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [April 7, 2025, 4:02pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/7 "2025-04-07T16:02:46Z")
Submitted an issue here:
>
>
> \`reinit!(integrator, u0)\` currently allocates a lot due to \`reinit\_dae\`. It woul…d be great if these allocations could be reduced.
>
> In this MWE, \`reinit\_dae\` causes around 1000 allocations:
> \`\`\`
> using ModelingToolkit, OrdinaryDiffEq
> using ModelingToolkit: D\_nounits as D, t\_nounits as t, setu, setp, getu
>
> Ts = 0.1
> @mtkmodel Pendulum begin
> @parameters begin
> g = 9.8
> L = 0.4
> K = 1.2
> m = 0.3
> τ = 0.0 # input
> end
> @variables begin
> θ(t) = 0.0 # state
> ω(t) = 0.0 # state
> y(t) # output
> end
> @equations begin
> D(θ) ~ ω
> D(ω) ~ -g/L\*sin(θ) - K/m\*ω + τ/m/L^2
> y ~ θ \* 180 / π
> end
> end
>
> @mtkbuild mtk\_model = Pendulum()
> prob = ODEProblem(mtk\_model, nothing, (0.0, Ts))
> integrator = OrdinaryDiffEq.init(prob, Tsit5(); dt=Ts, abstol=1e-8, reltol=1e-8, save\_on=false, save\_everystep=false)
>
> @time reinit!(integrator; reinit\_dae=true) # 1.05k allocs
> @time reinit!(integrator; reinit\_dae=false) # 6 allocs
> \`\`\`
> From \[this\](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771) discussion on discourse.
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [April 8, 2025, 9:32pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/8 "2025-04-08T21:32:52Z")
As I am still having problems with inconsistent results after resetting the integrator (probably because of some internal states being reset inconsistently), I wanted to try using `solve(prob)` instead. But now I am getting a new problem: when using `setu` to update the problem state, this new state is simply ignored by `solve`. Using `setp` to update the parameters works fine. What could be the problem here, and is there a workaround? @ChrisRackauckas
MWE:
```julia
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: D_nounits as D, t_nounits as t, setu, setp, getu, getp
Ts = 0.1
@mtkmodel Pendulum begin
@parameters begin
g = 9.8
L = 0.4
K = 1.2
m = 0.3
τ = 0.0 # input
end
@variables begin
θ(t) = 0.0 # state
ω(t) = 0.0 # state
y(t) # output
end
@equations begin
D(θ) ~ ω
D(ω) ~ -g/L*sin(θ) - K/m*ω + τ/m/L^2
y ~ θ * 180 / π
end
end
@mtkbuild mtk_model = Pendulum()
prob = ODEProblem(mtk_model, nothing, (0.0, Ts))
set_x = setu(mtk_model, unknowns(mtk_model))
get_x = getu(mtk_model, unknowns(mtk_model))
set_u = setp(mtk_model, [mtk_model.τ])
get_u = getp(mtk_model, [mtk_model.τ])
get_h = getu(mtk_model, [mtk_model.y])
p = (prob, set_x, set_u, get_h)
function f!(xnext, x, u, _, p)
(prob, set_x, set_u, _) = p
set_x(prob, x)
set_u(prob, u)
sol = solve(prob, Tsit5(); save_on=false, save_start=false)
xnext .= sol.u[end] # sol.u is the state, called x in the function
nothing
end
xnext = zeros(2)
println("Initial")
f!(xnext, zeros(2), 0.0, nothing, p)
@show xnext
println("Just changing state x - doesn't work, xnext stays at zero")
f!(xnext, ones(2), 0.0, nothing, p)
@show xnext
println("Just changing input u - does work, xnext is non-zero")
f!(xnext, zeros(2), 1.0, nothing, p)
@show xnext
nothing
```
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [April 9, 2025, 2:08pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/9 "2025-04-09T14:08:59Z")
The inconsistent results with the same inputs were due to the adaptive dt in the solver. This issue is solved by resetting the integrator time and dt before stepping:
```julia
function f!(xnext, x, u, _, p)
(integrator, set_x, set_u, _) = p
set_t!(integrator, 0.0)
set_proposed_dt!(integrator, Ts)
set_x(integrator, x)
set_u(integrator, u)
step!(integrator, Ts)
xnext .= integrator.u # sol.u is the state, called x in the function
nothing
end
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 9, 2025, 4:06pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/10 "2025-04-09T16:06:01Z")
> [@Bart\_van\_de\_Lint](#):
>
> when using `setu` to update the problem state, this new state is simply ignored by `solve`. Using `setp` to update the parameters works fine. What could be the problem here, and is there a workaround? @ChrisRackauckas
Open an issue with that in ModelingToolkit.jl. `setu` needs a similar treatment to what we did to `reinit!`, in that for MTK-built models it needs a symbolic hook.
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [April 9, 2025, 6:37pm UTC](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/11 "2025-04-09T18:37:47Z")
>
>
> When using setu to update the problem state, this new state is simply ignored by… solve. Using setp to update the parameters works fine.
>
> MWE:
> \`\`\`julia
> using ModelingToolkit, OrdinaryDiffEq
> using ModelingToolkit: D\_nounits as D, t\_nounits as t, setu, setp, getu, getp
>
> Ts = 0.1
> @mtkmodel Pendulum begin
> @parameters begin
> g = 9.8
> L = 0.4
> K = 1.2
> m = 0.3
> τ = 0.0 # input
> end
> @variables begin
> θ(t) = 0.0 # state
> ω(t) = 0.0 # state
> y(t) # output
> end
> @equations begin
> D(θ) ~ ω
> D(ω) ~ -g/L\*sin(θ) - K/m\*ω + τ/m/L^2
> y ~ θ \* 180 / π
> end
> end
>
> @mtkbuild mtk\_model = Pendulum()
> prob = ODEProblem(mtk\_model, nothing, (0.0, Ts))
> set\_x = setu(mtk\_model, unknowns(mtk\_model))
> get\_x = getu(mtk\_model, unknowns(mtk\_model))
> set\_u = setp(mtk\_model, \[mtk\_model.τ\])
> get\_u = getp(mtk\_model, \[mtk\_model.τ\])
> get\_h = getu(mtk\_model, \[mtk\_model.y\])
> p = (prob, set\_x, set\_u, get\_h)
>
> function f!(xnext, x, u, \_, p)
> (prob, set\_x, set\_u, \_) = p
> set\_x(prob, x)
> set\_u(prob, u)
> sol = solve(prob, Tsit5(); save\_on=false, save\_start=false)
> xnext .= sol.u\[end\] # sol.u is the state, called x in the function
> nothing
> end
>
> xnext = zeros(2)
>
> println("Initial")
> f!(xnext, zeros(2), 0.0, nothing, p)
> @show xnext
>
> println("Just changing state x - doesn't work, xnext stays at zero")
> f!(xnext, ones(2), 0.0, nothing, p)
> @show xnext
>
> println("Just changing input u - does work, xnext is non-zero")
> f!(xnext, zeros(2), 1.0, nothing, p)
> @show xnext
> nothing
> \`\`\`
>
> From this \[discussion\](https://discourse.julialang.org/t/how-to-reinit-an-integrator-in-a-non-allocating-and-consistent-way/127771/8?u=bart\_van\_de\_lint)