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 203 forks source link

Initial conditions set by mapping `unknowns(sys)` do not override defaults #3025

Closed hersle closed 1 week ago

hersle commented 1 week ago

This should work:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
@variables x(t) y(t)

# system 1 should solve to x = 1
ics1 = [x => 1]
sys1 = ODESystem([D(x) ~ 0], t; defaults = ics1, name = :sys1) |> structural_simplify

# system 2 should solve to x = y = 2
sys2 = extend(sys1, ODESystem([D(y) ~ 0], t; initialization_eqs = [y ~ 2], name = :sys2)) |> structural_simplify
ics2 = unknowns(sys1) .=> 2 # but "ics2 = [x => 2]" works
prob2 = ODEProblem(sys2, ics2, (0.0, 1.0), []; fully_determined = true)

In sys2, ics2 should override ics1, so x = 2 instead of 1. But the initialization fails with

ERROR: LoadError: ExtraEquationsSystemException: The system is unbalanced. There are 2 highest order derivative variables and 3 equations.
More equations than variables, here are the potential extra equation(s):
 0 ~ 1 - x(t)
Note that the process of determining extra equations is a best-effort heuristic. The true extra equations are dependent on the model and may not be in this list.
Stacktrace:
  [1] error_reporting(state::TearingState{…}, bad_idxs::Vector{…}, n_highest_vars::Int64, iseqs::Bool, orig_inputs::Set{…})
    @ ModelingToolkit.StructuralTransformations C:\Users\herma\.julia\dev\ModelingToolkit\src\structural_transformation\utils.jl:44
  [2] check_consistency(state::TearingState{NonlinearSystem}, orig_inputs::Set{Any})
    @ ModelingToolkit.StructuralTransformations C:\Users\herma\.julia\dev\ModelingToolkit\src\structural_transformation\utils.jl:85
  [3] _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:692
  [4] _structural_simplify!
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systemstructure.jl:675 [inlined]      
  [5] #structural_simplify!#1197
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systemstructure.jl:668 [inlined]      
  [6] __structural_simplify(sys::NonlinearSystem, io::Nothing; simplify::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systems.jl:85
  [7] __structural_simplify
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systems.jl:66 [inlined]
  [8] 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
  [9] structural_simplify (repeats 2 times)
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\systems.jl:20 [inlined]
 [10] ModelingToolkit.InitializationProblem{…}(sys::ODESystem, t::Float64, u0map::Vector{…}, 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:1445
 [11] (ModelingToolkit.InitializationProblem{…})(::ODESystem, ::Float64, ::Vararg{…}; kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1404
 [12] ModelingToolkit.InitializationProblem(::ODESystem, ::Float64, ::Vararg{…}; kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1392
 [13] 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
 [14] process_DEProblem
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:766 [inlined]
 [15] (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
 [16] (ODEProblem{true})(::ODESystem, ::Vector{Pair{…}}, ::Vararg{Any}; kwargs::@Kwargs{fully_determined::Bool})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:967
 [17] ODEProblem(::ODESystem, ::Vector{Pair{…}}, ::Vararg{Any}; kwargs::@Kwargs{fully_determined::Bool})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:956
hersle commented 1 week ago

The reason is that unknowns(sys) returns unwrapped symbolic variables, which are handled differently than wrapped symbolic variables when setting up the initialization system. This is shown by this simple example:

sys2 = ODESystem([D(x) ~ 0, D(y) ~ 0], t; defaults = [x => 1], initialization_eqs = [y ~ 2], name = :sys2) |> structural_simplify
prob2 = ODEProblem(sys2, [ModelingToolkit.unwrap(x) => 2], (0.0, 1.0), []; fully_determined = true)

It fails with the same error. But it runs if x is not unwrapped.

ChrisRackauckas commented 1 week ago

Interesting. Yeah that's a subtle bug that would need some digging.

ChrisRackauckas commented 1 week ago

Interesting. Yeah that's a subtle bug that would need some digging.