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 (
tankinstantiation above) into an “ODE-type” function liketank!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.