search

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

Initialization with derivative condition on an observed variable #2573

Open ChrisRackauckas opened 5 months ago

ChrisRackauckas commented 5 months ago

Take this case which comes from a periodic BC discretization of a PDE:

using ModelingToolkit
using OrdinaryDiffEq

using ModelingToolkit: Differential

using Plots

# Parameters, variables, and derivatives
@parameters t
@variables (u(..))[1:6]
D = Differential(t)

eqs = [0.9(u(t))[6] + D(D((u(t))[2])) ~ 0,
    - 0.9(u(t))[4] + D(D((u(t))[3])) ~ 0,
    - 0.9(u(t))[5] + D(D((u(t))[4])) ~ 0,
    0.9(u(t))[6] + D(D((u(t))[5])) ~ 0,
    - 0.9(u(t))[5] + D(D((u(t))[6])) ~ 0,
    (u(t))[1] ~ (u(t))[6]]

defaults = Dict(u(t) => zeros(6), D(u(t)) => ones(6))

ps = Pair[]

sys = ODESystem(eqs, t, [u(t)], ps, defaults=defaults, name = :sys)

simpsys = structural_simplify(sys)

prob = ODEProblem(simpsys, [], (0, 2pi); build_initializeprob=false)

sol = solve(prob, Tsit5(), saveat=pi / 40)

sol[u(t)]

It works if we skip the initialization problem. But the initialization does something funny. In particular, if you use prob = ODEProblem(simpsys, [D(u[i](t)) => 1.0 for i in 2:6], (0, 2pi)), then you get for the isys:

julia> unknowns(isys)
11-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 var"u(t)[2]ˍt"
 var"u(t)[3]ˍt"
 var"u(t)[4]ˍt"
 var"u(t)[5]ˍt"
 var"u(t)[6]ˍt"
 (u(t))[2]
 (u(t))[3]
 (u(t))[4]
 (u(t))[5]
 (u(t))[6]
 (u(t))[1]

julia> equations(isys)
12-element Vector{Equation}:
 0 ~ 1.0 - var"u(t)[2]ˍt"
 0 ~ 1.0 - var"u(t)[3]ˍt"
 0 ~ 1.0 - var"u(t)[4]ˍt"
 0 ~ 1.0 - var"u(t)[5]ˍt"
 0 ~ 1.0 - var"u(t)[6]ˍt"
 0 ~ -(u(t))[2]
 0 ~ -(u(t))[3]
 0 ~ -(u(t))[4]
 0 ~ -(u(t))[5]
 0 ~ -(u(t))[6]
 0 ~ -(u(t))[1]
 0 ~ -(u(t))[1] + (u(t))[6]

julia> simpisys = structural_simplify(isys; fully_determined = false)
Model sys with 5 equations
Unknowns (0):
Parameters (1):
  t

julia> equations(simpisys)
5-element Vector{Equation}:
 0 ~ 1.0 - var"u(t)[2]ˍt"
 0 ~ 1.0 - var"u(t)[3]ˍt"
 0 ~ 1.0 - var"u(t)[4]ˍt"
 0 ~ 1.0 - var"u(t)[5]ˍt"
 0 ~ 1.0 - var"u(t)[6]ˍt"

julia> unknowns(simpisys)
Any[]

This structural simplification bug is https://github.com/SciML/ModelingToolkit.jl/issues/2515. But ignoring that, if we do what we'd want to do, i.e. 

@parameters t
@variables (u(..))[1:6]
D = Differential(t)

eqs = [0.9(u(t))[6] + D(D((u(t))[2])) ~ 0,
    - 0.9(u(t))[4] + D(D((u(t))[3])) ~ 0,
    - 0.9(u(t))[5] + D(D((u(t))[4])) ~ 0,
    0.9(u(t))[6] + D(D((u(t))[5])) ~ 0,
    - 0.9(u(t))[5] + D(D((u(t))[6])) ~ 0,
    (u(t))[1] ~ (u(t))[6]]

defaults = Dict(vcat([u(t)[i] => 0.0 for i in 1:6], [D(u(t)[i]) => 1.0 for i in 1:6]))

ps = Pair[]

sys = ODESystem(eqs, t, [u(t)[i] for i in 1:6], ps, defaults=defaults, name = :sys)

simpsys = structural_simplify(sys)

prob = ODEProblem(simpsys, [D(u(t)[i]) => 1.0 for i in 1:6], (0, 2pi))

then you get:

julia> prob = ODEProblem(simpsys, [D(u(t)[i]) => 1.0 for i in 1:6], (0, 2pi))
ERROR: Initialization expression Differential(t)((u(t))[1]) is currently not supported. If its a higher order derivative expression, then only the dummy derivative expressions are supported.
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] generate_initializesystem(sys::ODESystem; u0map::Vector{…}, name::Symbol, guesses::Dict{…}, check_defguess::Bool, default_dd_value::Float64, kwargs::@Kwargs{})
    @ ModelingToolkit c:\Users\accou\.julia\dev\ModelingToolkit\src\systems\nonlinear\initializesystem.jl:54
  [3] generate_initializesystem
    @ c:\Users\accou\.julia\dev\ModelingToolkit\src\systems\nonlinear\initializesystem.jl:6 [inlined]
  [4] ModelingToolkit.InitializationProblem{…}(sys::ODESystem, t::Int64, u0map::Vector{…}, parammap::SciMLBase.NullParameters; guesses::Dict{…}, check_length::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit c:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1576
  [5] InitializationProblem
    @ c:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1560 [inlined]
  [6] #_#762
    @ C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1538 [inlined]
  [7] #InitializationProblem#760
    @ C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1526 [inlined]
  [8] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{…}, 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, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), guesses::Dict{…}, t::Int64, warn_initialize_determined::Bool, build_initializeprob::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit c:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:889
  [9] process_DEProblem
    @ c:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:845 [inlined]
 [10] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1048
 [11] ODEProblem
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1038 [inlined]
 [12] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Pair{…}}, tspan::Tuple{Int64, Float64})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1038
 [13] (ODEProblem{true})(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1025
 [14] (ODEProblem{true})(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1024
 [15] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1014
 [16] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1013
 [17] top-level scope
    @ REPL[30]:1
Some type information was truncated. Use `show(err)` to see complete types.

The reason that you get this is because u(t)[1] is eliminated from the expression, and so D(u(t)[1]) is not in the diffmap to know how to make it correspond to a differentiation equation.

ChrisRackauckas commented 5 months ago

https://github.com/SciML/ModelingToolkit.jl/pull/2574 partially fixes this in that you get the right isys, but simplification still fails because of https://github.com/SciML/ModelingToolkit.jl/issues/2515