# MTK linearization: Change of use, or bug?

**URL:** https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921
**Category:** Modelling & Simulations
**Created:** [August 7, 2024, 2:49pm UTC](https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921 "2024-08-07T14:49:16Z")
**Posts on this page:** 9
**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: [August 7, 2024, 2:49pm UTC](https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921/1 "2024-08-07T14:49:16Z")

</div>

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.]

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

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

```

---

<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: [August 7, 2024, 2:51pm UTC](https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921/2 "2024-08-07T14:51:30Z")

</div>

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

---

<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: [August 7, 2024, 2:52pm UTC](https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921/3 "2024-08-07T14:52:46Z")

</div>

> Using `complete` doesn’t change the response.

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

---

<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: [August 7, 2024, 3:08pm UTC](https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921/4 "2024-08-07T15:08:25Z")

</div>

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:

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

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

```

followed by:

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

---

<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: [August 7, 2024, 3:14pm UTC](https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921/5 "2024-08-07T15:14:09Z")

</div>

> [@BLI](#):
>
> This leads to an error message.

What does it say?

> [@BLI](#):
>
> I _assume_ the problem is that the above strategy attempts to set `tank.m=>m_ss` instead of `tank_noi.m=>m_ss)` ???

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).

> [@BLI](#):
>
> 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`…

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

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

```

---

<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: [August 7, 2024, 3:17pm UTC](https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921/6 "2024-08-07T15:17:16Z")

</div>

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.

---

<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: [August 8, 2024, 9:10am UTC](https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921/7 "2024-08-08T09:10:21Z")

</div>

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

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

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

---

<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: [August 8, 2024, 9:22am UTC](https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921/8 "2024-08-08T09:22:54Z")

</div>

> [@BLI](#):
>
> `[h]`

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](https://ieeexplore.ieee.org/abstract/document/8816763)

---

<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: [August 8, 2024, 10:07am UTC](https://discourse.julialang.org/t/mtk-linearization-change-of-use-or-bug/117921/9 "2024-08-08T10:07:04Z")

</div>

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
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.
