# ModelingToolkit: linearization failure

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.ṁ_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?

Likely a bug in MTK, please open an issue

1 Like