Accessing variables in a ModelingToolkit system

Specific Domains Modelling & Simulations
1

First off, the new functionality for connecting systems in ModelingToolkit is really powerful, awesome job!

I have a few issues with it that I can’t solve, however…

Suppose I have a ModelingToolkit ODESystem

od::ODESystem

To interactively access a state V(t) or a state W(t) in the subsystem sub, I need only write the following code

od.V
od.sub.W

How should I non-interactively access these variables? IE I want to build functions of the form

return_state(od, "V")
## gives od.V(t)

return substate(od, "sub", "W")
## gives od.sub.W(t) and not sub.W(t).

The difference between od.sub.W(t) and sub.W(t) seems to be important when connecting subsubsystems to other ODESystems.

Struggling to build such functions.

Thanks a lot in advance!

Also, to be a bit more informative as to why this is useful, I have a system with states x_1(t) … x_N(t). These are a subset of all the system states, and their position in the system.states array is variable. I need a function which takes in an integer i, and returns od.x_i(t)

1 Like
2

You can use the very new @nonamespace

using ModelingToolkit
@parameters v, t
@named d = ODESystem(Equation[],t,[],[v])
@named c = ODESystem(Equation[],t,[],[v],systems=[d])
@named b = ODESystem(Equation[],t,[],[v],systems=[c])
@named a = ODESystem(Equation[],t,[],[v],systems=[b])

a.b.c.d.v # a₊b₊c₊d₊v

@nonamespace a.b.c.d.v # v
@nonamespace tmp = a.b.c.d # Model d with 0 equations
tmp.v # d₊v

@nonamespace tmp = a.b # Model b with 0 equations
tmp.c.d.v # b.c.d.v

or getproperty(…; namepsace=false)

@parameters V, W
@named sub = ODESystem(Equation[],t,[],[V,W])
@named od = ODESystem(Equation[],t,[],[V,W],systems=[sub])

substate(sys,sub,state) = getproperty(getproperty(sys,Symbol(sub),namespace=false), Symbol(state))
substate(od, "sub", "W") # sub₊W

@variables V[1:3], W[1:3]
@named sub = ODESystem(Equation[],t,[V...,W...],[])
@named od = ODESystem(Equation[],t,[V...,W...],[],systems=[sub])
substate(sys,sub,state,i) = getproperty(getproperty(sys,Symbol(sub),namespace=false), Symbol("$(state)$(join('₀'+d for d in reverse(digits(i))))"))
substate(od,"sub","W",1) # sub₊W₁
6 Likes
3

We should probably document the programmatic ways to generate the namespaced names. Open an issue?

4

Maybe @pepijndevos has some ideas here/

5

works perfectly. thanks @lungd!

6

will do!

7

The nonamespace solution is ok, but I think there might be an underlying behavior that could be improved.

I think part of this is due to how namespacing of the toplevel system works currently.
Consider a flat system, if you do equations(sys) you get the toplevel equations without any namespacing, but if you do sys.var you get a namespaced variable. This then causes issues when connecting variables at the top level.

I’m not sure which way the solution would go, but if the namespacing on states(sys) and sys.var match, it’d avoid a lot of these challenges I think.

2 Likes
8

I’m also strugglig with this. Consider trying to supply initial conditions for the state P.x₁ in the following system

julia> fb
Model fb with 13 equations
States (13):
  u(t)
  loopgain₊y(t)
  loopgain₊u(t)
  y(t)
  loopgain₊C₊u(t)
  loopgain₊C₊y(t)
  loopgain₊nonlinear_P₊u(t)
  loopgain₊nonlinear_P₊y(t)
  loopgain₊C₊x₁(t) [defaults to 0.0]
  loopgain₊C₊x₂(t) [defaults to 0.0]
  loopgain₊nonlinear_P₊P₊u(t)
  loopgain₊nonlinear_P₊P₊y(t)
⋮
Parameters (0):

x0 = Pair[
    # fb.loopgain.nonlinear_P.P.x₁ => 1 # does not work
    loopgain.nonlinear_P.P.x₁ => 1      # works 
    # P.x₁ => 1                         # does not work
]

It’s not clear to me why I cannot use the last one, i.e., P.x₁ => 1. If I inspect the available states after structural simplification, I learn that the only remaining states are

julia> simplified_sys
Model fb with 3 equations
States (3):
  loopgain₊C₊x₁(t) [defaults to 0.0]
  loopgain₊C₊x₂(t) [defaults to 0.0]
  loopgain₊nonlinear_P₊P₊x₁(t) [defaults to 0.0]
Parameters (0):
Incidence matrix:
⠼⠎⠢

I wonder if it could be made smoother to supply the initial state for P₊x₁(t) here. If I have passed a system P into a function that returns something where I can only access the state as loopgain₊nonlinear_P₊P₊x₁(t) it feels like I can easily get lost.

9

@pepijndevos, Can you show a concrete example where you would face such an issue?
Maybe I already got used to the connector system you can see in the tutorials where you would define “flat” as a subsystem of a “full” system. By doing so, you connect the var with e.g. flat.var ~ 1 and there you need the namespacing.

@baggepinnen, I think you want to define default values while creating the subsystems.
P = ..System(..., defaults=[x₁ => 1]) or something like PSys(u0_x₁=1) = ..System(..., defaults=[x₁ => u0_x₁]).
If you add P to nonlinear_P’s subsystems you can later use defaults(simplified_sys) to get the namespaced initial states like loopgain₊nonlinear_P₊P₊x₁(t) => 1
If you want to change the initial state of the subsystem after simplifying “fb”, you must use the full name including namespaces.
Your system only contains one “P” sys. But how to select the correct state if there is another subsystem with P.x₁?

10

I think my confusion comes partly from the fact that I do not consider P the name of a system here, I have a Julia object P that happens to also have the name (in MTK sense) P. If I refer to P.x, my intuition says that it should be the same as loopgain₊nonlinear_P₊P₊x since the innermost P there came from the P I’m holding in my hand. To make it more concrete, I expected something like

julia> struct Bar a end

julia> a = [0]
1-element Vector{Int64}:
 0

julia> bar = Bar(a)
Bar([0])

julia> bar.a === a
true

which does not appear to hold for MTK systems.