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:
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…):
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
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.]
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).
