# MTK: looping through variables in plot?

I’m struggling to figure out how to loop through numbered variables in ModelingToolkit [MTK].

As a simple, yet concrete example, consider a model of a telegraph line, which consists in a telegraph line module component

which is cascaded + an input voltage and a terminal load model:

# Importing packages
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as Dt
using DifferentialEquations
using Plots, LaTeXStrings
#
# Some plot constants
LW1 = 2.5
#
LS1 = :solid
LS2 = :dash
LS3 = :dot
#
LC1 = :blue
LC2 = :red
LC3 = :green
#
# Telegraph Line module
@mtkmodel Telegraph_Module begin
# Model parameters
@parameters begin
# Default parameters X_, X ∈ {R,L,G,C} taken from https://en.wikipedia.org/wiki/Telegrapher%27s_equations
R_=172.28,      [description = "Unit resistance, Ω/km"]
L_=612.5e-6,    [description = "Unit inductance, H/km"]
G_=0.072e-6,    [description = "Unit conductance, S/km"]
C_=51.57e-9,    [description = "Unit capacitance, F/km"]
Δx=1,           [description = "Module length, km"]
R=R_*Δx,        [description = "Resistance, Ω"]
L=L_*Δx,        [description = "Inductance, H"]
G=G_*Δx,        [description = "Conductance, S"]
C=C_*Δx,        [description = "Capacitance, F"]
end
#
# Model variables, with initial values needed
@variables begin
# Inputs
u(t),           [description = "Voltage at input, V"]
# Variables
u₊(t)=0.5,      [description = "Voltage at output, V"]
i(t)=1e-3,      [description = "Current at input, A"]
i₊(t),          [description = "Current at output, A"]
end
#
# Telegraph line module equations
#
@equations begin
u ~ R*i + L*Dt(i) + u₊
i ~ G*u₊ + C*Dt(u₊) + i₊
end
end
# Telegraph Line system: cascaded modules with input voltage + terminal load
@mtkmodel Telegraph_Line begin
#
@parameters begin
R_L=600,    [description = "Load resistance, Ω"]
end
#
@structural_parameters begin
ℓ = 10 # km
N = 2
end
# Components used
@components begin
tm = [Telegraph_Module(; Δx=ℓ/N) for i in 1:N]
end
# Equations for connecting components
@equations begin
# Distribution manifold splitting
tm[1].u ~ sin(2π*1_000*t)                 # Input at 1 kHz
[tm[j].u ~ tm[j-1].u₊ for j in 2:N]...
[tm[j].i ~ tm[j-1].i₊ for j in 2:N]...
tm[N].i₊ ~ tm[N].u₊/R_L                      # Terminal load
end
end


One would then instantiate the model, and solve the model, e.g.,

ℓ = 18    # length in km
N = 3    # number of cascaded modules
# Instantiation
@mtkbuild tL = Telegraph_Line(; ℓ=ℓ, N=N)
#
# Simulation set-up
tspan = (0, 1e-2)
#
prob = ODEProblem(tL, [], tspan)
sol = solve(prob)


One can then plot the solution, e.g., voltages:

plot(sol, idxs=[tL.tm_1.u₊, tL.tm_2.u₊, tL.tm_3.u₊], lw=LW1, lc=[LC1 LC2 LC3], ls=[LS1 LS2 LS3],
xlabel=L"t \mathrm{\ [s]}", ylabel=L"u \mathrm{\ [V]}",
title="Telegraph line voltage",
frame=:box,
ylims=(-0.5,0.75))


Suppose I now want to use more cascaded modules, say N=10 or something? And that I want to plot all voltages? Then it becomes very impractical to explicitly write every voltage variable as:

plot(sol, idxs=[tL.tm_1.u₊, tL.tm_2.u₊, tL.tm_3.u₊, tL.tm_4.u₊, ....]


Question: How can I construct the vector of indices [tL.tm_1.u₊, tL.tm_2.u₊, tL.tm_3.u₊, tL.tm_4.u₊, ....] using comprehension, or similar techniques?

findfirst or other things on unknowns(sys), getfield.(observed(sys),:lhs). nameof and String doing some string matching to find all voltages or whatnot.

2 Likes

The unexported getname might be another useful function here:

import ModelingToolkit as MTK

filter(v -> contains(string(MTK.getname(v)), "₊u"), unknowns(tL))

1 Like

@hexaeder : OK – manually setting up the plot works (as before):

plot(sol, idxs=[tL.tm_1.u₊, tL.tm_2.u₊, tL.tm_3.u₊])


produces the plot. But here, N=3, and such manual set-up becomes untenable for large N or experimenting with different values of N.

Following your suggestion (assuming MTK as been defined), I get:

julia> u_vars = filter(v -> contains(string(MTK.getname(v)), "₊u"), unknowns(tL))
3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
tm_1₊u₊(t)
tm_2₊u₊(t)
tm_3₊u₊(t)


There are two problems with this vector:

• the time dependence ((t)) part of the vector elements can not be included when doing plotting, i.e., plot(sol, idxs=u_vars) does not work
• the vector elements must be prepended with tL. in order to work in the plotting

I had, of course, tried with simple string manipulation to set up the vector prior to posing the question, e.g.,

julia> [Symbol("tL.tm_($j).u₊") for j in 1:3] [Symbol("tL.tm_($j).u₊") for j in 1:3]


which can then be converted to a variable via mapping eval, but that doesn’t work because the variables are not variables in the global namespace, but rather symbolic names (I guess).

So I’m still stuck.

Which version of MTK are you using? I just tried it with N=8 and it worked for me:

Those are my versions

(jl_GORUmQ) pkg> st
Status /tmp/jl_GORUmQ/Project.toml
[b964fa9f] LaTeXStrings v1.3.1
[961ee093] ModelingToolkit v9.26.0
[1dea7af3] OrdinaryDiffEq v6.86.0
[91a5bcdd] Plots v1.40.5


Here is the full code I used including setup of the temp environment:

using Pkg
pkg"activate --temp"
# Importing packages
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as Dt, getname
using OrdinaryDiffEq
using Plots, LaTeXStrings
#
# Some plot constants
LW1 = 2.5
#
LS1 = :solid
LS2 = :dash
LS3 = :dot
#
LC1 = :blue
LC2 = :red
LC3 = :green
#
# Telegraph Line module
@mtkmodel Telegraph_Module begin
# Model parameters
@parameters begin
# Default parameters X_, X ∈ {R,L,G,C} taken from https://en.wikipedia.org/wiki/Telegrapher%27s_equations
R_=172.28,      [description = "Unit resistance, Ω/km"]
L_=612.5e-6,    [description = "Unit inductance, H/km"]
G_=0.072e-6,    [description = "Unit conductance, S/km"]
C_=51.57e-9,    [description = "Unit capacitance, F/km"]
Δx=1,           [description = "Module length, km"]
R=R_*Δx,        [description = "Resistance, Ω"]
L=L_*Δx,        [description = "Inductance, H"]
G=G_*Δx,        [description = "Conductance, S"]
C=C_*Δx,        [description = "Capacitance, F"]
end
#
# Model variables, with initial values needed
@variables begin
# Inputs
u(t),           [description = "Voltage at input, V"]
# Variables
u₊(t)=0.5,      [description = "Voltage at output, V"]
i(t)=1e-3,      [description = "Current at input, A"]
i₊(t),          [description = "Current at output, A"]
end
#
# Telegraph line module equations
#
@equations begin
u ~ R*i + L*Dt(i) + u₊
i ~ G*u₊ + C*Dt(u₊) + i₊
end
end
# Telegraph Line system: cascaded modules with input voltage + terminal load
@mtkmodel Telegraph_Line begin
#
@parameters begin
R_L=600,    [description = "Load resistance, Ω"]
end
#
@structural_parameters begin
ℓ = 10 # km
N = 2
end
# Components used
@components begin
tm = [Telegraph_Module(; Δx=ℓ/N) for i in 1:N]
end
# Equations for connecting components
@equations begin
# Distribution manifold splitting
tm[1].u ~ sin(2π*1_000*t)                 # Input at 1 kHz
[tm[j].u ~ tm[j-1].u₊ for j in 2:N]...
[tm[j].i ~ tm[j-1].i₊ for j in 2:N]...
tm[N].i₊ ~ tm[N].u₊/R_L                      # Terminal load
end
end
ℓ = 18    # length in km
N = 8    # number of cascaded modules
# Instantiation
@mtkbuild tL = Telegraph_Line(; ℓ=ℓ, N=N)
#
# Simulation set-up
tspan = (0, 1e-2)
#
prob = ODEProblem(tL, [], tspan)
sol = solve(prob)

u_vars = filter(v -> contains(string(getname(v)), "₊u"), unknowns(tL))
plot(sol, idxs=u_vars, lw=LW1, lc=[LC1 LC2 LC3], ls=[LS1 LS2 LS3],
xlabel=L"t \mathrm{\ [s]}", ylabel=L"u \mathrm{\ [V]}",
title="Telegraph line voltage",
frame=:box,
ylims=(-0.5,0.75))


Argh… I used the same version numbers as you, except v9.25.0 for MTK.

And your code works also on that version. I guess I got confused and didn’t check the eltype of the vector u_vars… I’m used to specifying plotting indices (idxs = ...) as Julia identifiers, and in that case, I need to prepend the variable names by the “object” (tL. in my case) + I cannot have explicit time dependency.

But with your proposed solution, the plotting indices are instead of eltype SymbolicUtils.BasicSymbolic{Real}, which I guess holds the required information.

Thanks for the help!

• Perhaps my choice of variable names is poor? I used u to indicate u_x and u₊ to inducate u_{x+\Delta x}… Maybe unfortunate in that MTK seems to use + instead of . when pretty-printing equations?
• It is somewhat confusing to me that you use contains(string(getname(v)), "₊u"). Would contains(string(getname(v)), ".u") also have worked? I can, of course check this myself, but I think this is a potential source of confusion.

It would perhaps be useful with some utility function for picking out variables? Something like:

function get_idxs(model, pattern, vartype=:unknowns)
if vartype == :unknowns
filter(v -> contains(string(getname(v)), pattern), unknowns(model))
elseif vartype == :observed
filter(v -> contains(string(getname(v)), pattern), getfield.(observed(model),:lhs))
end
end


??

Example of use:

u_vars = get_idxs(tL,"u₊")
i_vars = get_idxs(tL, "i₊", :observed)
plot(sol, idxs=[u_vars..., i_vars...])


I’m sure MTK developers would be able to improve on such a function, e.g., better name, etc., etc.

2 Likes