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