Add Interpolation & ParametrizedInterpolation

Project.toml

@@ -9,33 +9,42 @@ 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"
-DataInterpolations = "4.6"
+DataInterpolations = "4.6, 5, 6"
 DiffEqBase = "6.152"
 IfElse = "0.1"
 LinearAlgebra = "1.10"
-ModelingToolkit = "9.32"
+ModelingToolkit = "9.46"
+Optimization = "3.27"
 OrdinaryDiffEq = "6.87"
 OrdinaryDiffEqDefault = "1.1"
+PreallocationTools = "0.4.23"
 SafeTestsets = "0.1"
-Symbolics = "5.36, 6"
+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"
 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", "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",
@@ -92,8 +92,8 @@ plot(p1, p2, layout = (2, 1))
 When implementing and tuning a control system in simulation, it is a good practice to analyze the closed-loop properties and verify robustness of the closed-loop with respect to, e.g., modeling errors. To facilitate this, we added two analysis points to the set of connections above, more specifically, we added the analysis points named `:y` and `:u` to the connections (for more details on analysis points, see [Linear Analysis](@ref))
 
 ```julia
-connect(sys.speed_sensor.w, :y, feedback.input2)
-connect(sys.pi_controller.ctr_output, :u, source.V)
+connect(sys.speed_sensor.w, :y, sys.feedback.input2)
+connect(sys.pi_controller.ctr_output, :u, sys.source.V)
 ```
 
 one at the plant output (`:y`) and one at the plant input (`:u`). We may use these analysis points to calculate, e.g., sensitivity functions, illustrated below. Here, we calculate the sensitivity function $S(s)$ and the complimentary sensitivity function $T(s) = I - S(s)$, defined as

docs/src/tutorials/input_component.md

@@ -1,57 +1,162 @@
 # Running Models with Discrete Data
 
-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 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
+
+    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
 
-function System(f; name)
-    src = TimeVaryingFunction(f)
+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)=0 x(t)=0 dx(t)=0 ddx(t)=0
     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()
 
-f = LinearInterpolation(data, time)
+    eqs = [
+        connect(model.input, src.output)
+        connect(src.input, clk.output)
+    ]
 
-@named system = System(f)
+    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, 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 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.
+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(data))])
+sol2 = solve(prob2)
+plot(sol2)
+```
+
+!!! 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
 const rdata = Ref{Vector{Float64}}()
 
 # Data Sets
@@ -105,9 +210,9 @@ 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
+```@example sampled_data_component
 function System(; name)
-    src = SampledData(Float64)
+    @named 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
@@ -121,9 +226,9 @@ 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)
+prob = ODEProblem(sys, [], (0, time[end]); tofloat = false, use_union=true)
 defs = ModelingToolkit.defaults(sys)
 
 function get_prob(data)

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