ModelingToolkit 9.39 is broken, 9.40.1 works well

Specific Domains Modelling & Simulations
,
1

Steps to reproduce:

git clone https://github.com/ufechner7/Tethers.jl/
cd Tethers.jl
git checkout broken
julia --project

and then in Julia

using Pkg
Pkg.update()
include("src/Tether_01.jl")

This gives the following stack trace:

ERROR: LoadError: MethodError: no method matching iterate(::Nothing)

Closest candidates are:
  iterate(::Combinatorics.Combinations)
   @ Combinatorics ~/.julia/packages/Combinatorics/Udg6X/src/combinations.jl:13
  iterate(::Combinatorics.Combinations, ::Any)
   @ Combinatorics ~/.julia/packages/Combinatorics/Udg6X/src/combinations.jl:13
  iterate(::Combinatorics.CoolLexCombinations)
   @ Combinatorics ~/.julia/packages/Combinatorics/Udg6X/src/combinations.jl:87
  ...

Stacktrace:
  [1] isempty(itr::Nothing)
    @ Base ./essentials.jl:957
  [2] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Nothing, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, eval_module::Module, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), guesses::Dict{…}, t::Float64, warn_initialize_determined::Bool, build_initializeprob::Bool, initialization_eqs::Vector{…}, fully_determined::Bool, check_units::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/NWBfd/src/systems/diffeqs/abstractodesystem.jl:822
  [3] process_DEProblem
    @ ~/.julia/packages/ModelingToolkit/NWBfd/src/systems/diffeqs/abstractodesystem.jl:766 [inlined]
  [4] (ODEProblem{…})(sys::ODESystem, u0map::Nothing, tspan::Tuple{…}, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, eval_expression::Bool, eval_module::Module, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/NWBfd/src/systems/diffeqs/abstractodesystem.jl:995
  [5] ODEProblem
    @ ~/.julia/packages/ModelingToolkit/NWBfd/src/systems/diffeqs/abstractodesystem.jl:983 [inlined]
  [6] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Nothing, tspan::Tuple{Float64, Float64})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/NWBfd/src/systems/diffeqs/abstractodesystem.jl:983
  [7] (ODEProblem{true})(::ODESystem, ::Nothing, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/NWBfd/src/systems/diffeqs/abstractodesystem.jl:970
  [8] (ODEProblem{true})(::ODESystem, ::Nothing, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/NWBfd/src/systems/diffeqs/abstractodesystem.jl:969
  [9] ODEProblem(::ODESystem, ::Nothing, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/NWBfd/src/systems/diffeqs/abstractodesystem.jl:959
 [10] ODEProblem(::ODESystem, ::Nothing, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/NWBfd/src/systems/diffeqs/abstractodesystem.jl:958
 [11] top-level scope
    @ ~/repos/Tethers.jl/src/Tether_01.jl:25
 [12] include(fname::String)
    @ Base.MainInclude ./client.jl:489
 [13] top-level scope
    @ REPL[4]:1
in expression starting at /home/ufechner/repos/Tethers.jl/src/Tether_01.jl:25
Some type information was truncated. Use `show(err)` to see complete types.

The error is reported for this line:

prob = ODEProblem(simple_sys, nothing, (0.0, duration))

All works fine with ModelingToolkit 9.38 or older. It also fails with ModelingToolkit 9.40. Shall I create a bugreport?

2

@cryptic.ax didn’t you just handle something about this?

3

Yes, this looks like it should be fixed on MTK master

1 Like
4

Oh I’ll release, sorry didn’t realized I missed that part.

5

I updated the branch broken of the repo Tethers.jl to use MTK 9.40.1, and I can run the example now. But I get a warning that I did NOT get with the older MTK version:

julia> include("src/Tether_01.jl")
┌ Warning: Initialization system is overdetermined. 1 equations for 0 unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/AS70O/src/systems/diffeqs/abstractodesystem.jl:1465
  3.125912 seconds (4.50 M allocations: 295.207 MiB, 2.24% gc time, 99.99% compilation time)
OK

Any idea why?

My code:

# Example one: Falling mass.
using ModelingToolkit, OrdinaryDiffEq, ControlPlots
using ModelingToolkit: t_nounits as t, D_nounits as D

G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]

# definiting the model
@variables   pos(t)[1:3] = [0.0, 0.0,  0.0]
@variables   vel(t)[1:3] = [0.0, 0.0, 50.0] 
@variables   acc(t)[1:3] = [0.0, 0.0, -9.81] 

eqs = vcat(D(pos) ~ vel,
           D(vel) ~ acc,
           acc    ~ G_EARTH)

@named sys = ODESystem(eqs, t)
simple_sys = structural_simplify(sys)

# running the simulation
duration = 10.0
dt       = 0.02
tol      = 1e-6
ts       = 0:dt:duration

prob = ODEProblem(simple_sys, nothing, (0.0, duration))
@time sol = solve(prob, Rodas5(), dt=dt, abstol=tol, reltol=tol, saveat=ts)

# plotting the result
X = sol.t
POS_Z = stack(sol[pos], dims=1)[:,3]
VEL_Z = stack(sol[vel], dims=1)[:,3]

p = plot(X, POS_Z, VEL_Z; xlabel="time [s]", ylabels=["pos_z [m]", "vel_z [m/s]"], 
         labels=["pos_z [m]", "vel_z [m/s]"], fig="falling mass")
display(p)
6

It’s overdetermined because of the extra initial condition on the acceleration. That’s a computed value.

7

And why was there no warning with MTK 9.38 ?

8

@ChrisRackauckas Thanks a lot for the quick fix!

And the new warning is very useful, I could simplify my code a lot by remove redundant initializations, e.g.:

Old:

    POS0, VEL0, ACC0, SEGMENTS0, UNIT_VECTORS0 = calc_initial_state(se)
    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 pos(t)[1:3, 1:se.segments+1]  = POS0
    @variables vel(t)[1:3, 1:se.segments+1]  = VEL0
    @variables acc(t)[1:3, 1:se.segments+1]  = ACC0
    @variables segment(t)[1:3, 1:se.segments]  = SEGMENTS0
    @variables unit_vector(t)[1:3, 1:se.segments]  = UNIT_VECTORS0
    @variables length(t) = se.l0
    @variables c_spring(t) = c_spring0
    @variables damping(t) = se.damping  / l_seg
    @variables m_tether_particle(t) = mass_per_meter * l_seg
    @variables norm1(t)[1:se.segments] = l_seg * ones(se.segments)
    @variables rel_vel(t)[1:3, 1:se.segments]  = zeros(3, se.segments)
    @variables spring_vel(t)[1:se.segments] = zeros(se.segments)
    @variables c_spr(t)[1:se.segments] = c_spring0 * ones(se.segments)
    @variables spring_force(t)[1:3, 1:se.segments] = zeros(3, se.segments)
    @variables v_apparent(t)[1:3, 1:se.segments] = zeros(3, se.segments)
    @variables v_app_perp(t)[1:3, 1:se.segments] = zeros(3, se.segments)
    @variables norm_v_app(t)[1:se.segments] = ones(se.segments)
    @variables half_drag_force(t)[1:3, 1:se.segments] = zeros(3, se.segments)
    @variables total_force(t)[1:3, 1:se.segments+1] = zeros(3, se.segments+1)

New:

    POS0, VEL0 = calc_initial_state(se)
    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 pos(t)[1:3, 1:se.segments+1]  = POS0
    @variables vel(t)[1:3, 1:se.segments+1]  = VEL0
    @variables acc(t)[1:3, 1:se.segments+1]
    @variables segment(t)[1:3, 1:se.segments]
    @variables unit_vector(t)[1:3, 1:se.segments]
    @variables length(t)
    @variables c_spring(t)
    @variables damping(t)
    @variables m_tether_particle(t)
    @variables norm1(t)[1:se.segments]
    @variables rel_vel(t)[1:3, 1:se.segments]
    @variables spring_vel(t)[1:se.segments]
    @variables c_spr(t)[1:se.segments]
    @variables spring_force(t)[1:3, 1:se.segments]
    @variables v_apparent(t)[1:3, 1:se.segments]
    @variables v_app_perp(t)[1:3, 1:se.segments]
    @variables norm_v_app(t)[1:se.segments]
    @variables half_drag_force(t)[1:3, 1:se.segments]
    @variables total_force(t)[1:3, 1:se.segments+1]
1 Like