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

Unknowns in the `initialization_eqs` should be added to the variable list #2787

Closed ven-k closed 3 months ago

ven-k commented 3 months ago

Describe the bug 🐞

The unknowns in the initialization_eqs aren't automatically added to the variable map/list. ATM, it should additionally be added to the guesses. See the example (from docs) below.

Expected behavior

Unknowns that appear in the lhs should be just added to the list.

Minimal Reproducible Example 👇

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

@parameters g
@variables x(t) y(t) [state_priority = 10] λ(t)
eqs = [D(D(x)) ~ λ * x
       D(D(y)) ~ λ * y - g
       x^2 + y^2 ~ 1]
@mtkbuild pend = ODESystem(eqs, t)

prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = [λ => 0],
    initialization_eqs = [y ~ 1])

(The example is modified from here: https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/ prob in the docs: prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1], initialization_eqs = [y ~ 1]))

Error & Stacktrace ⚠️

┌ Warning: Initialization system is underdetermined. 3 equations for 4 unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:1629
ERROR: ArgumentError: Any[y(0.0)] are either missing from the variable map or missing from the system's unknowns/parameters list.
Stacktrace:
  [1] throw_missingvars_in_sys(vars::Vector{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/utils.jl:621
  [2] promote_to_concrete(vs::Vector{Any}; tofloat::Bool, use_union::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/utils.jl:638
  [3] varmap_to_vars(varmap::Dict{…}, varlist::Vector{…}; defaults::Dict{…}, check::Bool, toterm::Function, promotetoconcrete::Nothing, tofloat::Bool, use_union::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/variables.jl:163
  [4] get_u0(sys::NonlinearSystem, u0map::Dict{…}, parammap::Dict{…}; symbolic_u0::Bool, toterm::Function)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:830
  [5] get_u0
    @ ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:817 [inlined]
  [6] process_NonlinearProblem(constructor::Type, sys::NonlinearSystem, u0map::Dict{…}, parammap::Dict{…}; version::Nothing, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/nonlinear/nonlinearsystem.jl:398
  [7] process_NonlinearProblem
    @ ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/nonlinear/nonlinearsystem.jl:384 [inlined]
  [8] (NonlinearLeastSquaresProblem{…})(sys::NonlinearSystem, u0map::Dict{…}, parammap::Dict{…}; check_length::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/nonlinear/nonlinearsystem.jl:464
  [9] (NonlinearLeastSquaresProblem{…})(sys::NonlinearSystem, u0map::Dict{…}, parammap::Dict{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/nonlinear/nonlinearsystem.jl:458
 [10] NonlinearLeastSquaresProblem(::NonlinearSystem, ::Dict{Num, Int64}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/nonlinear/nonlinearsystem.jl:455
 [11] NonlinearLeastSquaresProblem(::NonlinearSystem, ::Dict{Num, Int64}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/nonlinear/nonlinearsystem.jl:454
 [12] ModelingToolkit.InitializationProblem{…}(sys::ODESystem, t::Float64, u0map::Vector{…}, parammap::Dict{…}; guesses::Vector{…}, check_length::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:1645
 [13] (ModelingToolkit.InitializationProblem{…})(::ODESystem, ::Float64, ::Vararg{…}; kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:1573
 [14] ModelingToolkit.InitializationProblem(::ODESystem, ::Float64, ::Vararg{…}; kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:1561
 [15] 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, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), guesses::Vector{…}, t::Float64, warn_initialize_determined::Bool, build_initializeprob::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:922
 [16] process_DEProblem
    @ ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:835 [inlined]
 [17] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::Vector{…}; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:1086
 [18] (ODEProblem{…})(::ODESystem, ::Vector{…}, ::Vararg{…}; kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:1063
 [19] ODEProblem(::ODESystem, ::Vector{…}, ::Vararg{…}; kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/diffeqs/abstractodesystem.jl:1052
 [20] top-level scope
    @ REPL[12]:1
Some type information was truncated. Use `show(err)` to see complete types

Environment (please complete the following information):

(mtkmsl) pkg> status
Status `~/.julia/dev/mtkmsl/Project.toml`
  [0c46a032] DifferentialEquations v7.13.0
  [961ee093] ModelingToolkit v9.16.0
  [16a59e39] ModelingToolkitStandardLibrary v2.7.1
  [1dea7af3] OrdinaryDiffEq v6.82.0
Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × AMD Ryzen 9 5900HS with Radeon Graphics
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 16 default, 0 interactive, 8 GC (on 16 virtual cores

Additional context

One might have to account for the Differential operator in intializtion_eqs of D(w) ~ 0 sort.

ChrisRackauckas commented 3 months ago

y is missing a guess. The error message needs to be handled but that is trying to throw an error about missing a guess on y which would be required to solve the nonlinear system.

SebastianM-C commented 3 months ago

does y need to be a guess even if it's specified in intializtion_eqs?

ven-k commented 3 months ago

does y need to be a guess even if it's specified in intializtion_eqs?

This is my doubt too. (and I wonder if is possible to avoid that)

SebastianM-C commented 3 months ago

As @ChrisRackauckas mentioned, in this case the initialization equations are not satisfiable (x=>1, y=>1, x^2 + y^2 ~ 1).

ChrisRackauckas commented 3 months ago

Yes, the equations are not satisfiable. Note that it's not (x=>1, y=>1, x^2 + y^2 ~ 1), it's (x=>1, y ~ 1, x^2 + y^2 ~ 1). That's different. initialization_eqs can have any nonlinear function. You cannot assume it's always analytically solvable. There is an issue of asking why this one isn't actually simplified out, which may be related to the unsatisfiabiltiy, but you should expect that if you can put in a nonlinear system then you need a guess.

But also, x = 1, y = 1, x^2 + y^2 != 1, so this test case is erroring for good reasons.

ven-k commented 3 months ago

Ok. That makes sense. Thanks

ChrisRackauckas commented 3 months ago

You also invented a new feature, which is easier to just support. https://github.com/SciML/ModelingToolkit.jl/pull/2789