Skip to content

Commit

Permalink
Merge pull request #165 from SciML/bgc/input_fix
Browse files Browse the repository at this point in the history
Input changed to SampledData
  • Loading branch information
YingboMa authored May 4, 2023
2 parents a3213be + e62e80f commit 4f90e60
Show file tree
Hide file tree
Showing 14 changed files with 386 additions and 206 deletions.
2 changes: 1 addition & 1 deletion docs/pages.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +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",
"SampledData Component" => "tutorials/input_component.md",
],
"About Acausal Connections" => "connectors/connections.md",
"API" => [
Expand Down
32 changes: 17 additions & 15 deletions docs/src/tutorials/input_component.md
Original file line number Diff line number Diff line change
@@ -1,16 +1,16 @@
# Running Models with Discrete Data

There are 2 ways to include data as part of a model.
There are 3 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`
2. using a custom component with external data
3. using `ModelingToolkitStandardLibrary.Blocks.SampledData`

This tutorial demonstrate each case and explain the pros and cons of each.

## TimeVaryingFunction Component
## `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.
The `ModelingToolkitStandardLibrary.Blocks.TimeVaryingFunction` 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 sampled discrete data.

```julia
using ModelingToolkit
Expand Down Expand Up @@ -47,11 +47,11 @@ 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.
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 how to do 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.
The below code shows how to include data using a `Ref` and registered `get_sampled_data` 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}}()
Expand All @@ -60,20 +60,20 @@ const rdata = Ref{Vector{Float64}}()
data1 = sin.(2 * pi * time * 100)
data2 = cos.(2 * pi * time * 50)

function input(t)
function get_sampled_data(t)
i = floor(Int, t / dt) + 1
x = rdata[][i]

return x
end

Symbolics.@register_symbolic input(t)
Symbolics.@register_symbolic get_sampled_data(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)
eqs = [f ~ get_sampled_data(t)
ddx * 10 ~ k * x + d * dx + f
D(x) ~ dx
D(dx) ~ ddx]
Expand All @@ -94,20 +94,22 @@ 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.
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`.

```julia
# the following test will fail
@test all(ddx1 .== sol1[sys.ddx]) #returns false
```

## Input Component
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.

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.
## `SampledData` Component

To resolve the issues presented above, the `ModelingToolkitStandardLibrary.Blocks.SampledData` 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 call to `solve()`.

```julia
function System(; name)
src = Input(Float64)
src = SampledData(Float64)

vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
pars = @parameters m=10 k=1000 d=1
Expand Down Expand Up @@ -144,4 +146,4 @@ sol2 = Ref{ODESolution}()
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`.
Note, in the above example, we can build the system with an empty `SampledData` component, only setting the expected data type: `@named src = SampledData(Float64)`. It's also possible to initialize the component with real sampled data: `@named src = SampledData(data, dt)`. Additionally note that before running an `ODEProblem` using the `SampledData` component, one must be careful about the parameter vector Type. The `SampledData` component contains a `buffer` parameter of type `Parameter`, therefore 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 a uniform `Parameter` type by running `p = Parameter.(p)`. This will wrap all the single values in a `Parameter` which will be mathmatically equivalent to a `Number`.
2 changes: 1 addition & 1 deletion src/Blocks/Blocks.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ export Log, Log10
include("math.jl")

export Constant, TimeVaryingFunction, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine,
Square, Triangular, Parameter, Input
Square, Triangular, Parameter, SampledData
include("sources.jl")

export Limiter, DeadZone, SlewRateLimiter
Expand Down
33 changes: 15 additions & 18 deletions src/Blocks/sources.jl
Original file line number Diff line number Diff line change
Expand Up @@ -488,18 +488,14 @@ function Parameter(x::T; tofloat = true) where {T <: Real}
end
Parameter(x::Vector{T}, dt::T) where {T <: Real} = Parameter(x, dt, length(x))

function input(t, memory::Parameter{T}) where {T}
function get_sampled_data(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
if T <: AbstractFloat
return T(NaN)
else
return zero(T)
end
Expand Down Expand Up @@ -528,27 +524,27 @@ end
get_sample_time(memory::Parameter) = memory.ref
Symbolics.@register_symbolic get_sample_time(memory)

Symbolics.@register_symbolic input(t, memory)
Symbolics.@register_symbolic get_sampled_data(t, memory)

function first_order_backwards_difference(t, memory)
Δt = get_sample_time(memory)
x1 = input(t, memory)
x0 = input(t - Δt, memory)
x1 = get_sampled_data(t, memory)
x0 = get_sampled_data(t - Δt, memory)

return (x1 - x0) / Δt
end

function Symbolics.derivative(::typeof(input), args::NTuple{2, Any}, ::Val{1})
function Symbolics.derivative(::typeof(get_sampled_data), 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))
SampledData(T::Type; name) = SampledData(T[], zero(T); name)
function SampledData(data::Vector{T}, dt::T; name) where {T <: Real}
SampledData(; name, buffer = Parameter(data, dt))
end

"""
Input(; name, buffer)
SampledData(; name, buffer)
data input component.
Expand All @@ -558,13 +554,14 @@ data input component.
# Connectors:
- `output`
"""
@component function Input(; name, buffer)
@component function SampledData(; name, buffer)
pars = @parameters begin buffer = buffer end
vars = []
systems = @named begin output = RealOutput() end
eqs = [
output.u ~ input(t, buffer),
output.u ~ get_sampled_data(t, buffer),
]
return ODESystem(eqs, t, vars, pars; name, systems,
defaults = [output.u => input(0.0, buffer)]) #TODO: get initial value from buffer
defaults = [output.u => get_sampled_data(0.0, buffer)])
end
@deprecate Input SampledData
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,10 @@ D = Differential(t)
export HydraulicPort, HydraulicFluid
include("utils.jl")

export Source, InputSource, Cap, Tube, FixedVolume, DynamicVolume
export Cap, Tube, FixedVolume, DynamicVolume
include("components.jl")

export MassFlow, Pressure, FixedPressure
include("sources.jl")

end
Loading

0 comments on commit 4f90e60

Please sign in to comment.