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.42k stars 206 forks source link

Incorrect initialization systems for high order systems #2520

Open YingboMa opened 7 months ago

YingboMa commented 7 months ago
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D
@parameters g L
@variables q₁(t) q₂(t) λ(t) θ(t) [state_priority = 10]

eqs = [D(D(q₁)) ~ -λ * q₁,
    D(D(q₂)) ~ -λ * q₂ - g,
    q₁ ~ L * sin(θ),
    q₂ ~ L * cos(θ)]

@named pend = ODESystem(eqs, t)
sys = structural_simplify(pend)
guesses = [θ => pi / 4, D(θ) => 0.0, λ => 1.0, D(D(θ)) => 1.0]
prob_mtk = ODEProblem(sys, [], (0.0, 10.0), [L => 1, g => 1]; guesses)
#=
ERROR: Initial condition underdefined. Some are missing from the variable map.
Please provide a default (`u0`), initialization equation, or guess
for the following variables:

SymbolicUtils.BasicSymbolic{Real}[q₁(t)]
=#

but we don't even have q₁ in unknowns

julia> unknowns(sys)
4-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 θˍt(t)
 θ(t)
 λ(t)
 θˍtt(t)
ChrisRackauckas commented 7 months ago

You do have it in the unknowns of the initialization problem:

julia> isys = generate_initializesystem(sys)
Model pend with 8 equations
Unknowns (10):
  θˍt(t) [defaults to 0.0]
  θ(t)
⋮
Parameters (3):
  g
  L
⋮

julia> isys = structural_simplify(isys; fully_determined = false)
Model pend with 4 equations
Unknowns (6):
  λ(t)
  θˍtt(t) [defaults to 0.0]
⋮
Parameters (3):
  g
  L
⋮

julia> equations(isys)
4-element Vector{Equation}:
 0 ~ -q₁ˍtt(t) + L*θˍtt(t)*cos(θ(t)) - L*sin(θ(t))*(θˍt(t)^2)
 0 ~ -q₁(t) + L*sin(θ(t))
 0 ~ -q₂ˍtt(t) - L*sin(θ(t))*θˍtt(t) - L*cos(θ(t))*(θˍt(t)^2)
 0 ~ -q₁ˍt(t) + L*cos(θ(t))*θˍt(t)

julia> unknowns(isys)
6-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 λ(t)
 θˍtt(t)
 θ(t)
 θˍt(t)
 q₁(t)
 q₁ˍt(t)

Why it doesn't solve 0 ~ -q₁(t) + L*sin(θ(t)) is a good question though since it does for the determined form:

julia> observed(sys)
6-element Vector{Equation}:
 q₁(t) ~ L*sin(θ(t))
 q₂(t) ~ L*cos(θ(t))
 q₂ˍtt(t) ~ -L*sin(θ(t))*θˍtt(t) - L*cos(θ(t))*(θˍt(t)^2)
 q₁ˍt(t) ~ L*cos(θ(t))*θˍt(t)
 q₂ˍt(t) ~ -L*sin(θ(t))*θˍt(t)
 q₁ˍtt(t) ~ -q₁(t)*λ(t)

But that just means you need to supply a guess. If the tearing can be made eager in the non-determined form that would be nice.