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.38k stars 196 forks source link

Better over-determined tearing #2793

Closed YingboMa closed 2 weeks ago

YingboMa commented 3 weeks ago

~I am not convinced that this is fully correct yet.~

I think it's ready.

Fixes https://github.com/SciML/ModelingToolkit.jl/issues/2788

ChrisRackauckas commented 2 weeks ago

Okay I thought this might've had a bug but it handled a difficult case nicely. https://github.com/SciML/ModelingToolkit.jl/issues/2788 has the actually work example:

using ModelingToolkitStandardLibrary.Blocks: Constant, SecondOrder
using ModelingToolkit: t_nounits as t, D_nounits as D
using ModelingToolkit

@component function SecondOrderTest(; name)
    systems = @named begin
        c = Constant(k=1)
        pt2 = SecondOrder(k=1, w=1, d=0.5)
    end
    eqs = [c.output.u ~ pt2.u]

    return ODESystem(eqs, t; systems, name)
end

@mtkbuild model = SecondOrderTest()

Now the case of the user doing "the right thing" works great:

isys = generate_initializesystem(model; u0map = [pt2.xd => 0.0])
sys = structural_simplify(isys; fully_determined = false)
julia> isys = generate_initializesystem(model; u0map = [pt2.xd => 0.0])
Model model with 8 equations
Unknowns (7):
  pt2₊x(t) [defaults to 0.0]: State of SecondOrder filter
  pt2₊xd(t) [defaults to 0.0]: Derivative state of SecondOrder filter
  c₊output₊u(t) [defaults to 0.0]: Inner variable in RealOutput output
  pt2₊y(t) [defaults to pt2₊y_start]: Output of SISO system
⋮
Parameters (8):
  u_start [defaults to 0.0]
  c₊k [defaults to 1]: Constant output value of block
  pt2₊u_start [defaults to 0.0]
  pt2₊y_start [defaults to 0.0]
⋮

julia> sys = structural_simplify(isys; fully_determined = false)
Model model with 1 equations
Unknowns (0):
Parameters (8):
  u_start [defaults to 0.0]
  c₊k [defaults to 1]: Constant output value of block
  pt2₊u_start [defaults to 0.0]
  pt2₊y_start [defaults to 0.0]

It truly has an extra equation because of a default that should've been a guess. But interestingly following case is seen. In this case, the user accidentally has 7 equations for 7 unknowns, but one of the variables is omitted from all equations:

isys = generate_initializesystem(model)
sys = structural_simplify(isys; fully_determined = false)

It simplifies 7 equations with 7 unknowns to 1 equation with 0 unknowns:

julia> isys = generate_initializesystem(model)
Model model with 7 equations
Unknowns (7):
  pt2₊x(t) [defaults to 0.0]: State of SecondOrder filter
  pt2₊xd(t) [defaults to 0.0]: Derivative state of SecondOrder filter
  c₊output₊u(t) [defaults to 0.0]: Inner variable in RealOutput output
  pt2₊y(t) [defaults to pt2₊y_start]: Output of SISO system
  pt2₊u(t) [defaults to u_start]: Input of SISO system
  pt2₊output₊u(t) [defaults to pt2₊y_start]: Inner variable in RealOutput output
  pt2₊input₊u(t) [defaults to pt2₊u_start]: Inner variable in RealInput input
Parameters (8):
  u_start [defaults to 0.0]
  c₊k [defaults to 1]: Constant output value of block
  pt2₊u_start [defaults to 0.0]
  pt2₊y_start [defaults to 0.0]
  pt2₊k [defaults to 1]: Gain
  pt2₊w [defaults to 1]: Bandwidth (angular frequency)
  pt2₊d [defaults to 0.5]: Relative damping
  t

julia> equations(isys)
7-element Vector{Equation}:
 0 ~ pt2₊y_start - pt2₊y(t)
 0 ~ u_start - pt2₊u(t)
 0 ~ c₊k - c₊output₊u(t)
 0 ~ -pt2₊y(t) + pt2₊x(t)
 0 ~ c₊output₊u(t) - pt2₊u(t)
 0 ~ -pt2₊output₊u(t) + pt2₊y(t)
 0 ~ -pt2₊input₊u(t) + pt2₊u(t)

julia> sys = structural_simplify(isys; fully_determined = false)
Model model with 1 equations
Unknowns (0):
Parameters (8):
  u_start [defaults to 0.0]
  c₊k [defaults to 1]: Constant output value of block
  pt2₊u_start [defaults to 0.0]
  pt2₊y_start [defaults to 0.0]
  pt2₊k [defaults to 1]: Gain
  pt2₊w [defaults to 1]: Bandwidth (angular frequency)
  pt2₊d [defaults to 0.5]: Relative damping
  t

julia> equations(sys)
1-element Vector{Equation}:
 0 ~ c₊output₊u(t) - pt2₊u(t)

julia> observed(sys)
6-element Vector{Equation}:
 pt2₊y(t) ~ pt2₊y_start
 pt2₊u(t) ~ u_start
 c₊output₊u(t) ~ c₊k
 pt2₊x(t) ~ pt2₊y(t)
 pt2₊output₊u(t) ~ pt2₊y(t)
 pt2₊input₊u(t) ~ pt2₊u(t)

You see pt2₊xd(t) can match with no equations since it's not in any equations. So this system has an overdetermined set but also doesn't determine pt2₊xd(t) at all. However, structural simplification drops it, presumably because that unknown is not matched to any equation?

But in theory I guess it's not really unknown if it has no constraint equations, so this may be a choice on our end.

How this manifests in the initialization system is the following error:

julia> prob = ODEProblem(model, [], (0, 10))
ERROR: Initialization incomplete. Not all of the state variables of the
DAE system can be determined by the initialization. Missing
variables:

Any[pt2₊xd(t)]

Stacktrace:
  [1] ModelingToolkit.InitializationProblem{…}(sys::ODESystem, t::Int64, u0map::Dict{…}, parammap::SciMLBase.NullParameters; guesses::Dict{…}, check_length::Bool, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1628
  [2] InitializationProblem
    @ C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1602 [inlined]
  [3] #_#830
    @ C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1580 [inlined]
  [4] InitializationProblem
    @ C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1579 [inlined]
  [5] #InitializationProblem#828
    @ C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1568 [inlined]
  [6] InitializationProblem
    @ C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1567 [inlined]
  [7] 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, initialization_eqs::Vector{…}, kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:928
  [8] process_DEProblem
    @ C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:840 [inlined]
  [9] (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:1092
 [10] ODEProblem
    @ C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1082 [inlined]
 [11] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Any}, tspan::Tuple{Int64, Int64})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1082
 [12] (ODEProblem{true})(::ODESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1069
 [13] (ODEProblem{true})(::ODESystem, ::Vector{Any}, ::Vararg{Any})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1068
 [14] ODEProblem(::ODESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1058
 [15] ODEProblem(::ODESystem, ::Vector{Any}, ::Vararg{Any})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1057
 [16] top-level scope
    @ REPL[19]:1
Some type information was truncated. Use `show(err)` to see complete types.

Since pt2₊xd(t) is dropped, it cannot get a value from the observed map and so it needs to error, and it correctly errors that it's simply not determined by anything. Maybe that's the right behavior?

mtiller-jh commented 2 weeks ago

Can you point me to the documentation that explains the semantics of default for a variable (not a parameter)?

I have the following mental model of these terms and I don't understand what role default plays given the existence of these other concepts:

So can some one point me to documentation to where default fits into this?

YingboMa commented 2 weeks ago

Merging because the test failures are unrelated and this is a strict improvement.

ChrisRackauckas commented 2 weeks ago

Default is just an initial condition if it's a number, or an initialization equation to determine variables initial condition if it's an equation. It has nothing to do with guesses.