Unkown connection number for ModelingToolkit system

Hello everyone,

I’ve run into a challenge with using ModelingToolkit, where I have a lots of (similar) systems, each with lots of inputs from a number of other systems connected to it.
Because of the large number of systems and connections involved, I need to generate these somehow procedurally.
Problem is, I cannot specify how many connections each subsystem has manually.

Say I have a system with variable X(t) and I want to add Y(t) to it from another system.
This works:

using ModelingToolkit, OrdinaryDiffEq

@variables t
D = Differential(t)

sys = ODESystem(Equation[], t; name=:System)

function subsystem(;name)
    states = @variables X(t)=0. In(t)=0.
    eqs = [X ~ In]

    ODESystem(eqs, t, states, []; name=name)
end

@named sub = subsystem(;name=:Sub)

sys = compose(sys, sub; name=:System)

function input_factory(;name)
    states = @variables Y(t)
    eq = [
        D(Y) ~ -Y
    ]
    ODESystem(eq, t, states, []; name=name)
end

function add_input(system, i)
    input = input_factory(name=Symbol(:input, i))
    sub = getproperty(sys, :Sub; namespace=false)

    eq = [
        sub.In ~ input.Y
    ]

    @named connection = ODESystem(eq, t)
    connected_input = compose(connection, input; name=Symbol(:input, i))
    compose(connected_input, system.systems; name=system.name)
end

julia> sys = add_input(sys, 1)
Model System with 3 equations
States (3):
  Sub₊In(t) [defaults to 0.0]
  input1₊Y(t)
  Sub₊X(t) [defaults to 0.0]
Parameters (0):

julia> equations(sys)
3-element Vector{Equation}:
 Sub₊In(t) ~ input1₊Y(t)
 Differential(t)(input1₊Y(t)) ~ -input1₊Y(t)
 Sub₊X(t) ~ Sub₊In(t)

But how to handle the case where I’d like to say X(t) = Y_1(t) + Y_2(t), for example, where each Y_i(t) originates from a different system?

Using add_input again fails, presumably because the In variable of sys gets overridden:

julia> sys = add_input(sys, 2)

julia> states(sys)
4-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 Sub₊In(t)
 input2₊Y(t)
 input1₊Y(t)
 Sub₊X(t)

julia> equations(sys)
4-element Vector{Equation}:
 Sub₊In(t) ~ input2₊Y(t)
 Differential(t)(input2₊Y(t)) ~ -input2₊Y(t)
 Differential(t)(input1₊Y(t)) ~ -input1₊Y(t)
 Sub₊X(t) ~ Sub₊In(t)

Any advice on how to circumvent this problem?

Thanks a lot in advance!

I don’t really understand the question. What’s the expected output of add_input(sys, 2)? Also, have you seen the connector example in Acausal Component-Based Modeling the RC Circuit · ModelingToolkit.jl?

Apologies for not explaining it clearly enough.
I did see the connector example, but I will have another look, thank you.

I will try to rephrase my question in a better way.
Say I have a subsystem A (with a variable X), receiving inputs from N subsystems B_i (with a variable Y).

What would be the best way (if any) to achieve X(t) = sum(B_1.Y(t) + B_2.Y(t) + ... + B_N.Y(t)),
given that I do not want to specify N manually for every subsystem.

To connect this to the example above, I would like to obtain


julia> states(sys)
4-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 Sub₊In(t)
 input2₊Y(t)
 input1₊Y(t)
 Sub₊X(t)

julia> equations(sys)
4-element Vector{Equation}:
 Sub₊In(t) ~ input1₊Y(t) + input2₊Y(t) 
 Differential(t)(input2₊Y(t)) ~ -input2₊Y(t)
 Differential(t)(input1₊Y(t)) ~ -input1₊Y(t)
 Sub₊X(t) ~ Sub₊In(t)

but without manually specifying all these equations a priori.

The reason for this is that I have many identical subsystems (i.e. their internal dynamics are the same), where the only difference to their overall dynamics is the number of inputs they receive.
So I would rather not have to build 100s of different subsystems if the underlying internal dynamics are the same!

Any insight would be greatly appreciated.

Do you mean something like this?

using ModelingToolkit

@variables t
D = Differential(t)

@named inner_sys = ODESystem(Equation[], t)

function subsystem(;name)
    states = @variables X(t)=0. In(t)=0.
    eqs = [X ~ In]

    ODESystem(eqs, t, states, []; name=name)
end

@named sub = subsystem()
@named sys = compose(inner_sys, sub)

function input_factory(;name)
    states = @variables Y(t)
    eq = [
        D(Y) ~ -Y
    ]
    ODESystem(eq, t, states, []; name=name)
end

@named inputs[1:2] = input_factory();
@named tmp = ODESystem([sub.In ~ sum(input->input.Y, inputs)])
@named sys_full = compose(extend(sys, tmp), inputs)

it gives

julia> equations(sys_full)
4-element Vector{Equation}:
 sub₊In(t) ~ inputs_1₊Y(t) + inputs_2₊Y(t)
 sub₊X(t) ~ sub₊In(t)
 Differential(t)(inputs_1₊Y(t)) ~ -inputs_1₊Y(t)
 Differential(t)(inputs_2₊Y(t)) ~ -inputs_2₊Y(t)

Thank you, that already looks like a nice step forward!

I need to see if I can shape this in a way that works for me, for it seems I still need to know all the inputs at the same time (where you specify @named inputs[1:2] = input_factory(); for instance).

What I had in mind is having one function, something like connect_B_to_A, that generates the ODESystem containing the equations necessary for B to provide input to A.
Then I would loop over all pairs of systems, and keep extend-ing the original system with these connections.
In summary, I would ideally like to provide the inputs one by one rather than all simultaneously.

But maybe that is not possible, since that would imply doing something along the lines of

function connect_B_to_A(sub, input)
    sub.In += input.Y
end

(or, I suppose, sub.In +~ input.Y in the Symbolics language)?

Mutating symbolic expressions might be an anti-pattern. Maybe you could save the operations in a transaction vector, and construct the final expression from that.

Agreed, thanks for the help!