From ddcb4f9679354e3c6cf5c4b0db66382796de1834 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 24 Apr 2023 19:42:51 -0400 Subject: [PATCH 01/16] Actuator system --- src/Blocks/Blocks.jl | 2 +- src/Blocks/sources.jl | 91 +++++++ .../IsothermalCompressible.jl | 4 +- .../IsothermalCompressible/components.jl | 234 +++++++++++++++--- src/Hydraulic/IsothermalCompressible/utils.jl | 25 +- src/Mechanical/Translational/sources.jl | 25 +- test/Hydraulic/isothermal_compressible.jl | 76 +++++- 7 files changed, 411 insertions(+), 46 deletions(-) diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl index 7997034b7..36a7a008c 100644 --- a/src/Blocks/Blocks.jl +++ b/src/Blocks/Blocks.jl @@ -17,7 +17,7 @@ export Log, Log10 include("math.jl") export Constant, TimeVaryingFunction, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine, - Square, Triangular + Square, Triangular, InputMemory, Input include("sources.jl") export Limiter, DeadZone, SlewRateLimiter diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 1cefb5ec1..a2273c48f 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -411,3 +411,94 @@ end # - Pulse Generate pulse signal of type Real # - SawTooth Generate saw tooth signal # - Trapezoid Generate trapezoidal signal of type Real + +function linear_interpolation(x1::T, x2::T, t1::T, t2::T, t::T) where {T<:Real} + if t1 != t2 + slope = (x2 - x1) / (t2 - t1) + intercept = x1 - slope * t1 + + return slope * t + intercept + else + @assert x1==x2 "x1 ($x1) and x2 ($x2) should be equal if t1 == t2" + + return x2 + end +end + +struct InputMemory{T<:Real} + data::Vector{T} + sample_time::T + n::Int +end + +Base.copy(x::InputMemory{T}) where {T} = InputMemory{T}(copy(x.data), x.sample_time, x.n) + +function InputMemory(data::Vector{T}, sample_time::T) where {T<:Real} + InputMemory{T}(data, sample_time, length(data)) +end + +function input(t, memory::InputMemory) + if t < 0 + t = zero(t) + end + + i1 = floor(Int, t / memory.sample_time) + 1 #expensive + i2 = i1 + 1 + + t1 = i1 * memory.sample_time + t2 = i2 * memory.sample_time + + #println("input: t=$t, i1=$i1, i2=$i2, t1=$t1, t2=$t2, Δt = $(memory.sample_time)") + + x1 = memory.data[i1 > memory.n ? memory.n : i1] + x2 = memory.data[i2 > memory.n ? memory.n : i2] + + return linear_interpolation(x1, x2, t1, t2, t) +end + +get_sample_time(memory::InputMemory) = memory.sample_time +Symbolics.@register_symbolic get_sample_time(memory) + +Symbolics.@register_symbolic input(t, memory) + +function first_order_backwards_difference(t, memory) + Δt = get_sample_time(memory) + x1 = input(t, memory) + x0 = input(t - Δt, memory) + + return (x1 - x0) / Δt +end + +function Symbolics.derivative(::typeof(input), args::NTuple{2, Any}, ::Val{1}) + first_order_backwards_difference(args[1], args[2]) +end + + + +Input(T::Type; name) = Input(T[], zero(T); name) +Input(data::Vector{T}, dt::T; name) where T<:Real = Input(; name, buffer = InputMemory(data, dt)) + +""" + Input(; name, buffer) + +data input component. + +# Parameters: + - `buffer`: an `InputMemory` parameter which holds the data and sample time + +# Connectors: + - `output` +""" +@component function Input(; name, buffer) + pars = @parameters begin + buffer = buffer + end + vars = [] + systems = @named begin + output = RealOutput() + end + eqs = [ + output.u ~ input(t, buffer) + ] + return ODESystem(eqs, t, vars, pars; name, systems) +end diff --git a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl index f2f7eb2e8..6f58bed1d 100644 --- a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl +++ b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl @@ -5,7 +5,7 @@ module IsothermalCompressible using ModelingToolkit, Symbolics, IfElse using ...Blocks: RealInput, RealOutput -using ...Mechanical.Translational: MechanicalPort +using ...Mechanical.Translational: MechanicalPort, Mass @parameters t D = Differential(t) @@ -13,7 +13,7 @@ D = Differential(t) export HydraulicPort, HydraulicFluid include("utils.jl") -export Source, InputSource, Cap, Pipe, FixedVolume, DynamicVolume +export Source, InputSource, Cap, Tube, FixedVolume, DynamicVolume include("components.jl") end diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index 98c90a370..d5d2bf043 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -96,7 +96,7 @@ Fixed fluid volume. systems = @named begin port = HydraulicPort(; p_int) end vars = @variables begin - rho(t) = density(port, p_int) + rho(t) = liquid_density(port) drho(t) = 0 end @@ -104,29 +104,29 @@ Fixed fluid volume. dm = port.dm eqs = [D(rho) ~ drho - rho ~ density(port, port.p) + rho ~ full_density(port) dm ~ drho * vol] ODESystem(eqs, t, vars, pars; name, systems) end """ - PipeBase(; p_int, area, length, perimeter=2*sqrt(area*pi), shape_factor=64, name) + TubeBase(; p_int, area, length, perimeter=2*sqrt(area*pi), shape_factor=64, name) -Pipe segement which models purely the fully developed flow friction, ignoring any compressibility. +Internal flow model of the fully developed flow friction, ignoring any compressibility. # Parameters: -- `p_int`: [Pa] initial pressure (set by `p_int` argument) -- `area`: [m^2] tube cross sectional area (set by `area` argument) -- `length`: [m] length of the pipe (set by `length` argument) -- `perimeter`: [m] perimeter of the pipe cross section (set by optional `perimeter` argument, needed only for non-circular pipes) +- `p_int`: [Pa] initial pressure +- `area`: [m^2] tube cross sectional area +- `length`: [m] length of the pipe +- `perimeter`: [m] perimeter of the pipe cross section (needed only for non-circular pipes) - `Φ`: shape factor, see `friction_factor` function (set by optional `shape_factor` argument, needed only for non-circular pipes). # Connectors: - `port_a`: hydraulic port - `port_b`: hydraulic port """ -@component function PipeBase(; p_int, area, length, perimeter = 2 * sqrt(area * pi), +@component function TubeBase(; p_int, area, length, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name) pars = @parameters begin p_int = p_int @@ -149,7 +149,7 @@ Pipe segement which models purely the fully developed flow friction, ignoring an d_h = 4 * area / perimeter - ρ = (density(port_a, port_a.p) + density(port_b, port_b.p)) / 2 + ρ = (full_density(port_a) + full_density(port_b)) / 2 μ = viscosity(port_a) f = friction_factor(dm, area, d_h, ρ, μ, Φ) @@ -162,30 +162,33 @@ Pipe segement which models purely the fully developed flow friction, ignoring an end """ - Pipe(N; p_int, area, length, perimeter=2*sqrt(area*pi), shape_factor=64, name) + Tube(N; p_int, area, length, effective_length=length, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name) -Pipe modeled with `N` segements which models the fully developed flow friction and compressibility. +Tube modeled with `N` segements which models the fully developed flow friction and compressibility. # Parameters: -- `p_int`: [Pa] initial pressure (set by `p_int` argument) -- `area`: [m^2] tube cross sectional area (set by `area` argument) -- `length`: [m] length of the pipe (set by `length` argument) -- `perimeter`: [m] perimeter of the pipe cross section (set by optional `perimeter` argument, needed only for non-circular pipes) -- `Φ`: shape factor, see `friction_factor` function (set by optional `shape_factor` argument, needed only for non-circular pipes). +- `p_int`: [Pa] initial pressure +- `area`: [m^2] tube cross sectional area +- `length`: [m] real length of the tube +- `effective_length`: [m] tube length to account for flow development and other restrictions, used to set the `TubeBase` length which calculates the flow resistance +- `perimeter`: [m] perimeter of the tube cross section (set by optional `perimeter` argument, needed only for non-circular tubes) +- `Φ`: shape factor, see `friction_factor` function (set by optional `shape_factor` argument, needed only for non-circular tubes) # Connectors: - `port_a`: hydraulic port - `port_b`: hydraulic port """ -@component function Pipe(N; p_int, area, length, perimeter = 2 * sqrt(area * pi), +@component function Tube(N; p_int, area, length, effective_length=length, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name) @assert(N>1, - "the pipe component must be defined with more than 1 segment (i.e. N>1), found N=$N") + "the Tube component must be defined with more than 1 segment (i.e. N>1), found N=$N") + #TODO: How to set an assert effective_length >= length ?? pars = @parameters begin p_int = p_int area = area length = length + effective_length=effective_length perimeter = perimeter Φ = shape_factor end @@ -199,9 +202,9 @@ Pipe modeled with `N` segements which models the fully developed flow friction a pipe_bases = [] for i in 1:(N - 1) - x = PipeBase(; name = Symbol("p$i"), shape_factor = ParentScope(Φ), + x = TubeBase(; name = Symbol("p$i"), shape_factor = ParentScope(Φ), p_int = ParentScope(p_int), area = ParentScope(area), - length = ParentScope(length) / (N - 1), + length = ParentScope(effective_length) / (N - 1), perimeter = ParentScope(perimeter)) push!(pipe_bases, x) end @@ -225,6 +228,44 @@ Pipe modeled with `N` segements which models the fully developed flow friction a ODESystem(eqs, t, vars, pars; name, systems = [ports; pipe_bases; volumes]) end +""" + FlowDivider(;p_int, n, name) + +Reduces the flow from `port_a` to `port_b` by `n`. Useful for modeling parallel tubes efficiently by placing a `FlowDivider` on each end of a tube. + +# Parameters: +- `p_int`: [Pa] initial pressure +- `n`: divide flow from `port_a` to `port_b` by `n` + +# Connectors: +- `port_a`: full flow hydraulic port +- `port_b`: part flow hydraulic port +""" +@component function FlowDivider(;p_int, n, name) + + #TODO: assert n >= 1 + + pars = @parameters begin + n = n + p_int = p_int + end + + vars = [] + + systems = @named begin + port_a = HydraulicPort(; p_int) + port_b = HydraulicPort(; p_int) + end + + eqs = [ + port_b.dm ~ port_a.dm/n + port_b.p ~ port_a.p + ] + + ODESystem(eqs, t, vars, pars; name, systems) +end + + """ DynamicVolume(; p_int, x_int=0, area, dead_volume=0, direction=+1, minimum_volume=0, name) @@ -240,6 +281,13 @@ dm ────► dead volume │ │ area │ └─► x (= flange.v * direction) ``` +# States: +- `x(t)`: [m] moving wall position +- `dx(t)`: [m/s] moving wall velocity +- `rho(t)`: [kg/m^3] density +- `drho(t)`: [kg/s-m^3] density derivative +- `vol(t)`: [m^3] volume +- `p(t)`: [Pa] dynamic pressure # Parameters: - `p_int`: [Pa] initial pressure @@ -271,20 +319,24 @@ dm ────► dead volume │ │ area vars = @variables begin x(t) = x_int dx(t) = 0 - rho(t) = density(port, p_int) + rho(t) = liquid_density(port) drho(t) = 0 vol(t) = dead_volume + area * x_int - p(t) = p_int + p(t) = p_int # p represents the dynamic pressure end - eqs = [0 ~ IfElse.ifelse(vol >= minimum_volume, p - port.p, port.dm) + # let + dm = port.dm + u = dm/(rho*area) + + eqs = [0 ~ IfElse.ifelse(vol >= minimum_volume, ( p ) - ( port.p - (1/2)*rho*u^2 ), dm) vol ~ dead_volume + area * x D(x) ~ dx D(rho) ~ drho dx ~ flange.v * direction - rho ~ density(port, p) - port.dm ~ drho * vol + rho * area * dx - flange.f ~ -p * area * direction] #TODO: update to dynamic pressure + rho ~ full_density(port) + dm ~ drho * vol + rho * area * dx + flange.f ~ -p * area * direction] ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0]) end @@ -305,7 +357,7 @@ Valve with area input and discharge coefficient `Cd` defined by https://en.wikip - `port_b`: hydraulic port - `input`: real input setting the valve `area`. Note: absolute value taken """ -@component function Valve(; p_a_int, p_b_int, area_int, Cd, name) +@component function Valve(reversible=false; p_a_int, p_b_int, area_int, Cd, name) pars = @parameters begin p_a_int = p_a_int p_b_int = p_b_int @@ -322,13 +374,135 @@ Valve with area input and discharge coefficient `Cd` defined by https://en.wikip vars = [] # let - ρ = (density(port_a, port_a.p) + density(port_b, port_b.p)) / 2 + ρ = (full_density(port_a) + full_density(port_b)) / 2 Δp = port_a.p - port_b.p dm = port_a.dm - area = abs(input.u) + area = if reversible + input.u + else + input.u*(input.u > 0) + end eqs = [sign(Δp) * dm ~ sqrt(2 * abs(Δp) * ρ / Cd) * area 0 ~ port_a.dm + port_b.dm] ODESystem(eqs, t, vars, pars; name, systems, defaults = [input.u => area_int]) end + + +@component function SpoolValve(reversible=false;p_a_int, p_b_int, x_int, Cd, d, name) + + pars = @parameters begin + p_a_int = p_a_int + p_b_int = p_b_int + d = d + x_int = x_int + Cd = Cd + end + + systems = @named begin + port_a = HydraulicPort(; p_int = p_a_int) + port_b = HydraulicPort(; p_int = p_b_int) + flange = MechanicalPort() + valve = Valve(reversible; p_a_int, p_b_int, area_int=ParentScope(x_int)*2π*ParentScope(d), Cd) + end + + vars = @variables begin + x(t) = x_int + dx(t) = 0 + end + + eqs = [ + D(x) ~ dx + flange.v ~ dx + flange.f ~ 0 #TODO: model flow force + connect(valve.port_a, port_a) + connect(valve.port_b, port_b) + valve.input.u ~ x*2π*d + ] + + ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0]) +end + +@component function SpoolValve2Way(reversible=false; p_s_int, p_a_int, p_b_int, p_r_int, m, g, x_int, Cd, d, name) + + pars = @parameters begin + p_s_int = p_s_int + p_a_int = p_a_int + p_b_int = p_b_int + p_r_int = p_r_int + + m=m + g=g + + x_int = x_int + + d=d + end + + vars = [] + + systems = @named begin + vSA = SpoolValve(reversible; p_a_int=p_s_int, p_b_int=p_a_int, x_int, Cd, d) + vBR = SpoolValve(reversible; p_a_int=p_b_int, p_b_int=p_r_int, x_int, Cd, d) + + port_s = HydraulicPort(; p_int = p_s_int) + port_a = HydraulicPort(; p_int = p_a_int) + port_b = HydraulicPort(; p_int = p_b_int) + port_r = HydraulicPort(; p_int = p_r_int) + + mass = Mass(; m=m, g=g) + + flange = MechanicalPort() + end + + eqs = [ + connect(vSA.port_a, port_s) + connect(vSA.port_b, port_a) + connect(vBR.port_a, port_b) + connect(vBR.port_b, port_r) + + connect(vSA.flange, vBR.flange, mass.flange, flange) + ] + + ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0]) +end + +@component function Actuator(; p_a_int, p_b_int, area_a, area_b, length_a_int, length_b_int, m, g, x_int = 0, minimum_volume_a = 0, minimum_volume_b = 0, name) + + pars = @parameters begin + p_a_int = p_a_int + p_b_int = p_b_int + area_a = area_a + area_b = area_b + x_int = x_int + length_a_int = length_a_int + length_b_int = length_b_int + minimum_volume_a = minimum_volume_a + minimum_volume_b = minimum_volume_b + end + + vars = @variables begin + x(t) = x_int + dx(t) = 0 + end + + systems = @named begin + vol_a = DynamicVolume(; p_int = p_a_int, x_int = +x_int, area=area_a, dead_volume = length_a_int*area_a, minimum_volume=minimum_volume_a, direction = +1) + vol_b = DynamicVolume(; p_int = p_b_int, x_int = -x_int, area=area_b, dead_volume = length_b_int*area_b, minimum_volume=minimum_volume_b, direction = -1) + mass = Mass(; m, g) + port_a = HydraulicPort(; p_int = p_a_int) + port_b = HydraulicPort(; p_int = p_b_int) + flange = MechanicalPort() + end + + eqs = [ + connect(vol_a.port, port_a) + connect(vol_b.port, port_b) + connect(vol_a.flange, vol_b.flange, mass.flange, flange) + x ~ vol_a.x + dx ~ vol_a.dx + ] + + ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0]) +end \ No newline at end of file diff --git a/src/Hydraulic/IsothermalCompressible/utils.jl b/src/Hydraulic/IsothermalCompressible/utils.jl index 4d112e27d..603d53431 100644 --- a/src/Hydraulic/IsothermalCompressible/utils.jl +++ b/src/Hydraulic/IsothermalCompressible/utils.jl @@ -8,7 +8,7 @@ Connector port for hydraulic components. - `p_int`: [Pa] initial gauge pressure # States: -- `p`: [Pa] gauge pressure +- `p`: [Pa] gauge total pressure - `dm`: [kg/s] mass flow """ @connector function HydraulicPort(; p_int, name) @@ -16,6 +16,9 @@ Connector port for hydraulic components. ρ β μ + n + ρ_gas + p_gas end vars = @variables begin @@ -29,20 +32,26 @@ end """ HydraulicFluid(; density=997, bulk_modulus=2.09e9, viscosity=0.0010016, name) -Fluid parameter setter for isothermal compressible fluid domain. Defaults given for water at 20°C and 1atm reference pressure. +Fluid parameter setter for isothermal compressible fluid domain. Defaults given for water at 20°C and 0Pa gage (1atm absolute) reference pressure. Density is modeled using the Tait equation of state. For pressures below the reference pressure, density is linearly interpolated to the gas state, this helps prevent pressures from going below the reference pressure. # Parameters: -- `ρ`: [kg/m^3] fluid density at reference pressure (set by `density` argument) +- `ρ`: [kg/m^3] fluid density at 0Pa reference gage pressure (set by `density` argument) - `Β`: [Pa] fluid bulk modulus describing the compressibility (set by `bulk_modulus` argument) - `μ`: [Pa*s] or [kg/m-s] fluid dynamic viscosity (set by `viscosity` argument) +- `n`: density exponent +- `ρ_gas`: [kg/m^3] density of fluid in gas state at reference gage pressure `p_gas` (set by `gas_density` argument) +- `p_gas`: [Pa] reference pressure (set by `gas_pressure` argument) """ @connector function HydraulicFluid(; density = 997, bulk_modulus = 2.09e9, - viscosity = 0.0010016, name) + viscosity = 0.0010016, gas_density = 0.0073955, gas_pressure = -1000, n = 1, name) pars = @parameters begin ρ = density β = bulk_modulus μ = viscosity + n = n + ρ_gas = gas_density + p_gas = gas_pressure end vars = @variables begin dm(t), [connect = Flow] end @@ -102,7 +111,11 @@ function transition(x1, x2, y1, y2, x) end end -density(port) = port.ρ +density_ref(port) = port.ρ +gas_density_ref(port) = port.ρ_gas +gas_pressure_ref(port) = port.p_gas bulk_modulus(port) = port.β viscosity(port) = port.μ -density(port, p) = density(port) * (1 + p / bulk_modulus(port)) +liquid_density(port) = density_ref(port) * (1 + port.p / bulk_modulus(port)) +gas_density(port) = density_ref(port) - port.p*( density_ref(port) - gas_density_ref(port) ) / gas_pressure_ref(port) +full_density(port) = liquid_density(port) #ifelse( port.p > 0, liquid_density(port), gas_density(port) ) diff --git a/src/Mechanical/Translational/sources.jl b/src/Mechanical/Translational/sources.jl index 73225e7c6..b3643f327 100644 --- a/src/Mechanical/Translational/sources.jl +++ b/src/Mechanical/Translational/sources.jl @@ -2,12 +2,31 @@ using ModelingToolkitStandardLibrary.Blocks @component function Force(; name) @named flange = MechanicalPort() - @named f = RealInput() + @named input = RealInput() vars = pars = [] eqs = [ - flange.f ~ -f.u, + flange.f ~ -input.u, ] compose(ODESystem(eqs, t, vars, pars; name = name, defaults = [flange.v => 0]), - [flange, f]) + [flange, input]) end + +@component function Position(;x_int=0, name) + @named flange = MechanicalPort() + @named input = RealInput() + pars = @parameters begin + x_int = x_int + end + vars = @variables begin + x(t) = x_int + dx(t) = 0 + end + eqs = [ + x ~ input.u + D(x) ~ dx + flange.v ~ dx + ] + compose(ODESystem(eqs, t, vars, pars; name = name, defaults = [flange.v => 0]), + [flange, input]) +end \ No newline at end of file diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index f657b48f2..2d95e2f06 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -7,7 +7,7 @@ import ModelingToolkitStandardLibrary.Mechanical.Translational as T D = Differential(t) # ------------------------------------------------------------ -# Test fluid domain and Pipe +# Test fluid domain and Tube function System(N; bulk_modulus, name) pars = @parameters begin bulk_modulus = bulk_modulus end @@ -20,9 +20,9 @@ function System(N; bulk_modulus, name) end if N == 1 - @named res = IC.PipeBase(; p_int = 0, area = 0.01, length = 500.0) + @named res = IC.TubeBase(; p_int = 0, area = 0.01, length = 500.0) else - @named res = IC.Pipe(N; p_int = 0, area = 0.01, length = 500.0) + @named res = IC.Tube(N; p_int = 0, area = 0.01, length = 500.0) end push!(systems, res) @@ -41,7 +41,7 @@ end NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax = 4 // 10) syss = structural_simplify.([sys1_2, sys1_1, sys5_1]) -probs = [ODEProblem(sys, [], (0, 0.2)) for sys in syss]; +probs = [ODEProblem(sys, ModelingToolkit.dummy_derivative_defaults(sys), (0, 0.2)) for sys in syss]; sols = [solve(prob, ImplicitEuler(nlsolve = NEWTON); initializealg = NoInit()) for prob in probs]; @@ -149,3 +149,71 @@ sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4, @test round(sol[s.vol2.x][1]; digits = 2) == 0.1 @test round(sol[s.vol1.x][end]; digits = 2) == 0.1 @test round(sol[s.vol2.x][end]; digits = 2) == 0.0 + + + +# Actuator System +function System(;name) + + @parameters t + + pars = @parameters begin + p_s = 285e5 + p_r = 5e5 + p_1 = 17e5 + p_2 = 42.9 + A_1 = 908e-4 + A_2 = 360e-4 + l_1 = 0.7 + l_2 = 2.7 + m_f = 276 + g = 0 + x_f_int = 0 + + d = 230e-3 + + + Cd= 70e5*2*876/(876*50/(25e-3*2π*230e-3))^2 #50_000 lpm + + m_piston = 880 + m_sled = 1500 + end + + vars = [] + + systems = @named begin + src = IC.Source(;p=p_s) + valve = IC.SpoolValve2Way(; p_s_int=p_s, p_a_int=p_1, p_b_int=p_2, p_r_int=p_r, g, m=m_f, x_int=x_f_int, d, Cd) + piston = IC.Actuator(;p_a_int=p_1, p_b_int=p_2, area_a=A_1, area_b=A_2, length_a_int=l_1, length_b_int=l_2, m=m_piston, g=0, x_int = 0, minimum_volume_a = 0, minimum_volume_b = 0) + sled = T.Mass(; m=m_sled) + pipe = IC.Tube(5; p_int=p_2, area = A_2, length = 2.0) + snk = IC.Source(;p=p_r) + input = B.Input(Float64) + pos = T.Position(; x_int = x_f_int) + + m1 = IC.FlowDivider(; p_int=p_2, n=3) + m2 = IC.FlowDivider(; p_int=p_2, n=3) + + fluid = IC.HydraulicFluid() + end + + eqs = [ + connect(input.output, pos.input) + connect(valve.flange, pos.flange) + connect(valve.port_a, piston.port_a) + connect(piston.flange, sled.flange) + connect(piston.port_b, m1.port_a) + connect(m1.port_b, pipe.port_b) + connect(pipe.port_a, m2.port_b) + connect(m2.port_a, valve.port_b) + connect(src.port, valve.port_s) + connect(snk.port, valve.port_r) + connect(fluid, src.port) + ] + + ODESystem(eqs, t, vars, pars; name, systems) +end + +@named system = System() + +sys = expand_connections(system) \ No newline at end of file From 9a6e92c6cd338894608ec4b1d9e03b8b16148c83 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 24 Apr 2023 19:52:19 -0400 Subject: [PATCH 02/16] added tests --- test/Hydraulic/isothermal_compressible.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index 2d95e2f06..c825aa3b6 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -216,4 +216,13 @@ end @named system = System() -sys = expand_connections(system) \ No newline at end of file +sys = expand_connections(system) +defs = ModelingToolkit.defaults(sys) +s = complete(system) + +@test Symbol(defs[s.src.port.ρ]) == Symbol(s.fluid.ρ) +@test Symbol(defs[s.valve.port_s.ρ]) == Symbol(s.fluid.ρ) +@test Symbol(defs[s.valve.port_a.ρ]) == Symbol(s.fluid.ρ) +@test Symbol(defs[s.valve.port_b.ρ]) == Symbol(s.fluid.ρ) +@test Symbol(defs[s.valve.port_r.ρ]) == Symbol(s.fluid.ρ) +@test Symbol(defs[s.snk.port.ρ]) == Symbol(s.fluid.ρ) From 59db8551949c44bf135a5da29285d9ad90281285 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Wed, 26 Apr 2023 19:14:37 -0400 Subject: [PATCH 03/16] fixed system test --- src/Blocks/sources.jl | 22 ++-- .../IsothermalCompressible/components.jl | 112 +++++++++--------- src/Hydraulic/IsothermalCompressible/utils.jl | 8 +- src/Mechanical/Translational/sources.jl | 16 +-- test/Hydraulic/isothermal_compressible.jl | 61 +++++----- 5 files changed, 104 insertions(+), 115 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index a2273c48f..ee35a4c40 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -412,7 +412,7 @@ end # - SawTooth Generate saw tooth signal # - Trapezoid Generate trapezoidal signal of type Real -function linear_interpolation(x1::T, x2::T, t1::T, t2::T, t::T) where {T<:Real} +function linear_interpolation(x1::T, x2::T, t1::T, t2::T, t::T) where {T <: Real} if t1 != t2 slope = (x2 - x1) / (t2 - t1) intercept = x1 - slope * t1 @@ -425,7 +425,7 @@ function linear_interpolation(x1::T, x2::T, t1::T, t2::T, t::T) where {T<:Real} end end -struct InputMemory{T<:Real} +struct InputMemory{T <: Real} data::Vector{T} sample_time::T n::Int @@ -433,7 +433,7 @@ end Base.copy(x::InputMemory{T}) where {T} = InputMemory{T}(copy(x.data), x.sample_time, x.n) -function InputMemory(data::Vector{T}, sample_time::T) where {T<:Real} +function InputMemory(data::Vector{T}, sample_time::T) where {T <: Real} InputMemory{T}(data, sample_time, length(data)) end @@ -473,10 +473,10 @@ function Symbolics.derivative(::typeof(input), args::NTuple{2, Any}, ::Val{1}) first_order_backwards_difference(args[1], args[2]) end - - Input(T::Type; name) = Input(T[], zero(T); name) -Input(data::Vector{T}, dt::T; name) where T<:Real = Input(; name, buffer = InputMemory(data, dt)) +function Input(data::Vector{T}, dt::T; name) where {T <: Real} + Input(; name, buffer = InputMemory(data, dt)) +end """ Input(; name, buffer) @@ -490,15 +490,11 @@ data input component. - `output` """ @component function Input(; name, buffer) - pars = @parameters begin - buffer = buffer - end + pars = @parameters begin buffer = buffer end vars = [] - systems = @named begin - output = RealOutput() - end + systems = @named begin output = RealOutput() end eqs = [ - output.u ~ input(t, buffer) + output.u ~ input(t, buffer), ] return ODESystem(eqs, t, vars, pars; name, systems) end diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index d5d2bf043..8977e4ae4 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -178,7 +178,8 @@ Tube modeled with `N` segements which models the fully developed flow friction a - `port_a`: hydraulic port - `port_b`: hydraulic port """ -@component function Tube(N; p_int, area, length, effective_length=length, perimeter = 2 * sqrt(area * pi), +@component function Tube(N; p_int, area, length, effective_length = length, + perimeter = 2 * sqrt(area * pi), shape_factor = 64, name) @assert(N>1, "the Tube component must be defined with more than 1 segment (i.e. N>1), found N=$N") @@ -188,7 +189,7 @@ Tube modeled with `N` segements which models the fully developed flow friction a p_int = p_int area = area length = length - effective_length=effective_length + effective_length = effective_length perimeter = perimeter Φ = shape_factor end @@ -241,7 +242,7 @@ Reduces the flow from `port_a` to `port_b` by `n`. Useful for modeling parallel - `port_a`: full flow hydraulic port - `port_b`: part flow hydraulic port """ -@component function FlowDivider(;p_int, n, name) +@component function FlowDivider(; p_int, n, name) #TODO: assert n >= 1 @@ -257,15 +258,12 @@ Reduces the flow from `port_a` to `port_b` by `n`. Useful for modeling parallel port_b = HydraulicPort(; p_int) end - eqs = [ - port_b.dm ~ port_a.dm/n - port_b.p ~ port_a.p - ] + eqs = [port_b.dm ~ port_a.dm / n + port_b.p ~ port_a.p] ODESystem(eqs, t, vars, pars; name, systems) end - """ DynamicVolume(; p_int, x_int=0, area, dead_volume=0, direction=+1, minimum_volume=0, name) @@ -327,16 +325,17 @@ dm ────► dead volume │ │ area # let dm = port.dm - u = dm/(rho*area) + u = dm / (rho * area) - eqs = [0 ~ IfElse.ifelse(vol >= minimum_volume, ( p ) - ( port.p - (1/2)*rho*u^2 ), dm) + eqs = [0 ~ IfElse.ifelse(vol >= minimum_volume, (p) - (port.p - (1 / 2) * rho * u^2), + dm) vol ~ dead_volume + area * x D(x) ~ dx D(rho) ~ drho dx ~ flange.v * direction rho ~ full_density(port) dm ~ drho * vol + rho * area * dx - flange.f ~ -p * area * direction] + flange.f ~ -p * area * direction] ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0]) end @@ -357,7 +356,7 @@ Valve with area input and discharge coefficient `Cd` defined by https://en.wikip - `port_b`: hydraulic port - `input`: real input setting the valve `area`. Note: absolute value taken """ -@component function Valve(reversible=false; p_a_int, p_b_int, area_int, Cd, name) +@component function Valve(reversible = false; p_a_int, p_b_int, area_int, Cd, name) pars = @parameters begin p_a_int = p_a_int p_b_int = p_b_int @@ -380,7 +379,7 @@ Valve with area input and discharge coefficient `Cd` defined by https://en.wikip area = if reversible input.u else - input.u*(input.u > 0) + input.u * (input.u > 0) end eqs = [sign(Δp) * dm ~ sqrt(2 * abs(Δp) * ρ / Cd) * area @@ -389,9 +388,7 @@ Valve with area input and discharge coefficient `Cd` defined by https://en.wikip ODESystem(eqs, t, vars, pars; name, systems, defaults = [input.u => area_int]) end - -@component function SpoolValve(reversible=false;p_a_int, p_b_int, x_int, Cd, d, name) - +@component function SpoolValve(reversible = false; p_a_int, p_b_int, x_int, Cd, d, name) pars = @parameters begin p_a_int = p_a_int p_b_int = p_b_int @@ -399,77 +396,74 @@ end x_int = x_int Cd = Cd end - + systems = @named begin port_a = HydraulicPort(; p_int = p_a_int) port_b = HydraulicPort(; p_int = p_b_int) flange = MechanicalPort() - valve = Valve(reversible; p_a_int, p_b_int, area_int=ParentScope(x_int)*2π*ParentScope(d), Cd) - end - + valve = Valve(reversible; p_a_int, p_b_int, + area_int = ParentScope(x_int) * 2π * ParentScope(d), Cd) + end + vars = @variables begin x(t) = x_int dx(t) = 0 end - - eqs = [ - D(x) ~ dx - flange.v ~ dx - flange.f ~ 0 #TODO: model flow force - connect(valve.port_a, port_a) - connect(valve.port_b, port_b) - valve.input.u ~ x*2π*d - ] + + eqs = [D(x) ~ dx + flange.v ~ dx + flange.f ~ 0 #TODO: model flow force + connect(valve.port_a, port_a) + connect(valve.port_b, port_b) + valve.input.u ~ x * 2π * d] ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0]) end -@component function SpoolValve2Way(reversible=false; p_s_int, p_a_int, p_b_int, p_r_int, m, g, x_int, Cd, d, name) - +@component function SpoolValve2Way(reversible = false; p_s_int, p_a_int, p_b_int, p_r_int, + m, g, x_int, Cd, d, name) pars = @parameters begin p_s_int = p_s_int p_a_int = p_a_int p_b_int = p_b_int p_r_int = p_r_int - m=m - g=g + m = m + g = g x_int = x_int - d=d + d = d end vars = [] systems = @named begin - vSA = SpoolValve(reversible; p_a_int=p_s_int, p_b_int=p_a_int, x_int, Cd, d) - vBR = SpoolValve(reversible; p_a_int=p_b_int, p_b_int=p_r_int, x_int, Cd, d) + vSA = SpoolValve(reversible; p_a_int = p_s_int, p_b_int = p_a_int, x_int, Cd, d) + vBR = SpoolValve(reversible; p_a_int = p_b_int, p_b_int = p_r_int, x_int, Cd, d) port_s = HydraulicPort(; p_int = p_s_int) port_a = HydraulicPort(; p_int = p_a_int) port_b = HydraulicPort(; p_int = p_b_int) port_r = HydraulicPort(; p_int = p_r_int) - mass = Mass(; m=m, g=g) + mass = Mass(; m = m, g = g) flange = MechanicalPort() end - eqs = [ - connect(vSA.port_a, port_s) - connect(vSA.port_b, port_a) - connect(vBR.port_a, port_b) - connect(vBR.port_b, port_r) - - connect(vSA.flange, vBR.flange, mass.flange, flange) - ] + eqs = [connect(vSA.port_a, port_s) + connect(vSA.port_b, port_a) + connect(vBR.port_a, port_b) + connect(vBR.port_b, port_r) + connect(vSA.flange, vBR.flange, mass.flange, flange)] ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0]) end -@component function Actuator(; p_a_int, p_b_int, area_a, area_b, length_a_int, length_b_int, m, g, x_int = 0, minimum_volume_a = 0, minimum_volume_b = 0, name) - +@component function Actuator(; p_a_int, p_b_int, area_a, area_b, length_a_int, length_b_int, + m, g, x_int = 0, minimum_volume_a = 0, minimum_volume_b = 0, + name) pars = @parameters begin p_a_int = p_a_int p_b_int = p_b_int @@ -478,7 +472,7 @@ end x_int = x_int length_a_int = length_a_int length_b_int = length_b_int - minimum_volume_a = minimum_volume_a + minimum_volume_a = minimum_volume_a minimum_volume_b = minimum_volume_b end @@ -488,21 +482,23 @@ end end systems = @named begin - vol_a = DynamicVolume(; p_int = p_a_int, x_int = +x_int, area=area_a, dead_volume = length_a_int*area_a, minimum_volume=minimum_volume_a, direction = +1) - vol_b = DynamicVolume(; p_int = p_b_int, x_int = -x_int, area=area_b, dead_volume = length_b_int*area_b, minimum_volume=minimum_volume_b, direction = -1) + vol_a = DynamicVolume(; p_int = p_a_int, x_int = +x_int, area = area_a, + dead_volume = length_a_int * area_a, + minimum_volume = minimum_volume_a, direction = +1) + vol_b = DynamicVolume(; p_int = p_b_int, x_int = -x_int, area = area_b, + dead_volume = length_b_int * area_b, + minimum_volume = minimum_volume_b, direction = -1) mass = Mass(; m, g) port_a = HydraulicPort(; p_int = p_a_int) port_b = HydraulicPort(; p_int = p_b_int) flange = MechanicalPort() end - eqs = [ - connect(vol_a.port, port_a) - connect(vol_b.port, port_b) - connect(vol_a.flange, vol_b.flange, mass.flange, flange) - x ~ vol_a.x - dx ~ vol_a.dx - ] + eqs = [connect(vol_a.port, port_a) + connect(vol_b.port, port_b) + connect(vol_a.flange, vol_b.flange, mass.flange, flange) + x ~ vol_a.x + dx ~ vol_a.dx] ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0]) -end \ No newline at end of file +end diff --git a/src/Hydraulic/IsothermalCompressible/utils.jl b/src/Hydraulic/IsothermalCompressible/utils.jl index 603d53431..e8816017c 100644 --- a/src/Hydraulic/IsothermalCompressible/utils.jl +++ b/src/Hydraulic/IsothermalCompressible/utils.jl @@ -44,7 +44,8 @@ Fluid parameter setter for isothermal compressible fluid domain. Defaults given - `p_gas`: [Pa] reference pressure (set by `gas_pressure` argument) """ @connector function HydraulicFluid(; density = 997, bulk_modulus = 2.09e9, - viscosity = 0.0010016, gas_density = 0.0073955, gas_pressure = -1000, n = 1, name) + viscosity = 0.0010016, gas_density = 0.0073955, + gas_pressure = -1000, n = 1, name) pars = @parameters begin ρ = density β = bulk_modulus @@ -117,5 +118,8 @@ gas_pressure_ref(port) = port.p_gas bulk_modulus(port) = port.β viscosity(port) = port.μ liquid_density(port) = density_ref(port) * (1 + port.p / bulk_modulus(port)) -gas_density(port) = density_ref(port) - port.p*( density_ref(port) - gas_density_ref(port) ) / gas_pressure_ref(port) +function gas_density(port) + density_ref(port) - + port.p * (density_ref(port) - gas_density_ref(port)) / gas_pressure_ref(port) +end full_density(port) = liquid_density(port) #ifelse( port.p > 0, liquid_density(port), gas_density(port) ) diff --git a/src/Mechanical/Translational/sources.jl b/src/Mechanical/Translational/sources.jl index b3643f327..2e60580c6 100644 --- a/src/Mechanical/Translational/sources.jl +++ b/src/Mechanical/Translational/sources.jl @@ -12,21 +12,17 @@ using ModelingToolkitStandardLibrary.Blocks [flange, input]) end -@component function Position(;x_int=0, name) +@component function Position(; x_int = 0, name) @named flange = MechanicalPort() @named input = RealInput() - pars = @parameters begin - x_int = x_int - end + pars = @parameters begin x_int = x_int end vars = @variables begin x(t) = x_int dx(t) = 0 end - eqs = [ - x ~ input.u - D(x) ~ dx - flange.v ~ dx - ] + eqs = [x ~ input.u + D(x) ~ dx + flange.v ~ dx] compose(ODESystem(eqs, t, vars, pars; name = name, defaults = [flange.v => 0]), [flange, input]) -end \ No newline at end of file +end diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index c825aa3b6..41f3888b3 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -41,7 +41,7 @@ end NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax = 4 // 10) syss = structural_simplify.([sys1_2, sys1_1, sys5_1]) -probs = [ODEProblem(sys, ModelingToolkit.dummy_derivative_defaults(sys), (0, 0.2)) for sys in syss]; +probs = [ODEProblem(sys, [], (0, 0.2)) for sys in syss]; #ModelingToolkit.missing_variable_defaults(sys) sols = [solve(prob, ImplicitEuler(nlsolve = NEWTON); initializealg = NoInit()) for prob in probs]; @@ -150,11 +150,8 @@ sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4, @test round(sol[s.vol1.x][end]; digits = 2) == 0.1 @test round(sol[s.vol2.x][end]; digits = 2) == 0.0 - - # Actuator System -function System(;name) - +function System(; name) @parameters t pars = @parameters begin @@ -169,54 +166,54 @@ function System(;name) m_f = 276 g = 0 x_f_int = 0 - + d = 230e-3 - - - Cd= 70e5*2*876/(876*50/(25e-3*2π*230e-3))^2 #50_000 lpm + + Cd = 70e5 * 2 * 876 / (876 * 50 / (25e-3 * 2π * 230e-3))^2 #50_000 lpm m_piston = 880 m_sled = 1500 end - + vars = [] systems = @named begin - src = IC.Source(;p=p_s) - valve = IC.SpoolValve2Way(; p_s_int=p_s, p_a_int=p_1, p_b_int=p_2, p_r_int=p_r, g, m=m_f, x_int=x_f_int, d, Cd) - piston = IC.Actuator(;p_a_int=p_1, p_b_int=p_2, area_a=A_1, area_b=A_2, length_a_int=l_1, length_b_int=l_2, m=m_piston, g=0, x_int = 0, minimum_volume_a = 0, minimum_volume_b = 0) - sled = T.Mass(; m=m_sled) - pipe = IC.Tube(5; p_int=p_2, area = A_2, length = 2.0) - snk = IC.Source(;p=p_r) + src = IC.Source(; p = p_s) + valve = IC.SpoolValve2Way(; p_s_int = p_s, p_a_int = p_1, p_b_int = p_2, + p_r_int = p_r, g, m = m_f, x_int = x_f_int, d, Cd) + piston = IC.Actuator(; p_a_int = p_1, p_b_int = p_2, area_a = A_1, area_b = A_2, + length_a_int = l_1, length_b_int = l_2, m = m_piston, g = 0, + x_int = 0, minimum_volume_a = 0, minimum_volume_b = 0) + sled = T.Mass(; m = m_sled) + pipe = IC.Tube(5; p_int = p_2, area = A_2, length = 2.0) + snk = IC.Source(; p = p_r) input = B.Input(Float64) pos = T.Position(; x_int = x_f_int) - m1 = IC.FlowDivider(; p_int=p_2, n=3) - m2 = IC.FlowDivider(; p_int=p_2, n=3) + m1 = IC.FlowDivider(; p_int = p_2, n = 3) + m2 = IC.FlowDivider(; p_int = p_2, n = 3) fluid = IC.HydraulicFluid() end - eqs = [ - connect(input.output, pos.input) - connect(valve.flange, pos.flange) - connect(valve.port_a, piston.port_a) - connect(piston.flange, sled.flange) - connect(piston.port_b, m1.port_a) - connect(m1.port_b, pipe.port_b) - connect(pipe.port_a, m2.port_b) - connect(m2.port_a, valve.port_b) - connect(src.port, valve.port_s) - connect(snk.port, valve.port_r) - connect(fluid, src.port) - ] + eqs = [connect(input.output, pos.input) + connect(valve.flange, pos.flange) + connect(valve.port_a, piston.port_a) + connect(piston.flange, sled.flange) + connect(piston.port_b, m1.port_a) + connect(m1.port_b, pipe.port_b) + connect(pipe.port_a, m2.port_b) + connect(m2.port_a, valve.port_b) + connect(src.port, valve.port_s) + connect(snk.port, valve.port_r) + connect(fluid, src.port, snk.port)] ODESystem(eqs, t, vars, pars; name, systems) end @named system = System() -sys = expand_connections(system) +sys = structural_simplify(system) defs = ModelingToolkit.defaults(sys) s = complete(system) From e0ca1c7d9048fe4a9e4010f01cbbba2e8715918a Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Fri, 28 Apr 2023 14:33:16 -0400 Subject: [PATCH 04/16] fixed FlowDivider --- src/Blocks/sources.jl | 64 ++++++++++++++----- .../IsothermalCompressible/components.jl | 39 ++++++++++- test/Hydraulic/isothermal_compressible.jl | 51 +++++++++++++-- 3 files changed, 131 insertions(+), 23 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index ee35a4c40..a81fb3a1e 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -425,30 +425,62 @@ function linear_interpolation(x1::T, x2::T, t1::T, t2::T, t::T) where {T <: Real end end -struct InputMemory{T <: Real} +struct Parameter{T <: Real} data::Vector{T} - sample_time::T + ref::T n::Int end -Base.copy(x::InputMemory{T}) where {T} = InputMemory{T}(copy(x.data), x.sample_time, x.n) +Base.:*(x::Number, y::Parameter) = x*y.ref +Base.:*(y::Parameter, x::Number) = Base.:*(x,y) +Base.:*(x::Parameter, y::Parameter) = x.ref*y.ref -function InputMemory(data::Vector{T}, sample_time::T) where {T <: Real} - InputMemory{T}(data, sample_time, length(data)) +Base.:/(x::Number, y::Parameter) = x/y.ref +Base.:/(y::Parameter, x::Number) = y.ref/x +Base.:/(x::Parameter, y::Parameter) = x.ref/y.ref + +Base.:+(x::Number, y::Parameter) = x + y.ref +Base.:+(y::Parameter, x::Number) = Base.:+(x,y) +Base.:+(x::Parameter, y::Parameter) = x.ref + y.ref + +Base.:-(x::Number, y::Parameter) = x - y.ref +Base.:-(y::Parameter, x::Number) = y.ref - x +Base.:-(x::Parameter, y::Parameter) = x.ref - y.ref + +Base.:^(x::Number, y::Parameter) = Base.power_by_squaring(x, y.ref) +Base.:^(y::Parameter, x::Number) = Base.power_by_squaring(y.ref, x) +Base.:^(x::Parameter, y::Parameter) = Base.power_by_squaring(x.ref, y.ref) + +Base.isless(x::Parameter, y::Number) = Base.isless(x.ref, y) +Base.isless(y::Number, x::Parameter) = Base.isless(y, x.ref) + + +Base.copy(x::Parameter{T}) where {T} = Parameter{T}(copy(x.data), x.ref, x.n) + +function Base.show(io::IO, m::MIME"text/plain", p::Parameter) + if !isempty(p.data) + print(io, p.data) + else + print(io, p.ref) + end end -function input(t, memory::InputMemory) +Parameter(x::Parameter) = x +Parameter(x::T) where T <: Real = Parameter(T[], x, 0) +Parameter(x::Vector{T}, dt::T) where T <: Real = Parameter(x, dt, length(x)) + +function input(t, memory::Parameter) if t < 0 t = zero(t) end - i1 = floor(Int, t / memory.sample_time) + 1 #expensive + i1 = floor(Int, t / memory.ref) + 1 #expensive i2 = i1 + 1 - t1 = i1 * memory.sample_time - t2 = i2 * memory.sample_time + t1 = i1 * memory.ref + t2 = i2 * memory.ref - #println("input: t=$t, i1=$i1, i2=$i2, t1=$t1, t2=$t2, Δt = $(memory.sample_time)") + #println("input: t=$t, i1=$i1, i2=$i2, t1=$t1, t2=$t2, Δt = $(memory.ref)") x1 = memory.data[i1 > memory.n ? memory.n : i1] x2 = memory.data[i2 > memory.n ? memory.n : i2] @@ -456,7 +488,7 @@ function input(t, memory::InputMemory) return linear_interpolation(x1, x2, t1, t2, t) end -get_sample_time(memory::InputMemory) = memory.sample_time +get_sample_time(memory::Parameter) = memory.ref Symbolics.@register_symbolic get_sample_time(memory) Symbolics.@register_symbolic input(t, memory) @@ -475,7 +507,7 @@ end Input(T::Type; name) = Input(T[], zero(T); name) function Input(data::Vector{T}, dt::T; name) where {T <: Real} - Input(; name, buffer = InputMemory(data, dt)) + Input(; name, buffer = Parameter(data, dt)) end """ @@ -484,17 +516,17 @@ end data input component. # Parameters: - - `buffer`: an `InputMemory` parameter which holds the data and sample time + - `buffer`: a `Parameter` type which holds the data and sample time # Connectors: - `output` """ @component function Input(; name, buffer) pars = @parameters begin buffer = buffer end - vars = [] + vars = [] systems = @named begin output = RealOutput() end eqs = [ - output.u ~ input(t, buffer), + output.u ~ input(t, buffer) ] - return ODESystem(eqs, t, vars, pars; name, systems) + return ODESystem(eqs, t, vars, pars; name, systems, defaults=[output.u => 0.0]) #TODO: get initial value from buffer end diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index 8977e4ae4..9aceebd26 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -75,6 +75,25 @@ Caps a hydrualic port to prevent mass flow in or out. ODESystem(eqs, t, vars, pars; name, systems) end + +@component function Open(; p_int, name) + pars = @parameters p_int = p_int + + vars = @variables begin + p(t) = p_int + dm(t) = 0 + end + + + systems = @named begin port = HydraulicPort(; p_int = p_int) end + + eqs = [port.p ~ p + port.dm ~ dm] + + ODESystem(eqs, t, vars, pars; name, systems) +end + + """ FixedVolume(; vol, p_int, name) @@ -251,15 +270,25 @@ Reduces the flow from `port_a` to `port_b` by `n`. Useful for modeling parallel p_int = p_int end - vars = [] + vars = @variables begin + dm_a(t) = 0 + dm_b(t) = 0 + end systems = @named begin port_a = HydraulicPort(; p_int) port_b = HydraulicPort(; p_int) + open = Open(; p_int) + end - eqs = [port_b.dm ~ port_a.dm / n - port_b.p ~ port_a.p] + eqs = [ + connect(port_a, port_b, open.port) + dm_a ~ port_a.dm + dm_b ~ dm_a/n + open.dm ~ dm_a - dm_b # extra flow dumps into an open port + # port_b.dm ~ dm_b # divided flow goes to port_b + ] ODESystem(eqs, t, vars, pars; name, systems) end @@ -434,6 +463,8 @@ end x_int = x_int d = d + + Cd=Cd end vars = [] @@ -474,6 +505,8 @@ end length_b_int = length_b_int minimum_volume_a = minimum_volume_a minimum_volume_b = minimum_volume_b + m=m + g=g end vars = @variables begin diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index 41f3888b3..71c7fdf5b 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -157,10 +157,13 @@ function System(; name) pars = @parameters begin p_s = 285e5 p_r = 5e5 - p_1 = 17e5 - p_2 = 42.9 + A_1 = 908e-4 A_2 = 360e-4 + + p_1 = 17e5 + p_2 = p_1*A_1/A_2 + l_1 = 0.7 l_2 = 2.7 m_f = 276 @@ -175,7 +178,9 @@ function System(; name) m_sled = 1500 end - vars = [] + vars = @variables begin + ddx(t) = 0 + end systems = @named begin src = IC.Source(; p = p_s) @@ -206,7 +211,8 @@ function System(; name) connect(m2.port_a, valve.port_b) connect(src.port, valve.port_s) connect(snk.port, valve.port_r) - connect(fluid, src.port, snk.port)] + connect(fluid, src.port, snk.port) + D(sled.v) ~ ddx] ODESystem(eqs, t, vars, pars; name, systems) end @@ -223,3 +229,40 @@ s = complete(system) @test Symbol(defs[s.valve.port_b.ρ]) == Symbol(s.fluid.ρ) @test Symbol(defs[s.valve.port_r.ρ]) == Symbol(s.fluid.ρ) @test Symbol(defs[s.snk.port.ρ]) == Symbol(s.fluid.ρ) + +prob = ODEProblem(sys, [], (0, 0.1); tofloat=false) + +# Generate Optimization Data + +dt = 1e-4 +defs[s.Cd] = 0.005 + + +# Data Set #1 +time = 0:dt:0.055 +x = @. 0.9*(time>0.015)*(time - 0.015)^2 - 25*(time>0.02)*(time - 0.02)^3 + +defs[s.input.buffer] = Parameter(x, dt) +defs[s.m_sled] = 1500 + +p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false)) +prob = remake(prob; p, tspan=(0, time[end])) + +@time sol1 = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); + + + +# Data Set #2 +time = 0:dt:0.1 +x = @. 5*((time>0.015)*(time - 0.015))^(3) - 650.0 * ((time>0.02)*(time - 0.02))^5 + +defs[s.input.buffer] = Parameter(x, dt) +defs[s.m_sled] = 500 + +p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false)) +prob = remake(prob; p, tspan=(0, time[end])) + +@time sol2 = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); + + + From 53eab0b903ece9e86d1891e5fd6bbfc64852cc3e Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Fri, 28 Apr 2023 16:54:45 -0400 Subject: [PATCH 05/16] update --- src/Blocks/sources.jl | 9 +++++++-- test/Hydraulic/isothermal_compressible.jl | 17 ++++++++++++++++- 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index a81fb3a1e..7265e4645 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -477,13 +477,18 @@ function input(t, memory::Parameter) i1 = floor(Int, t / memory.ref) + 1 #expensive i2 = i1 + 1 + if i2 > memory.n + i2 = memory.n + i1 = i2-1 + end + t1 = i1 * memory.ref t2 = i2 * memory.ref #println("input: t=$t, i1=$i1, i2=$i2, t1=$t1, t2=$t2, Δt = $(memory.ref)") - x1 = memory.data[i1 > memory.n ? memory.n : i1] - x2 = memory.data[i2 > memory.n ? memory.n : i2] + x1 = memory.data[i1] + x2 = memory.data[i2] return linear_interpolation(x1, x2, t1, t2, t) end diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index 71c7fdf5b..22cdb857b 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -3,6 +3,8 @@ import ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible as IC import ModelingToolkitStandardLibrary.Blocks as B import ModelingToolkitStandardLibrary.Mechanical.Translational as T +using ModelingToolkitStandardLibrary.Blocks: Parameter + @parameters t D = Differential(t) @@ -248,7 +250,7 @@ defs[s.m_sled] = 1500 p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false)) prob = remake(prob; p, tspan=(0, time[end])) -@time sol1 = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); +@time sol1 = solve(prob, ImplicitEuler(nlsolve = NEWTON) ; adaptive = false, dt, initializealg = NoInit()); @@ -266,3 +268,16 @@ prob = remake(prob; p, tspan=(0, time[end])) +# #= +# julia> @time sol2 = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); +# 0.644283 seconds (8.60 M allocations: 535.142 MiB, 14.96% gc time) +# =# + + +# # ModelOptimizer +# # In my real life optimization, I only know the sled acceleration s.ddx +# trial1 = Experiment(DataFrame(s.ddx => sol1[s.ddx]), sys, tspan = (0,0.055)) # How does this experiment know that s.m_sled = 1500, and s.input.buffer = input curve #1?? +# trial2 = Experiment(DataFrame(s.ddx => sol2[s.ddx]), sys, tspan = (0,0.1)) # How does this experiment know that s.m_sled = 500, and s.input.buffer = input curve #2?? + +# # How can I change the guess value of Cd to 0.008?? +# invprob = InverseProblem([trial1, trial2], sys, [s.Cd => (0.001, 0.1)]) # this should find that Cd = 0.005 \ No newline at end of file From 01b2e134a051d65e6a8400122993a8341a443bcb Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Sat, 29 Apr 2023 07:52:30 -0400 Subject: [PATCH 06/16] added time comparison using TimeVaryingFunction --- test/Hydraulic/isothermal_compressible.jl | 26 +++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index 22cdb857b..efbfba533 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -153,7 +153,7 @@ sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4, @test round(sol[s.vol2.x][end]; digits = 2) == 0.0 # Actuator System -function System(; name) +function System(use_input, f; name) @parameters t pars = @parameters begin @@ -194,7 +194,6 @@ function System(; name) sled = T.Mass(; m = m_sled) pipe = IC.Tube(5; p_int = p_2, area = A_2, length = 2.0) snk = IC.Source(; p = p_r) - input = B.Input(Float64) pos = T.Position(; x_int = x_f_int) m1 = IC.FlowDivider(; p_int = p_2, n = 3) @@ -203,6 +202,14 @@ function System(; name) fluid = IC.HydraulicFluid() end + if use_input + @named input = B.Input(Float64) + else + @named input = B.TimeVaryingFunction(f;t) + end + + push!(systems, input) + eqs = [connect(input.output, pos.input) connect(valve.flange, pos.flange) connect(valve.port_a, piston.port_a) @@ -219,7 +226,7 @@ function System(; name) ODESystem(eqs, t, vars, pars; name, systems) end -@named system = System() +@named system = System(true, nothing) sys = structural_simplify(system) defs = ModelingToolkit.defaults(sys) @@ -242,7 +249,8 @@ defs[s.Cd] = 0.005 # Data Set #1 time = 0:dt:0.055 -x = @. 0.9*(time>0.015)*(time - 0.015)^2 - 25*(time>0.02)*(time - 0.02)^3 +xf(time) = 0.9*(time>0.015)*(time - 0.015)^2 - 25*(time>0.02)*(time - 0.02)^3 +x = xf.(time) defs[s.input.buffer] = Parameter(x, dt) defs[s.m_sled] = 1500 @@ -253,6 +261,16 @@ prob = remake(prob; p, tspan=(0, time[end])) @time sol1 = solve(prob, ImplicitEuler(nlsolve = NEWTON) ; adaptive = false, dt, initializealg = NoInit()); +@register_symbolic xf(time) +Symbolics.derivative(::typeof(xf), args::NTuple{1, Any}, ::Val{1}) = 0 +@named system_ = System(false, xf) +s_ = complete(system_) +sys_ = structural_simplify(system_) + +prob_ = ODEProblem(sys_, [], (0, 0.055)) +@time sol_ = solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); + + # Data Set #2 time = 0:dt:0.1 From daf8d06aa8ab3cb904e56573fe370ab09661ed8e Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Sat, 29 Apr 2023 10:10:06 -0400 Subject: [PATCH 07/16] added Fast Solve --- test/Hydraulic/isothermal_compressible.jl | 35 +++++++---------------- 1 file changed, 11 insertions(+), 24 deletions(-) diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index efbfba533..af6ff3c7e 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -257,21 +257,9 @@ defs[s.m_sled] = 1500 p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false)) prob = remake(prob; p, tspan=(0, time[end])) - @time sol1 = solve(prob, ImplicitEuler(nlsolve = NEWTON) ; adaptive = false, dt, initializealg = NoInit()); -@register_symbolic xf(time) -Symbolics.derivative(::typeof(xf), args::NTuple{1, Any}, ::Val{1}) = 0 -@named system_ = System(false, xf) -s_ = complete(system_) -sys_ = structural_simplify(system_) - -prob_ = ODEProblem(sys_, [], (0, 0.055)) -@time sol_ = solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); - - - # Data Set #2 time = 0:dt:0.1 x = @. 5*((time>0.015)*(time - 0.015))^(3) - 650.0 * ((time>0.02)*(time - 0.02))^5 @@ -285,17 +273,16 @@ prob = remake(prob; p, tspan=(0, time[end])) @time sol2 = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); +# Fast Solve ------------ +@register_symbolic xf(time) +Symbolics.derivative(::typeof(xf), args::NTuple{1, Any}, ::Val{1}) = 0 +@named system_ = System(false, xf) +s_ = complete(system_) +sys_ = structural_simplify(system_) -# #= -# julia> @time sol2 = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); -# 0.644283 seconds (8.60 M allocations: 535.142 MiB, 14.96% gc time) -# =# - - -# # ModelOptimizer -# # In my real life optimization, I only know the sled acceleration s.ddx -# trial1 = Experiment(DataFrame(s.ddx => sol1[s.ddx]), sys, tspan = (0,0.055)) # How does this experiment know that s.m_sled = 1500, and s.input.buffer = input curve #1?? -# trial2 = Experiment(DataFrame(s.ddx => sol2[s.ddx]), sys, tspan = (0,0.1)) # How does this experiment know that s.m_sled = 500, and s.input.buffer = input curve #2?? +prob_ = ODEProblem(sys_, [], (0, 0.055)) +@time sol_ = solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); -# # How can I change the guess value of Cd to 0.008?? -# invprob = InverseProblem([trial1, trial2], sys, [s.Cd => (0.001, 0.1)]) # this should find that Cd = 0.005 \ No newline at end of file +p_ = Parameter.(prob_.p) +prob_ = remake(prob_; p=p_) +@time sol_ = solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); \ No newline at end of file From a24b11dfe13b36a25f9f9af53b8e68d45d57eaa9 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Sat, 29 Apr 2023 20:17:47 -0400 Subject: [PATCH 08/16] added isequal --- src/Blocks/sources.jl | 19 +++++++++++++++---- test/Hydraulic/isothermal_compressible.jl | 2 +- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 7265e4645..ccfeedbcf 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -431,6 +431,17 @@ struct Parameter{T <: Real} n::Int end +function Base.isequal(x::Parameter, y::Parameter) + b0 = x.n == y.n + if b0 + b1 = all(x.data .== y.data) + b2 = x.ref == y.ref + return b1 & b2 + else + return false + end +end + Base.:*(x::Number, y::Parameter) = x*y.ref Base.:*(y::Parameter, x::Number) = Base.:*(x,y) Base.:*(x::Parameter, y::Parameter) = x.ref*y.ref @@ -455,6 +466,8 @@ Base.isless(x::Parameter, y::Number) = Base.isless(x.ref, y) Base.isless(y::Number, x::Parameter) = Base.isless(y, x.ref) + + Base.copy(x::Parameter{T}) where {T} = Parameter{T}(copy(x.data), x.ref, x.n) function Base.show(io::IO, m::MIME"text/plain", p::Parameter) @@ -485,10 +498,8 @@ function input(t, memory::Parameter) t1 = i1 * memory.ref t2 = i2 * memory.ref - #println("input: t=$t, i1=$i1, i2=$i2, t1=$t1, t2=$t2, Δt = $(memory.ref)") - - x1 = memory.data[i1] - x2 = memory.data[i2] + x1 = @inbounds getindex(memory.data, i1) + x2 = @inbounds getindex(memory.data, i2) return linear_interpolation(x1, x2, t1, t2, t) end diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index af6ff3c7e..42603c0b7 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -280,7 +280,7 @@ Symbolics.derivative(::typeof(xf), args::NTuple{1, Any}, ::Val{1}) = 0 s_ = complete(system_) sys_ = structural_simplify(system_) -prob_ = ODEProblem(sys_, [], (0, 0.055)) +prob_ = ODEProblem(sys_, [], (0, 0.055), [s.Cd => 0.005]) @time sol_ = solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); p_ = Parameter.(prob_.p) From 705b440d49fb0cbf5f45ec3cfe26ed341d202884 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Sun, 30 Apr 2023 10:12:56 -0400 Subject: [PATCH 09/16] fixed speed with Parameter type stability --- src/Blocks/sources.jl | 25 +++++++++++++++---- .../IsothermalCompressible/components.jl | 2 +- src/Hydraulic/IsothermalCompressible/utils.jl | 2 ++ test/Hydraulic/isothermal_compressible.jl | 25 +++++++++++-------- 4 files changed, 37 insertions(+), 17 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index ccfeedbcf..f77c837ba 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -479,9 +479,19 @@ function Base.show(io::IO, m::MIME"text/plain", p::Parameter) end Parameter(x::Parameter) = x -Parameter(x::T) where T <: Real = Parameter(T[], x, 0) +function Parameter(x::T; tofloat=true) where T <: Real + if tofloat + x = float(x) + P = typeof(x) + else + P = T + end + + return Parameter(P[], x, 0) +end Parameter(x::Vector{T}, dt::T) where T <: Real = Parameter(x, dt, length(x)) + function input(t, memory::Parameter) if t < 0 t = zero(t) @@ -496,12 +506,17 @@ function input(t, memory::Parameter) end t1 = i1 * memory.ref - t2 = i2 * memory.ref - x1 = @inbounds getindex(memory.data, i1) - x2 = @inbounds getindex(memory.data, i2) - return linear_interpolation(x1, x2, t1, t2, t) + if t == t1 + return x1 + else + + t2 = i2 * memory.ref + x2 = @inbounds getindex(memory.data, i2) + return linear_interpolation(x1, x2, t1, t2, t) + end + end get_sample_time(memory::Parameter) = memory.ref diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index 9aceebd26..78055a6f1 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -408,7 +408,7 @@ Valve with area input and discharge coefficient `Cd` defined by https://en.wikip area = if reversible input.u else - input.u * (input.u > 0) + ifelse(input.u > 0, input.u , 0) end eqs = [sign(Δp) * dm ~ sqrt(2 * abs(Δp) * ρ / Cd) * area diff --git a/src/Hydraulic/IsothermalCompressible/utils.jl b/src/Hydraulic/IsothermalCompressible/utils.jl index e8816017c..06b2b1f87 100644 --- a/src/Hydraulic/IsothermalCompressible/utils.jl +++ b/src/Hydraulic/IsothermalCompressible/utils.jl @@ -99,6 +99,8 @@ function friction_factor(dm, area, d_h, density, viscosity, shape_factor) return f end @register_symbolic friction_factor(dm, area, d_h, density, viscosity, shape_factor) +Symbolics.derivative(::typeof(friction_factor), args, ::Val{1}) = 0 +Symbolics.derivative(::typeof(friction_factor), args, ::Val{4}) = 0 function transition(x1, x2, y1, y2, x) if x < x1 diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index 42603c0b7..244f9ca2a 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -231,6 +231,7 @@ end sys = structural_simplify(system) defs = ModelingToolkit.defaults(sys) s = complete(system) +prob = ODEProblem(sys, [], (0, 1.0); tofloat=false, jac=true) @test Symbol(defs[s.src.port.ρ]) == Symbol(s.fluid.ρ) @test Symbol(defs[s.valve.port_s.ρ]) == Symbol(s.fluid.ρ) @@ -239,7 +240,7 @@ s = complete(system) @test Symbol(defs[s.valve.port_r.ρ]) == Symbol(s.fluid.ρ) @test Symbol(defs[s.snk.port.ρ]) == Symbol(s.fluid.ρ) -prob = ODEProblem(sys, [], (0, 0.1); tofloat=false) + # Generate Optimization Data @@ -257,8 +258,11 @@ defs[s.m_sled] = 1500 p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false)) prob = remake(prob; p, tspan=(0, time[end])) -@time sol1 = solve(prob, ImplicitEuler(nlsolve = NEWTON) ; adaptive = false, dt, initializealg = NoInit()); - +solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); # precompile +@time "p is Vector{Parameter{Float64}}" sol1 = solve(prob, ImplicitEuler(nlsolve = NEWTON) ; adaptive = false, dt, initializealg = NoInit()); +u = copy(prob.u0) +prob.f.f(u,prob.u0, prob.p, 0.0) +@time "f()" prob.f.f(u,prob.u0, prob.p, 0.0) # Data Set #2 time = 0:dt:0.1 @@ -270,7 +274,7 @@ defs[s.m_sled] = 500 p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false)) prob = remake(prob; p, tspan=(0, time[end])) -@time sol2 = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); +sol2 = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); # Fast Solve ------------ @@ -279,10 +283,9 @@ Symbolics.derivative(::typeof(xf), args::NTuple{1, Any}, ::Val{1}) = 0 @named system_ = System(false, xf) s_ = complete(system_) sys_ = structural_simplify(system_) - -prob_ = ODEProblem(sys_, [], (0, 0.055), [s.Cd => 0.005]) -@time sol_ = solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); - -p_ = Parameter.(prob_.p) -prob_ = remake(prob_; p=p_) -@time sol_ = solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); \ No newline at end of file +prob_ = ODEProblem(sys_, [], (0, 0.055), [s_.Cd => 0.005]) +solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); #precompile +@time "p is Vector{Float64}" sol_ = solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); +u = copy(prob_.u0) +prob_.f.f(u,prob_.u0, prob_.p, 0.0) +@time "f()" prob_.f.f(u,prob_.u0, prob_.p, 0.0) From 07ee92176d9487d6321c7f4e18444092b9c10419 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 1 May 2023 10:41:32 -0400 Subject: [PATCH 10/16] publish update --- Project.toml | 3 +- docs/src/index.md | 1 + docs/src/tutorials/input_component.md | 147 ++++++ src/Blocks/Blocks.jl | 2 +- src/Blocks/sources.jl | 63 +-- .../IsothermalCompressible/components.jl | 25 +- test/Blocks/sources.jl | 31 +- test/Hydraulic/isothermal_compressible.jl | 467 ++++++++---------- 8 files changed, 442 insertions(+), 297 deletions(-) create mode 100644 docs/src/tutorials/input_component.md diff --git a/Project.toml b/Project.toml index 8896d1590..d9718057e 100644 --- a/Project.toml +++ b/Project.toml @@ -17,9 +17,10 @@ julia = "1.6" [extras] ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["OrdinaryDiffEq", "SafeTestsets", "Test", "ControlSystemsBase"] +test = ["OrdinaryDiffEq", "SafeTestsets", "Test", "ControlSystemsBase", "DataInterpolations"] diff --git a/docs/src/index.md b/docs/src/index.md index 296c41a1b..4d7f19e30 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -18,6 +18,7 @@ Pkg.add("ModelingToolkitStandardLibrary") - [Custom Component](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/tutorials/custom_component/) - [Thermal Model](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/tutorials/thermal_model/) - [DC Motor with PI-controller](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/tutorials/dc_motor_pi/) + - [Input Component](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/tutorials/input_component/) ## Libraries diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md new file mode 100644 index 000000000..8f40a9575 --- /dev/null +++ b/docs/src/tutorials/input_component.md @@ -0,0 +1,147 @@ +# Running Models with Discrete Data + +There are 2 ways to include data as part of a model. + + 1. using `ModelingToolkitStandardLibrary.Blocks.TimeVaryingFunction` + 2. using a custom component with external data or sensor + 3. using `ModelingToolkitStandardLibrary.Blocks.Input` + +This tutorial demonstrate each case and explain the pros and cons of each. + +## TimeVaryingFunction Component + +This component is easy to use and is performative. However the data is locked to the `ODESystem` and can only be changed by building a new `ODESystem`. Therefore, running a batch of data would not be efficient. Below is an example of how to use the `TimeVaryingFunction` with `DataInterpolations` to build the function from discrete data. + +```julia +using ModelingToolkit +using ModelingToolkitStandardLibrary.Blocks +using DataInterpolations +using OrdinaryDiffEq + +@parameters t +D = Differential(t) + +function System(f; name) + src = TimeVaryingFunction(f) + + vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0 + pars = @parameters m=10 k=1000 d=1 + + eqs = [f ~ src.output.u + ddx * 10 ~ k * x + d * dx + f + D(x) ~ dx + D(dx) ~ ddx] + + ODESystem(eqs, t, vars, pars; systems = [src], name) +end + +dt = 4e-4 +time = 0:dt:0.1 +data = sin.(2 * pi * time * 100) # example data + +f = LinearInterpolation(data, time) + +@named system = System(f) +sys = structural_simplify(system) +prob = ODEProblem(sys, [], (0, time[end])) +sol = solve(prob, ImplicitEuler()) +``` + +If we want to run a new data set, this requires building a new `LinearInterpolation` and `ODESystem` followed by running `structural_simplify`, all of which takes time. Therefore, to run serveral pieces of data it's better to re-use an `ODESystem`. The next couple methods will demonstrate this. + +## Custom Component with External Data + +The below code shows how to include data using a `Ref` and registered `input` function. This example uses a very basic function which requires non-adaptive solving and sampled data. As can be seen, the data can easily be set and changed before solving. + +```julia +const rdata = Ref{Vector{Float64}}() + +# Data Sets +data1 = sin.(2 * pi * time * 100) +data2 = cos.(2 * pi * time * 50) + +function input(t) + i = floor(Int, t / dt) + 1 + x = rdata[][i] + + return x +end + +Symbolics.@register_symbolic input(t) + +function System(; name) + vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0 + pars = @parameters m=10 k=1000 d=1 + + eqs = [f ~ input(t) + ddx * 10 ~ k * x + d * dx + f + D(x) ~ dx + D(dx) ~ ddx] + + ODESystem(eqs, t, vars, pars; name) +end + +@named system = System() +sys = structural_simplify(system) +prob = ODEProblem(sys, [], (0, time[end])) + +rdata[] = data1 +sol1 = solve(prob, ImplicitEuler(); dt, adaptive = false); +ddx1 = sol1[sys.ddx] + +rdata[] = data2 +sol2 = solve(prob, ImplicitEuler(); dt, adaptive = false) +ddx2 = sol2[sys.ddx] +``` + +The drawback of this method is that the solution observables can be linked to the data `Ref`, which means that if the data changes then the observables are no longer valid. In this case `ddx` is an observable that is derived directly from the data. Therefore, `sol1[sys.ddx]` is no longer correct after the data is changed for `sol2`. Additional code could be added to resolve this issue, for example by using a `Ref{Dict}` that could link a parameter of the model to the data source. This would also be necessary for parallel processing. + +```julia +# the following test will fail +@test all(ddx1 .== sol1[sys.ddx]) #returns false +``` + +## Input Component + +To resolve the issues presented above, the `Input` component can be used which allows for a resusable `ODESystem` and self contained data which ensures a solution which remains valid for it's lifetime. Now it's possible to also parallize the solving. + +```julia +function System(; name) + src = Input(Float64) + + vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0 + pars = @parameters m=10 k=1000 d=1 + + eqs = [f ~ src.output.u + ddx * 10 ~ k * x + d * dx + f + D(x) ~ dx + D(dx) ~ ddx] + + ODESystem(eqs, t, vars, pars; systems = [src], name) +end + +@named system = System() +sys = structural_simplify(system) +s = complete(system) +prob = ODEProblem(sys, [], (0, time[end]); tofloat = false) +defs = ModelingToolkit.defaults(sys) + +function get_prob(data) + defs[s.src.buffer] = Parameter(data, dt) + # ensure p is a uniform type of Vector{Parameter{Float64}} (converting from Vector{Any}) + p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false)) + remake(prob; p) +end + +prob1 = get_prob(data1) +prob2 = get_prob(data2) + +sol1 = Ref{ODESolution}() +sol2 = Ref{ODESolution}() +@sync begin + @async sol1[] = solve(prob1, ImplicitEuler()) + @async sol2[] = solve(prob2, ImplicitEuler()) +end +``` + +Note, in the above example, we can build the system with an empty `Input` component, only setting the expected data type: `@named src = Input(Float64)`. It's also possible to initialize the component with real data: `@named src = Input(data, dt)`. Additionally note that before running an `ODEProblem` using the `Input` component that the parameter vector should be a uniform type so that Julia is not slowed down by type instability. Because the `Input` component contains the `buffer` parameter of type `Parameter`, we must generate the problem using `tofloat=false`. This will initially give a parameter vector of type `Vector{Any}` with a mix of numbers and `Parameter` type. We can convert the vector to all `Parameter` type by running `p = Parameter.(p)`. This will wrap all the single values in a `Parameter` type which will be mathmatically equivalent to a `Number`. diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl index 36a7a008c..dab81e5f6 100644 --- a/src/Blocks/Blocks.jl +++ b/src/Blocks/Blocks.jl @@ -17,7 +17,7 @@ export Log, Log10 include("math.jl") export Constant, TimeVaryingFunction, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine, - Square, Triangular, InputMemory, Input + Square, Triangular, Parameter, Input include("sources.jl") export Limiter, DeadZone, SlewRateLimiter diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index f77c837ba..8ec900db1 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -412,7 +412,7 @@ end # - SawTooth Generate saw tooth signal # - Trapezoid Generate trapezoidal signal of type Real -function linear_interpolation(x1::T, x2::T, t1::T, t2::T, t::T) where {T <: Real} +function linear_interpolation(x1::T, x2::T, t1::T, t2::T, t) where {T <: Real} if t1 != t2 slope = (x2 - x1) / (t2 - t1) intercept = x1 - slope * t1 @@ -431,7 +431,7 @@ struct Parameter{T <: Real} n::Int end -function Base.isequal(x::Parameter, y::Parameter) +function Base.isequal(x::Parameter, y::Parameter) b0 = x.n == y.n if b0 b1 = all(x.data .== y.data) @@ -442,16 +442,16 @@ function Base.isequal(x::Parameter, y::Parameter) end end -Base.:*(x::Number, y::Parameter) = x*y.ref -Base.:*(y::Parameter, x::Number) = Base.:*(x,y) -Base.:*(x::Parameter, y::Parameter) = x.ref*y.ref +Base.:*(x::Number, y::Parameter) = x * y.ref +Base.:*(y::Parameter, x::Number) = Base.:*(x, y) +Base.:*(x::Parameter, y::Parameter) = x.ref * y.ref -Base.:/(x::Number, y::Parameter) = x/y.ref -Base.:/(y::Parameter, x::Number) = y.ref/x -Base.:/(x::Parameter, y::Parameter) = x.ref/y.ref +Base.:/(x::Number, y::Parameter) = x / y.ref +Base.:/(y::Parameter, x::Number) = y.ref / x +Base.:/(x::Parameter, y::Parameter) = x.ref / y.ref Base.:+(x::Number, y::Parameter) = x + y.ref -Base.:+(y::Parameter, x::Number) = Base.:+(x,y) +Base.:+(y::Parameter, x::Number) = Base.:+(x, y) Base.:+(x::Parameter, y::Parameter) = x.ref + y.ref Base.:-(x::Number, y::Parameter) = x - y.ref @@ -465,21 +465,18 @@ Base.:^(x::Parameter, y::Parameter) = Base.power_by_squaring(x.ref, y.ref) Base.isless(x::Parameter, y::Number) = Base.isless(x.ref, y) Base.isless(y::Number, x::Parameter) = Base.isless(y, x.ref) - - - Base.copy(x::Parameter{T}) where {T} = Parameter{T}(copy(x.data), x.ref, x.n) function Base.show(io::IO, m::MIME"text/plain", p::Parameter) if !isempty(p.data) - print(io, p.data) + print(io, p.data) else print(io, p.ref) end end Parameter(x::Parameter) = x -function Parameter(x::T; tofloat=true) where T <: Real +function Parameter(x::T; tofloat = true) where {T <: Real} if tofloat x = float(x) P = typeof(x) @@ -489,34 +486,43 @@ function Parameter(x::T; tofloat=true) where T <: Real return Parameter(P[], x, 0) end -Parameter(x::Vector{T}, dt::T) where T <: Real = Parameter(x, dt, length(x)) +Parameter(x::Vector{T}, dt::T) where {T <: Real} = Parameter(x, dt, length(x)) - -function input(t, memory::Parameter) +function input(t, memory::Parameter{T}) where {T} if t < 0 t = zero(t) end + if isempty(memory.data) + if T isa Float16 + return NaN16 + elseif T isa Float32 + return NaN32 + elseif T isa Float64 + return NaN64 + else + return zero(T) + end + end + i1 = floor(Int, t / memory.ref) + 1 #expensive i2 = i1 + 1 - if i2 > memory.n - i2 = memory.n - i1 = i2-1 - end - - t1 = i1 * memory.ref + t1 = (i1 - 1) * memory.ref x1 = @inbounds getindex(memory.data, i1) if t == t1 return x1 else + if i2 > memory.n + i2 = memory.n + i1 = i2 - 1 + end - t2 = i2 * memory.ref + t2 = (i2 - 1) * memory.ref x2 = @inbounds getindex(memory.data, i2) return linear_interpolation(x1, x2, t1, t2, t) end - end get_sample_time(memory::Parameter) = memory.ref @@ -554,10 +560,11 @@ data input component. """ @component function Input(; name, buffer) pars = @parameters begin buffer = buffer end - vars = [] + vars = [] systems = @named begin output = RealOutput() end eqs = [ - output.u ~ input(t, buffer) + output.u ~ input(t, buffer), ] - return ODESystem(eqs, t, vars, pars; name, systems, defaults=[output.u => 0.0]) #TODO: get initial value from buffer + return ODESystem(eqs, t, vars, pars; name, systems, + defaults = [output.u => input(0.0, buffer)]) #TODO: get initial value from buffer end diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index 78055a6f1..67a4831ff 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -75,7 +75,6 @@ Caps a hydrualic port to prevent mass flow in or out. ODESystem(eqs, t, vars, pars; name, systems) end - @component function Open(; p_int, name) pars = @parameters p_int = p_int @@ -83,7 +82,6 @@ end p(t) = p_int dm(t) = 0 end - systems = @named begin port = HydraulicPort(; p_int = p_int) end @@ -93,7 +91,6 @@ end ODESystem(eqs, t, vars, pars; name, systems) end - """ FixedVolume(; vol, p_int, name) @@ -279,16 +276,14 @@ Reduces the flow from `port_a` to `port_b` by `n`. Useful for modeling parallel port_a = HydraulicPort(; p_int) port_b = HydraulicPort(; p_int) open = Open(; p_int) - end - eqs = [ - connect(port_a, port_b, open.port) - dm_a ~ port_a.dm - dm_b ~ dm_a/n - open.dm ~ dm_a - dm_b # extra flow dumps into an open port - # port_b.dm ~ dm_b # divided flow goes to port_b - ] + eqs = [connect(port_a, port_b, open.port) + dm_a ~ port_a.dm + dm_b ~ dm_a / n + open.dm ~ dm_a - dm_b # extra flow dumps into an open port + # port_b.dm ~ dm_b # divided flow goes to port_b + ] ODESystem(eqs, t, vars, pars; name, systems) end @@ -408,7 +403,7 @@ Valve with area input and discharge coefficient `Cd` defined by https://en.wikip area = if reversible input.u else - ifelse(input.u > 0, input.u , 0) + ifelse(input.u > 0, input.u, 0) end eqs = [sign(Δp) * dm ~ sqrt(2 * abs(Δp) * ρ / Cd) * area @@ -464,7 +459,7 @@ end d = d - Cd=Cd + Cd = Cd end vars = [] @@ -505,8 +500,8 @@ end length_b_int = length_b_int minimum_volume_a = minimum_volume_a minimum_volume_b = minimum_volume_b - m=m - g=g + m = m + g = g end vars = @variables begin diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index 0e54a0f03..1c87f36d3 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -25,7 +25,7 @@ using OrdinaryDiffEq: ReturnCode.Success end @testset "TimeVaryingFunction" begin - f(t) = t^2 + f(t) = t^2 + 1 @named src = TimeVaryingFunction(f) @named int = Integrator() @@ -41,7 +41,7 @@ end sol = solve(prob, Rodas4()) @test sol.retcode == Success @test sol[src.output.u]≈f.(sol.t) atol=1e-3 - @test sol[int.output.u][end]≈1 / 3 * 10^3 atol=1e-3 # closed-form solution to integral + @test sol[int.output.u][end]≈1 / 3 * 10^3 + 10 atol=1e-3 # closed-form solution to integral end @testset "Sine" begin @@ -406,3 +406,30 @@ end amplitude, damping, phase, offset, start_time) atol=1e-3 end + +@testset "Input" begin + dt = 4e-4 + t_end = 10.0 + time = 0:dt:t_end + x = @. time^2 + 1.0 + + @named src = Input(Float64) + @named int = Integrator() + @named iosys = ODESystem([ + connect(src.output, int.input), + ], + t, + systems = [int, src]) + sys = structural_simplify(iosys) + s = complete(iosys) + prob = ODEProblem(sys, [], (0.0, t_end), [s.src.buffer => Parameter(x, dt)]; + tofloat = false) + prob = remake(prob; p = Parameter.(prob.p)) + + sol = solve(prob, Rodas4(); initializealg = NoInit()) + @test sol.retcode == Success + @test sol[src.output.u][1] == 1.0 #check correct initial condition + + @test sol(time)[src.output.u]≈x atol=1e-3 + @test sol[int.output.u][end]≈1 / 3 * 10^3 + 10.0 atol=1e-3 # closed-form solution to integral +end diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index 244f9ca2a..c8523b852 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -8,284 +8,251 @@ using ModelingToolkitStandardLibrary.Blocks: Parameter @parameters t D = Differential(t) -# ------------------------------------------------------------ -# Test fluid domain and Tube -function System(N; bulk_modulus, name) - pars = @parameters begin bulk_modulus = bulk_modulus end - - systems = @named begin - fluid = IC.HydraulicFluid(; bulk_modulus) - stp = B.Step(; height = 10e5, offset = 0, start_time = 0.005, duration = Inf, - smooth = true) - src = IC.InputSource(; p_int = 0) - vol = IC.FixedVolume(; p_int = 0, vol = 10.0) - end - - if N == 1 - @named res = IC.TubeBase(; p_int = 0, area = 0.01, length = 500.0) - else - @named res = IC.Tube(N; p_int = 0, area = 0.01, length = 500.0) - end - push!(systems, res) - - eqs = [connect(stp.output, src.input) - connect(fluid, src.port) - connect(src.port, res.port_a) - connect(res.port_b, vol.port)] - - ODESystem(eqs, t, [], pars; name, systems) -end - -@named sys1_2 = System(1; bulk_modulus = 2e9) -@named sys1_1 = System(1; bulk_modulus = 1e9) -@named sys5_1 = System(5; bulk_modulus = 1e9) - NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax = 4 // 10) -syss = structural_simplify.([sys1_2, sys1_1, sys5_1]) -probs = [ODEProblem(sys, [], (0, 0.2)) for sys in syss]; #ModelingToolkit.missing_variable_defaults(sys) -sols = [solve(prob, ImplicitEuler(nlsolve = NEWTON); initializealg = NoInit()) - for prob in probs]; - -s1_2 = complete(sys1_2) -s1_1 = complete(sys1_1) -s5_1 = complete(sys5_1) - -# higher stiffness should compress more quickly and give a higher pressure -@test sols[1][s1_2.vol.port.p][end] > sols[2][s1_1.vol.port.p][end] - -# N=5 pipe is compressible, will pressurize more slowly -@test sols[2][s1_1.vol.port.p][end] > sols[3][s5_1.vol.port.p][end] - -# ------------------------------------------------------------ -# Test the Valve -function System(; name) - pars = [] - - systems = @named begin - fluid = IC.HydraulicFluid() - sink = IC.Source(; p = 10e5) - vol = IC.FixedVolume(; vol = 0.1, p_int = 100e5) - valve = IC.Valve(; p_a_int = 10e5, p_b_int = 100e5, area_int = 0, Cd = 1e6) - ramp = B.Ramp(; height = 1, duration = 0.001, offset = 0, start_time = 0.001, - smooth = true) +@testset "Fluid Domain and Tube" begin + function System(N; bulk_modulus, name) + pars = @parameters begin bulk_modulus = bulk_modulus end + + systems = @named begin + fluid = IC.HydraulicFluid(; bulk_modulus) + stp = B.Step(; height = 10e5, offset = 0, start_time = 0.005, duration = Inf, + smooth = true) + src = IC.InputSource(; p_int = 0) + vol = IC.FixedVolume(; p_int = 0, vol = 10.0) + end + + if N == 1 + @named res = IC.TubeBase(; p_int = 0, area = 0.01, length = 500.0) + else + @named res = IC.Tube(N; p_int = 0, area = 0.01, length = 500.0) + end + push!(systems, res) + + eqs = [connect(stp.output, src.input) + connect(fluid, src.port) + connect(src.port, res.port_a) + connect(res.port_b, vol.port)] + + ODESystem(eqs, t, [], pars; name, systems) end - eqs = [connect(fluid, sink.port) - connect(sink.port, valve.port_a) - connect(valve.port_b, vol.port) - connect(valve.input, ramp.output)] - - ODESystem(eqs, t, [], pars; name, systems) -end - -@named valve_system = System() -s = complete(valve_system) -sys = structural_simplify(valve_system) -prob = ODEProblem(sys, [], (0, 0.01)) -sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4, - initializealg = NoInit()) - -# the volume should discharge to 10bar -@test sol[s.vol.port.p][end] ≈ 10e5 - -# ------------------------------------------------------------ -# Test DynmaicVolume and minimum_volume feature -function System(; name) - pars = [] + @named sys1_2 = System(1; bulk_modulus = 2e9) + @named sys1_1 = System(1; bulk_modulus = 1e9) + @named sys5_1 = System(5; bulk_modulus = 1e9) - # DynamicVolume values - area = 0.01 - length = 0.1 + syss = structural_simplify.([sys1_2, sys1_1, sys5_1]) + probs = [ODEProblem(sys, [], (0, 0.2)) for sys in syss] #ModelingToolkit.missing_variable_defaults(sys) + sols = [solve(prob, ImplicitEuler(nlsolve = NEWTON); initializealg = NoInit()) + for prob in probs] - systems = @named begin - fluid = IC.HydraulicFluid() - src1 = IC.InputSource(; p_int = 10e5) - src2 = IC.InputSource(; p_int = 10e5) + s1_2 = complete(sys1_2) + s1_1 = complete(sys1_1) + s5_1 = complete(sys5_1) - vol1 = IC.DynamicVolume(; p_int = 10e5, area, direction = +1, dead_volume = 2e-4, - minimum_volume = 2e-4) - vol2 = IC.DynamicVolume(; p_int = 10e5, area, direction = -1, dead_volume = 2e-4, - minimum_volume = 2e-4, x_int = length) + # higher stiffness should compress more quickly and give a higher pressure + @test sols[1][s1_2.vol.port.p][end] > sols[2][s1_1.vol.port.p][end] - mass = T.Mass(; m = 100) - - sin1 = B.Sine(; frequency = 0.5, amplitude = +1e5, offset = 10e5) - sin2 = B.Sine(; frequency = 0.5, amplitude = -1e5, offset = 10e5) - end - - eqs = [connect(fluid, src1.port, src2.port) - connect(src1.port, vol1.port) - connect(src2.port, vol2.port) - connect(vol1.flange, mass.flange, vol2.flange) - connect(src1.input, sin1.output) - connect(src2.input, sin2.output)] - - ODESystem(eqs, t, [], pars; name, systems) + # N=5 pipe is compressible, will pressurize more slowly + @test sols[2][s1_1.vol.port.p][end] > sols[3][s5_1.vol.port.p][end] end -@named min_vol_system = System() -s = complete(min_vol_system) -sys = structural_simplify(min_vol_system) -prob = ODEProblem(sys, [], (0, 5.0)) - -sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4, - initializealg = NoInit()) - -# begin -# fig = Figure() - -# ax = Axis(fig[1,1], ylabel="position [m]", xlabel="time [s]") -# lines!(ax, sol.t, sol[s.vol1.x]; label="vol1") -# lines!(ax, sol.t, sol[s.vol2.x]; label="vol2") - -# ax = Axis(fig[2,1], ylabel="pressure [bar]", xlabel="time [s]") -# lines!(ax, sol.t, sol[s.vol1.port.p]/1e5; label="vol1") -# lines!(ax, sol.t, sol[s.vol2.port.p]/1e5; label="vol2") -# Legend(fig[2,2], ax) -# fig -# end - -# volume/mass should stop moving at opposite ends -@test round(sol[s.vol1.x][1]; digits = 2) == 0.0 -@test round(sol[s.vol2.x][1]; digits = 2) == 0.1 -@test round(sol[s.vol1.x][end]; digits = 2) == 0.1 -@test round(sol[s.vol2.x][end]; digits = 2) == 0.0 - -# Actuator System -function System(use_input, f; name) - @parameters t - - pars = @parameters begin - p_s = 285e5 - p_r = 5e5 - - A_1 = 908e-4 - A_2 = 360e-4 - - p_1 = 17e5 - p_2 = p_1*A_1/A_2 - - l_1 = 0.7 - l_2 = 2.7 - m_f = 276 - g = 0 - x_f_int = 0 - - d = 230e-3 - - Cd = 70e5 * 2 * 876 / (876 * 50 / (25e-3 * 2π * 230e-3))^2 #50_000 lpm - - m_piston = 880 - m_sled = 1500 +@testset "Valve" begin + function System(; name) + pars = [] + + systems = @named begin + fluid = IC.HydraulicFluid() + sink = IC.Source(; p = 10e5) + vol = IC.FixedVolume(; vol = 0.1, p_int = 100e5) + valve = IC.Valve(; p_a_int = 10e5, p_b_int = 100e5, area_int = 0, Cd = 1e6) + ramp = B.Ramp(; height = 1, duration = 0.001, offset = 0, start_time = 0.001, + smooth = true) + end + + eqs = [connect(fluid, sink.port) + connect(sink.port, valve.port_a) + connect(valve.port_b, vol.port) + connect(valve.input, ramp.output)] + + ODESystem(eqs, t, [], pars; name, systems) end - vars = @variables begin - ddx(t) = 0 - end - - systems = @named begin - src = IC.Source(; p = p_s) - valve = IC.SpoolValve2Way(; p_s_int = p_s, p_a_int = p_1, p_b_int = p_2, - p_r_int = p_r, g, m = m_f, x_int = x_f_int, d, Cd) - piston = IC.Actuator(; p_a_int = p_1, p_b_int = p_2, area_a = A_1, area_b = A_2, - length_a_int = l_1, length_b_int = l_2, m = m_piston, g = 0, - x_int = 0, minimum_volume_a = 0, minimum_volume_b = 0) - sled = T.Mass(; m = m_sled) - pipe = IC.Tube(5; p_int = p_2, area = A_2, length = 2.0) - snk = IC.Source(; p = p_r) - pos = T.Position(; x_int = x_f_int) - - m1 = IC.FlowDivider(; p_int = p_2, n = 3) - m2 = IC.FlowDivider(; p_int = p_2, n = 3) - - fluid = IC.HydraulicFluid() - end + @named valve_system = System() + s = complete(valve_system) + sys = structural_simplify(valve_system) + prob = ODEProblem(sys, [], (0, 0.01)) + sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4, + initializealg = NoInit()) - if use_input - @named input = B.Input(Float64) - else - @named input = B.TimeVaryingFunction(f;t) - end - - push!(systems, input) - - eqs = [connect(input.output, pos.input) - connect(valve.flange, pos.flange) - connect(valve.port_a, piston.port_a) - connect(piston.flange, sled.flange) - connect(piston.port_b, m1.port_a) - connect(m1.port_b, pipe.port_b) - connect(pipe.port_a, m2.port_b) - connect(m2.port_a, valve.port_b) - connect(src.port, valve.port_s) - connect(snk.port, valve.port_r) - connect(fluid, src.port, snk.port) - D(sled.v) ~ ddx] - - ODESystem(eqs, t, vars, pars; name, systems) + # the volume should discharge to 10bar + @test sol[s.vol.port.p][end] ≈ 10e5 end -@named system = System(true, nothing) +@testset "DynamicVolume and minimum_volume feature" begin + function System(; name) + pars = [] -sys = structural_simplify(system) -defs = ModelingToolkit.defaults(sys) -s = complete(system) -prob = ODEProblem(sys, [], (0, 1.0); tofloat=false, jac=true) + # DynamicVolume values + area = 0.01 + length = 0.1 -@test Symbol(defs[s.src.port.ρ]) == Symbol(s.fluid.ρ) -@test Symbol(defs[s.valve.port_s.ρ]) == Symbol(s.fluid.ρ) -@test Symbol(defs[s.valve.port_a.ρ]) == Symbol(s.fluid.ρ) -@test Symbol(defs[s.valve.port_b.ρ]) == Symbol(s.fluid.ρ) -@test Symbol(defs[s.valve.port_r.ρ]) == Symbol(s.fluid.ρ) -@test Symbol(defs[s.snk.port.ρ]) == Symbol(s.fluid.ρ) + systems = @named begin + fluid = IC.HydraulicFluid() + src1 = IC.InputSource(; p_int = 10e5) + src2 = IC.InputSource(; p_int = 10e5) + vol1 = IC.DynamicVolume(; p_int = 10e5, area, direction = +1, + dead_volume = 2e-4, + minimum_volume = 2e-4) + vol2 = IC.DynamicVolume(; p_int = 10e5, area, direction = -1, + dead_volume = 2e-4, + minimum_volume = 2e-4, x_int = length) + mass = T.Mass(; m = 100) -# Generate Optimization Data + sin1 = B.Sine(; frequency = 0.5, amplitude = +1e5, offset = 10e5) + sin2 = B.Sine(; frequency = 0.5, amplitude = -1e5, offset = 10e5) + end -dt = 1e-4 -defs[s.Cd] = 0.005 + eqs = [connect(fluid, src1.port, src2.port) + connect(src1.port, vol1.port) + connect(src2.port, vol2.port) + connect(vol1.flange, mass.flange, vol2.flange) + connect(src1.input, sin1.output) + connect(src2.input, sin2.output)] + ODESystem(eqs, t, [], pars; name, systems) + end -# Data Set #1 -time = 0:dt:0.055 -xf(time) = 0.9*(time>0.015)*(time - 0.015)^2 - 25*(time>0.02)*(time - 0.02)^3 -x = xf.(time) + @named min_vol_system = System() + s = complete(min_vol_system) + sys = structural_simplify(min_vol_system) + prob = ODEProblem(sys, [], (0, 5.0)) + + sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4, + initializealg = NoInit()) + + # begin + # fig = Figure() + + # ax = Axis(fig[1,1], ylabel="position [m]", xlabel="time [s]") + # lines!(ax, sol.t, sol[s.vol1.x]; label="vol1") + # lines!(ax, sol.t, sol[s.vol2.x]; label="vol2") + + # ax = Axis(fig[2,1], ylabel="pressure [bar]", xlabel="time [s]") + # lines!(ax, sol.t, sol[s.vol1.port.p]/1e5; label="vol1") + # lines!(ax, sol.t, sol[s.vol2.port.p]/1e5; label="vol2") + # Legend(fig[2,2], ax) + # fig + # end + + # volume/mass should stop moving at opposite ends + @test round(sol[s.vol1.x][1]; digits = 2) == 0.0 + @test round(sol[s.vol2.x][1]; digits = 2) == 0.1 + @test round(sol[s.vol1.x][end]; digits = 2) == 0.1 + @test round(sol[s.vol2.x][end]; digits = 2) == 0.0 +end -defs[s.input.buffer] = Parameter(x, dt) -defs[s.m_sled] = 1500 +@testset "Actuator System" begin + function System(use_input, f; name) + @parameters t + + pars = @parameters begin + p_s = 285e5 + p_r = 5e5 + + A_1 = 908e-4 + A_2 = 360e-4 + + p_1 = 17e5 + p_2 = p_1 * A_1 / A_2 + + l_1 = 0.7 + l_2 = 2.7 + m_f = 276 + g = 0 + x_f_int = 0 + + d = 230e-3 + + Cd = 70e5 * 2 * 876 / (876 * 50 / (25e-3 * 2π * 230e-3))^2 #50_000 lpm + + m_piston = 880 + m_body = 1500 + end + + vars = @variables begin ddx(t) = 0 end + + systems = @named begin + src = IC.Source(; p = p_s) + valve = IC.SpoolValve2Way(; p_s_int = p_s, p_a_int = p_1, p_b_int = p_2, + p_r_int = p_r, g, m = m_f, x_int = x_f_int, d, Cd) + piston = IC.Actuator(; p_a_int = p_1, p_b_int = p_2, area_a = A_1, area_b = A_2, + length_a_int = l_1, length_b_int = l_2, m = m_piston, + g = 0, + x_int = 0, minimum_volume_a = 0, minimum_volume_b = 0) + body = T.Mass(; m = m_body) + pipe = IC.Tube(5; p_int = p_2, area = A_2, length = 2.0) + snk = IC.Source(; p = p_r) + pos = T.Position(; x_int = x_f_int) + + m1 = IC.FlowDivider(; p_int = p_2, n = 3) + m2 = IC.FlowDivider(; p_int = p_2, n = 3) + + fluid = IC.HydraulicFluid() + end + + if use_input + @named input = B.Input(Float64) + else + @named input = B.TimeVaryingFunction(f; t) + end + + push!(systems, input) + + eqs = [connect(input.output, pos.input) + connect(valve.flange, pos.flange) + connect(valve.port_a, piston.port_a) + connect(piston.flange, body.flange) + connect(piston.port_b, m1.port_a) + connect(m1.port_b, pipe.port_b) + connect(pipe.port_a, m2.port_b) + connect(m2.port_a, valve.port_b) + connect(src.port, valve.port_s) + connect(snk.port, valve.port_r) + connect(fluid, src.port, snk.port) + D(body.v) ~ ddx] + + ODESystem(eqs, t, vars, pars; name, systems) + end -p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false)) -prob = remake(prob; p, tspan=(0, time[end])) -solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); # precompile -@time "p is Vector{Parameter{Float64}}" sol1 = solve(prob, ImplicitEuler(nlsolve = NEWTON) ; adaptive = false, dt, initializealg = NoInit()); -u = copy(prob.u0) -prob.f.f(u,prob.u0, prob.p, 0.0) -@time "f()" prob.f.f(u,prob.u0, prob.p, 0.0) + @named system = System(true, nothing) -# Data Set #2 -time = 0:dt:0.1 -x = @. 5*((time>0.015)*(time - 0.015))^(3) - 650.0 * ((time>0.02)*(time - 0.02))^5 + sys = structural_simplify(system) + defs = ModelingToolkit.defaults(sys) + s = complete(system) + prob = ODEProblem(sys, [], (0, 1.0); tofloat = false, jac = true) -defs[s.input.buffer] = Parameter(x, dt) -defs[s.m_sled] = 500 + # check the fluid domain + @test Symbol(defs[s.src.port.ρ]) == Symbol(s.fluid.ρ) + @test Symbol(defs[s.valve.port_s.ρ]) == Symbol(s.fluid.ρ) + @test Symbol(defs[s.valve.port_a.ρ]) == Symbol(s.fluid.ρ) + @test Symbol(defs[s.valve.port_b.ρ]) == Symbol(s.fluid.ρ) + @test Symbol(defs[s.valve.port_r.ρ]) == Symbol(s.fluid.ρ) + @test Symbol(defs[s.snk.port.ρ]) == Symbol(s.fluid.ρ) -p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false)) -prob = remake(prob; p, tspan=(0, time[end])) + time = 0:dt:0.055 + x = @. 0.9 * (time > 0.015) * (time - 0.015)^2 - 25 * (time > 0.02) * (time - 0.02)^3 -sol2 = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); + defs[s.input.buffer] = Parameter(x, dt) + p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false)) + prob = remake(prob; p, tspan = (0, time[end])) + sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, + initializealg = NoInit()) -# Fast Solve ------------ -@register_symbolic xf(time) -Symbolics.derivative(::typeof(xf), args::NTuple{1, Any}, ::Val{1}) = 0 -@named system_ = System(false, xf) -s_ = complete(system_) -sys_ = structural_simplify(system_) -prob_ = ODEProblem(sys_, [], (0, 0.055), [s_.Cd => 0.005]) -solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); #precompile -@time "p is Vector{Float64}" sol_ = solve(prob_, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit()); -u = copy(prob_.u0) -prob_.f.f(u,prob_.u0, prob_.p, 0.0) -@time "f()" prob_.f.f(u,prob_.u0, prob_.p, 0.0) + @test sol[sys.ddx][1] == 0.0 + @test maximum(sol[sys.ddx]) > 500 + @test sol[sys.ddx][end] < 100 +end From 2606160237ea2ef58f21ee6939d81510745c6046 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 1 May 2023 11:08:31 -0400 Subject: [PATCH 11/16] fixed docs --- docs/pages.jl | 1 + docs/src/index.md | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/pages.jl b/docs/pages.jl index 24ded55eb..05b09ec83 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -5,6 +5,7 @@ pages = [ "Custom Components" => "tutorials/custom_component.md", "Thermal Conduction Model" => "tutorials/thermal_model.md", "DC Motor with Speed Controller" => "tutorials/dc_motor_pi.md", + "Input Component" => "tutorials/input_component.md", ], "About Acausal Connections" => "connectors/connections.md", "API" => [ diff --git a/docs/src/index.md b/docs/src/index.md index 4d7f19e30..296c41a1b 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -18,7 +18,6 @@ Pkg.add("ModelingToolkitStandardLibrary") - [Custom Component](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/tutorials/custom_component/) - [Thermal Model](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/tutorials/thermal_model/) - [DC Motor with PI-controller](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/tutorials/dc_motor_pi/) - - [Input Component](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/tutorials/input_component/) ## Libraries From dc76969a22fe6a4271d899018f070cd11e656c61 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 1 May 2023 12:14:43 -0400 Subject: [PATCH 12/16] fixed tests --- test/Hydraulic/isothermal_compressible.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index c8523b852..5d2075f05 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -242,6 +242,7 @@ end @test Symbol(defs[s.valve.port_r.ρ]) == Symbol(s.fluid.ρ) @test Symbol(defs[s.snk.port.ρ]) == Symbol(s.fluid.ρ) + dt = 1e-4 time = 0:dt:0.055 x = @. 0.9 * (time > 0.015) * (time - 0.015)^2 - 25 * (time > 0.02) * (time - 0.02)^3 From 57058b64a5b9800569e0ce6e49da81964f7eeb59 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 1 May 2023 13:12:15 -0400 Subject: [PATCH 13/16] fixed Mechanical test and Force component --- src/Mechanical/TranslationalPosition/sources.jl | 7 ++++--- test/Mechanical/translational.jl | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/Mechanical/TranslationalPosition/sources.jl b/src/Mechanical/TranslationalPosition/sources.jl index 02ad284df..0e11b4d20 100644 --- a/src/Mechanical/TranslationalPosition/sources.jl +++ b/src/Mechanical/TranslationalPosition/sources.jl @@ -6,7 +6,8 @@ Input signal acting as external force on a flange @component function Force(; name, use_support = false) @named partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support) @unpack flange = partial_element - @named f = RealInput() # Accelerating force acting at flange (= -flange.tau) - eqs = [flange.f ~ -f.u] - return extend(ODESystem(eqs, t, [], []; name = name, systems = [f]), partial_element) + @named input = RealInput() # Accelerating force acting at flange (= -flange.tau) + eqs = [flange.f ~ -input.u] + return extend(ODESystem(eqs, t, [], []; name = name, systems = [input]), + partial_element) end diff --git a/test/Mechanical/translational.jl b/test/Mechanical/translational.jl index b9729752a..b2e3bb3e8 100644 --- a/test/Mechanical/translational.jl +++ b/test/Mechanical/translational.jl @@ -65,7 +65,7 @@ end @named source = Sine(frequency = 3, amplitude = 2) function simplify_and_solve(damping, spring, body, ground, f, source) - eqs = [connect(f.f, source.output) + eqs = [connect(f.input, source.output) connect(f.flange, body.flange) connect(spring.flange_a, body.flange, damping.flange_a) connect(spring.flange_b, damping.flange_b, ground.flange)] From a48af53ce06ba838e32c8eaba4366b54519c56fc Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 1 May 2023 13:59:05 -0400 Subject: [PATCH 14/16] IfElse fix --- .../IsothermalCompressible/IsothermalCompressible.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl index 6f58bed1d..8ac4ffd0e 100644 --- a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl +++ b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl @@ -3,10 +3,13 @@ Library to model iso-thermal compressible liquid fluid flow """ module IsothermalCompressible -using ModelingToolkit, Symbolics, IfElse +using ModelingToolkit, Symbolics + using ...Blocks: RealInput, RealOutput using ...Mechanical.Translational: MechanicalPort, Mass +using IfElse: ifelse + @parameters t D = Differential(t) From 8ae44276aa1946bf19d439884ee21ebe208f8369 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 1 May 2023 14:31:34 -0400 Subject: [PATCH 15/16] fixed IfElse issue --- src/Hydraulic/IsothermalCompressible/components.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index 67a4831ff..a0c45a2e1 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -351,7 +351,7 @@ dm ────► dead volume │ │ area dm = port.dm u = dm / (rho * area) - eqs = [0 ~ IfElse.ifelse(vol >= minimum_volume, (p) - (port.p - (1 / 2) * rho * u^2), + eqs = [0 ~ ifelse(vol >= minimum_volume, (p) - (port.p - (1 / 2) * rho * u^2), dm) vol ~ dead_volume + area * x D(x) ~ dx From 7ca521c1c332cd2f1ced22df3f29c4564139911f Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 1 May 2023 14:35:05 -0400 Subject: [PATCH 16/16] format --- src/Hydraulic/IsothermalCompressible/components.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index a0c45a2e1..da787dae1 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -352,7 +352,7 @@ dm ────► dead volume │ │ area u = dm / (rho * area) eqs = [0 ~ ifelse(vol >= minimum_volume, (p) - (port.p - (1 / 2) * rho * u^2), - dm) + dm) vol ~ dead_volume + area * x D(x) ~ dx D(rho) ~ drho