SciML / ModelingToolkit.jl

An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning and automated transformations of differential equations
https://mtk.sciml.ai/dev/
Other
1.41k stars 204 forks source link

Cannot initialize system with D(x)-like initialization equations #3029

Closed hersle closed 2 weeks ago

hersle commented 2 weeks ago

This should work:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

@variables x(t)
sys = ODESystem([D(D(x)) ~ 0], t; initialization_eqs = [x ~ 0, D(x) ~ 1], name = :sys) |> structural_simplify
prob = ODEProblem(sys, [], (0.0, 1.0), [])

It fails with

ERROR: MethodError: objects of type Nothing are not callable
Stacktrace:
  [1] tearing_reassemble(state::TearingState{…}, var_eq_matching::ModelingToolkit.BipartiteGraphs.Matching{…}, full_var_eq_matching::ModelingToolkit.BipartiteGraphs.Matching{…}; simplify::Bool, mm::ModelingToolkit.SparseMatrixCLIL{…})
    @ ModelingToolkit.StructuralTransformations C:\Users\herma\.julia\dev\ModelingToolkit\src\structural_transformation\symbolics_tearing.jl:457
  [2] tearing_reassemble
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\structural_transformation\symbolics_tearing.jl:229 [inlined]
  [3] tearing(sys::NonlinearSystem, state::TearingState{…}; mm::ModelingToolkit.SparseMatrixCLIL{…}, simplify::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit.StructuralTransformations C:\Users\herma\.julia\dev\ModelingToolkit\src\structural_transformation\symbolics_tearing.jl:637
  [4] _structural_simplify!(state::TearingState{…}, io::Nothing; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, dummy_derivative::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systemstructure.jl:705
  [5] _structural_simplify!
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systemstructure.jl:675 [inlined]
  [6] #structural_simplify!#1198
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systemstructure.jl:668 [inlined]
  [7] __structural_simplify(sys::NonlinearSystem, io::Nothing; simplify::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systems.jl:85
  [8] __structural_simplify
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systems.jl:66 [inlined]
  [9] structural_simplify(sys::NonlinearSystem, io::Nothing; simplify::Bool, split::Bool, allow_symbolic::Bool, allow_parameter::Bool, conservative::Bool,
 fully_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systems.jl:24
 [10] structural_simplify (repeats 2 times)
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systems.jl:20 [inlined]
 [11] ModelingToolkit.InitializationProblem{…}(sys::ODESystem, t::Float64, u0map::Dict{…}, parammap::SciMLBase.NullParameters; guesses::Dict{…}, check_length::Bool, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, fully_determined::Bool, check_units::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1442
 [12] (ModelingToolkit.InitializationProblem{…})(::ODESystem, ::Float64, ::Vararg{…}; kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1404
 [13] InitializationProblem
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1403 [inlined]
 [14] #InitializationProblem#876
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1392 [inlined]
 [15] InitializationProblem
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1391 [inlined]
 [16] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{…}, parammap::Vector{…}; 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 C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:822
 [17] process_DEProblem
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:766 [inlined]
 [18] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::Vector{…}; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, eval_expression::Bool, eval_module::Module, kwargs::@Kwargs{})  
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:992
 [19] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::Vector{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:980
 [20] (ODEProblem{true})(::ODESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:967
 [21] (ODEProblem{true})(::ODESystem, ::Vector{Any}, ::Vararg{Any})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:966
 [22] ODEProblem(::ODESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:956
 [23] ODEProblem(::ODESystem, ::Vector{Any}, ::Vararg{Any})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:955

A workaround is to define an intermediate observed variable for D(x):

@variables x(t) Dx(t)
sys = ODESystem([D(D(x)) ~ 0, Dx ~ D(x)], t; initialization_eqs = [x ~ 0, Dx ~ 1], name = :sys) |> structural_simplify
prob = ODEProblem(sys, [], (0.0, 1.0), [])

But this should not be necessary. It also works with if x and D(x) are provided through defaults in ODESystem or u0 in ODEProblem.

hersle commented 2 weeks ago

The error happens because tearing_reassemble() is called with the initialization system, which is a NonlinearSystem without an associated independent variable, so the differential operator is set to D = nothing. I'm not sure if this is what is causing the bug itself, though.

ChrisRackauckas commented 2 weeks ago

This should be handled by the dummy derivative substitutions here:

https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/nonlinear/initializesystem.jl#L36-L38

I can play with this and see what's going on.

ChrisRackauckas commented 2 weeks ago

Got it

hersle commented 2 weeks ago

Nice, thanks!