I’m trying to linearize a “complex” (for me) ModelingToolkit [MTK] model consisting of several layers of components that I have written. My components do not use ports (I’m in the process of learning this), and I try to do linearization without using the standard library, so no Analysis points at the moment.

First – let me illustrate the *procedure* I use for linearization using a simple model consisting of a liquid tank.

Packages:

```
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as Dt
using DifferentialEquations
using ControlSystems
```

I define a mass flow rate `md`

[math notation \dot{m}, `m`

with a dot, `md`

]

```
md_const(t) = 2
md_step(t) = t < 1 ? 2 : 1.5
#
@register_symbolic md(t)
md(t) = md_const(t)
```

The balanced tank model:

```
# Below, macro `mtkmodel` creates a model instantiator/"class" with name `Tank`
#
@mtkmodel Tank begin
# Model parameters
@parameters begin
ρ=1, [description = "Liquid density"]
A=5, [description = "Cross sectional tank area"]
K=5, [description = "Effluent valve constant"]
h_ς=3, [description = "Scaling level in valve model"]
end
# Model variables, with initial values needed
@variables begin
m(t)=1.5*ρ*A, [description = "Liquid mass"]
md_i(t), [description = "Influent mass flow rate"]
md_e(t), [description = "Effluent mass flow rate"]
V(t), [description = "Liquid volume"]
h(t), [description = "level"]
end
# Providing model equations
@equations begin
Dt(m) ~ md_i-md_e
m ~ ρ*V
V ~ A*h
md_e ~ K*sqrt(h/h_ς)
md_i ~ md(t)
end
end
#
# Instantiation
@mtkbuild tank = Tank()
;
```

Next, I define the same model instantiator *without input*, with the name `Tank_noi`

meaning `Tank`

with no/undefined input – observe that `Tank_noi`

is identical to `Tank`

*except* that I have commented out the equation connecting the inflow `md_i`

to the mass flow function `md`

.

```
# Below, macro `mtkmodel` creates a model instantiator/"class" with name `Tank`
#
@mtkmodel Tank_noi begin
# Model parameters
@parameters begin
ρ=1, [description = "Liquid density"]
A=5, [description = "Cross sectional tank area"]
K=5, [description = "Effluent valve constant"]
h_ς=3, [description = "Scaling level in valve model"]
end
# Model variables, with initial values needed
@variables begin
m(t)=1.5*ρ*A, [description = "Liquid mass"]
md_i(t), [description = "Influent mass flow rate"]
md_e(t), [description = "Effluent mass flow rate"]
V(t), [description = "Liquid volume"]
h(t), [description = "level"]
end
# Providing model equations
@equations begin
Dt(m) ~ md_i-md_e
m ~ ρ*V
V ~ A*h
md_e ~ K*sqrt(h/h_ς)
# md_i ~ md(t)
end
end
;
```

When I instantiate `Tank_noi`

, the result is an *imbalanced* model that can not be simplified nor can be simulated – but it can be linearized:

```
@named tank_noi = Tank_noi()
tank_noi = complete(tank_noi)
```

The balanced instance `tank`

can now be simulated “for a long time” to find the steady state mass `m_ss`

[which is the unknown in the model]:

```
tspan = (0,1e3)
#
prob = ODEProblem(tank, [], tspan)
sol = solve(prob)
m_ss = sol(tspan[2])[1]
```

Finally, I can substitute `m_ss`

into the *imbalanced* model and do (numeric) linearization, where I specify `md_i`

as input, and want to use level `h`

as output:

```
julia> mats, tank_ = linearize(tank_noi, [md_i], [h]; op = Dict(m=>m_ss, md_i=>md(0)))
julia> ss(mats...)
ControlSystemsBase.StateSpace{Continuous, Float64}
A =
-0.4166666666386296
B =
1.0
C =
0.2
D =
-0.0
Continuous-time state-space model
```

So far, so good: the procedure works for this simple example.

Next, I try to apply exactly the same procedure to a more complex model with several layers of components. Then I get an error message…

```
julia> @named of_noi = OilField()
julia> of_noi = complete(of_noi)
julia> unknowns(of_noi)
302-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
dm₊ṁ_i(t)
dm₊ṁ_e(t)
dm₊m(t)
dm₊p(t)
dm₊ρ(t)
dm₊Ṽ(t)
dm₊z(t)
w_1₊ṁ_i(t)
w_1₊p_i(t)
w_1₊ρ_i(t)
⋮
```

When I attempt to do the linearization, it fails even if I try to use the default operating point:

```
julia> linearize(of_noi, [of_noi.dm.ṁ_i], [of_noi.dm.p])
BoundsError: attempt to access 302-element Vector{Vector{Int64}} at index [303]
Stacktrace:
[1] getindex
@ .\essentials.jl:13 [inlined]
...
```

The input I have used is element 1 of the list of unknowns of `of_noi`

, and the output is element 5 of this list.

- Why does the error message say that I attempt to use element 303???
- Any ideas why my linearization attempt fails?