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

Parameters with default values that depend on variables' initial conditions no longer supported #2556

Closed TorkelE closed 5 months ago

TorkelE commented 5 months ago

Needed for e.g. conservation laws in Catalyst.

MWE:

using ModelingToolkit
import ModelingToolkit.t_nounits as t
import ModelingToolkit.D_nounits as D
@variables X(t) Y(t)
@parameters P=Y
eq = D(X) ~ P - X
@mtkbuild osys = ODESystem([eq], t, [X, Y], [P])
u0 = [X => 1.0, Y => 2.0]
ps = []
ODEProblem(osys, u0, (0.0, 1.0), ps)

gives a

RROR: MethodError: Cannot `convert` an object of type SymbolicUtils.BasicSymbolic{Real} to an object of type Float64

Closest candidates are:
  convert(::Type{T}, ::Unitful.Quantity) where T<:Real
   @ Unitful ~/.julia/packages/Unitful/R4J37/src/conversion.jl:139
  convert(::Type{T}, ::Unitful.Level) where T<:Real
   @ Unitful ~/.julia/packages/Unitful/R4J37/src/logarithm.jl:22
  convert(::Type{T}, ::Unitful.Gain) where T<:Real
   @ Unitful ~/.julia/packages/Unitful/R4J37/src/logarithm.jl:62
  ...

Stacktrace:
  [1] symconvert(::Type{Float64}, x::SymbolicUtils.BasicSymbolic{Real})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/parameter_buffer.jl:2
  [2] ModelingToolkit.MTKParameters(sys::ODESystem, p::Vector{Any}; tofloat::Bool, use_union::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/parameter_buffer.jl:86
  [3] MTKParameters
    @ ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/parameter_buffer.jl:13 [inlined]
  [4] 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::Dict{…}, t::Float64, warn_initialize_determined::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/diffeqs/abstractodesystem.jl:889
  [5] process_DEProblem
    @ ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/diffeqs/abstractodesystem.jl:831 [inlined]
  [6] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::Vector{…}; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/diffeqs/abstractodesystem.jl:1033
  [7] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::Vector{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/diffeqs/abstractodesystem.jl:1023
  [8] (ODEProblem{true})(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/diffeqs/abstractodesystem.jl:1010
  [9] (ODEProblem{true})(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/diffeqs/abstractodesystem.jl:1009
 [10] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/diffeqs/abstractodesystem.jl:999
 [11] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Mxj1Q/src/systems/diffeqs/abstractodesystem.jl:998
 [12] top-level scope
    @ ~/Desktop/Julia Playground/Environment - Catalyst test/playground.jl:92
Some type information was truncated. Use `show(err)` to see complete types.

@AayushSabharwal

ChrisRackauckas commented 5 months ago

Was this ever supported? I thought there was always a data flow of parameters to values. You can solve this in a modeling sense by defining a parameter that is then used as the initial condition?

isaacsas commented 5 months ago

Yes, in MTK8 if you define a parameter as a function of states then it is calculated based on their initial values, and Catalyst relied on this to determine the value of conserved constants.

isaacsas commented 5 months ago

Defining a parameter that gives an initial condition wouldn't work, we want to allow someone to specify initial values for all species and then automatically calculate the value of the conserved quantity from it to use in the reduced model as a parameter.

isaacsas commented 5 months ago

Would it make sense to create a decorator to indicate we want the initial value? i.e. p ~ initialvalue(A) + initialvalue(B) for @species A(t) B(t)?