Getting correct namespacing of subsystems in ModelingToolkit.jl

Specific Domains Modelling & Simulations
1

Hello,

I really like ModelingToolkit, but I am running into some issues I hope someone can help me with.

My usecase is building lots of similar / identical systems, and connecting them in a specific way (where maybe the connections also have dynamics).
I’d like to have a nice namespace for each subsystem.

Take the following example:

using ModelingToolkit, OrdinaryDiffEq

@variables t
D = Differential(t)

function sysA(;name, τ=2)
    states = @variables X(t)
    params = @parameters τ=τ
    eqs = [
        D(X) ~ -X / τ + sin(t)
    ]
    ODESystem(eqs, t, states, params; name=name)
end

function sysB(;name, τ=1)
    states = @variables Y(t) c(t)
    params = @parameters τ=τ
    eqs = [
        D(Y) ~ -Y \ τ + c
    ]
    ODESystem(eqs, t, states, params; name=name)
end

function connection(;name, ω=5)
    states = @variables A(t) C_in(t) C_out(t)
    params = @variables ω=ω
    eqs = [
        D(A) ~ (C_in - C_out) / ω
    ]
    ODESystem(eqs, t, states, params; name=name)
end

function connect_systems(sys_out, sys_in, connection; name)
    connection = ModelingToolkit.rename(connection, Symbol("coupling"))

    eqs = [
        connection.C_in ~ sys_out.X
        connection.C_out ~ (connection.A - sys_in.Y)
        sys_in.c ~ connection.C_out
    ]

   temp_connection = ODESystem(eqs, t; name=:temp) 
   full_connection = extend(connection, temp_connection)

   compose(sys_in, sys_out, full_connection; name=name)
end

@named sys1 = sysA()
@named sys2 = sysB()
@named conn = connection()

@named connected_sys = connect_systems(sys1, sys2, conn)

I am trying to achieve states as:

 sys2₊Y(t)
 sys2₊c(t)
 sys1₊X(t)
 coupling₊C_in(t)
 coupling₊C_out(t)
 coupling₊A(t)

with all these variables connected via the equations:

julia> equations(connected_sys)
6-element Vector{Equation}:
 Differential(t)(sys2₊Y(t)) ~ τ / (-sys2₊Y(t)) + sys2₊c(t)
 Differential(t)(sys1₊X(t)) ~ (-sys1₊X(t)) / sys1₊τ + sin(t)
 coupling₊C_in(t) ~ sys1₊X(t)
 coupling₊C_out(t) ~ coupling₊A(t) - sys2₊Y(t)
 sys2₊c(t) ~ coupling₊C_out(t)
 Differential(t)(coupling₊A(t)) ~ (coupling₊C_in(t) - coupling₊C_out(t)) / coupling₊ω

Instead I get for the states:

 Y(t)
 c(t)
 sys1₊X(t)
 coupling₊coupling₊C_in(t)
 coupling₊sys1₊X(t)
 coupling₊coupling₊C_out(t)
 coupling₊sys2₊Y(t)
 coupling₊coupling₊A(t)
 coupling₊sys2₊c(t)
 coupling₊A(t)
 coupling₊C_in(t)
 coupling₊C_out(t)

So lots of extra states, where the equations are not correct either:

julia> equations(connected_sys)
6-element Vector{Equation}:
 Differential(t)(Y(t)) ~ τ / (-Y(t)) + c(t)
 Differential(t)(sys1₊X(t)) ~ (-sys1₊X(t)) / sys1₊τ + sin(t)
 coupling₊coupling₊C_in(t) ~ coupling₊sys1₊X(t)
 coupling₊coupling₊C_out(t) ~ coupling₊coupling₊A(t) - coupling₊sys2₊Y(t)
 coupling₊sys2₊c(t) ~ coupling₊coupling₊C_out(t)
 Differential(t)(coupling₊A(t)) ~ (coupling₊C_in(t) - coupling₊C_out(t)) / coupling₊ω

Note that e.g. coupling₊sys1₊X(t) is undetermined in this system.

For this toy example, losing the namespace and simply using extend is not a big deal, but since I want to deal with lots of subsystems that loses its convenience very quickly.

I am struggling to achieve this at the moment.
Is there a good way of doing what I am trying to do with ModelingToolkit?

Many thanks in advance!

2

I think I have narrowed down where my understanding of ModelingToolkit is lacking.

Say we take sys2 from the example above, and do:

extension_eq = [
    sys2.c ~ cos(t)
]
@named extension = ODESystem(extension_eq, t)
extended_sys2 = compose(sys2, extension)

This produces:

julia> extended_sys2 = compose(sys2, extension)
Model sys2 with 2 equations
States (3):
  Y(t)
  c(t)
  extension₊sys2₊c(t)
Parameters (1):
  τ [defaults to 1]

julia> equations(extended_sys2)
2-element Vector{Equation}:
 Differential(t)(Y(t)) ~ τ / (-Y(t)) + c(t)
 extension₊sys2₊c(t) ~ cos(t)

and so extending sys2 actually fails.
(Using extend leads to the same result.)
To me this behaviour is counterintuitive.
Is the ‘top-level’ system in ModelingToolkit somehow privileged?
And is there a way of easily circumventing this behaviour?

Edit:

It appears compose is not a commutative function, because this seems to work just fine:

julia> compose(extension, sys2)
Model extension with 2 equations
States (2):
  sys2₊c(t)
  sys2₊Y(t)
Parameters (1):
  sys2₊τ [defaults to 1]

julia> equations(compose(extension, sys2))
2-element Vector{Equation}:
 sys2₊c(t) ~ cos(t)
 Differential(t)(sys2₊Y(t)) ~ sys2₊τ / (-sys2₊Y(t)) + sys2₊c(t)

I find this all slightly puzzling, but that could just be my unfamiliarity with ModelingToolkit (which otherwise I find wonderful).