Merge pull request #163 from SciML/bgc/input

Input component and Hydraulic changes

Author: YingboMa

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

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" => [

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

src/Blocks/Blocks.jl

@@ -18,7 +18,7 @@ export Log, Log10
 include("math.jl")
 
 export Constant, TimeVaryingFunction, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine,
-       Square, Triangular
+       Square, Triangular, Parameter, Input
 include("sources.jl")
 
 export Limiter, DeadZone, SlewRateLimiter

src/Blocks/sources.jl

@@ -411,3 +411,160 @@ 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) 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 Parameter{T <: Real}
+    data::Vector{T}
+    ref::T
+    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
+
+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
+
+Parameter(x::Parameter) = x
+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{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
+
+    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 - 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
+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)
+function Input(data::Vector{T}, dt::T; name) where {T <: Real}
+    Input(; name, buffer = Parameter(data, dt))
+end
+
+"""
+    Input(; name, buffer)
+
+data input component.  
+
+# Parameters:
+  - `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 = []
+    systems = @named begin output = RealOutput() end
+    eqs = [
+        output.u ~ input(t, buffer),
+    ]
+    return ODESystem(eqs, t, vars, pars; name, systems,
+                     defaults = [output.u => input(0.0, buffer)]) #TODO: get initial value from buffer
+end

src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl

@@ -3,17 +3,20 @@ 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
+using ...Mechanical.Translational: MechanicalPort, Mass
+
+using IfElse: ifelse
 
 @parameters t
 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