Structural_simplify not working on tiered model

[Snowd1n]

For context: I am building a model using ModellingToolkit that has layers of connections as each component within the overall ODE will be built at runtime from symbolics equations etc. I have built this model with the minimal requirements and when I go to structural_simply the final ODE it doesn’t work and produces the error I have attached beneath the MWE. This is being created as an issue from a Julia discourse: https://discourse.julialang.org/t/structural-simplify-not-working-on-tiered-model/114428

Here is the MWE:

using ModelingToolkit,DifferentialEquations,Symbolics,Latexify,LaTeXStrings
using ModelingToolkit: t_nounits as t, D_nounits as D 

function dX(;name)
    @parameters c
    @variables x(t) y(t) s(t)

    eqs = D(x) ~ -c*(x-y) + s

    ODESystem(eqs,t;name)
end

function dY(;name)
    @parameters c
    @variables x(t) y(t)

    eqs = D(y) ~ c*(x-y)

    ODESystem(eqs,t;name)
end

function symboliceq()
    @parameters a b
    s = sqrt(4*log(2)/pi)*a*exp(-(t-b))
    return s
end

function x_factory(;name)
    @named DX = dX()
    symbolic=symboliceq()
    connections=[DX.s~symbolic]
    compose(ODESystem(connections,t;name),DX)
end

function main()
    @named x_eq = x_factory()
    @named y_eq = dY()
    connections = [x_eq.DX.y ~ y_eq.y,
                   y_eq.x ~ x_eq.DX.x,
                   x_eq.DX.c ~ y_eq.c]
    connected = compose(ODESystem(connections,t,name=:connected),x_eq,y_eq)
    connected_simp = structural_simplify(connected)
end

main()

And the error message

ERROR: BoundsError: attempt to access 7-element Vector{Vector{Int64}} at index [8] Stacktrace: [1] getindex @ .\essentials.jl:13 [inlined] [2] 𝑑neighbors @ C:\Users\u5522838\.julia\packages\ModelingToolkit\w72sG\src\bipartite_graph.jl:369 [inlined] [3] 𝑑neighbors @ C:\Users\u5522838\.julia\packages\ModelingToolkit\w72sG\src\bipartite_graph.jl:368 [inlined] [4] check_consistency(state::TearingState{ODESystem}, orig_inputs::Set{Any}) @ ModelingToolkit.StructuralTransformations C:\Users\u5522838\.julia\packages\ModelingToolkit\w72sG\src\structural_transformation\utils.jl:96 [5] _structural_simplify!(state::TearingState{…}, io::Nothing; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, dummy_derivative::Bool, kwargs::@Kwargs{}) @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\w72sG\src\systems\systemstructure.jl:689 [6] structural_simplify!(state::TearingState{…}, io::Nothing; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{}) @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\w72sG\src\systems\systemstructure.jl:635 [7] __structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, kwargs::@Kwargs{}) @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\w72sG\src\systems\systems.jl:74 [8] __structural_simplify @ C:\Users\u5522838\.julia\packages\ModelingToolkit\w72sG\src\systems\systems.jl:55 [inlined] [9] structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, split::Bool, kwargs::@Kwargs{}) @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\w72sG\src\systems\systems.jl:22 [10] structural_simplify (repeats 2 times) @ C:\Users\u5522838\.julia\packages\ModelingToolkit\w72sG\src\systems\systems.jl:19 [inlined] [11] main() @ Main .\Untitled-1:42 [12] top-level scope @ Untitled-1:45 Some type information was truncated. Useshow(err)to see complete types.

[hersle]

The problem here is just the equation x_eq.DX.c ~ y_eq.c. This creates an equation relating two parameters, which makes no sense to MTK inside an ODE because parameters are always specified numbers and thus known. If you want one parameter to take the value of another parameter (by default; unless both are specified), you can use defaults. This works:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D 

function dX(;name)
    @parameters c
    @variables x(t) y(t) s(t)

    eqs = D(x) ~ -c*(x-y) + s

    ODESystem(eqs,t;name)
end

function dY(;name)
    @parameters c
    @variables x(t) y(t)

    eqs = D(y) ~ c*(x-y)

    ODESystem(eqs,t;name)
end

function symboliceq()
    @parameters a b
    s = sqrt(4*log(2)/pi)*a*exp(-(t-b))
    return s
end

function x_factory(;name)
    @named DX = dX()
    symbolic=symboliceq()
    connections=[DX.s~symbolic]
    compose(ODESystem(connections,t;name),DX)
end

function main()
    @named x_eq = x_factory()
    @named y_eq = dY()
    connections = [x_eq.DX.y ~ y_eq.y,
                   y_eq.x ~ x_eq.DX.x]
    defaults = [x_eq.DX.c => y_eq.c]                                                          # <-- only
    connected = compose(ODESystem(connections,t;name=:connected,defaults=defaults),x_eq,y_eq) # <-- change
    connected_simp = structural_simplify(connected)
end

main()

[ChrisRackauckas]

Thanks for figuring out the root cause. I agree that a proper error message is a good solution here.