# MTK linearization: Change of use, or bug?

I recently asked about a linearization problem I had – the problem persists. However, here I raise another question: ModelingToolkit linearization fails for a simple problem that used to work.

Question: What is it I do incorrectly below?
[I have used Julia v10.4 and ModelingToolkit v9.30 + DifferentialEquations v7.13.0, i.e., updated just now.]

``````# IMPORTING PACKAGES
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as Dt
using DifferentialEquations
using ControlSystems
#
# BALANCED SIMULATION 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
#
# INPUT FUNCTION
md_const(t) = 2
md(t) = md_const(t)
#
# UNBALANCED MODEL FOR LINEARIZATION
#=
The following model Tank_noi [noi: no input] is the same as Tank, except
that the equation connecting variable md_i to the input function md(t) has
been commented out
=#
#
@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
#
# INSTANTIATING SIMULATION MODEL + LINEARIZATION MODEL
@mtkbuild tank = Tank()
@named tank_noi = Tank_noi()
#
# SIMULATING SYSTEM TO FIND STEADY STATE
tspan = (0,1e3)
#
prob = ODEProblem(tank, [], tspan)
sol = solve(prob)
m_ss = sol(tspan[2])[1]
#
# ATTEMPTING TO LINEARIZE SYSTEM
mats_modern, tank_ = linearize(tank_noi, [tank_noi.md_i], [tank_noi.h]; op = Dict(tank_noi.m=>m_ss, tank_noi.md_i=>md(0)))
#ss(mats_modern...)
``````

The error message I get is:

``````Some specified inputs were not found in system. The following variables were not found Any[tank_noi₊md_i(t)]
``````

you forgot to call `tank_noi = complete(tank_noi)`?

1 Like

Using `complete` doesn’t change the response.

Correction: it does. Hm, I was positive I tested this before, with no success.

Follow-up question. In this case, I had few variables. I have another case (more complicated, with many variables).

How can I avoid manually listing all variables in the `op` dictionary? Just to indicate the problem… I need to specify the numeric values of the “unknowns” of the original model (`tank`) + the input `md_i`. Here, it is simple to write:

``````...op = Dict(tank_noi.m=>m_ss, tank_noi.md_i=>md(0))
``````

For large models, this is “unfeasable”. What I have done in the other model (which fails) is akin to:

``````unk = unknowns(tank)
op_unk = Dict(unk.=>sol(0, idxs=unk))
``````

followed by:

``````... op = Dict(op_unk, tank_noi.md_i=>md(0))
``````

This leads to an error message. I assume the problem is that the above strategy attempts to set `tank.m=>m_ss` instead of `tank_noi.m=>m_ss` ???

So … if my assumption above is correct, how can I pull out the subset of unknowns from `tank_noi` which corresponds to the unknowns in `tank_noi`… and in the correct order… so that I can associate the correct unknowns of `tank_noi` to the numeric value `m_ss`

… if my question makes sense??

Or is there a better way to do it?

What does it say?

This does not have to be a problem, as long as you have correctly called `complete` on both models (or used `@mtkbuild` for the simulation model).

your approach above should work, otherwise you can index the solution object with the unknown variables from `tank_noi`, something like this

``````op = [var => sol(sol.t[end], idxs = var) for var in unknowns(tank_noi)]
``````
1 Like

OK – I’ll check again tomorrow… I left my files open at work PC, and don’t want to mess them up by changing them on my home PC.

If there is still a problem, I’ll post the code as a Github issue instead.

I’m confused… I have my above `Tank` and `Tank_noi` models in two different Jupyter notebooks. In one notebook (the original one), doing:

``````...
@named tank_noi = Tank_noi()
@mtkbuild tank = Tank()
...
linearize(tank_noi, [md_i], [h]; op = Dict(m=>m_ss, md_i=>md(0)))
``````

works, while in the other notebook, I have to do:

``````...
@named tank_noi = Tank_noi()
tank_noi = complete(tank_noi)
@mtkbuild tank = Tank()
...
linearize(tank_noi, [tank_noi.md_i], [tank_noi.h]; op = Dict(tank_noi.m=>m_ss, tank_noi.md_i=>md(0)))
``````

to make it work.

In the first case (original notebook), I have loads of other code and models, some where I develop the models from scratch (i.e., without using `mtkmodel`/functions). I’m re-using the words `tank`, and possibly also `tank_noi`.

In the second case, I only have the model as listed above.

Any idea of what may cause this different behavior?

you probably have this defined in this notebook and it happens to have the correct namespace. Your problem is likely that you have too many global variables and some things work by coincidence because you happened to have a suitable variable defined. I recommend that you always use the latter approach with `complete` to avoid this.

Jupyter notebooks are notorious for leading to messy code in global namespace that fails to reproduce after restart etc. [1]

1 Like

Could be a coincidence, yes.

In my more complex problem, when I do the same things as above and do `complete` the unbalanced model, I get a weird error message:

``````julia> linearize(hps_noi, [hps_noi.u_v], [hps_noi.w], op=Dict(hps_noi.u_v=>0.89))
BoundsError: attempt to access ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa148094c, 0xba54127a, 0xe432ae56, 0x70bc53cc, 0xa7efd9f2), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a3c5c5b, 0x93608ad1, 0xbcc33a9c, 0x1971196f, 0xe3aae881), Nothing}} at index [2]

Stacktrace:
``````

I’ll post that as an issue in GitHub instead.