# ModelingToolkit: linearization failure

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-linearization-failure/116899>\
**Category:** Modelling & Simulations\
**Created:** [July 11, 2024, 10:42am UTC](https://discourse.julialang.org/t/modelingtoolkit-linearization-failure/116899 "2024-07-11T10:42:21Z")\
**Posts on this page:** 2\
**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 11, 2024, 10:42am UTC](https://discourse.julialang.org/t/modelingtoolkit-linearization-failure/116899/1 "2024-07-11T10:42:21Z")

</div>

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:

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

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

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

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

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

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 11, 2024, 12:33pm UTC](https://discourse.julialang.org/t/modelingtoolkit-linearization-failure/116899/2 "2024-07-11T12:33:08Z")

</div>

Likely a bug in MTK, please open an issue 😊