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?
