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

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

Could you share the whole model?

I found the reason. I forgot one equation.

1 Like