# MTK: computing solutions of vectors of unknowns?

With ModelingToolkit, I have a number of variables of.w_1.m, of.w_2.m, …, of_w_5.m and I want to check the solution of these variables, e.g., at time = 0.5.

With solution structure sol, I can do this as follows:

sol(0.5, idxs=[of.w_1.m, of.w_2.m, ..., of.w_5.m])


Question: is it possible to create a loop in order to set up the vector:
[of.w_1.m, of.w_2.m, ..., of.w_5.m]?

I have tried with

sol(0.5, idxs=[eval(Symbol("of.w_$j.m")) for j in 1:5])  without success. I’d like to use the same idea for plots, i.e., replace: plot(sol, idxs = [of.w_1.m, of.w_2.m, ..., of.w_5.m])  with some loop structure. If you have an array variable you just just name that the index. I don’t understand. I created w as an array inside of an @mtkmodel macro, and I have figured out how to use w as an array within this model. But when I instantiate the model using @mtkbuild ..., I have to address vector w elements by w_1, w_2, etc. instead of w[1], etc. I have tried with… sol(0.5, idxs = of.w.m) sol(0.5, idxs = of.w[:].m) sol(0.5, idxs = [of.w_j.m for j in 1:5]) ...  No luck. I can try to create a simple example. How about this getvarbyname(sys, name) = mapfoldl(Symbol, getproperty, split(name, '.'); init=sys) # and then [getvarbyname(of, "w_$j.m") for j in 1:5]


Although you may want to use ₊ for the separator rather than . to stay with the MTK convention.

So here is a relatively simple example… a bank of N RC circuits, all driven by the same voltage v_\mathrm{i} = b + \sin(\omega\cdot t), with initial values for the capacitor voltage v_\mathrm{c} generated by a random number generator, and resistance R also given random values.

So the N RC circuits all may have different currents i and different capacitor voltages v_\mathrm{c}. I also compute the total capacitor power P.

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations
using Plots
#
@mtkmodel RC begin
@parameters begin
R
C = 10e-3
end
@variables begin
v_i(t)
i(t)
v_c(t)=randn()
end
@equations begin
i ~ C*D(v_c)
v_i ~ R*i + v_c
end
end
#
@mtkmodel RCBank begin
@parameters begin
ω = 2pi*50/3
b = 1
end
@structural_parameters begin
N = 5
end
@variables begin
P(t)
end
@components begin
bank = [RC(; R=10 + 10*rand()) for j in 1:N]
end
@equations begin
[bank[j].v_i ~ b + sin(ω*t) for j in 1:N]...
P ~ sum([bank[j].v_c*bank[j].i for j in 1:N])
end
end
#
@mtkbuild sys = RCBank()

tspan=(0,1)

prob = ODEProblem(sys, [], tspan)
sol = solve(prob)

plot(sol)


The result is (e.g., since the system is random…):

My question is really: how can I plot all the currents in an elegant and simple way?

plot(sol, idxs=[sys.bank_1.i, sys.bank_2.i, sys.bank_3.i, sys.bank_4.i, sys.bank_5.i])


works, but it seems clumsy.

This seems kind of like what you want but still has 2 problems. N can’t be used in variable definitions (LoadError: promotion of types Int64 and Symbol failed to change any arguments, maybe a bug?). And the final plot just uses I as the legend for all lines, the index is lost.

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations
using Plots
#
@mtkmodel RC begin
@parameters begin
R
C = 10e-3
end
@variables begin
v_i(t)
i(t)
v_c(t)=randn()
end
@equations begin
i ~ C*D(v_c)
v_i ~ R*i + v_c
end
end
#
@mtkmodel RCBank begin
@parameters begin
ω = 2pi*50/3
b = 1
end
@structural_parameters begin
N = 5
end
@variables begin
P(t)
(I(t))[1:5]
end
@components begin
bank = [RC(; R=10 + 10*rand()) for j in 1:N]
end
@equations begin
[bank[j].v_i ~ b + sin(ω*t) for j in 1:N]...
P ~ sum([bank[j].v_c*bank[j].i for j in 1:N])
I ~ [bank[j].i for j in 1:N]
end
end

function testrcbank()
@mtkbuild sys = RCBank()

tspan=(0,1)

prob = ODEProblem(sys, [], tspan)
sol = solve(prob)

plot(sol; idxs=sys.I)
end

1 Like

Thanks for suggestion. Let me try to summarize: essentially, you suggest one solution, and point to two possible bugs:

A. Your solution is to copy the currents that are not “unknowns” (states in this case) up from the RC model to the RCBank model. [By introducing new variable I(t) in the RCBank model and copying currents bank[j].i into the new variables.]

• This may appear as an unnecessary introduction of extra variables.
• However, it seems like MTK/the plot recipe can handle arrays at the end of the variable name, e.g., sys.I (where I is an array), but can not handle the situation when the array is inside of the variable name, e.g., sys.bank.i (where bank is the array).

B: Possible bug?

• It is a little bit strange that the structural parameter N can not be used to define the size of variable I, i.e., I(t)[1:5] works but I(t)[1:N] does not work.
• In the legend of the plot, the plot recipe strips off sys from sys.I and only uses the I in the labels. Ideally, perhaps it should have used I_1, I_2, etc.

[OK – normally, I manually set the label anyway, so not a big deal that all plots get the same label.]

I’m trying to use your getvarbyname function on the RC circuit. It doesn’t quite work:

getvarbyname(sys, name) = mapfoldl(Symbol, getproperty, split(name, '.'); init=sys)
# application
plotvar = [getvarbyname(sys, "bank_$j.i") for j in 1:5]  produces the vector: Num[bank_1₊i(t), bank_2₊i(t), bank_3₊i(t), bank_4₊i(t), bank_5₊i(t)]  The vector elements need to be prepended by sys. in order to work. It is not quite clear to me how to achieve that. I think that is just a printing convention. In the context of your example script, this works for me: getvarbyname(sys, name) = mapfoldl(Symbol, getproperty, split(name, '.'); init=sys) function testbyname() @mtkbuild sys = RCBank() tspan=(0,1) prob = ODEProblem(sys, [], tspan) sol = solve(prob) plotvar = [getvarbyname(sys, "bank_$j.i") for j in 1:5]

plot(sol; idxs=plotvar)
end


Arghh. Yes, it works… I did a silly mistake (plot(sys; ...) instead of plot(sol; ...)).

MTK is very good at simplifying these out, this will have no cost during solving and only be evaluated during plotting.

This is actually introducing a symbolic array, which is different from a vector of systems (which is what bank is).

1 Like

Yeah, I guess the I will be made into an “observed” variable – it should be trivial since it is just an assignment…

But I’m curious as to why there was a problem with N.

Anyways, perhaps it is better to “promote” the extra variables (I) to the top level of the system as a “symbolic array” just to avoid introducing a new function (getvarbyname).