# MTK: computing solutions of vectors of unknowns?

**URL:** <https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614>\
**Category:** Modelling & Simulations\
**Created:** [April 29, 2024, 3:30pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614 "2024-04-29T15:30:50Z")\
**Posts on this page:** 12\
**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:** [April 29, 2024, 3:30pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/1 "2024-04-29T15:30:50Z")

</div>

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:

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

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

```julia
plot(sol, idxs = [of.w_1.m, of.w_2.m, ..., of.w_5.m])

```

with some loop structure.

---

<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:** [April 29, 2024, 4:20pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/2 "2024-04-29T16:20:47Z")

</div>

If you have an array variable you just just name that the index.

---

<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:** [April 29, 2024, 4:53pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/3 "2024-04-29T16:53:42Z")

</div>

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…

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

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [April 29, 2024, 5:30pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/4 "2024-04-29T17:30:04Z")

</div>

How about this

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

---

<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:** [April 29, 2024, 6:55pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/5 "2024-04-29T18:55:41Z")

</div>

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.

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

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/d/0d0223a39331186fdba28246e0d498910a5ff05b.png)

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

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

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [April 29, 2024, 7:50pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/6 "2024-04-29T19:50:18Z")

</div>

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.

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

```

---

<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:** [April 30, 2024, 7:11am UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/7 "2024-04-30T07:11:36Z")

</div>

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

---

<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:** [April 30, 2024, 2:53pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/8 "2024-04-30T14:53:47Z")

</div>

> [@contradict](#):
>
> ```julia
> getvarbyname(sys, name) = mapfoldl(Symbol, getproperty, split(name, '.'); init=sys)
> # and then
> [getvarbyname(of, "w_$j.m") for j in 1:5]
> 
> ```

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

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

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

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [April 30, 2024, 3:05pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/9 "2024-04-30T15:05:12Z")

</div>

I think that is just a printing convention. In the context of your example script, this works for me:

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

```

---

<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:** [April 30, 2024, 3:14pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/10 "2024-04-30T15:14:51Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [April 30, 2024, 3:21pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/11 "2024-04-30T15:21:17Z")

</div>

> [@BLI](#):
>
> - This may appear as an unnecessary introduction of extra variables.

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

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

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

---

<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:** [April 30, 2024, 3:29pm UTC](https://discourse.julialang.org/t/mtk-computing-solutions-of-vectors-of-unknowns/113614/12 "2024-04-30T15:29:08Z")

</div>

> [@contradict](#):
>
> MTK is very good at simplifying these out, this will have no cost during solving and only be evaluated during plotting.

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