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.43k stars 208 forks source link

ODAEProblem's observed function is not always topologically sorted #1726

Closed LeFiLF closed 8 months ago

LeFiLF commented 2 years ago

Hi,

The input/output causality specifiers do not seem to work properly. Or it is rather an equation system simplification issue given the time when the error occurs. The following test model is a serial RC-circuit with variable resistance. Commenting out the input and output specifiers of the RealInput and RealOuput declarations seems to fix it.

The error message is

ERROR: LoadError: UndefVarError: resistor₊R₊u(t) not defined
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/ModelingToolkit/9ppm9/src/structural_transformation/codegen.jl:199 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qulQp/src/code.jl:351 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:129 [inlined]
  [4] macro expansion
    @ ./none:0 [inlined]
  [5] generated_callfunc
    @ ./none:0 [inlined]
  [6] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#3445568178391039249"), Symbol("##arg#15332867572972984708"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x222941a6, 0xb1dff6e4, 0x7f05668a, 0x65a1b59d, 0x16c7ae4e)})(::SubArray{Float64, 1, RecursiveArrayTools.DiffEqArray{Float64, 2, Vector{Vector{Float64}}, Vector{Float64}, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Vector{Float64}}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, false}, ::Vector{Float64}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:117
  [7] (::ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}})(obsvar::Num, u::SubArray{Float64, 1, RecursiveArrayTools.DiffEqArray{Float64, 2, Vector{Vector{Float64}}, Vector{Float64}, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Vector{Float64}}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, false}, p::Vector{Float64}, t::Float64)
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/9ppm9/src/structural_transformation/codegen.jl:328
  [8] _broadcast_getindex_evalf
    @ ./broadcast.jl:670 [inlined]
  [9] _broadcast_getindex
    @ ./broadcast.jl:643 [inlined]
 [10] getindex
    @ ./broadcast.jl:597 [inlined]
 [11] copy
    @ ./broadcast.jl:899 [inlined]
 [12] materialize
    @ ./broadcast.jl:860 [inlined]
 [13] u_n(timeseries::RecursiveArrayTools.DiffEqArray{Float64, 2, Vector{Vector{Float64}}, Vector{Float64}, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Vector{Float64}}, sym::Num, sol::ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#295"), Symbol("##arg#3445568178391039249"), Symbol("##arg#15332867572972984708"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0970482c, 0x4943e995, 0xa92f8bc9, 0x2c0935d6, 0xcb950406)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#295"), Symbol("##arg#3445568178391039249"), Symbol("##arg#15332867572972984708"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0970482c, 0x4943e995, 0xa92f8bc9, 0x2c0935d6, 0xcb950406)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, plott::Vector{Float64}, plot_timeseries::RecursiveArrayTools.DiffEqArray{Float64, 2, Vector{Vector{Float64}}, Vector{Float64}, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Vector{Float64}})
    @ SciMLBase ~/.julia/packages/SciMLBase/YasGG/src/solutions/solution_interface.jl:617
 [14] solplot_vecs_and_labels(dims::Int64, vars::Vector{Tuple}, plot_timeseries::RecursiveArrayTools.DiffEqArray{Float64, 2, Vector{Vector{Float64}}, Vector{Float64}, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Vector{Float64}}, plott::Vector{Float64}, sol::ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#295"), Symbol("##arg#3445568178391039249"), Symbol("##arg#15332867572972984708"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0970482c, 0x4943e995, 0xa92f8bc9, 0x2c0935d6, 0xcb950406)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#295"), Symbol("##arg#3445568178391039249"), Symbol("##arg#15332867572972984708"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0970482c, 0x4943e995, 0xa92f8bc9, 0x2c0935d6, 0xcb950406)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, plot_analytic::Bool, plot_analytic_timeseries::Nothing, strs::Vector{String})
    @ SciMLBase ~/.julia/packages/SciMLBase/YasGG/src/solutions/solution_interface.jl:628
 [15] diffeq_to_arrays(sol::ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#295"), Symbol("##arg#3445568178391039249"), Symbol("##arg#15332867572972984708"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0970482c, 0x4943e995, 0xa92f8bc9, 0x2c0935d6, 0xcb950406)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#295"), Symbol("##arg#3445568178391039249"), Symbol("##arg#15332867572972984708"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0970482c, 0x4943e995, 0xa92f8bc9, 0x2c0935d6, 0xcb950406)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#219"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, plot_analytic::Bool, denseplot::Bool, plotdensity::Int64, tspan::Nothing, axis_safety::Float64, vars::Vector{Num}, int_vars::Vector{Tuple}, tscale::Symbol, strs::Vector{String})
    @ SciMLBase ~/.julia/packages/SciMLBase/YasGG/src/solutions/solution_interface.jl:415
 [16] macro expansion
    @ ~/.julia/packages/SciMLBase/YasGG/src/solutions/solution_interface.jl:191 [inlined]
 [17] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, sol::SciMLBase.AbstractTimeseriesSolution)
    @ SciMLBase ~/.julia/packages/RecipesBase/qpxEX/src/RecipesBase.jl:289
 [18] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/OXGmH/src/user_recipe.jl:36
 [19] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/OXGmH/src/RecipesPipeline.jl:70
 [20] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
    @ Plots ~/.julia/packages/Plots/SkUg1/src/plot.jl:209
 [21] plot(args::Any; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})
    @ Plots ~/.julia/packages/Plots/SkUg1/src/plot.jl:91

The script is

using ModelingToolkit, OrdinaryDiffEq, IfElse, Plots
import ModelingToolkitStandardLibrary.Electrical:
    OnePort,
    Capacitor,
    Voltage,
    Ground
import ModelingToolkitStandardLibrary.Blocks:
    Constant,
    Sine,
    RealInput

@variables t
D = Differential(t)

function VariableResistor(; name, Rmin=0.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    @named R = RealInput()

    pars = @parameters Rmin = Rmin

    eqs = [
        v ~ IfElse.ifelse(
            R.u > Rmin,
            i*R.u,
            i*Rmin,
        )
    ]

    return extend(ODESystem(eqs, t, [], pars; name=name, systems=[R]), oneport)
end

@named resistor = VariableResistor(Rmin=1)
@named r_signal = Sine(frequency=5, amplitude=1000)
@named capacitor = Capacitor(C=1)
@named controlled_voltage = Voltage()
@named voltage_signal = Constant(k=1)
@named ground = Ground()

eqs = [
    connect(controlled_voltage.p, resistor.p)
    connect(resistor.n, capacitor.p)
    connect(resistor.R, r_signal.output)
    connect(capacitor.n, controlled_voltage.n)
    connect(controlled_voltage.V, voltage_signal.output)
    connect(controlled_voltage.n, ground.g)
]

@named os = ODESystem(eqs, t, systems=[
    resistor,
    capacitor,
    controlled_voltage,
    voltage_signal,
    r_signal,
    ground
])
ss = structural_simplify(os)
prob = ODAEProblem(ss, Pair[], (0.0, 10.0))
sol = solve(prob, Tsit5())
plot(sol, vars=[capacitor.v, resistor.i])
savefig("plot.png")

My environment is:

[6e4b80f9] BenchmarkTools v1.3.1
[e2ed5e7c] Bijections v0.1.4
[336ed68f] CSV v0.10.4
[a6e380b2] ControlSystems v1.2.0
[a93c6f00] DataFrames v1.3.4
[abce61dc] Decimals v0.4.1
[0c46a032] DifferentialEquations v7.2.0
[31c24e10] Distributions v0.25.66
[e30172f5] Documenter v0.27.22
[59287772] Formatting v0.4.2
[ff7dd447] FromFile v0.1.5
[86223c79] Graphs v1.7.1
[615f187c] IfElse v0.1.1
[a98d9a8b] Interpolations v0.14.3
[682c06a0] JSON v0.21.3
[98e50ef6] JuliaFormatter v1.0.8
[093fc24a] LightGraphs v1.3.5
[33e6dc65] MKL v0.5.0
[626554b9] MetaGraphs v0.7.1
[961ee093] ModelingToolkit v8.18.4
[16a59e39] ModelingToolkitStandardLibrary v1.5.0
[1dea7af3] OrdinaryDiffEq v6.19.3
[9b87118b] PackageCompiler v2.0.7
[91a5bcdd] Plots v1.31.5
[c3572dad] Sundials v4.9.4
[0c5d862f] Symbolics v4.10.3
[1986cc42] Unitful v1.11.0
[c2297ded] ZMQ v1.2.1
[37e2e46d] LinearAlgebra
[56ddb016] Logging
[9a3f8284] Random
[cf7118a7] UUIDs
YingboMa commented 2 years ago

It's not related to input/outputs. It's just the observed function generation isn't robust. I recommend to use ODEProblem when you encounter issues with the ODAEProblem.

ChrisRackauckas commented 8 months ago

ODAEProblem removed.