# ModelingToolkit.jl hydraulic model

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084>\
**Category:** Modelling & Simulations\
**Created:** [March 13, 2021, 4:39pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084 "2021-03-13T16:39:33Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [March 13, 2021, 4:39pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084/1 "2021-03-13T16:39:33Z")

</div>

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

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

```julia
@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](https://www.maplesoft.com/documentation_center/online_manuals/modelica/Modelica_Fluid_Utilities.html#Modelica.Fluid.Utilities.regRoot) of the `sign(x) * sqrt(abs(x))` term:

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

```julia
@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](https://mtk.sciml.ai/dev/tutorials/acausal_components/#Simplifying-and-Solving-this-System-1) 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.

```julia
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.  
And thx in advance!

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [March 13, 2021, 5:04pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084/2 "2021-03-13T17:04:40Z")

</div>

Hey, thanks for trying MTK out. Could you update your packages and try the second example again? I fixed that problem yesterday afternoon. [https://github.com/SciML/ModelingToolkit.jl/pull/880](https://github.com/SciML/ModelingToolkit.jl/pull/880)

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. [https://github.com/SciML/ModelingToolkit.jl/issues/38](https://github.com/SciML/ModelingToolkit.jl/issues/38) So at the moment, I still would recommend to use a smooth approximation if that’s suitable in your model.

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [March 13, 2021, 6:29pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084/3 "2021-03-13T18:29:50Z")

</div>

Thanks!  
I tried that and now I get:

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

```

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [March 13, 2021, 6:37pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084/4 "2021-03-13T18:37:14Z")

</div>

I did get it to work with:

```julia
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())

```

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [March 13, 2021, 6:58pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084/5 "2021-03-13T18:58:38Z")

</div>

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

```

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [March 13, 2021, 7:15pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084/6 "2021-03-13T19:15:11Z")

</div>

Thank you very much!

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [March 14, 2021, 10:43am UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084/7 "2021-03-14T10:43:18Z")

</div>

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

```julia
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))
ERROR: KeyError: key 5 not found
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

```

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [March 14, 2021, 8:15pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084/8 "2021-03-14T20:15:14Z")

</div>

Could you share the whole model?

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [March 15, 2021, 3:36pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-hydraulic-model/57084/9 "2021-03-15T15:36:49Z")

</div>

I found the reason. I forgot one equation.