Merge pull request #314 from SebastianM-C/smc/interp

Add `Interpolation` & `ParametrizedInterpolation`

Project.toml

@@ -9,33 +9,45 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 Aqua = "0.8"
 ChainRulesCore = "1.24"
 ControlSystemsBase = "1.4"
+DataFrames = "1.7"
 DataInterpolations = "6"
+ForwardDiff = "0.10"
 DiffEqBase = "6.152"
 IfElse = "0.1"
 LinearAlgebra = "1.10"
-ModelingToolkit = "9.46.1"
+Optimization = "4"
+ModelingToolkit = "9.47"
 OrdinaryDiffEq = "6.87"
 OrdinaryDiffEqDefault = "1.1"
+PreallocationTools = "0.4.23"
 SafeTestsets = "0.1"
+SciMLStructures = "1.4.2"
+SymbolicIndexingInterface = "0.3.28"
 Symbolics = "6.14"
 Test = "1"
 julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "LinearAlgebra", "OrdinaryDiffEq", "OrdinaryDiffEqDefault", "SafeTestsets", "Test", "ControlSystemsBase", "DataInterpolations"]
+test = ["Aqua", "LinearAlgebra", "OrdinaryDiffEqDefault", "OrdinaryDiffEq", "Optimization", "SafeTestsets", "Test", "ControlSystemsBase", "DataFrames", "DataInterpolations", "SciMLStructures", "SymbolicIndexingInterface", "ForwardDiff"]

docs/Project.toml

@@ -1,5 +1,7 @@
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
@@ -10,6 +12,8 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [compat]
 ControlSystemsBase = "1.1"
+DataFrames = "1.7"
+DataInterpolations = "6.4"
 DifferentialEquations = "7.6"
 Documenter = "1"
 IfElse = "0.1"

docs/src/tutorials/dc_motor_pi.md

@@ -77,7 +77,7 @@ so that it can be represented as a system of `ODEs` (ordinary differential equat
 
 ```@example dc_motor_pi
 sys = structural_simplify(model)
-prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0, 6.0))
+prob = ODEProblem(sys, [sys.L1.i => 0.0], (0, 6.0))
 sol = solve(prob)
 
 p1 = plot(sol.t, sol[sys.inertia.w], ylabel = "Angular Vel. in rad/s",

docs/src/tutorials/input_component.md

@@ -1,59 +1,171 @@
-# Running Models with Discrete Data
+# Building Models with Discrete Data, Interpolations, and Lookup Tables
 
-There are 3 ways to include data as part of a model.
+There are 4 ways to include data as part of a model.
 
- 1. using `ModelingToolkitStandardLibrary.Blocks.TimeVaryingFunction`
- 2. using a custom component with external data
- 3. using `ModelingToolkitStandardLibrary.Blocks.SampledData`
+ 1. using `ModelingToolkitStandardLibrary.Blocks.Interpolation`
+ 2. using `ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolation`
+ 3. using a custom component with external data (not recommended)
+ 4. using `ModelingToolkitStandardLibrary.Blocks.SampledData` (legacy)
 
 This tutorial demonstrate each case and explain the pros and cons of each.
 
-## `TimeVaryingFunction` Component
+## `Interpolation` Block
 
-The `ModelingToolkitStandardLibrary.Blocks.TimeVaryingFunction` component is easy to use and is performant.  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.
+The `ModelingToolkitStandardLibrary.Blocks.Interpolation` component is easy to use and is performant.
+It is similar to using callable parameters, but it provides a block interface with `RealInput` and `RealOutput` connectors.
+The `Interpolation` is compatible with interpolation types from `DataInterpolation`.
 
-```julia
+```@docs
+ModelingToolkitStandardLibrary.Blocks.Interpolation
+```
+
+Here is an example on how to use it. Let's consider a mass-spring-damper system, where
+we have an external force as an input. We then generate some example data in a `DataFrame`
+that would represent a measurement of the input. In a more realistic case, this `DataFrame`
+would be read from a file.
+
+```@example interpolation_block
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using ModelingToolkitStandardLibrary.Blocks
 using DataInterpolations
 using OrdinaryDiffEq
+using DataFrames
+using Plots
 
-function System(f; name)
-    src = TimeVaryingFunction(f)
+function MassSpringDamper(; name)
+    @named input = RealInput()
+    @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
+    @parameters m=10 k=1000 d=1
 
-    vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
+    eqs = [
+           f ~ input.u
+           ddx * 10 ~ k * x + d * dx + f
+           D(x) ~ dx
+           D(dx) ~ ddx]
+
+    ODESystem(eqs, t; name, systems = [input])
+end
+
+function MassSpringDamperSystem(data, time; name)
+    @named src = Interpolation(LinearInterpolation, data, time)
+    @named clk = ContinuousClock()
+    @named model = MassSpringDamper()
+
+    eqs = [
+        connect(src.input, clk.output)
+        connect(src.output, model.input)
+    ]
+
+    ODESystem(eqs, t, [], []; name, systems = [src, clk, model])
+end
+
+function generate_data()
+    dt = 4e-4
+    time = 0:dt:0.1
+    data = sin.(2 * pi * time * 100)
+
+    return DataFrame(; time, data)
+end
+
+df = generate_data() # example data
+
+@named system = MassSpringDamperSystem(df.data, df.time)
+sys = structural_simplify(system)
+prob = ODEProblem(sys, [], (0, df.time[end]))
+sol = solve(prob)
+plot(sol)
+```
+
+## `ParametrizedInterpolation` Block
+
+The `ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolation` component is similar to `Interpolation`, but as the name suggests, it is parametrized by the data, allowing one to change the underlying data without rebuilding the model as the data is represented via vector parameters.
+The main advantage of this block over the [`Interpolation`](@ref) one is that one can use it for optimization problems. Currently, this supports forward mode AD via ForwardDiff, but due to the increased flexibility of the types in the component, this is not as fast as the `Interpolation` block,
+so it is recommended to use only when the added flexibility is required.
+
+```@docs
+ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolation
+```
+
+Here is an example on how to use it
+
+```@example parametrized_interpolation
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using ModelingToolkitStandardLibrary.Blocks
+using DataInterpolations
+using OrdinaryDiffEq
+using DataFrames
+using Plots
+
+function MassSpringDamper(; name)
+    @named input = RealInput()
+    vars = @variables f(t) x(t)=0 dx(t) [guess=0] ddx(t)
     pars = @parameters m=10 k=1000 d=1
 
-    eqs = [f ~ src.output.u
+    eqs = [
+           f ~ input.u
            ddx * 10 ~ k * x + d * dx + f
            D(x) ~ dx
            D(dx) ~ ddx]
 
-    ODESystem(eqs, t, vars, pars; systems = [src], name)
+    ODESystem(eqs, t, vars, pars; name, systems = [input])
 end
 
-dt = 4e-4
-time = 0:dt:0.1
-data = sin.(2 * pi * time * 100) # example data
+function MassSpringDamperSystem(data, time; name)
+    @named src = ParametrizedInterpolation(LinearInterpolation, data, time)
+    @named clk = ContinuousClock()
+    @named model = MassSpringDamper()
+
+    eqs = [
+        connect(model.input, src.output)
+        connect(src.input, clk.output)
+    ]
 
-f = LinearInterpolation(data, time)
+    ODESystem(eqs, t; name, systems = [src, clk, model])
+end
+
+function generate_data()
+    dt = 4e-4
+    time = 0:dt:0.1
+    data = sin.(2 * pi * time * 100)
+
+    return DataFrame(; time, data)
+end
 
-@named system = System(f)
+df = generate_data() # example data
+
+@named system = MassSpringDamperSystem(df.data, df.time)
 sys = structural_simplify(system)
-prob = ODEProblem(sys, [], (0, time[end]))
-sol = solve(prob, ImplicitEuler())
+prob = ODEProblem(sys, [], (0, df.time[end]))
+sol = solve(prob)
+plot(sol)
+```
+
+If we want to run a new data set, this requires only remaking the problem and solving again
+```@example parametrized_interpolation
+prob2 = remake(prob, p = [sys.src.data => ones(length(df.data))])
+sol2 = solve(prob2)
+plot(sol2)
 ```
 
-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 several pieces of data it's better to re-use an `ODESystem`.  The next couple methods will demonstrate how to do this.
+!!! note
+    Note that when changing the data, the length of the new data must be the same as the length of the original data.
 
 ## Custom Component with External Data
 
 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
+```@example custom_component_external_data
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using ModelingToolkitStandardLibrary.Blocks
+using OrdinaryDiffEq
+
 const rdata = Ref{Vector{Float64}}()
 
+dt = 4e-4
+time = 0:dt:0.1
 # Data Sets
 data1 = sin.(2 * pi * time * 100)
 data2 = cos.(2 * pi * time * 50)
@@ -106,8 +218,13 @@ Additional code could be added to resolve this issue, for example by using a `Re
 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 parallelize the call to `solve()`.
 
 ```julia
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using ModelingToolkitStandardLibrary.Blocks
+using OrdinaryDiffEq
+
 function System(; name)
-    src = SampledData(Float64)
+    src = SampledData(Float64, name=:src)
 
     vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
     pars = @parameters m=10 k=1000 d=1
@@ -121,16 +238,22 @@ function System(; name)
 end
 
 @named system = System()
-sys = structural_simplify(system)
+sys = structural_simplify(system, split=false)
 s = complete(system)
-prob = ODEProblem(sys, [], (0, time[end]); tofloat = false)
+
+dt = 4e-4
+time = 0:dt:0.1
+data1 = sin.(2 * pi * time * 100)
+data2 = cos.(2 * pi * time * 50)
+
+prob = ODEProblem(sys, [], (0, time[end]); split=false, tofloat = false, use_union=true)
 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)
+    remake(prob; p, build_initializeprob=false)
 end
 
 prob1 = get_prob(data1)

src/Blocks/Blocks.jl

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