# ModelingToolkit -- possible to do "numeric" linearization?

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-possible-to-do-numeric-linearization/118933>\
**Category:** Modelling & Simulations\
**Created:** [September 2, 2024, 8:59am UTC](https://discourse.julialang.org/t/modelingtoolkit-possible-to-do-numeric-linearization/118933 "2024-09-02T08:59:45Z")\
**Posts on this page:** 7\
**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:** [September 2, 2024, 8:59am UTC](https://discourse.julialang.org/t/modelingtoolkit-possible-to-do-numeric-linearization/118933/1 "2024-09-02T08:59:45Z")

</div>

I’ve had problems doing linearization with ModelingToolkit [MTK] for “complicated” models since May (due to bugs?). I’m wondering whether it is possible to “massage” the MTK model into a form that allows for “old-fashioned” linearization by finite differences??

* * *

**Example model** :  
Consider a simple water tank with mass m, influent mass flow rate \dot{m}\_\mathrm{i}, effluent mass flow rate \dot{m}\_\mathrm{e},

\frac{\mathrm{d}m}{\mathrm{d}t} = \dot{m}\_\mathrm{i} - \dot{m}\_\mathrm{e}

Water has density \rho, leading to a water volume V contained in the tank which has constant cross sectional surface area A. The resulting water level is h,

m = \rho\cdot V \\ V = A\cdot h

Assuming effluent driven by gravity,

\dot{m}\_\mathrm{e} = K\cdot \sqrt{\frac{h}{h^\varsigma}}

As an example:

\rho = 1 \\ A = 5 \\ K = 5 \\ h^\varsigma = 3

* * *

Non-MTK (“classical”) model implementation:

```julia
function tank!(u, p, t)
    m = u
    p = ρ, A, K, h_ς, md(t)
    V = m/ρ 
    h = V/A 
    md_e = K*sqrt(h/h_ς)
    dmdt = md(t) - md_e   
end
#
ρ, A, K, h_ς = 1, 5, 5, 3
#
u0 = 1.5*ρ*S
md(t) = 2
par = [ρ, A, K, h_ς, md(t)]
#
t_fin = 20
tspan = (0.0, t_fin)
prob = ODEProblem(tank!, u0, tspan, par)
sol_c = solve(prob)

```

Alternatively, I can do an MTK implementation:

```julia
@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
#
@mtkbuild tank = Tank()
#
tspan = (0,t_fin)
#
prob = ODEProblem(tank, [], tspan)
sol = solve(prob)

plot(sol,lw=LW1, lc=LC1, label="m, MTK")
plot!(sol_c, lw=LW1, lc=LC2, ls=LS2, label="m, ODE")
plot!(title="Tank model")

```

The results are, of course, the same:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/d/cddc5ed6ceb33c61f69c0e67b11fdfa094d77031.png)

Here, steady state is `m_ss = 2.4`.

* * *

**Finite difference linearization** :

```julia
julia> m_ss = 2.4
julia> dxdt_0 = tank!(m_ss,0,par)
julia> deps = 1e-5 
julia> A = (tank!(m_ss+deps, 0, par) - dxdt_0)/deps
-0.4166662326277048

```

**Alternative: complex theory linearization:**

```julia
julia> deps0 = eps(Float64)
julia> A = tank!(m_ss + deps0*im,0,par)/deps0 |> imag
-0.41666666663862956

```

**Question** :

- Is there a way to _convert_ the MTK-model (`tank` instantiation above) into an “ODE-type” function like `tank!` which I can do finite difference or complex theory linearization on?
- What additional problems would there be if the MTK model contains algebraic equations in addition to differential equations?

NOTE: for this simple model (`tank`), the MTK function `linearize` works. But I get error messages for more complicated models.

---

<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:** [September 2, 2024, 9:44am UTC](https://discourse.julialang.org/t/modelingtoolkit-possible-to-do-numeric-linearization/118933/2 "2024-09-02T09:44:14Z")

</div>

When you have created an `ODEProblem`, the function f that computes the right-hand side  
\dot x = f(dx, x, p, t) is available as `prob.f`.

> What additional problems would there be if the MTK model contains algebraic equations in addition to differential equations?

Your finite perturbation could (probably would) violate the algebraic equations.

> [@BLI](#):
>
> NOTE: for this simple model (`tank`), the MTK function `linearize` works. But I get error messages for more complicated models.

How about trying to solve those problems rather than using finite differences?

---

<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:** [September 2, 2024, 10:01am UTC](https://discourse.julialang.org/t/modelingtoolkit-possible-to-do-numeric-linearization/118933/3 "2024-09-02T10:01:43Z")

</div>

> [@baggepinnen](#):
>
> How about trying to solve those problems rather than using finite differences?

That would be the best solution, of course. The reported bug ( #2941) is still pending, I think. The error message may have changed slightly, though. [Tested with all packages I used updated today.]

```julia
BoundsError: attempt to access Tuple{SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, NonlinearLeastSquaresProblem{Nothing, true, MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}}, NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#636"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x2c6790d7, 0x190a4316, 0x8cb7cea9, 0xf28b1af7, 0x85751e05), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8c32180b, 0xaafd1b53, 0xcfdb1533, 0xa7a00f2d, 0xd64d3dc5), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{NonlinearSystem}, Nothing, NonlinearSystem, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing}} at index [2]

```

---

<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:** [September 2, 2024, 11:16am UTC](https://discourse.julialang.org/t/modelingtoolkit-possible-to-do-numeric-linearization/118933/4 "2024-09-02T11:16:02Z")

</div>

Hopefully fixed on this PR branch

> <https://github.com/SciML/ModelingToolkit.jl/pull/3006>
>
> Close #2941
> 
> \## Checklist
> 
> \- \[\] Appropriate tests were added
> \- \[\] Any cod…e changes were done in a way that does not break public API
> \- \[\] All documentation related to code changes were updated
> \- \[\] The new code follows the
> \[contributor guidelines\](https://github.com/SciML/.github/blob/master/CONTRIBUTING.md), in particular the \[SciML Style Guide\](https://github.com/SciML/SciMLStyle) and
> \[COLPRAC\](https://github.com/SciML/COLPRAC).
> \- \[\] Any new documentation only uses public API
>   
> \## Additional context
> 
> Add any other context about the problem here.

---

<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:** [September 2, 2024, 11:24am UTC](https://discourse.julialang.org/t/modelingtoolkit-possible-to-do-numeric-linearization/118933/5 "2024-09-02T11:24:24Z")

</div>

Will check when it is pushed out. Hopefully, that will also fix the previous problem in #2888 [I haven’t checked that in a while].

---

<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:** [September 2, 2024, 11:29am UTC](https://discourse.julialang.org/t/modelingtoolkit-possible-to-do-numeric-linearization/118933/6 "2024-09-02T11:29:25Z")

</div>

This fix will probably not fix the other issue, which appears to be in the structural simplification rather than the code generation. #2888 does not have a reproducing example making it hard to fix.

---

<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:** [September 2, 2024, 12:00pm UTC](https://discourse.julialang.org/t/modelingtoolkit-possible-to-do-numeric-linearization/118933/7 "2024-09-02T12:00:22Z")

</div>

I tried to make a reproducible example for #2888, but my attempt “failed” in the sense that my slightly simplified example worked…

I can share the code that doesn’t work (as reported in #2888) to make it simpler to fix, but would prefer to not share it in the open – I can share it with you, and you can share it with the one who is tasked to look into it.

Would that be ok?