Reinitializing the InitializationProb of an ODEProblem

Specific Domains Modelling & Simulations
1

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

    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:

┌ 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

2

@cryptic.ax

3

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

4

@cryptic.ax You can run this:

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
5

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.

6

feat: propagate `ODEProblem` guesses to `remake` by AayushSabharwal · Pull Request #3226 · SciML/ModelingToolkit.jl · GitHub fixes this issue

1 Like
7

thanks!