# MTK: looping through variables in plot?

**URL:** <https://discourse.julialang.org/t/mtk-looping-through-variables-in-plot/117339>\
**Category:** Modelling & Simulations\
**Created:** [July 22, 2024, 12:31pm UTC](https://discourse.julialang.org/t/mtk-looping-through-variables-in-plot/117339 "2024-07-22T12:31:58Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [July 22, 2024, 12:31pm UTC](https://discourse.julialang.org/t/mtk-looping-through-variables-in-plot/117339/1 "2024-07-22T12:31:58Z")

</div>

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

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/1/5121e509225d080e2610a22c6d7832df6f20b5ef.png)  
which is cascaded + an input voltage and a terminal load model:

```julia
# 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.,

```julia
ℓ = 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:

```julia
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))

```

leading to (example):

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/7/a7ba858c5df556ad99365e6ebbec12cb2f6c6264.png)

* * *

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:

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 22, 2024, 1:09pm UTC](https://discourse.julialang.org/t/mtk-looping-through-variables-in-plot/117339/2 "2024-07-22T13:09:05Z")

</div>

> [@BLI](#):
>
> **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.

---

<div class="post-metadata">

**Author:** ![hexaeder](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hexaeder/32/24403_2.png) [@hexaeder](https://discourse.julialang.org/u/hexaeder)\
**Post date:** [July 22, 2024, 1:12pm UTC](https://discourse.julialang.org/t/mtk-looping-through-variables-in-plot/117339/3 "2024-07-22T13:12:25Z")

</div>

The unexported `getname` might be another useful function here:

```julia
import ModelingToolkit as MTK

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

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [July 23, 2024, 10:08am UTC](https://discourse.julialang.org/t/mtk-looping-through-variables-in-plot/117339/4 "2024-07-23T10:08:14Z")

</div>

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

```julia
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
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
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.

---

<div class="post-metadata">

**Author:** ![hexaeder](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hexaeder/32/24403_2.png) [@hexaeder](https://discourse.julialang.org/u/hexaeder)\
**Post date:** [July 24, 2024, 8:11am UTC](https://discourse.julialang.org/t/mtk-looping-through-variables-in-plot/117339/5 "2024-07-24T08:11:26Z")

</div>

Which version of MTK are you using? I just tried it with `N=8` and it worked for me:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/3/03b243c6c00f67722d44361d68e92b056328e02b.jpeg)

Those are my versions

```julia
(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:

```julia
using Pkg
pkg"activate --temp"
pkg"add ModelingToolkit, OrdinaryDiffEq, Plots, LaTeXStrings"
# 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))

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [July 24, 2024, 10:23am UTC](https://discourse.julialang.org/t/mtk-looping-through-variables-in-plot/117339/6 "2024-07-24T10:23:13Z")

</div>

> [@hexaeder](#):
>
> `u_vars = filter(v -> contains(string(getname(v)), "₊u"), unknowns(tL))`

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!

* * *

Two final comments…

- 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.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [July 24, 2024, 11:03am UTC](https://discourse.julialang.org/t/mtk-looping-through-variables-in-plot/117339/7 "2024-07-24T11:03:26Z")

</div>

> [@hexaeder](#):
>
> Which version of MTK are you using?

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

```julia
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:

```julia
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.