# Hierarchical Component Models in ModelingToolkit, How to Connect

I am trying to build a hierarchical model of groups of Lorenz oscillators. I am following the tutorial, but instead of solving the pair of coupled oscillators, I’d like to solve a set of many pairs of coupled oscillators, i.e. use the connected system as a component in the next level system.

I am having trouble figuring out how to assemble the pairs into an ODESystem.

``````# I am following the tutorial https://mtk.sciml.ai/stable/tutorials/ode_modeling/
using ModelingToolkit

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
D(y) ~ x*(ρ-z)-y,
D(z) ~ x*y - β*z]

lorenz1 = ODESystem(eqs,name=:lorenz1)
lorenz2 = ODESystem(eqs,name=:lorenz2)

@variables a(t)
@parameters γ
connections = [0 ~ lorenz1.x + lorenz2.y + a*γ]

# Now deviate from tutorial by introducing a 2nd level of connected blocks
# Other issue: How does scoping work? Is it ok to re-use components named lorenz1 and lorenz2, or will the hierarchy be flattened eventually
# connected1 and 2 are pairs of Lorenz oscillators. I'd like to solve multiple pairs together.
connected1 = ODESystem(connections,t,[a],[γ],systems=[lorenz1,lorenz2], name=:connected1)
connected2 = ODESystem(connections,t,[a],[γ],systems=[lorenz1,lorenz2], name=:connected2)

connections2ndLevel = [] # Pairs are not connected
variables2ndLevel = [] # No extra variables
parameters2ndLevel = [] # No extra parameters

totalSystem = ODESystem(connections2ndLevel,t,variables2ndLevel, parameters2ndLevel, systems = [connected1, connected2], name=:totalSystem)
# This does not work:

ERROR: MethodError: no method matching ODESystem(::Array{Any,1}, ::Operation, ::Array{Any,1}, ::Array{Any,1}; systems=ODESystem[ODESystem(Equation[Equation(ModelingToolkit.Constant(0), (lorenz1₊x(t) + lorenz2₊y(t)) + a(t) * γ)], t, Variable[a], Variable[γ], Variable[], Equation[], Base.RefValue{Array{Expression,1}}(Expression[]), Base.RefValue{Any}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), :connected1, ODESystem[ODESystem(Equation[Equation(derivative(x(t), t), σ * (y(t) - x(t))), Equation(derivative(y(t), t), x(t) * (ρ - z(t)) - y(t)), Equation(derivative(z(t), t), x(t) * y(t) - β * z(t))], t, Variable[x, y, z], Variable[β, σ, ρ], Variable[], Equation[], Base.RefValue{Array{Expression,1}}(Expression[]), Base.RefValue{Any}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), :lorenz1, ODESystem[]), ODESystem(Equation[Equation(derivative(x(t), t), σ * (y(t) - x(t))), Equation(derivative(y(t), t), x(t) * (ρ - z(t)) - y(t)), Equation(derivative(z(t), t), x(t) * y(t) - β * z(t))], t, Variable[x, y, z], Variable[β, σ, ρ], Variable[], Equation[], Base.RefValue{Array{Expression,1}}(Expression[]), Base.RefValue{Any}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), :lorenz2, ODESystem[])]), ODESystem(Equation[Equation(ModelingToolkit.Constant(0), (lorenz1₊x(t) + lorenz2₊y(t)) + a(t) * γ)], t, Variable[a], Variable[γ], Variable[], Equation[], Base.RefValue{Array{Expression,1}}(Expression[]), Base.RefValue{Any}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), :connected2, ODESystem[ODESystem(Equation[Equation(derivative(x(t), t), σ * (y(t) - x(t))), Equation(derivative(y(t), t), x(t) * (ρ - z(t)) - y(t)), Equation(derivative(z(t), t), x(t) * y(t) - β * z(t))], t, Variable[x, y, z], Variable[β, σ, ρ], Variable[], Equation[], Base.RefValue{Array{Expression,1}}(Expression[]), Base.RefValue{Any}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), :lorenz1, ODESystem[]), ODESystem(Equation[Equation(derivative(x(t), t), σ * (y(t) - x(t))), Equation(derivative(y(t), t), x(t) * (ρ - z(t)) - y(t)), Equation(derivative(z(t), t), x(t) * y(t) - β * z(t))], t, Variable[x, y, z], Variable[β, σ, ρ], Variable[], Equation[], Base.RefValue{Array{Expression,1}}(Expression[]), Base.RefValue{Any}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), Base.RefValue{Array{Expression,2}}(Array{Expression}(undef,0,0)), :lorenz2, ODESystem[])])], name=:totalSystem)
...
``````

How do I re-use the oscillator pairs as further components?

How does scoping work? Do all of the components eventually form a hierarchical namespace, or are all components flattened? In this case I would need to give each pair an individual name.

``````# I couldn't address individual components either:
connections2ndLevel = [0 ~ connected1.lorenz1.x + connected2.lorenz1.x]
ERROR: MethodError: no method matching ODESystem(::Array{Equation,1}, ::Variable{ModelingToolkit.Parameter{Number}}, ::Array{Variable,1}, ::Array{Variable,1}, ::Base.RefValue{Array{Expression,1}}, ::Base.RefValue{Any}, ::Base.RefValue{Array{Expression,2}}, ::Base.RefValue{Array{Expression,2}}, ::Symbol, ::Array{ODESystem,1})
``````

If nothing is connected then what’s happening at the higher level? Anyways, weirdness aside, the error message says you need to have the first array be `<:Equation`, i.e.

``````connections2ndLevel = Equation[] # Pairs are not connected
variables2ndLevel = [] # No extra variables
parameters2ndLevel = [] # No extra parameters

totalSystem = ODESystem(connections2ndLevel,t,variables2ndLevel, parameters2ndLevel, systems = [connected1, connected2],
name=:totalSystem)
``````

works fine.

This is now fixed in https://github.com/SciML/ModelingToolkit.jl/pull/600 and you example was added as a test.

1 Like

I am trying so solve a system where each coupled system will have a different time-dependent forcing term, otherwise they are the same. The systems are not coupled to each other. I am not yet sure how to best introduce these forcing terms, I could either add the term to each individual oscillator when I define `eqs` (but then all systems look different, even though they just have a different RHS and are the same otherwise), or add RHS ports to `eqs` and then use the connections at the highest level to hook up each system to the proper external function. This is why I left the connections empty for now.

Thank you for the explanation how to interpret the error message, very helpful.

Many thanks for the fix as well. Works!