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. A Large-Scale Study About Quality and Reproducibility of Jupyter Notebooks | IEEE Conference Publication | IEEE Xplore ↩︎

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.