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

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:

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:

@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:

Here, steady state is m_ss = 2.4.

Finite difference linearization:

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

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.

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

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]


Hopefully fixed on this PR branch

1 Like

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

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.

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?