# Reinitializing the InitializationProb of an ODEProblem

**URL:** <https://discourse.julialang.org/t/reinitializing-the-initializationprob-of-an-odeproblem/122778>\
**Category:** Modelling & Simulations\
**Created:** [November 18, 2024, 2:39pm UTC](https://discourse.julialang.org/t/reinitializing-the-initializationprob-of-an-odeproblem/122778 "2024-11-18T14:39:20Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [November 18, 2024, 2:39pm UTC](https://discourse.julialang.org/t/reinitializing-the-initializationprob-of-an-odeproblem/122778/1 "2024-11-18T14:39:20Z")

</div>

I am trying to re-initialize an ODEProblem with a new x-position on the last point `pos[1, end]`. A bit like in the example here: [Initialization of ODESystems · ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/dev/tutorials/initialization/#Numerical-Isolation:-InitializationProblem)

```julia
    dt = 0.05
    tol = 1e-6
    tspan = (0.0, se.duration)
    ts = 0:dt:se.duration
    p2 = [-40, 0, -47]
    p1 = [0, 0, 0]
    u0map = vcat(
        [simple_sys.acc[j, i] => 0 for j in 1:3 for i in 1:se.segments],
        [simple_sys.vel[j, i] => 0 for j in 1:3 for i in 1:se.segments+1],
        [simple_sys.pos[j, end] => p2[j] for j in 1:3],
        [simple_sys.pos[j, 1] => p1[j] for j in 1:3],
    )
    guesses = vcat(
        [simple_sys.segment[j, i] => 1.0 for j in 1:3 for i in 1:se.segments],
        [simple_sys.total_force[j, i] => 1.0 for j in 1:3 for i in 1:se.segments+1],
    )

    prob = ODEProblem(simple_sys, u0map, tspan; guesses, fully_determined=true) # how to remake iprob with new parameters
    integ = init(prob, FBDF(autodiff=true); dt, abstol=tol, reltol=tol, saveat=ts)

    prob = remake(prob; u0=[simple_sys.pos[1, end] => 40])
    integ = init(prob, FBDF(autodiff=true); dt, abstol=tol, reltol=tol, saveat=ts)

```

But the line `prob = remake(prob; u0=[simple_sys.pos[1, end] => 40])` is causing the error:

```julia
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/GYV9b/src/variable.jl:546
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/GYV9b/src/variable.jl:546
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/GYV9b/src/variable.jl:546
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/GYV9b/src/variable.jl:546
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/GYV9b/src/variable.jl:546
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/GYV9b/src/variable.jl:546
┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics ~/.julia/packages/Symbolics/GYV9b/src/variable.jl:546
ERROR: LoadError: Found symbolic value -(pos(t))[3, 1] + (pos(t))[3, 2] for variable (segment(t))[3, 1]. You may be missing an initial condition or have cyclic initial conditions. If this is intended, pass `symbolic_u0 = true`. In case the initial conditions are not cyclic but require more substitutions to resolve, increase `substitution_limit`. To report cycles in initial conditions of unknowns/parameters, pass `warn_cyclic_dependency = true`. If the cycles are still not reported, you may need to pass a larger value for `circular_dependency_max_cycle_length` or `circular_dependency_max_cycles`.

Stacktrace:
  [1] better_varmap_to_vars(varmap::Dict{…}, vars::Vector{…}; tofloat::Bool, use_union::Bool, container_type::Type, toterm::Function, promotetoconcrete::Nothing, check::Bool, allow_symbolic::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/89pF2/src/systems/problem_utils.jl:303
  [2] process_SciMLProblem(constructor::Type, sys::NonlinearSystem, u0map::Dict{…}, pmap::Dict{…}; build_initializeprob::Bool, implicit_dae::Bool, t::Nothing, guesses::Dict{…}, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, eval_expression::Bool, eval_module::Module, fully_determined::Bool, check_initialization_units::Bool, tofloat::Bool, use_union::Bool, u0_constructor::typeof(identity), du0map::Nothing, check_length::Bool, symbolic_u0::Bool, warn_cyclic_dependency::Bool, circular_dependency_max_cycle_length::Int64, circular_dependency_max_cycles::Int64, substitution_limit::Int64, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/89pF2/src/systems/problem_utils.jl:610
  [3] (NonlinearProblem{true})(sys::NonlinearSystem, u0map::Dict{Any, Int64}, parammap::Dict{SymbolicUtils.BasicSymbolic{Real}, Float64}; check_length::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/89pF2/src/systems/nonlinear/nonlinearsystem.jl:499
  [4] (NonlinearProblem{true})(sys::NonlinearSystem, u0map::Dict{Any, Int64}, parammap::Dict{SymbolicUtils.BasicSymbolic{Real}, Float64})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/89pF2/src/systems/nonlinear/nonlinearsystem.jl:493
  [5] NonlinearProblem(::NonlinearSystem, ::Dict{Any, Int64}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/89pF2/src/systems/nonlinear/nonlinearsystem.jl:490
  [6] NonlinearProblem(::NonlinearSystem, ::Dict{Any, Int64}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/89pF2/src/systems/nonlinear/nonlinearsystem.jl:489
  [7] ModelingToolkit.InitializationProblem{…}(sys::ODESystem, t::Float64, u0map::Dict{…}, parammap::Dict{…}; guesses::Vector{…}, check_length::Bool, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, fully_determined::Bool, check_units::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/89pF2/src/systems/diffeqs/abstractodesystem.jl:1330
  [8] ModelingToolkit.InitializationProblem{true, SciMLBase.AutoSpecialize}(sys::ODESystem, t::Float64, u0map::Dict{Any, Any}, parammap::Dict{Any, Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/89pF2/src/systems/diffeqs/abstractodesystem.jl:1254
  [9] (ModelingToolkit.InitializationProblem{true})(::ODESystem, ::Float64, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/89pF2/src/systems/diffeqs/abstractodesystem.jl:1232
 [10] InitializationProblem
    @ ~/.julia/packages/ModelingToolkit/89pF2/src/systems/diffeqs/abstractodesystem.jl:1231 [inlined]
 [11] #InitializationProblem#991
    @ ~/.julia/packages/ModelingToolkit/89pF2/src/systems/diffeqs/abstractodesystem.jl:1220 [inlined]
 [12] InitializationProblem
    @ ~/.julia/packages/ModelingToolkit/89pF2/src/systems/diffeqs/abstractodesystem.jl:1219 [inlined]
 [13] remake_initializeprob(sys::ODESystem, odefn::Function, u0::Vector{Pair{Num, Int64}}, t0::Float64, p::Missing)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/89pF2/src/systems/nonlinear/initializesystem.jl:322
 [14] remake_initialization_data
    @ ~/.julia/packages/SciMLBase/gZRU3/src/remake.jl:216 [inlined]
 [15] remake(prob::ODEProblem{…}; f::Missing, u0::Vector{…}, tspan::Missing, p::Missing, kwargs::Missing, interpret_symbolicmap::Bool, build_initializeprob::Bool, use_defaults::Bool, _kwargs::@Kwargs{})
    @ SciMLBase ~/.julia/packages/SciMLBase/gZRU3/src/remake.jl:128
 [16] macro expansion
    @ ./timing.jl:279 [inlined]
 [17] simulate(se::Settings3, simple_sys::ODESystem, p1::Vector{Int64}, p2::Vector{Int64})
    @ Main ~/Code/Tethers.jl/src/Tether_10.jl:185
 [18] main(; p1::Vector{Int64}, p2::Vector{Int64}, fix_p1::Bool, fix_p2::Bool)
    @ Main ~/Code/Tethers.jl/src/Tether_10.jl:228
 [19] top-level scope
    @ ~/Code/Tethers.jl/src/Tether_10.jl:236
 [20] include(fname::String)
    @ Base.MainInclude ./client.jl:494
 [21] top-level scope
    @ REPL[21]:1
in expression starting at /home/bart/Code/Tethers.jl/src/Tether_10.jl:236
Some type information was truncated. Use `show(err)` to see complete types.

```

What could be the issue here? I am on ModelingToolkit v. 9.51.0

---

<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:** [November 18, 2024, 4:58pm UTC](https://discourse.julialang.org/t/reinitializing-the-initializationprob-of-an-odeproblem/122778/2 "2024-11-18T16:58:02Z")

</div>

@cryptic.ax

---

<div class="post-metadata">

**Author:** ![cryptic.ax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cryptic.ax/32/220540_2.png) [@cryptic.ax](https://discourse.julialang.org/u/cryptic.ax)\
**Post date:** [November 18, 2024, 5:05pm UTC](https://discourse.julialang.org/t/reinitializing-the-initializationprob-of-an-odeproblem/122778/3 "2024-11-18T17:05:46Z")

</div>

@Bart_van_de_Lint could you share a runnable example? The code snippet does not include the model definition.

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [November 18, 2024, 6:07pm UTC](https://discourse.julialang.org/t/reinitializing-the-initializationprob-of-an-odeproblem/122778/4 "2024-11-18T18:07:09Z")

</div>

@cryptic.ax You can run this:

```julia
using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Parameters
using ModelingToolkit: t_nounits as t, D_nounits as D

@with_kw mutable struct Settings3 @deftype Float64
    g_earth::Vector{Float64} = [0.0, 0.0, -9.81] # gravitational acceleration [m/s²]
    v_wind_tether::Vector{Float64} = [2, 0.0, 0.0]
    rho = 1.225
    cd_tether = 0.958
    l0 = 70 # initial tether length [m]
    v_ro = 0.3 # reel-out speed [m/s]
    d_tether = 4 # tether diameter [mm]
    rho_tether = 724 # density of Dyneema [kg/m³]
    c_spring = 614600 # unit spring constant [N]
    rel_compression_stiffness = 0.01 # relative compression stiffness [-]
    damping = 473 # unit damping constant [Ns]
    segments::Int64 = 10 # number of tether segments [-]
    α0 = π/10 # initial tether angle [rad]
    duration = 5 # duration of the simulation [s]
    save::Bool = false # save png files in folder video
end

function set_tether_diameter!(se, d; c_spring_4mm = 614600, damping_4mm = 473)
    se.d_tether = d
    se.c_spring = c_spring_4mm * (d/4.0)^2
    se.damping = damping_4mm * (d/4.0)^2
end
                              
function model(se, p1, p2, fix_p1, fix_p2)
    mass_per_meter = se.rho_tether * π * (se.d_tether/2000.0)^2
    @parameters c_spring0=se.c_spring/(se.l0/se.segments) l_seg=se.l0/se.segments
    @parameters rel_compression_stiffness = se.rel_compression_stiffness
    @variables begin 
        pos(t)[1:3, 1:se.segments+1] # = POS0
        vel(t)[1:3, 1:se.segments+1] # = VEL0
        acc(t)[1:3, 1:se.segments+1]
        segment(t)[1:3, 1:se.segments]
        unit_vector(t)[1:3, 1:se.segments]
        l_spring(t), c_spring(t), damping(t), m_tether_particle(t)
        len(t)[1:se.segments]
        rel_vel(t)[1:3, 1:se.segments]
        spring_vel(t)[1:se.segments]
        c_spr(t)[1:se.segments]
        spring_force(t)[1:3, 1:se.segments]
        v_apparent(t)[1:3, 1:se.segments]
        v_app_perp(t)[1:3, 1:se.segments]
        norm_v_app(t)[1:se.segments]
        half_drag_force(t)[1:3, 1:se.segments]
        total_force(t)[1:3, 1:se.segments+1]
    end

    # basic differential equations
    eqs1 = vcat(D.(pos) .~ vel,
                D.(vel) .~ acc)
    eqs2 = vcat(eqs1...)
    # loop over all segments to calculate the spring and drag forces
    for i in 1:se.segments
        eqs = [segment[:, i] ~ pos[:, i+1] - pos[:, i],
               len[i] ~ norm(segment[:, i]),
               unit_vector[:, i] ~ -segment[:, i]/len[i],
               rel_vel[:, i] ~ vel[:, i+1] - vel[:, i],
               spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i],
               c_spr[i] ~ c_spring / (1+rel_compression_stiffness) 
                                     * (rel_compression_stiffness+(len[i] > l_spring)),
               spring_force[:, i] ~ (c_spr[i] * (len[i] - l_spring) 
                                     + damping * spring_vel[i]) * unit_vector[:, i],
               v_apparent[:, i] ~ se.v_wind_tether .- (vel[:, i] + vel[:, i+1])/2,
               v_app_perp[:, i] ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i],
               norm_v_app[i] ~ norm(v_app_perp[:, i]),
               half_drag_force[:, i] ~ 0.25 * se.rho * se.cd_tether * norm_v_app[i] * (len[i]*se.d_tether/1000.0)
                                        * v_app_perp[:, i]]
        eqs2 = vcat(eqs2, reduce(vcat, eqs))
    end
    # loop over all tether particles to apply the forces and calculate the accelerations
    for i in 1:(se.segments+1)
        eqs = []
        if i == se.segments+1
            push!(eqs, total_force[:, i] ~ spring_force[:, i-1] + half_drag_force[:, i-1])
            if isnothing(p2) || ! fix_p2
                push!(eqs, acc[:, i] ~ se.g_earth + total_force[:, i] / (0.5 * m_tether_particle))
            else
                push!(eqs, acc[:, i] ~ zeros(3))
            end
        elseif i == 1
            push!(eqs, total_force[:, i] ~ spring_force[:, i] + half_drag_force[:, i])
            if isnothing(p1) || ! fix_p1
                push!(eqs, acc[:, i] ~ se.g_earth + total_force[:, i] / (0.5 * m_tether_particle))
            else
                push!(eqs, acc[:, i] ~ zeros(3))
            end
        else
            push!(eqs, total_force[:, i] ~ spring_force[:, i-1] - spring_force[:, i] 
                                           + half_drag_force[:, i-1] + half_drag_force[:, i])
            push!(eqs, acc[:, i] ~ se.g_earth + total_force[:, i] / m_tether_particle)
        end
        eqs2 = vcat(eqs2, reduce(vcat, eqs))
    end
    # scalar equations
    eqs = [l_spring ~ (se.l0 + se.v_ro*t) / se.segments,
    c_spring ~ se.c_spring / l_spring,
    m_tether_particle ~ mass_per_meter * l_spring,
    damping ~ se.damping / l_spring]
    eqs2 = vcat(eqs2, reduce(vcat, eqs))  
        
    @mtkbuild sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t)
    sys, pos, vel, len, c_spr
end

function simulate(se, simple_sys, p1, p2)
    global prob, u0map
    dt = 0.05
    tol = 1e-6
    tspan = (0.0, se.duration)
    ts = 0:dt:se.duration
    u0map = [
        [simple_sys.acc[j, i] => 0 for j in 1:3 for i in 1:se.segments]
        [simple_sys.vel[j, i] => 0 for j in 1:3 for i in 1:se.segments+1]
        [simple_sys.pos[j, end] => p2[j] for j in 1:3]
        [simple_sys.pos[j, 1] => p1[j] for j in 1:3]
    ]
    guesses = vcat(
        [simple_sys.segment[j, i] => 1.0 for j in 1:3 for i in 1:se.segments],
        [simple_sys.total_force[j, i] => 1.0 for j in 1:3 for i in 1:se.segments+1],
    )

    prob = ODEProblem(simple_sys, u0map, tspan; guesses, fully_determined=true) # how to remake iprob with new parameters
    integ = init(prob, FBDF(autodiff=true); dt, abstol=tol, reltol=tol, saveat=ts)

    # comment out to make it run
    prob = remake(prob; u0=[simple_sys.pos[1, end] => 40])
    integ = init(prob, FBDF(autodiff=true); dt, abstol=tol, reltol=tol, saveat=ts)

    elapsed_time = @elapsed step!(integ, se.duration, true)
    # elapsed_time = @elapsed sol = step!(integ, se.duration, true)
    integ.sol, elapsed_time
end

function main(; p1=[0,0,0], p2=nothing, fix_p1=true, fix_p2=false)
    global sol, pos, vel, len, c_spr
    se = Settings3()
    set_tether_diameter!(se, se.d_tether) # adapt spring and damping constants to tether diameter
    simple_sys, pos, vel, len, c_spr = model(se, p1, p2, fix_p1, fix_p2)
    sol, elapsed_time = simulate(se, simple_sys, p1, p2)
    if @isdefined __PC
        return sol, pos, vel, simple_sys
    end
    println("Elapsed time: $(elapsed_time) s, speed: $(round(se.duration/elapsed_time)) times real-time")
    println("Number of evaluations per step: ", round(sol.stats.nf/(se.duration/0.02), digits=1))
    sol, pos, vel, simple_sys
end

main(p2=[-40,0,-47], fix_p2=false);

nothing

```

---

<div class="post-metadata">

**Author:** ![cryptic.ax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cryptic.ax/32/220540_2.png) [@cryptic.ax](https://discourse.julialang.org/u/cryptic.ax)\
**Post date:** [November 20, 2024, 6:40am UTC](https://discourse.julialang.org/t/reinitializing-the-initializationprob-of-an-odeproblem/122778/5 "2024-11-20T06:40:39Z")

</div>

So the problem here is that the guesses are not part of the system, so in `remake` while re-creating the initialization problem, MTK doesn’t have access to them and errors. The fix isn’t very difficult.

---

<div class="post-metadata">

**Author:** ![cryptic.ax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cryptic.ax/32/220540_2.png) [@cryptic.ax](https://discourse.julialang.org/u/cryptic.ax)\
**Post date:** [November 20, 2024, 7:27am UTC](https://discourse.julialang.org/t/reinitializing-the-initializationprob-of-an-odeproblem/122778/6 "2024-11-20T07:27:21Z")

</div>

[feat: propagate `ODEProblem` guesses to `remake` by AayushSabharwal · Pull Request #3226 · SciML/ModelingToolkit.jl · GitHub](https://github.com/SciML/ModelingToolkit.jl/pull/3226) fixes this issue

---

<div class="post-metadata">

**Author:** ![Bart\_van\_de\_Lint](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bart_van_de_lint/32/212161_2.png) [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)\
**Post date:** [November 20, 2024, 7:27am UTC](https://discourse.julialang.org/t/reinitializing-the-initializationprob-of-an-odeproblem/122778/7 "2024-11-20T07:27:57Z")

</div>

thanks!