# ModelingToolkit.jl hydraulic model

@ChrisRackauckas Thx for `ModelingToolkit.jl` I really like it!

I tried to model a simple hydraulic system. However, I soon ran into some problems.

Consider the following example (adaptation of the RC example):

``````using ModelingToolkit, Plots, DifferentialEquations

@parameters t

function Pin(; name)
@variables p(t) mdot(t)
return ODESystem(Equation[], t, [p, mdot], []; name=name, default_u0=[p => 1e5, mdot => 0.0])
end

function Reservoir(; name, p=1e5)
p_val = p
@parameters p
@named a = Pin()
eqs = [
a.p ~ p
]
return ODESystem(eqs, t, [], [p]; systems=[a], default_p=Dict(p => p_val), name=name)
end

@register sign(x)

function FixedFlowResistance(; name, A=1e-7)
A_val = A
@parameters A # effective area of valve
@named a = Pin()
@named b = Pin()
@variables Î”p(t)

rho = 998.0 # Density of water

eqs = [
Î”p ~ b.p - a.p
a.mdot + b.mdot ~ 0.0 # Mass balance
b.mdot ~ A * sqrt(2.0 * rho) * sqrt(abs(Î”p)) * sign(Î”p) # Momentum balance
]

return ODESystem(eqs, t, [Î”p], [A]; systems=[a, b], default_p=Dict(A => A_val), name=name)
end

function connect_pins(pins...)
# mass balance
eqs = [
sum(pin -> pin.mdot, pins) ~ 0.0, # sum(mdot) = 0
]

# pressure balance
for pin in pins[2:end]
push!(eqs, pins[1].p ~ pin.p) # p_1 = p_2, p_1 = p_3,...
end

return eqs
end
``````

I than assemble a model with two reservoirs at different pressures and a fixed flow resistance in between. Since the pressure at the left side is higher than on the right, the mass flow rate should be into the valve (i.e. positive) at pin `a` and negative at pin `b`.
However, I do not get that far because of an error when calling `structural_simplify`.

``````@named reservoir_left = Reservoir(; p=2e3)
@named valve = FixedFlowResistance(; A=0.7 * 1e-4)
@named reservoir_right = Reservoir(; p=1e3)
connections = [
connect_pins(reservoir_left.a, valve.a)
connect_pins(valve.b, reservoir_right.a)
]
@named model = ODESystem(connections, t; systems=[reservoir_left, valve, reservoir_right])
sys = structural_simplify(model)

ERROR: TypeError: non-boolean (Term{Real,Nothing}) used in boolean context
Stacktrace:
[1] ifelse(::Term{Real,Nothing}, ::Int64, ::Int64) at .julia\packages\IfElse\u1l5W\src\IfElse.jl:3
[2] derivative(::typeof(abs), ::Tuple{Term{Real,Nothing}}, ::Val{1}) at .julia\packages\Symbolics\ZuF87\src\extra_functions.jl:7
[3] derivative_idx(::Term{Real,Nothing}, ::Int64) at .julia\packages\Symbolics\ZuF87\src\diff.jl:204
[4] expand_derivatives(::Term{Real,Nothing}, ::Bool; occurances::Term{Real,Nothing}) at .julia\packages\Symbolics\ZuF87\src\diff.jl:0
[5] expand_derivatives(::Term{Real,Nothing}, ::Bool; occurances::Term{Real,Nothing}) at .julia\packages\Symbolics\ZuF87\src\diff.jl:118 (repeats 3 times)
[6] expand_derivatives(::Term{Real,Nothing}, ::Bool) at .julia\packages\Symbolics\ZuF87\src\diff.jl:78
[7] find_linear_equations(::ODESystem) at .julia\packages\ModelingToolkit\0BgwD\src\systems\systemstructure.jl:194
[8] alias_elimination(::ODESystem) at .julia\packages\ModelingToolkit\0BgwD\src\systems\alias_elimination.jl:13
[9] structural_simplify(::ODESystem) at .julia\packages\ModelingToolkit\0BgwD\src\systems\abstractsystem.jl:483
[10] top-level scope at REPL[15]:1
``````

I guess `Symbolics.jl` does not like the `sign` function.
Is it possible to have `ifelse` in the equations of a system?

I also tried a smooth approximation of the `sign(x) * sqrt(abs(x))` term:

``````# Smooth approximation of sign(x) * sqrt(abs(x)), see: https://www.maplesoft.com/documentation_center/online_manuals/modelica/Modelica_Fluid_Utilities.html#Modelica.Fluid.Utilities.regRoot
regRoot(x, delta) = x / (x^2 + delta^2)^0.25
@register regRoot(x, delta)

function FixedFlowResistance(; name, A=1e-7)
A_val = A
@parameters A # effective area of valve
@named a = Pin()
@named b = Pin()
@variables Î”p(t) p_cr(t)

rho = 998.0 # Density of water
p_atm = 101325.0 # Atmospheric pressure
B_lam = 0.999 # Laminar flow pressure ratio

eqs = [
p_cr ~ (a.p + b.p + 2.0 * p_atm) / 2.0 * (1.0 - B_lam) # critical pressure
Î”p ~ b.p - a.p
a.mdot + b.mdot ~ 0.0 # Mass balance
b.mdot ~ A * sqrt(2.0 * rho) * regRoot(Î”p, p_cr) # Momentum balance
]

return ODESystem(eqs, t, [Î”p, p_cr], [A]; systems=[a, b], default_p=Dict(A => A_val), name=name)
end
``````

but now I get another error:

``````@named reservoir_left = Reservoir(; p=2e3)
@named valve = FixedFlowResistance(; A=0.7 * 1e-4)
@named reservoir_right = Reservoir(; p=1e3)
connections = [
connect_pins(reservoir_left.a, valve.a)
connect_pins(valve.b, reservoir_right.a)
]
@named model = ODESystem(connections, t; systems=[reservoir_left, valve, reservoir_right])
sys = structural_simplify(model)

ERROR: InexactError: Int64(0.0005000000000000004)
Stacktrace:
[1] Int64 at .\float.jl:710 [inlined]
[2] convert at .\number.jl:7 [inlined]
[3] push!(::Array{Int64,1}, ::Float64) at .\array.jl:934
[4] find_linear_equations(::ODESystem) at .julia\packages\ModelingToolkit\0BgwD\src\systems\systemstructure.jl:203
[5] alias_elimination(::ODESystem) at .julia\packages\ModelingToolkit\0BgwD\src\systems\alias_elimination.jl:13
[6] structural_simplify(::ODESystem) at
.julia\packages\ModelingToolkit\0BgwD\src\systems\abstractsystem.jl:483
[7] top-level scope at REPL[9]:1
``````

In the documentation you write that: â€śThis system could be solved directly as a DAE using one of the DAE solvers from DifferentialEquations.jlâ€ť. However, I am not sure how to obtain a `DAEProblem` from the `ODESystem`.

I am aware that this system has no dynamics but I tried a simple voltage divider circuit (which also has no dynamics) with the stuff from the RC example and it worked fine.

``````using ModelingToolkit, Plots, DifferentialEquations

@parameters t

function Pin(; name)
@variables v(t) i(t)
return ODESystem(Equation[], t, [v, i], []; name=name, default_u0=[v => 1.0, i => 1.0])
end

function Reference(; name)
@named p = Pin()
eqs = [p.v ~ 0]
return ODESystem(eqs, t, [], []; systems=[p], name=name)
end

function Resistor(; name, R=1.0)
R_val = R
@named p = Pin()
@named n = Pin()
@variables v(t)
@parameters R
eqs = [
v ~ p.v - n.v
0 ~ p.i + n.i
v ~ p.i * R
]
return ODESystem(eqs, t, [v], [R]; systems=[p, n], default_p=Dict(R => R_val), name=name)
end

function ConstantVoltage(; name, V=1.0)
V_val = V
@named p = Pin()
@named n = Pin()
@parameters V
eqs = [
V ~ p.v - n.v
0 ~ p.i + n.i
]
return ODESystem(eqs, t, [], [V]; systems=[p, n], default_p=Dict(V => V_val), name=name)
end

function connect_pins(ps...)
eqs = [
0 ~ sum(p -> p.i, ps), # KCL
]
# KVL
for i in 1:(length(ps) - 1)
push!(eqs, ps[i].v ~ ps[i + 1].v)
end

return eqs
end

@named R1 = Resistor(; R=1e3)
@named R2 = Resistor(; R=2e3)
@named source = ConstantVoltage(; V=10)
@named ref = Reference()

connections = [
connect_pins(source.p, R1.p)
connect_pins(R1.n, R2.p)
connect_pins(R2.n, source.n, ref.p)
]
@named model = ODESystem(connections, t; systems=[R1, R2, source, ref])

sys = structural_simplify(model)
u0 = [
R2.p.i => 0.0,
]
tspan = (0.0, 10.0)
prob = ODEProblem(sys, u0, tspan)
sol = solve(prob, Rodas4())

plot(sol, vars=[source.p.v, R1.v, R2.v], ylims=(-0.3, 10.3))
``````

Sry, for the lengthy question.

2 Likes

Hey, thanks for trying MTK out. Could you update your packages and try the second example again? I fixed that problem yesterday afternoon. Fix linear equation finding by YingboMa Â· Pull Request #880 Â· SciML/ModelingToolkit.jl Â· GitHub

As for the problem with `sign` and `abs`. I can fix that soon. However, do note that robust handling of non-smooth functions requires event handling. That feature is planed. Support for events Â· Issue #38 Â· SciML/ModelingToolkit.jl Â· GitHub So at the moment, I still would recommend to use a smooth approximation if thatâ€™s suitable in your model.

2 Likes

Thanks!
I tried that and now I get:

``````julia> sys = structural_simplify(model)
Model model with 0 equations
States (0):
Parameters (3):
reservoir_leftâ‚Šp [defaults to 2000.0]
valveâ‚ŠA [defaults to 7.0e-5]
reservoir_rightâ‚Šp [defaults to 1000.0]
Incidence matrix:0Ă—0 SparseArrays.SparseMatrixCSC{Num,Int64} with 0 stored entries
``````

I did get it to work with:

``````sys = alias_elimination(model)
u0 = [
valve.a.p => 2e5,
valve.b.mdot => 0.0,
reservoir_right.a.p => 1e5,
]
tspan = (0.0, 10.0)
prob = ODEProblem(sys, u0, tspan)
sol = solve(prob, Rodas4())
``````
1 Like

This means that `structural_simplify` is able to simply all the equations. You can check `observed` for the solved expressions. If you want to evaluate it, then you can use substitute. You can also obtain the observed variables by

``````julia> sys = structural_simplify(model);

julia> prob = ODEProblem(sys, Pair[], tspan);

julia> prob.f.observed(valve.a.p, prob.u0, prob.p, 0.0)
2000.0

julia> prob.f.observed(reservoir_right.a.p, prob.u0, prob.p, 0.0)
1000.0
``````

Thank you very much!

@YingboMa
I encountered another issue.
This time i was modeling a pipe with two states: temperature and massflow rate along the pipe.

``````julia> sys = structural_simplify(model)
Model model with 2 equations
States (3):
pipeâ‚Šmdot(t) [defaults to 0.0]
pipeâ‚ŠT(t) [defaults to 323.15]
pipeâ‚Šbâ‚Šmdot(t) [defaults to 0.0]
Parameters (4):
reservoir_leftâ‚Šp [defaults to 200000.0]
reservoir_leftâ‚ŠT [defaults to 343.15]
reservoir_rightâ‚Šp [defaults to 100000.0]
reservoir_rightâ‚ŠT [defaults to 323.15]
Incidence matrix:
[1, 1]  =  Ă—
[2, 2]  =  Ă—
[1, 3]  =  Ă—
[2, 4]  =  Ă—
[2, 5]  =  Ă—

julia> prob = ODAEProblem(sys, Pair[], (0, 10.0))
Stacktrace:
[1] getindex at .\dict.jl:467 [inlined]
[2] torn_system_jacobian_sparsity(::ODESystem) at .julia\packages\ModelingToolkit\GBGLe\src\structural_transformation\codegen.jl:81
[3] build_torn_function(::ODESystem; expression::Bool, jacobian_sparsity::Bool, checkbounds::Bool, kw::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at .julia\packages\ModelingToolkit\GBGLe\src\structural_transformation\codegen.jl:230
[4] build_torn_function at .julia\packages\ModelingToolkit\GBGLe\src\structural_transformation\codegen.jl:187 [inlined]
[5] ODAEProblem{true}(::ODESystem, ::Array{Pair,1}, ::Tuple{Int64,Float64}, ::SciMLBase.NullParameters; kw::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at .julia\packages\ModelingToolkit\GBGLe\src\structural_transformation\codegen.jl:358
[6] ODAEProblem at .julia\packages\ModelingToolkit\GBGLe\src\structural_transformation\codegen.jl:338 [inlined] (repeats 2 times)
[7] ODAEProblem(::ODESystem, ::Vararg{Any,N} where N; kw::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at .julia\packages\ModelingToolkit\GBGLe\src\structural_transformation\codegen.jl:330
[8] ODAEProblem(::ODESystem, ::Vararg{Any,N} where N) at .julia\packages\ModelingToolkit\GBGLe\src\structural_transformation\codegen.jl:330
[9] top-level scope at REPL[42]:1
``````

Could you share the whole model?

I found the reason. I forgot one equation.

1 Like