Hierarchical Component Models in ModelingToolkit, How to Connect

Specific Domains Modelling & Simulations
1

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})
2

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 fix nested hierarchy by ChrisRackauckas · Pull Request #600 · SciML/ModelingToolkit.jl · GitHub and you example was added as a test.

1 Like
3

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!