From 0ba28ad1a26e6ac9cc992eae8b8d2be8d34e4e0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 2 Aug 2024 03:07:55 +0300 Subject: [PATCH 01/44] feat: add `ParametrizedInterpolation` Co-authored-by: Fredrik Bagge Carlson --- Project.toml | 8 +++++- ...kitStandardLibraryDataInterpolationsExt.jl | 10 +++++++ src/Blocks/Blocks.jl | 2 +- src/Blocks/sources.jl | 26 +++++++++++++++++++ 4 files changed, 44 insertions(+), 2 deletions(-) create mode 100644 ext/ModelingToolkitStandardLibraryDataInterpolationsExt.jl diff --git a/Project.toml b/Project.toml index b604de635..6084eb19f 100644 --- a/Project.toml +++ b/Project.toml @@ -11,11 +11,17 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +[weakdeps] +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" + +[extensions] +ModelingToolkitStandardLibraryDataInterpolationsExt = "DataInterpolations" + [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" diff --git a/ext/ModelingToolkitStandardLibraryDataInterpolationsExt.jl b/ext/ModelingToolkitStandardLibraryDataInterpolationsExt.jl new file mode 100644 index 000000000..992ec99b4 --- /dev/null +++ b/ext/ModelingToolkitStandardLibraryDataInterpolationsExt.jl @@ -0,0 +1,10 @@ +module ModelingToolkitStandardLibraryDataInterpolationsExt + +using ModelingToolkitStandardLibrary +using ModelingToolkitStandardLibrary.Blocks.Symbolics +using DataInterpolations + +@register_symbolic ModelingToolkitStandardLibrary.Blocks.apply_interpolation( + i::DataInterpolations.AbstractInterpolation, t) + +end diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl index 3eceab637..0bf733c92 100644 --- a/src/Blocks/Blocks.jl +++ b/src/Blocks/Blocks.jl @@ -16,7 +16,7 @@ export Log, Log10 include("math.jl") export Constant, TimeVaryingFunction, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine, - Square, Triangular, Parameter, SampledData + Square, Triangular, Parameter, SampledData, ParametrizedInterpolation include("sources.jl") export Limiter, DeadZone, SlewRateLimiter diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 5930be1e9..938ef021d 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -730,3 +730,29 @@ end function SampledData(; name, buffer, sample_time, circular_buffer) SampledData(SampledDataType.vector_based; name, buffer, sample_time, circular_buffer) end + +# This needs to be extend for interpolation types +apply_interpolation(interp, t) = interp(t) + + +@register_symbolic build_interpolation( + interpolation_type::UnionAll, u::AbstractArray, x::AbstractArray) +build_interpolation(interpolation_type, u, x) = interpolation_type(u, x) + +function ParametrizedInterpolation(interp_type::T, u, x, t = t; name) where {T} + @parameters data[1:length(x)] = u + @parameters ts[1:length(x)] = x + @parameters interpolation_type::T=interp_type [tunable = false] + @parameters interpolator::interp_type + + @named output = RealOutput() + + eqs = [output.u ~ apply_interpolation(interpolator, t)] + + ODESystem(eqs, t, [], [u, x, interpolation_type, interpolator]; + parameter_dependencies = [ + interpolator => build_interpolation(interpolation_type, u, x) + ], + systems = [output], + name) +end From 1b3d83e2b27d76feee5019b2184608bd80c6fa5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 5 Aug 2024 21:58:24 +0300 Subject: [PATCH 02/44] refactor: add derivative pass-through --- src/Blocks/sources.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 938ef021d..d49abea85 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -734,6 +734,9 @@ end # This needs to be extend for interpolation types apply_interpolation(interp, t) = interp(t) +function Symbolics.derivative(::typeof(apply_interpolation), args::NTuple{2, Any}, ::Val{2}) + Symbolics.derivative(args[1], (args[2],), Val(1)) +end @register_symbolic build_interpolation( interpolation_type::UnionAll, u::AbstractArray, x::AbstractArray) From aed04814cfeeb816d14c503bdd61d7f05a5a25a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 5 Aug 2024 22:02:53 +0300 Subject: [PATCH 03/44] test: add test for ParametrizedInterpolation --- test/Blocks/sources.jl | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index f03a6de4a..efab52627 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -5,6 +5,7 @@ using ModelingToolkitStandardLibrary.Blocks: smooth_sin, smooth_cos, smooth_damp smooth_square, smooth_step, smooth_ramp, smooth_triangular, triangular, square using OrdinaryDiffEq: ReturnCode.Success +using DataInterpolations @testset "Constant" begin @named src = Constant(k = 2) @@ -412,8 +413,6 @@ end end @testset "SampledData" begin - using DataInterpolations - dt = 4e-4 t_end = 10.0 time = 0:dt:t_end @@ -478,3 +477,20 @@ end @test sol[ddy][end]≈2 atol=1e-3 end end + +@testset "ParametrizedInterpolation" begin + @variables y(t) = 0 + @parameters u[1:15] = rand(15) + @parameters x[1:15] = 0:14 + + @named i = ParametrizedInterpolation(LinearInterpolation, u, x) + eqs = [D(y) ~ i.output.u] + + @named model = ODESystem(eqs, t, systems = [i]) + sys = structural_simplify(model) + + prob = ODEProblem(sys, [], (0.0, 4)) + sol = solve(prob) + + @test SciMLBase.successful_retcode(sol) +end From 770ea8dc5815f25aa091f5927fb3f6b63a0b76e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 5 Aug 2024 23:15:39 +0300 Subject: [PATCH 04/44] refactor: add support for arbitrary args in the interpolation constructor also use the standard MTK time --- src/Blocks/sources.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index d49abea85..a277a7170 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -739,22 +739,23 @@ function Symbolics.derivative(::typeof(apply_interpolation), args::NTuple{2, Any end @register_symbolic build_interpolation( - interpolation_type::UnionAll, u::AbstractArray, x::AbstractArray) -build_interpolation(interpolation_type, u, x) = interpolation_type(u, x) + interpolation_type::UnionAll, u::AbstractArray, x::AbstractArray, args::Tuple) +build_interpolation(interpolation_type, u, x, args) = interpolation_type(u, x, args...) -function ParametrizedInterpolation(interp_type::T, u, x, t = t; name) where {T} +function ParametrizedInterpolation(interp_type::T, u, x, args...; name) where {T} @parameters data[1:length(x)] = u @parameters ts[1:length(x)] = x - @parameters interpolation_type::T=interp_type [tunable = false] + @parameters interpolation_type::T=interp_type [tunable = false] interpolation_args::Tuple=args [tunable = false] @parameters interpolator::interp_type @named output = RealOutput() eqs = [output.u ~ apply_interpolation(interpolator, t)] - ODESystem(eqs, t, [], [u, x, interpolation_type, interpolator]; + ODESystem(eqs, t, [], [u, x, interpolation_type, interpolator, interpolation_args]; parameter_dependencies = [ - interpolator => build_interpolation(interpolation_type, u, x) + interpolator => build_interpolation( + interpolation_type, u, x, interpolation_args) ], systems = [output], name) From 89cc33322c8830edf2b6c7dfb5cac18da59f6106 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 5 Aug 2024 23:16:39 +0300 Subject: [PATCH 05/44] docs: add docstring for `ParametrizedInterpolation` --- src/Blocks/sources.jl | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index a277a7170..de768abbd 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -742,6 +742,26 @@ end interpolation_type::UnionAll, u::AbstractArray, x::AbstractArray, args::Tuple) build_interpolation(interpolation_type, u, x, args) = interpolation_type(u, x, args...) +""" + ParametrizedInterpolation(interp_type, u, x, args...; name) + +Represent function interpolation symbolically as a block component. By default interpolation types +from [`DataInterpolations.jl`](https://github.com/SciML/DataInterpolations.jl) are supported. +# Arguments: + - `interp_type`: the type of the interpolation. For `DataInterpolations`, +these would be any of [the available interpolations](https://github.com/SciML/DataInterpolations.jl?tab=readme-ov-file#available-interpolations), +such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. + - `u`: the data used for interpolation. For `DataInterpolations` this will be an `AbstractVector` + - `x`: the values that each data points correspond to, usually the times corresponding to each value in `u`. + - `args`: any other arguments beeded to build the interpolation + +# Parameters: + - `data`: the symbolic representation of the data passed at construction time via `u`. + - `ts`: the symbolic representation of times corresponding to the data passed at construction time via `x`. + +# Connectors: + - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value +""" function ParametrizedInterpolation(interp_type::T, u, x, args...; name) where {T} @parameters data[1:length(x)] = u @parameters ts[1:length(x)] = x From 2b82a13967492d0643f182bf878f7c562f80e247 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 5 Aug 2024 23:17:00 +0300 Subject: [PATCH 06/44] test: add test for `BSplineInterpolation` --- test/Blocks/sources.jl | 32 ++++++++++++++++++++++++-------- 1 file changed, 24 insertions(+), 8 deletions(-) diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index efab52627..be321c77b 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -481,16 +481,32 @@ end @testset "ParametrizedInterpolation" begin @variables y(t) = 0 @parameters u[1:15] = rand(15) - @parameters x[1:15] = 0:14 + @parameters x[1:15] = 0:14.0 - @named i = ParametrizedInterpolation(LinearInterpolation, u, x) - eqs = [D(y) ~ i.output.u] + @testset "LinearInterpolation" begin + @named i = ParametrizedInterpolation(LinearInterpolation, u, x) + eqs = [D(y) ~ i.output.u] - @named model = ODESystem(eqs, t, systems = [i]) - sys = structural_simplify(model) + @named model = ODESystem(eqs, t, systems = [i]) + sys = structural_simplify(model) - prob = ODEProblem(sys, [], (0.0, 4)) - sol = solve(prob) + prob = ODEProblem(sys, [], (0.0, 4)) + sol = solve(prob) - @test SciMLBase.successful_retcode(sol) + @test SciMLBase.successful_retcode(sol) + end + + @testset "BSplineInterpolation" begin + @named i = ParametrizedInterpolation( + BSplineInterpolation, u, x, 3, :Uniform, :Uniform) + eqs = [D(y) ~ i.output.u] + + @named model = ODESystem(eqs, t, systems = [i]) + sys = structural_simplify(model) + + prob = ODEProblem(sys, [], (0.0, 4)) + sol = solve(prob) + + @test SciMLBase.successful_retcode(sol) + end end From fd51e97519b0f8fda0a9727d1c34cfa61190c99f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Tue, 6 Aug 2024 01:55:11 +0300 Subject: [PATCH 07/44] docs: fix namespacing in DC motor tutorial --- docs/src/tutorials/dc_motor_pi.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorials/dc_motor_pi.md b/docs/src/tutorials/dc_motor_pi.md index 420eb43d4..5b5f6982e 100644 --- a/docs/src/tutorials/dc_motor_pi.md +++ b/docs/src/tutorials/dc_motor_pi.md @@ -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 From 339db5a47475bc67f6d891c4fc2ca622005df981 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 16 Aug 2024 03:28:31 +0300 Subject: [PATCH 08/44] refactor: cache the interpolation manually As of MTK@9.33 the parameter dependencies are no longer cached by MTKParameters, so we need to cache the interpolation object so that it's not recreated on every function call --- src/Blocks/sources.jl | 37 ++++++++++++++++++++++++++++++------- 1 file changed, 30 insertions(+), 7 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index de768abbd..342b7e810 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -738,9 +738,30 @@ function Symbolics.derivative(::typeof(apply_interpolation), args::NTuple{2, Any Symbolics.derivative(args[1], (args[2],), Val(1)) end -@register_symbolic build_interpolation( - interpolation_type::UnionAll, u::AbstractArray, x::AbstractArray, args::Tuple) -build_interpolation(interpolation_type, u, x, args) = interpolation_type(u, x, args...) +function cached_interpolation(interpolation_type, u, x, args) + prev_u = Ref(collect(copy(u))) + prev_x = Ref(collect(copy(x))) + # MTKParameters use views, so we want to ensure that the type is the same + interp = Ref(interpolation_type(prev_u[], prev_x[], args...)) + + let prev_u = prev_u, + prev_x = prev_x, + interp = interp, + interpolation_type = interpolation_type + function build_interpolation(u, x, args) + if (u, x) ≠ (prev_u[], prev_x[]) + prev_u[] = collect(u) + prev_x[] = collect(x) + interp[] = interpolation_type(prev_u[], prev_x[], args...) + else + interp[] + end + end + end +end + +@register_symbolic interpolation_builder(f::Function, u::AbstractArray, x::AbstractArray, args::Tuple) +interpolation_builder(f, u, x, args) = f(u, x, args) """ ParametrizedInterpolation(interp_type, u, x, args...; name) @@ -762,20 +783,22 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. # Connectors: - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value """ -function ParametrizedInterpolation(interp_type::T, u, x, args...; name) where {T} +function ParametrizedInterpolation(interp_type::T, u::AbstractVector, x::AbstractVector, args...; name) where {T} @parameters data[1:length(x)] = u @parameters ts[1:length(x)] = x @parameters interpolation_type::T=interp_type [tunable = false] interpolation_args::Tuple=args [tunable = false] @parameters interpolator::interp_type + build_interpolation = cached_interpolation(interp_type, u, x, args) + @parameters memoized_builder::typeof(build_interpolation)=build_interpolation [tunable = false] + @named output = RealOutput() eqs = [output.u ~ apply_interpolation(interpolator, t)] - ODESystem(eqs, t, [], [u, x, interpolation_type, interpolator, interpolation_args]; + ODESystem(eqs, t, [], [data, ts, interpolation_type, interpolation_args, interpolator, memoized_builder]; parameter_dependencies = [ - interpolator => build_interpolation( - interpolation_type, u, x, interpolation_args) + interpolator => interpolation_builder(memoized_builder, data, ts, interpolation_args) ], systems = [output], name) From 29dac5ffea9e03f50f57e2f08e664a814fcdb10a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 16 Aug 2024 04:47:31 +0300 Subject: [PATCH 09/44] refactor: make cached_interpolation work with Duals --- Project.toml | 2 ++ src/Blocks/sources.jl | 25 ++++++++++++++++--------- 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 6084eb19f..c3eaa51ce 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ 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" [weakdeps] @@ -28,6 +29,7 @@ LinearAlgebra = "1.10" ModelingToolkit = "9.32" OrdinaryDiffEq = "6.87" OrdinaryDiffEqDefault = "1.1" +PreallocationTools = "0.4.23" SafeTestsets = "0.1" Symbolics = "5.36, 6" Test = "1" diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 342b7e810..538db6547 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -1,5 +1,6 @@ using DiffEqBase import ChainRulesCore +using PreallocationTools # Define and register smooth functions # These are "smooth" aka differentiable and avoid Gibbs effect @@ -739,22 +740,28 @@ function Symbolics.derivative(::typeof(apply_interpolation), args::NTuple{2, Any end function cached_interpolation(interpolation_type, u, x, args) - prev_u = Ref(collect(copy(u))) - prev_x = Ref(collect(copy(x))) - # MTKParameters use views, so we want to ensure that the type is the same - interp = Ref(interpolation_type(prev_u[], prev_x[], args...)) + prev_u = DiffCache(u) + # Interpolation points can be a range, but we want to be able + # to update the cache if needed (and setindex! is not defined on ranges) + # with a view from MTKParameters, so we collect to get a vector + prev_x = DiffCache(collect(x)) + interp = GeneralLazyBufferCache() do (u,x) + interpolation_type(get_tmp(prev_u, u), get_tmp(prev_x, x), args...) + end let prev_u = prev_u, prev_x = prev_x, interp = interp, interpolation_type = interpolation_type + function build_interpolation(u, x, args) - if (u, x) ≠ (prev_u[], prev_x[]) - prev_u[] = collect(u) - prev_x[] = collect(x) - interp[] = interpolation_type(prev_u[], prev_x[], args...) + if (u, x) ≠ (get_tmp(prev_u, u), get_tmp(prev_x, x)) + get_tmp(prev_u, u) .= u + get_tmp(prev_x, x) .= x + interp.bufs[(u,x)] = interpolation_type( + get_tmp(prev_u, u), get_tmp(prev_x, x), args...) else - interp[] + interp[(u,x)] end end end From 4110a3542ed57bf837d76febafa884374d57bd9f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 16 Aug 2024 04:54:30 +0300 Subject: [PATCH 10/44] test: add tests for changing data in ParametrizedInterpolation --- Project.toml | 5 ++++- test/Blocks/sources.jl | 42 ++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 44 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index c3eaa51ce..db829932c 100644 --- a/Project.toml +++ b/Project.toml @@ -40,10 +40,13 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" 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", "OrdinaryDiffEq", "OrdinaryDiffEqDefault", "Optimization", "SafeTestsets", "Test", "ControlSystemsBase", "DataInterpolations", "SciMLStructures", "SymbolicIndexingInterface"] diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index be321c77b..da378e083 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -6,6 +6,9 @@ using ModelingToolkitStandardLibrary.Blocks: smooth_sin, smooth_cos, smooth_damp smooth_triangular, triangular, square using OrdinaryDiffEq: ReturnCode.Success using DataInterpolations +using SymbolicIndexingInterface +using SciMLStructures: SciMLStructures, Tunable +using Optimization @testset "Constant" begin @named src = Constant(k = 2) @@ -480,8 +483,8 @@ end @testset "ParametrizedInterpolation" begin @variables y(t) = 0 - @parameters u[1:15] = rand(15) - @parameters x[1:15] = 0:14.0 + u = rand(15) + x = 0:14.0 @testset "LinearInterpolation" begin @named i = ParametrizedInterpolation(LinearInterpolation, u, x) @@ -494,6 +497,41 @@ end sol = solve(prob) @test SciMLBase.successful_retcode(sol) + + prob2 = remake(prob, p=[i.data => ones(15)]) + sol2 = solve(prob2) + + @test SciMLBase.successful_retcode(sol2) + @test all(only.(sol2.u) .≈ sol2.t) # the solution for y' = 1 is y(t) = t + + set_data! = setp(prob2, i.data) + set_data!(prob2, zeros(15)) + sol3 = solve(prob2) + @test SciMLBase.successful_retcode(sol3) + @test iszero(sol3) + + function loss(x, p) + prob0, set_data! = p + ps = parameter_values(prob0) + arr, repack, alias = SciMLStructures.canonicalize(Tunable(), ps) + T = promote_type(eltype(x), eltype(arr)) + promoted_ps = SciMLStructures.replace(Tunable(), ps, T.(arr)) + prob = remake(prob0; p = promoted_ps) + + set_data!(prob, x) + sol = solve(prob) + sum(abs2.(only.(sol.u) .- sol.t)) + end + + set_data! = setp(prob, i.data) + of = OptimizationFunction(loss, AutoForwardDiff()) + op = OptimizationProblem(of, u, (prob, set_data!), lb = zeros(15), ub = fill(2.0, 15)) + + # check that type changing works + @test length(ForwardDiff.gradient(x -> of(x, (prob, set_data!)), u)) == 15 + + r = solve(op, Optimization.LBFGS(), maxiters = 1000) + @test of(r.u, (prob, set_data!)) < of(u, (prob, set_data!)) end @testset "BSplineInterpolation" begin From df9020943f5094d93f120761e9b459948ffc6d5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Tue, 27 Aug 2024 16:00:48 +0300 Subject: [PATCH 11/44] build: add compat for test deps --- Project.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Project.toml b/Project.toml index db829932c..0124a1518 100644 --- a/Project.toml +++ b/Project.toml @@ -27,10 +27,13 @@ DiffEqBase = "6.152" IfElse = "0.1" LinearAlgebra = "1.10" ModelingToolkit = "9.32" +Optimization = "3.27" OrdinaryDiffEq = "6.87" OrdinaryDiffEqDefault = "1.1" PreallocationTools = "0.4.23" SafeTestsets = "0.1" +SciMLStructures = "1.4.2" +SymbolicIndexingInterface = "0.3.28" Symbolics = "5.36, 6" Test = "1" julia = "1.10" From 9eed4da868da022a613e5fb65ec4b3a37ecc1e34 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Tue, 27 Aug 2024 16:02:37 +0300 Subject: [PATCH 12/44] style: fix formatting --- src/Blocks/sources.jl | 18 +++++++++++------- .../IsothermalCompressible.jl | 2 +- test/Blocks/sources.jl | 5 +++-- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 538db6547..9b17e33d1 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -745,7 +745,7 @@ function cached_interpolation(interpolation_type, u, x, args) # to update the cache if needed (and setindex! is not defined on ranges) # with a view from MTKParameters, so we collect to get a vector prev_x = DiffCache(collect(x)) - interp = GeneralLazyBufferCache() do (u,x) + interp = GeneralLazyBufferCache() do (u, x) interpolation_type(get_tmp(prev_u, u), get_tmp(prev_x, x), args...) end @@ -758,16 +758,17 @@ function cached_interpolation(interpolation_type, u, x, args) if (u, x) ≠ (get_tmp(prev_u, u), get_tmp(prev_x, x)) get_tmp(prev_u, u) .= u get_tmp(prev_x, x) .= x - interp.bufs[(u,x)] = interpolation_type( + interp.bufs[(u, x)] = interpolation_type( get_tmp(prev_u, u), get_tmp(prev_x, x), args...) else - interp[(u,x)] + interp[(u, x)] end end end end -@register_symbolic interpolation_builder(f::Function, u::AbstractArray, x::AbstractArray, args::Tuple) +@register_symbolic interpolation_builder( + f::Function, u::AbstractArray, x::AbstractArray, args::Tuple) interpolation_builder(f, u, x, args) = f(u, x, args) """ @@ -790,7 +791,8 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. # Connectors: - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value """ -function ParametrizedInterpolation(interp_type::T, u::AbstractVector, x::AbstractVector, args...; name) where {T} +function ParametrizedInterpolation( + interp_type::T, u::AbstractVector, x::AbstractVector, args...; name) where {T} @parameters data[1:length(x)] = u @parameters ts[1:length(x)] = x @parameters interpolation_type::T=interp_type [tunable = false] interpolation_args::Tuple=args [tunable = false] @@ -803,9 +805,11 @@ function ParametrizedInterpolation(interp_type::T, u::AbstractVector, x::Abstrac eqs = [output.u ~ apply_interpolation(interpolator, t)] - ODESystem(eqs, t, [], [data, ts, interpolation_type, interpolation_args, interpolator, memoized_builder]; + ODESystem(eqs, t, [], + [data, ts, interpolation_type, interpolation_args, interpolator, memoized_builder]; parameter_dependencies = [ - interpolator => interpolation_builder(memoized_builder, data, ts, interpolation_args) + interpolator => interpolation_builder( + memoized_builder, data, ts, interpolation_args) ], systems = [output], name) diff --git a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl index ce5e26907..3403ee022 100644 --- a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl +++ b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl @@ -1,5 +1,5 @@ """ -Library to model iso-thermal compressible liquid fluid flow +Library to model iso-thermal compressible liquid fluid flow """ module IsothermalCompressible diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index da378e083..af23dbff8 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -498,7 +498,7 @@ end @test SciMLBase.successful_retcode(sol) - prob2 = remake(prob, p=[i.data => ones(15)]) + prob2 = remake(prob, p = [i.data => ones(15)]) sol2 = solve(prob2) @test SciMLBase.successful_retcode(sol2) @@ -525,7 +525,8 @@ end set_data! = setp(prob, i.data) of = OptimizationFunction(loss, AutoForwardDiff()) - op = OptimizationProblem(of, u, (prob, set_data!), lb = zeros(15), ub = fill(2.0, 15)) + op = OptimizationProblem( + of, u, (prob, set_data!), lb = zeros(15), ub = fill(2.0, 15)) # check that type changing works @test length(ForwardDiff.gradient(x -> of(x, (prob, set_data!)), u)) == 15 From 748444a9bf4bc712a671174930e4a122ff06b5c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Tue, 27 Aug 2024 18:26:09 +0300 Subject: [PATCH 13/44] test: fix test dep for ParametrizedInterpolation --- Project.toml | 3 ++- test/Blocks/sources.jl | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0124a1518..edc2dce28 100644 --- a/Project.toml +++ b/Project.toml @@ -42,6 +42,7 @@ julia = "1.10" 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" @@ -52,4 +53,4 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "LinearAlgebra", "OrdinaryDiffEq", "OrdinaryDiffEqDefault", "Optimization", "SafeTestsets", "Test", "ControlSystemsBase", "DataInterpolations", "SciMLStructures", "SymbolicIndexingInterface"] +test = ["Aqua", "LinearAlgebra", "OrdinaryDiffEqDefault", "OrdinaryDiffEq", "Optimization", "SafeTestsets", "Test", "ControlSystemsBase", "DataInterpolations", "SciMLStructures", "SymbolicIndexingInterface", "ForwardDiff"] diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index af23dbff8..fd9ca9617 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -9,6 +9,7 @@ using DataInterpolations using SymbolicIndexingInterface using SciMLStructures: SciMLStructures, Tunable using Optimization +using ForwardDiff @testset "Constant" begin @named src = Constant(k = 2) From de3ccbd099d193fd9d360b81946aa69f5c4a19f8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Tue, 27 Aug 2024 18:26:32 +0300 Subject: [PATCH 14/44] docs: update docs to ParametrizedInterpolation --- docs/src/tutorials/input_component.md | 47 +++++++++++++++++---------- 1 file changed, 29 insertions(+), 18 deletions(-) diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md index 414a22eed..e085fca03 100644 --- a/docs/src/tutorials/input_component.md +++ b/docs/src/tutorials/input_component.md @@ -2,25 +2,29 @@ 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 - 3. using `ModelingToolkitStandardLibrary.Blocks.SampledData` + 1. using `ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolation` & `DataInterpolations` + 2. using a custom component with external data (not recommended) + 3. using `ModelingToolkitStandardLibrary.Blocks.SampledData` (legacy) This tutorial demonstrate each case and explain the pros and cons of each. -## `TimeVaryingFunction` Component +## `ParametrizedInterpolation` Component -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.ParametrizedInterpolation` component is easy to use and is performant. +It allows one to change the underlying data without rebuilding the model as the data is represented via vector parameters. +The `ParametrizedInterpolation` is compatible with interpolation types from `DataInterpolation`. +Here is an example on how to use it -```julia +```@example parametrized_interpolation using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D using ModelingToolkitStandardLibrary.Blocks using DataInterpolations using OrdinaryDiffEq +using Plots -function System(f; name) - src = TimeVaryingFunction(f) +function System(data, time; name) + @named src = ParametrizedInterpolation(LinearInterpolation, data, time) vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0 pars = @parameters m=10 k=1000 d=1 @@ -37,21 +41,28 @@ dt = 4e-4 time = 0:dt:0.1 data = sin.(2 * pi * time * 100) # example data -f = LinearInterpolation(data, time) - -@named system = System(f) +@named system = System(data, time) sys = structural_simplify(system) prob = ODEProblem(sys, [], (0, time[end])) -sol = solve(prob, ImplicitEuler()) +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(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 lenght 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 +116,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 +132,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) From 40dcd284dae7d6bc9f4bea50eed762bf4e5f20d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 27 Sep 2024 00:56:48 +0300 Subject: [PATCH 15/44] refactor: make use of callable parameters to simplify the implementation --- src/Blocks/sources.jl | 121 ++++++++++++++++++++++++++++-------------- 1 file changed, 81 insertions(+), 40 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 9b17e33d1..d7f9f57e1 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -733,43 +733,86 @@ function SampledData(; name, buffer, sample_time, circular_buffer) end # This needs to be extend for interpolation types -apply_interpolation(interp, t) = interp(t) +# apply_interpolation(interp, t) = interp(t) + +# function Symbolics.derivative(::typeof(apply_interpolation), args::NTuple{2, Any}, ::Val{2}) +# Symbolics.derivative(args[1], (args[2],), Val(1)) +# end + +# function cached_interpolation(interpolation_type, u, x, args) +# prev_u = DiffCache(u) +# # Interpolation points can be a range, but we want to be able +# # to update the cache if needed (and setindex! is not defined on ranges) +# # with a view from MTKParameters, so we collect to get a vector +# prev_x = DiffCache(collect(x)) +# interp = GeneralLazyBufferCache() do (u, x) +# interpolation_type(get_tmp(prev_u, u), get_tmp(prev_x, x), args...) +# end + +# let prev_u = prev_u, +# prev_x = prev_x, +# interp = interp, +# interpolation_type = interpolation_type + +# function build_interpolation(u, x, args) +# if (u, x) ≠ (get_tmp(prev_u, u), get_tmp(prev_x, x)) +# get_tmp(prev_u, u) .= u +# get_tmp(prev_x, x) .= x +# interp.bufs[(u, x)] = interpolation_type( +# get_tmp(prev_u, u), get_tmp(prev_x, x), args...) +# else +# interp[(u, x)] +# end +# end +# end +# end + +struct CachedInterpolation{T,I,U,X,C} + interpolation_type::I + prev_u::U + prev_x::X + cache::C + + function CachedInterpolation(interpolation_type, u, x, args) + # we need to copy the inputs to avoid aliasing + prev_u = DiffCache(copy(u)) + # Interpolation points can be a range, but we want to be able + # to update the cache if needed (and setindex! is not defined on ranges) + # with a view from MTKParameters, so we collect to get a vector + prev_x = DiffCache(collect(copy(x))) + cache = GeneralLazyBufferCache() do (u, x) + interpolation_type(get_tmp(prev_u, u), get_tmp(prev_x, x), args...) + end + T = typeof(cache[(get_tmp(prev_u, u), get_tmp(prev_x, x))]) + I = typeof(interpolation_type) + U = typeof(prev_u) + X = typeof(prev_x) + C = typeof(cache) -function Symbolics.derivative(::typeof(apply_interpolation), args::NTuple{2, Any}, ::Val{2}) - Symbolics.derivative(args[1], (args[2],), Val(1)) + new{T,I,U,X,C}(interpolation_type, prev_u, prev_x, cache) + end end -function cached_interpolation(interpolation_type, u, x, args) - prev_u = DiffCache(u) - # Interpolation points can be a range, but we want to be able - # to update the cache if needed (and setindex! is not defined on ranges) - # with a view from MTKParameters, so we collect to get a vector - prev_x = DiffCache(collect(x)) - interp = GeneralLazyBufferCache() do (u, x) - interpolation_type(get_tmp(prev_u, u), get_tmp(prev_x, x), args...) - end +function (f::CachedInterpolation{T})(u, x, args) where T + (;prev_u, prev_x, cache, interpolation_type) = f - let prev_u = prev_u, - prev_x = prev_x, - interp = interp, - interpolation_type = interpolation_type - - function build_interpolation(u, x, args) - if (u, x) ≠ (get_tmp(prev_u, u), get_tmp(prev_x, x)) - get_tmp(prev_u, u) .= u - get_tmp(prev_x, x) .= x - interp.bufs[(u, x)] = interpolation_type( - get_tmp(prev_u, u), get_tmp(prev_x, x), args...) - else - interp[(u, x)] - end - end + interp = @inbounds if (u, x) ≠ (get_tmp(prev_u, u), get_tmp(prev_x, x)) + get_tmp(prev_u, u) .= u + get_tmp(prev_x, x) .= x + # @info "cache miss" + cache.bufs[(u, x)] = interpolation_type( + get_tmp(prev_u, u), get_tmp(prev_x, x), args...) + else + # @info "cache hit" + cache[(u, x)] end + + return interp end -@register_symbolic interpolation_builder( - f::Function, u::AbstractArray, x::AbstractArray, args::Tuple) -interpolation_builder(f, u, x, args) = f(u, x, args) +Base.nameof(::CachedInterpolation) = :CachedInterpolation + +@register_symbolic (f::CachedInterpolation)(u::AbstractArray, x::AbstractArray, args::Tuple) """ ParametrizedInterpolation(interp_type, u, x, args...; name) @@ -792,24 +835,22 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value """ function ParametrizedInterpolation( - interp_type::T, u::AbstractVector, x::AbstractVector, args...; name) where {T} + interp_type::T, u::AbstractVector, x::AbstractVector, args...; name) where T + @parameters data[1:length(x)] = u @parameters ts[1:length(x)] = x - @parameters interpolation_type::T=interp_type [tunable = false] interpolation_args::Tuple=args [tunable = false] - @parameters interpolator::interp_type - - build_interpolation = cached_interpolation(interp_type, u, x, args) - @parameters memoized_builder::typeof(build_interpolation)=build_interpolation [tunable = false] + @parameters interpolation_type::T=interp_type [tunable = false] + build_interpolation = CachedInterpolation(interp_type, u, x, args) + @parameters (interpolator::interp_type)(..)::eltype(u) @named output = RealOutput() - eqs = [output.u ~ apply_interpolation(interpolator, t)] + eqs = [output.u ~ interpolator(t)] ODESystem(eqs, t, [], - [data, ts, interpolation_type, interpolation_args, interpolator, memoized_builder]; + [data, ts, interpolation_type, interpolator]; parameter_dependencies = [ - interpolator => interpolation_builder( - memoized_builder, data, ts, interpolation_args) + interpolator ~ build_interpolation(data, ts, args) ], systems = [output], name) From 9bc6f1d1a9b35b722b7d6cf91e873493e92b3075 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 27 Sep 2024 00:57:10 +0300 Subject: [PATCH 16/44] test: refactor test --- test/Blocks/sources.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index fd9ca9617..3b02bbac5 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -494,8 +494,8 @@ end @named model = ODESystem(eqs, t, systems = [i]) sys = structural_simplify(model) - prob = ODEProblem(sys, [], (0.0, 4)) - sol = solve(prob) + prob = ODEProblem{true, SciMLBase.FullSpecialize}(sys, [], (0.0, 4)) + sol = solve(prob, Tsit5()) @test SciMLBase.successful_retcode(sol) From 50c953976944ad631a3805793db7d369442553f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 13:56:31 +0300 Subject: [PATCH 17/44] feat: add InterpolationBlock --- src/Blocks/Blocks.jl | 3 +- src/Blocks/sources.jl | 74 +++++++++++++++++++++++-------------------- 2 files changed, 42 insertions(+), 35 deletions(-) diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl index 0bf733c92..0d2b839a5 100644 --- a/src/Blocks/Blocks.jl +++ b/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, ParametrizedInterpolation + Square, Triangular, Parameter, SampledData, + InterpolationBlock, ParametrizedInterpolationBlock include("sources.jl") export Limiter, DeadZone, SlewRateLimiter diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index d7f9f57e1..dd13bd424 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -732,40 +732,46 @@ function SampledData(; name, buffer, sample_time, circular_buffer) SampledData(SampledDataType.vector_based; name, buffer, sample_time, circular_buffer) end -# This needs to be extend for interpolation types -# apply_interpolation(interp, t) = interp(t) - -# function Symbolics.derivative(::typeof(apply_interpolation), args::NTuple{2, Any}, ::Val{2}) -# Symbolics.derivative(args[1], (args[2],), Val(1)) -# end - -# function cached_interpolation(interpolation_type, u, x, args) -# prev_u = DiffCache(u) -# # Interpolation points can be a range, but we want to be able -# # to update the cache if needed (and setindex! is not defined on ranges) -# # with a view from MTKParameters, so we collect to get a vector -# prev_x = DiffCache(collect(x)) -# interp = GeneralLazyBufferCache() do (u, x) -# interpolation_type(get_tmp(prev_u, u), get_tmp(prev_x, x), args...) -# end - -# let prev_u = prev_u, -# prev_x = prev_x, -# interp = interp, -# interpolation_type = interpolation_type - -# function build_interpolation(u, x, args) -# if (u, x) ≠ (get_tmp(prev_u, u), get_tmp(prev_x, x)) -# get_tmp(prev_u, u) .= u -# get_tmp(prev_x, x) .= x -# interp.bufs[(u, x)] = interpolation_type( -# get_tmp(prev_u, u), get_tmp(prev_x, x), args...) -# else -# interp[(u, x)] -# end -# end -# end -# end +""" + InterpolationBlock(interp_type, u, x, args...; name, t = ModelingToolkit.t_nounits) + +Represent function interpolation symbolically as a block component. +By default interpolation types from [`DataInterpolations.jl`](https://github.com/SciML/DataInterpolations.jl) are supported, +but in general any callable type that builds the interpolation object via `itp = interpolation_type(u, x, args...)` and calls +the interpolation with `itp(t)` should work. This does not need to represent an interpolation, it can be any type that satisfies +the interface, such as lookup tables. +# Arguments: + - `interp_type`: the type of the interpolation. For `DataInterpolations`, +these would be any of [the available interpolations](https://github.com/SciML/DataInterpolations.jl?tab=readme-ov-file#available-interpolations), +such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. + - `u`: the data used for interpolation. For `DataInterpolations` this will be an `AbstractVector` + - `x`: the values that each data points correspond to, usually the times corresponding to each value in `u`. + - `args`: any other arguments beeded to build the interpolation +# Keyword arguments: + - `name`: the name of the component + - `t`: the interpolation parameter, this is the time (`ModelingToolkit.t_nounits`) by default + +# Parameters: + - `interpolator`: the symbolic + - `t`: the parameter used for interpolation + +# Connectors: + - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value +""" +function InterpolationBlock(interp_type, u, x, args...; name, t = ModelingToolkit.t_nounits) + itp = interp_type(u, x, args...) + InterpolationBlock(itp; name, t) +end + +function InterpolationBlock(itp; name, t = ModelingToolkit.t_nounits) + @parameters (interpolator::typeof(itp))(..) = itp + + @named output = RealOutput() + + eqs = [output.u ~ interpolator(t)] + + ODESystem(eqs, ModelingToolkit.t_nounits, [], [interpolator, t]; name, systems = [output]) +end struct CachedInterpolation{T,I,U,X,C} interpolation_type::I From 367a4c91f13472a2e557e47dd7e2c01d1c7e59de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 13:57:42 +0300 Subject: [PATCH 18/44] refactor: remove DataInterpolations extension --- Project.toml | 6 ------ ...elingToolkitStandardLibraryDataInterpolationsExt.jl | 10 ---------- 2 files changed, 16 deletions(-) delete mode 100644 ext/ModelingToolkitStandardLibraryDataInterpolationsExt.jl diff --git a/Project.toml b/Project.toml index edc2dce28..0110426cf 100644 --- a/Project.toml +++ b/Project.toml @@ -12,12 +12,6 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" -[weakdeps] -DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" - -[extensions] -ModelingToolkitStandardLibraryDataInterpolationsExt = "DataInterpolations" - [compat] Aqua = "0.8" ChainRulesCore = "1.24" diff --git a/ext/ModelingToolkitStandardLibraryDataInterpolationsExt.jl b/ext/ModelingToolkitStandardLibraryDataInterpolationsExt.jl deleted file mode 100644 index 992ec99b4..000000000 --- a/ext/ModelingToolkitStandardLibraryDataInterpolationsExt.jl +++ /dev/null @@ -1,10 +0,0 @@ -module ModelingToolkitStandardLibraryDataInterpolationsExt - -using ModelingToolkitStandardLibrary -using ModelingToolkitStandardLibrary.Blocks.Symbolics -using DataInterpolations - -@register_symbolic ModelingToolkitStandardLibrary.Blocks.apply_interpolation( - i::DataInterpolations.AbstractInterpolation, t) - -end From bd276ececdb65f5d72bf27c2fd10dd412df796ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 13:58:21 +0300 Subject: [PATCH 19/44] refactor: add custom `t` for ParametrizedInterpolationBlock --- src/Blocks/sources.jl | 29 +++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index dd13bd424..8b570e1d9 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -816,15 +816,22 @@ function (f::CachedInterpolation{T})(u, x, args) where T return interp end +# function Symbolics.derivative(::typeof(CachedInterpolation), args::NTuple{2, Any}, ::Val{2}) +# Symbolics.derivative(args[1], (args[2],), Val(1)) +# end + Base.nameof(::CachedInterpolation) = :CachedInterpolation @register_symbolic (f::CachedInterpolation)(u::AbstractArray, x::AbstractArray, args::Tuple) """ - ParametrizedInterpolation(interp_type, u, x, args...; name) + ParametrizedInterpolationBlock(interp_type, u, x, args...; name, t = ModelingToolkit.t_nounits) -Represent function interpolation symbolically as a block component. By default interpolation types -from [`DataInterpolations.jl`](https://github.com/SciML/DataInterpolations.jl) are supported. +Represent function interpolation symbolically as a block component, with the interpolation data represented parametrically. +By default interpolation types from [`DataInterpolations.jl`](https://github.com/SciML/DataInterpolations.jl) are supported, +but in general any callable type that builds the interpolation object via `itp = interpolation_type(u, x, args...)` and calls +the interpolation with `itp(t)` should work. This does not need to represent an interpolation, it can be any type that satisfies +the interface, such as lookup tables. # Arguments: - `interp_type`: the type of the interpolation. For `DataInterpolations`, these would be any of [the available interpolations](https://github.com/SciML/DataInterpolations.jl?tab=readme-ov-file#available-interpolations), @@ -832,29 +839,35 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. - `u`: the data used for interpolation. For `DataInterpolations` this will be an `AbstractVector` - `x`: the values that each data points correspond to, usually the times corresponding to each value in `u`. - `args`: any other arguments beeded to build the interpolation +# Keyword arguments: + - `name`: the name of the component + - `t`: the interpolation parameter, this is the time (`ModelingToolkit.t_nounits`) by default # Parameters: - `data`: the symbolic representation of the data passed at construction time via `u`. - `ts`: the symbolic representation of times corresponding to the data passed at construction time via `x`. + - `t`: the parameter used for interpolation # Connectors: - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value """ -function ParametrizedInterpolation( - interp_type::T, u::AbstractVector, x::AbstractVector, args...; name) where T +function ParametrizedInterpolationBlock( + interp_type::T, u::AbstractVector, x::AbstractVector, args...; + name, t = ModelingToolkit.t_nounits) where {T} + + build_interpolation = CachedInterpolation(interp_type, u, x, args) @parameters data[1:length(x)] = u @parameters ts[1:length(x)] = x @parameters interpolation_type::T=interp_type [tunable = false] - build_interpolation = CachedInterpolation(interp_type, u, x, args) @parameters (interpolator::interp_type)(..)::eltype(u) @named output = RealOutput() eqs = [output.u ~ interpolator(t)] - ODESystem(eqs, t, [], - [data, ts, interpolation_type, interpolator]; + ODESystem(eqs, ModelingToolkit.t_nounits, [], + [data, ts, interpolation_type, interpolator, t]; parameter_dependencies = [ interpolator ~ build_interpolation(data, ts, args) ], From 3cf6e6cf4d93e740e4dc02b84c74eb72276928ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 13:58:41 +0300 Subject: [PATCH 20/44] test: add tests for `InterpolationBlock` --- test/Blocks/sources.jl | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index 3b02bbac5..ca4ee1065 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -482,13 +482,30 @@ end end end -@testset "ParametrizedInterpolation" begin +@testset "InterpolationBlock" begin + @variables y(t) = 0 + u = rand(15) + x = 0:14.0 + + @named i = InterpolationBlock(LinearInterpolation, u, x) + eqs = [D(y) ~ i.output.u] + + @named model = ODESystem(eqs, t, systems = [i]) + sys = structural_simplify(model) + + prob = ODEProblem{true, SciMLBase.FullSpecialize}(sys, [], (0.0, 4)) + sol = solve(prob, Tsit5()) + + @test SciMLBase.successful_retcode(sol) +end + +@testset "ParametrizedInterpolationBlock" begin @variables y(t) = 0 u = rand(15) x = 0:14.0 @testset "LinearInterpolation" begin - @named i = ParametrizedInterpolation(LinearInterpolation, u, x) + @named i = ParametrizedInterpolationBlock(LinearInterpolation, u, x) eqs = [D(y) ~ i.output.u] @named model = ODESystem(eqs, t, systems = [i]) @@ -537,7 +554,7 @@ end end @testset "BSplineInterpolation" begin - @named i = ParametrizedInterpolation( + @named i = ParametrizedInterpolationBlock( BSplineInterpolation, u, x, 3, :Uniform, :Uniform) eqs = [D(y) ~ i.output.u] From 223826ff8db852dc6c85dee7aa2d0c2b515fef7a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 13:59:05 +0300 Subject: [PATCH 21/44] docs: add docs for int4erpolation blocks --- docs/src/tutorials/input_component.md | 56 +++++++++++++++++++++++---- 1 file changed, 48 insertions(+), 8 deletions(-) diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md index e085fca03..ba80f7fde 100644 --- a/docs/src/tutorials/input_component.md +++ b/docs/src/tutorials/input_component.md @@ -1,18 +1,58 @@ # 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.ParametrizedInterpolation` & `DataInterpolations` - 2. using a custom component with external data (not recommended) - 3. using `ModelingToolkitStandardLibrary.Blocks.SampledData` (legacy) + 1. using `ModelingToolkitStandardLibrary.Blocks.InterpolationBlock` + 2. using `ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolationBlock` + 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. -## `ParametrizedInterpolation` Component +## `InterpolationBlock` Component -The `ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolation` component is easy to use and is performant. -It allows one to change the underlying data without rebuilding the model as the data is represented via vector parameters. -The `ParametrizedInterpolation` is compatible with interpolation types from `DataInterpolation`. +The `ModelingToolkitStandardLibrary.Blocks.InterpolationBlock` component is easy to use and is performant. +It is simlar to using callable paramterers, but it provides a block interface and a `RealOutput` connector. +The `InterpolationBlock` is compatible with interpolation types from `DataInterpolation`. +Here is an example on how to use it + +```@example interpolation_block +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D +using ModelingToolkitStandardLibrary.Blocks +using DataInterpolations +using OrdinaryDiffEq +using Plots + +function System(data, time; name) + @named src = InterpolationBlock(LinearInterpolation, data, time) + + 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 + +@named system = System(data, time) +sys = structural_simplify(system) +prob = ODEProblem(sys, [], (0, time[end])) +sol = solve(prob) +plot(sol) +``` + +## `ParametrizedInterpolationBlock` Component + +The `ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolationBlock` component is similar to `InterpolationBlock`, 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 `ParametrizedInterpolationBlock` is compatible with interpolation types from `DataInterpolation`. Here is an example on how to use it ```@example parametrized_interpolation From 6a9707245a169fabba31307853053978b4f88f7f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 21:55:36 +0300 Subject: [PATCH 22/44] refactor: revert renaming for the interpolation components Co-authored-by: Fredrik Bagge Carlson --- src/Blocks/Blocks.jl | 2 +- src/Blocks/sources.jl | 11 +++++------ test/Blocks/sources.jl | 8 ++++---- 3 files changed, 10 insertions(+), 11 deletions(-) diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl index 0d2b839a5..adb5db6e5 100644 --- a/src/Blocks/Blocks.jl +++ b/src/Blocks/Blocks.jl @@ -17,7 +17,7 @@ include("math.jl") export Constant, TimeVaryingFunction, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine, Square, Triangular, Parameter, SampledData, - InterpolationBlock, ParametrizedInterpolationBlock + Interpolation, ParametrizedInterpolation include("sources.jl") export Limiter, DeadZone, SlewRateLimiter diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 8b570e1d9..feaa4154d 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -733,7 +733,7 @@ function SampledData(; name, buffer, sample_time, circular_buffer) end """ - InterpolationBlock(interp_type, u, x, args...; name, t = ModelingToolkit.t_nounits) + Interpolation(interp_type, u, x, args...; name) Represent function interpolation symbolically as a block component. By default interpolation types from [`DataInterpolations.jl`](https://github.com/SciML/DataInterpolations.jl) are supported, @@ -749,7 +749,6 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. - `args`: any other arguments beeded to build the interpolation # Keyword arguments: - `name`: the name of the component - - `t`: the interpolation parameter, this is the time (`ModelingToolkit.t_nounits`) by default # Parameters: - `interpolator`: the symbolic @@ -758,12 +757,12 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. # Connectors: - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value """ -function InterpolationBlock(interp_type, u, x, args...; name, t = ModelingToolkit.t_nounits) +function Interpolation(interp_type, u, x, args...; name) itp = interp_type(u, x, args...) InterpolationBlock(itp; name, t) end -function InterpolationBlock(itp; name, t = ModelingToolkit.t_nounits) +function Interpolation(itp; name) @parameters (interpolator::typeof(itp))(..) = itp @named output = RealOutput() @@ -825,7 +824,7 @@ Base.nameof(::CachedInterpolation) = :CachedInterpolation @register_symbolic (f::CachedInterpolation)(u::AbstractArray, x::AbstractArray, args::Tuple) """ - ParametrizedInterpolationBlock(interp_type, u, x, args...; name, t = ModelingToolkit.t_nounits) + ParametrizedInterpolation(interp_type, u, x, args...; name, t = ModelingToolkit.t_nounits) Represent function interpolation symbolically as a block component, with the interpolation data represented parametrically. By default interpolation types from [`DataInterpolations.jl`](https://github.com/SciML/DataInterpolations.jl) are supported, @@ -851,7 +850,7 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. # Connectors: - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value """ -function ParametrizedInterpolationBlock( +function ParametrizedInterpolation( interp_type::T, u::AbstractVector, x::AbstractVector, args...; name, t = ModelingToolkit.t_nounits) where {T} diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index ca4ee1065..1329523e1 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -482,12 +482,12 @@ end end end -@testset "InterpolationBlock" begin +@testset "Interpolation" begin @variables y(t) = 0 u = rand(15) x = 0:14.0 - @named i = InterpolationBlock(LinearInterpolation, u, x) + @named i = Interpolation(LinearInterpolation, u, x) eqs = [D(y) ~ i.output.u] @named model = ODESystem(eqs, t, systems = [i]) @@ -499,13 +499,13 @@ end @test SciMLBase.successful_retcode(sol) end -@testset "ParametrizedInterpolationBlock" begin +@testset "ParametrizedInterpolation" begin @variables y(t) = 0 u = rand(15) x = 0:14.0 @testset "LinearInterpolation" begin - @named i = ParametrizedInterpolationBlock(LinearInterpolation, u, x) + @named i = ParametrizedInterpolation(LinearInterpolation, u, x) eqs = [D(y) ~ i.output.u] @named model = ODESystem(eqs, t, systems = [i]) From 9cbae31837bfa7fa90e4a5180b0ef5a3b5454f30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 21:56:15 +0300 Subject: [PATCH 23/44] refactor: add inputs to the interpolation components Co-authored-by: Fredrik Bagge Carlson --- src/Blocks/sources.jl | 24 ++++++++++++------------ test/Blocks/sources.jl | 6 +++--- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index feaa4154d..e670b36ce 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -752,24 +752,25 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. # Parameters: - `interpolator`: the symbolic - - `t`: the parameter used for interpolation # Connectors: + - `input`: a [`RealInput`](@ref) connector corresponding to the independent variable - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value """ function Interpolation(interp_type, u, x, args...; name) itp = interp_type(u, x, args...) - InterpolationBlock(itp; name, t) + Interpolation(itp; name) end function Interpolation(itp; name) @parameters (interpolator::typeof(itp))(..) = itp - + @named input = RealInput() @named output = RealOutput() - eqs = [output.u ~ interpolator(t)] + eqs = [output.u ~ interpolator(input.u)] - ODESystem(eqs, ModelingToolkit.t_nounits, [], [interpolator, t]; name, systems = [output]) + ODESystem( + eqs, t, [], [interpolator]; name, systems = [input, output]) end struct CachedInterpolation{T,I,U,X,C} @@ -840,20 +841,18 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. - `args`: any other arguments beeded to build the interpolation # Keyword arguments: - `name`: the name of the component - - `t`: the interpolation parameter, this is the time (`ModelingToolkit.t_nounits`) by default # Parameters: - `data`: the symbolic representation of the data passed at construction time via `u`. - `ts`: the symbolic representation of times corresponding to the data passed at construction time via `x`. - - `t`: the parameter used for interpolation # Connectors: + - `input`: a [`RealInput`](@ref) connector corresponding to the independent variable - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value """ function ParametrizedInterpolation( interp_type::T, u::AbstractVector, x::AbstractVector, args...; - name, t = ModelingToolkit.t_nounits) where {T} - + name) where {T} build_interpolation = CachedInterpolation(interp_type, u, x, args) @parameters data[1:length(x)] = u @@ -861,15 +860,16 @@ function ParametrizedInterpolation( @parameters interpolation_type::T=interp_type [tunable = false] @parameters (interpolator::interp_type)(..)::eltype(u) + @named input = RealInput() @named output = RealOutput() - eqs = [output.u ~ interpolator(t)] + eqs = [output.u ~ interpolator(input.u)] ODESystem(eqs, ModelingToolkit.t_nounits, [], - [data, ts, interpolation_type, interpolator, t]; + [data, ts, interpolation_type, interpolator]; parameter_dependencies = [ interpolator ~ build_interpolation(data, ts, args) ], - systems = [output], + systems = [input, output], name) end diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index 1329523e1..4c1797265 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -488,7 +488,7 @@ end x = 0:14.0 @named i = Interpolation(LinearInterpolation, u, x) - eqs = [D(y) ~ i.output.u] + eqs = [i.input.u ~ t, D(y) ~ i.output.u] @named model = ODESystem(eqs, t, systems = [i]) sys = structural_simplify(model) @@ -506,7 +506,7 @@ end @testset "LinearInterpolation" begin @named i = ParametrizedInterpolation(LinearInterpolation, u, x) - eqs = [D(y) ~ i.output.u] + eqs = [i.input.u ~ t, D(y) ~ i.output.u] @named model = ODESystem(eqs, t, systems = [i]) sys = structural_simplify(model) @@ -556,7 +556,7 @@ end @testset "BSplineInterpolation" begin @named i = ParametrizedInterpolationBlock( BSplineInterpolation, u, x, 3, :Uniform, :Uniform) - eqs = [D(y) ~ i.output.u] + eqs = [i.input.u ~ t, D(y) ~ i.output.u] @named model = ODESystem(eqs, t, systems = [i]) sys = structural_simplify(model) From 59c62fc78c7d03426749409c851e1e8141fe22ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 21:57:10 +0300 Subject: [PATCH 24/44] docs: update docs for interpolations Co-authored-by: Fredrik Bagge Carlson --- docs/Project.toml | 4 + docs/src/tutorials/input_component.md | 118 +++++++++++++++++++------- 2 files changed, 90 insertions(+), 32 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index dffb4fce9..2c3701147 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md index ba80f7fde..779421bad 100644 --- a/docs/src/tutorials/input_component.md +++ b/docs/src/tutorials/input_component.md @@ -2,19 +2,27 @@ There are 4 ways to include data as part of a model. - 1. using `ModelingToolkitStandardLibrary.Blocks.InterpolationBlock` - 2. using `ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolationBlock` + 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. -## `InterpolationBlock` Component +## `Interpolation` Block -The `ModelingToolkitStandardLibrary.Blocks.InterpolationBlock` component is easy to use and is performant. -It is simlar to using callable paramterers, but it provides a block interface and a `RealOutput` connector. -The `InterpolationBlock` is compatible with interpolation types from `DataInterpolation`. -Here is an example on how to use it +The `ModelingToolkitStandardLibrary.Blocks.Interpolation` component is easy to use and is performant. +It is simlar to using callable paramterers, but it provides a block interface with `RealInput` and `RealOutput` connectors. +The `Interpolation` is compatible with interpolation types from `DataInterpolation`. + +```@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 @@ -22,37 +30,63 @@ using ModelingToolkit: t_nounits as t, D_nounits as D using ModelingToolkitStandardLibrary.Blocks using DataInterpolations using OrdinaryDiffEq +using DataFrames using Plots -function System(data, time; name) - @named src = InterpolationBlock(LinearInterpolation, data, time) +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 - 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, [], []; 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 = Interpolation(LinearInterpolation, data, time) + @named clk = ContinuousClock() + @named model = MassSpringDamper() -@named system = System(data, time) + 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, time[end])) +prob = ODEProblem(sys, [], (0, df.time[end])) sol = solve(prob) plot(sol) ``` -## `ParametrizedInterpolationBlock` Component +## `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 +``` -The `ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolationBlock` component is similar to `InterpolationBlock`, 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 `ParametrizedInterpolationBlock` is compatible with interpolation types from `DataInterpolation`. Here is an example on how to use it ```@example parametrized_interpolation @@ -61,29 +95,49 @@ using ModelingToolkit: t_nounits as t, D_nounits as D using ModelingToolkitStandardLibrary.Blocks using DataInterpolations using OrdinaryDiffEq +using DataFrames using Plots -function System(data, time; name) - @named src = ParametrizedInterpolation(LinearInterpolation, data, time) - +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 + +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) + ] + + ODESystem(eqs, t, [], []; name, systems = [src, clk, model]) end -dt = 4e-4 -time = 0:dt:0.1 -data = sin.(2 * pi * time * 100) # example data +function generate_data() + dt = 4e-4 + time = 0:dt:0.1 + data = sin.(2 * pi * time * 100) -@named system = System(data, time) + 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])) +prob = ODEProblem(sys, [], (0, df.time[end])) sol = solve(prob) plot(sol) ``` From c8abefea494ef8934330f31b72b98f37813fd891 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 21:57:32 +0300 Subject: [PATCH 25/44] docs: fix dc motor initialization --- docs/src/tutorials/dc_motor_pi.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/tutorials/dc_motor_pi.md b/docs/src/tutorials/dc_motor_pi.md index 5b5f6982e..864dfa93e 100644 --- a/docs/src/tutorials/dc_motor_pi.md +++ b/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", From 57a101ee88f9f1dc3c216baad28afeffe901e8e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 22:04:04 +0300 Subject: [PATCH 26/44] style: fix formatting --- src/Blocks/sources.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index e670b36ce..8659428e9 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -773,7 +773,7 @@ function Interpolation(itp; name) eqs, t, [], [interpolator]; name, systems = [input, output]) end -struct CachedInterpolation{T,I,U,X,C} +struct CachedInterpolation{T, I, U, X, C} interpolation_type::I prev_u::U prev_x::X @@ -795,12 +795,12 @@ struct CachedInterpolation{T,I,U,X,C} X = typeof(prev_x) C = typeof(cache) - new{T,I,U,X,C}(interpolation_type, prev_u, prev_x, cache) + new{T, I, U, X, C}(interpolation_type, prev_u, prev_x, cache) end end -function (f::CachedInterpolation{T})(u, x, args) where T - (;prev_u, prev_x, cache, interpolation_type) = f +function (f::CachedInterpolation{T})(u, x, args) where {T} + (; prev_u, prev_x, cache, interpolation_type) = f interp = @inbounds if (u, x) ≠ (get_tmp(prev_u, u), get_tmp(prev_x, x)) get_tmp(prev_u, u) .= u From d26935b09b2450a651ac23fd729dd05909f06a89 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 3 Oct 2024 22:04:12 +0300 Subject: [PATCH 27/44] chore: fix typos --- docs/src/tutorials/dc_motor_pi.md | 2 +- docs/src/tutorials/input_component.md | 4 ++-- test/Blocks/sources.jl | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/tutorials/dc_motor_pi.md b/docs/src/tutorials/dc_motor_pi.md index 864dfa93e..b4e84f38b 100644 --- a/docs/src/tutorials/dc_motor_pi.md +++ b/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, [sys.L1.i => 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", diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md index 779421bad..7aa242f07 100644 --- a/docs/src/tutorials/input_component.md +++ b/docs/src/tutorials/input_component.md @@ -12,7 +12,7 @@ This tutorial demonstrate each case and explain the pros and cons of each. ## `Interpolation` Block The `ModelingToolkitStandardLibrary.Blocks.Interpolation` component is easy to use and is performant. -It is simlar to using callable paramterers, but it provides a block interface with `RealInput` and `RealOutput` connectors. +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`. ```@docs @@ -150,7 +150,7 @@ plot(sol2) ``` !!! note - Note that when changing the data, the length of the new data must be the same as the lenght of the original data. + 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 diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index 4c1797265..017da8d5e 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -554,7 +554,7 @@ end end @testset "BSplineInterpolation" begin - @named i = ParametrizedInterpolationBlock( + @named i = ParametrizedInterpolation( BSplineInterpolation, u, x, 3, :Uniform, :Uniform) eqs = [i.input.u ~ t, D(y) ~ i.output.u] From 96040a192767ed8711bbc91467e5916a156d779f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Sun, 13 Oct 2024 01:17:08 +0300 Subject: [PATCH 28/44] build: bump MTK and Symbolics --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 0110426cf..94d402ea4 100644 --- a/Project.toml +++ b/Project.toml @@ -20,7 +20,7 @@ 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" @@ -28,7 +28,7 @@ PreallocationTools = "0.4.23" SafeTestsets = "0.1" SciMLStructures = "1.4.2" SymbolicIndexingInterface = "0.3.28" -Symbolics = "5.36, 6" +Symbolics = "6.14" Test = "1" julia = "1.10" From dc710175c4fb0d5a5ddf229786d5a51fa34b2eb6 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 16 Oct 2024 11:23:12 -0400 Subject: [PATCH 29/44] Update src/Blocks/sources.jl Co-authored-by: Aayush Sabharwal --- src/Blocks/sources.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 8659428e9..7bdd04230 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -746,7 +746,7 @@ these would be any of [the available interpolations](https://github.com/SciML/Da such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. - `u`: the data used for interpolation. For `DataInterpolations` this will be an `AbstractVector` - `x`: the values that each data points correspond to, usually the times corresponding to each value in `u`. - - `args`: any other arguments beeded to build the interpolation + - `args`: any other arguments needed to build the interpolation # Keyword arguments: - `name`: the name of the component From 1bb50dde4f789914240c698db2952c4e8425938f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 16 Oct 2024 11:23:19 -0400 Subject: [PATCH 30/44] Update src/Blocks/sources.jl Co-authored-by: Aayush Sabharwal --- src/Blocks/sources.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 7bdd04230..cf9928243 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -751,7 +751,7 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. - `name`: the name of the component # Parameters: - - `interpolator`: the symbolic + - `interpolator`: the symbolic representation of the interpolation object, callable as `interpolator(t)` # Connectors: - `input`: a [`RealInput`](@ref) connector corresponding to the independent variable From 46fdf8257f25a4d610536ae3e4a78b0f093e760f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 16 Oct 2024 11:25:07 -0400 Subject: [PATCH 31/44] Update src/Blocks/sources.jl Co-authored-by: Aayush Sabharwal --- src/Blocks/sources.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index cf9928243..de2e2c09b 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -805,11 +805,9 @@ function (f::CachedInterpolation{T})(u, x, args) where {T} interp = @inbounds if (u, x) ≠ (get_tmp(prev_u, u), get_tmp(prev_x, x)) get_tmp(prev_u, u) .= u get_tmp(prev_x, x) .= x - # @info "cache miss" cache.bufs[(u, x)] = interpolation_type( get_tmp(prev_u, u), get_tmp(prev_x, x), args...) else - # @info "cache hit" cache[(u, x)] end From 559014675f04c1994bb310bf9d70c7d2ab8abe33 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 18 Oct 2024 11:18:12 -0400 Subject: [PATCH 32/44] Update src/Blocks/sources.jl Co-authored-by: Aayush Sabharwal --- src/Blocks/sources.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index de2e2c09b..585fa5440 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -814,9 +814,6 @@ function (f::CachedInterpolation{T})(u, x, args) where {T} return interp end -# function Symbolics.derivative(::typeof(CachedInterpolation), args::NTuple{2, Any}, ::Val{2}) -# Symbolics.derivative(args[1], (args[2],), Val(1)) -# end Base.nameof(::CachedInterpolation) = :CachedInterpolation From 4054d641e6387afd936292ede1945c5057c4ce77 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 20 Oct 2024 03:35:30 -0400 Subject: [PATCH 33/44] Update src/Blocks/sources.jl --- src/Blocks/sources.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 585fa5440..dcdfa0f5e 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -754,7 +754,7 @@ such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`. - `interpolator`: the symbolic representation of the interpolation object, callable as `interpolator(t)` # Connectors: - - `input`: a [`RealInput`](@ref) connector corresponding to the independent variable + - `input`: a [`RealInput`](@ref) connector corresponding to the input variable - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value """ function Interpolation(interp_type, u, x, args...; name) From 762eb744534bfe5808e9ad0f783d7424cd1f934d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Sun, 20 Oct 2024 10:47:54 +0300 Subject: [PATCH 34/44] add compats --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 6ca7cd6c4..1bca6cf3e 100644 --- a/Project.toml +++ b/Project.toml @@ -17,9 +17,11 @@ Aqua = "0.8" ChainRulesCore = "1.24" ControlSystemsBase = "1.4" DataInterpolations = "6" +ForwardDiff = "0.10" DiffEqBase = "6.152" IfElse = "0.1" LinearAlgebra = "1.10" +Optimization = "4" ModelingToolkit = "9.46.1" OrdinaryDiffEq = "6.87" OrdinaryDiffEqDefault = "1.1" From 5818f52355e3b075153444040d704bcad8a16b50 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Sun, 20 Oct 2024 10:52:07 +0300 Subject: [PATCH 35/44] add docstring --- src/Blocks/sources.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index dcdfa0f5e..31a2a291f 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -773,6 +773,11 @@ function Interpolation(itp; name) eqs, t, [], [interpolator]; name, systems = [input, output]) end +""" + CachedInterpolation + +This callable struct caches the calls to an interpolation object via PreallocationTools. +""" struct CachedInterpolation{T, I, U, X, C} interpolation_type::I prev_u::U @@ -814,7 +819,6 @@ function (f::CachedInterpolation{T})(u, x, args) where {T} return interp end - Base.nameof(::CachedInterpolation) = :CachedInterpolation @register_symbolic (f::CachedInterpolation)(u::AbstractArray, x::AbstractArray, args::Tuple) From c91104d936613c9856f4d82153fea4ea174e2722 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 20 Oct 2024 03:52:50 -0400 Subject: [PATCH 36/44] Update input_component.md --- docs/src/tutorials/input_component.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md index 7aa242f07..9e2006c40 100644 --- a/docs/src/tutorials/input_component.md +++ b/docs/src/tutorials/input_component.md @@ -1,4 +1,4 @@ -# Running Models with Discrete Data +# Building Models with Discrete Data, Interpolations, and Lookup Tables There are 4 ways to include data as part of a model. From 4de943f831115bd0d616dd0b97c991745046ea54 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Sun, 20 Oct 2024 11:30:42 +0300 Subject: [PATCH 37/44] try to fix docs --- docs/src/tutorials/input_component.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md index 9e2006c40..c451949fc 100644 --- a/docs/src/tutorials/input_component.md +++ b/docs/src/tutorials/input_component.md @@ -212,7 +212,7 @@ To resolve the issues presented above, the `ModelingToolkitStandardLibrary.Block ```@example sampled_data_component function System(; name) - @named 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 From 2ee5ecdf274752ddd732ea787e2f9f7928615630 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Sun, 20 Oct 2024 11:59:27 +0300 Subject: [PATCH 38/44] fix sampled data docs --- docs/src/tutorials/input_component.md | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md index c451949fc..2ea926a9c 100644 --- a/docs/src/tutorials/input_component.md +++ b/docs/src/tutorials/input_component.md @@ -157,8 +157,15 @@ plot(sol2) 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. ```@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) @@ -211,6 +218,11 @@ 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()`. ```@example sampled_data_component +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D +using ModelingToolkitStandardLibrary.Blocks +using OrdinaryDiffEq + function System(; name) src = SampledData(Float64, name=:src) @@ -228,6 +240,12 @@ end @named system = System() sys = structural_simplify(system, split=false) s = complete(system) + +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]); tofloat = false, use_union=true) defs = ModelingToolkit.defaults(sys) From 95c24d8bec1f44ae6f761c80e51674b3e4ec200f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Sun, 20 Oct 2024 12:39:15 +0300 Subject: [PATCH 39/44] fix typos & avoid running SampledData for now --- docs/src/tutorials/input_component.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md index 2ea926a9c..4cbd03759 100644 --- a/docs/src/tutorials/input_component.md +++ b/docs/src/tutorials/input_component.md @@ -44,7 +44,7 @@ function MassSpringDamper(; name) D(x) ~ dx D(dx) ~ ddx] - ODESystem(eqs, t, [], []; name, systems = [input]) + ODESystem(eqs, t; name, systems = [input]) end function MassSpringDamperSystem(data, time; name) @@ -122,7 +122,7 @@ function MassSpringDamperSystem(data, time; name) connect(src.input, clk.output) ] - ODESystem(eqs, t, [], []; name, systems = [src, clk, model]) + ODESystem(eqs, t; name, systems = [src, clk, model]) end function generate_data() @@ -217,7 +217,7 @@ 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()`. -```@example sampled_data_component +```julia using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D using ModelingToolkitStandardLibrary.Blocks @@ -246,14 +246,14 @@ time = 0:dt:0.1 data1 = sin.(2 * pi * time * 100) data2 = cos.(2 * pi * time * 50) -prob = ODEProblem(sys, [], (0, time[end]); tofloat = false, use_union=true) +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) From 6f549b535290400764acb7ec82b047c6a40750e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Tue, 22 Oct 2024 19:22:09 +0300 Subject: [PATCH 40/44] test: add test making sure that the initialization works this is the example in the docs --- test/Blocks/sources.jl | 48 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index fb138dde7..8437f79df 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -566,4 +566,52 @@ end @test SciMLBase.successful_retcode(sol) end + + @testset "Initialization" begin + 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 ~ input.u + ddx * 10 ~ k * x + d * dx + f + D(x) ~ dx + D(dx) ~ ddx] + + ODESystem(eqs, t, vars, pars; name, systems = [input]) + end + + 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)] + + 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) + + @test SciMLBase.successful_retcode(sol) + + prob2 = remake(prob, p = [sys.src.data => ones(length(df.data))]) + sol2 = solve(prob2) + + @test SciMLBase.successful_retcode(sol2) + end end From 3c8f2da5bb2027dd058be8e148ade75a66585d8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Tue, 22 Oct 2024 19:22:29 +0300 Subject: [PATCH 41/44] docs: fix initialization --- docs/src/tutorials/input_component.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md index 4cbd03759..f48dd96b1 100644 --- a/docs/src/tutorials/input_component.md +++ b/docs/src/tutorials/input_component.md @@ -100,7 +100,7 @@ using Plots function MassSpringDamper(; name) @named input = RealInput() - vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0 + vars = @variables f(t) x(t)=0 dx(t) [guess=0] ddx(t) pars = @parameters m=10 k=1000 d=1 eqs = [ @@ -144,7 +144,7 @@ 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(data))]) +prob2 = remake(prob, p = [sys.src.data => ones(length(df.data))]) sol2 = solve(prob2) plot(sol2) ``` From e2f285920f64d9b271407ef389844a9788b4679c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 24 Oct 2024 13:32:06 +0300 Subject: [PATCH 42/44] build: bump MTK --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1bca6cf3e..f353c9815 100644 --- a/Project.toml +++ b/Project.toml @@ -22,7 +22,7 @@ DiffEqBase = "6.152" IfElse = "0.1" LinearAlgebra = "1.10" Optimization = "4" -ModelingToolkit = "9.46.1" +ModelingToolkit = "9.47" OrdinaryDiffEq = "6.87" OrdinaryDiffEqDefault = "1.1" PreallocationTools = "0.4.23" From 8dc639a6262a5a13653f0c3828cfdeb2aecdb587 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 24 Oct 2024 14:08:57 +0300 Subject: [PATCH 43/44] test: add DataFrames to test deps --- Project.toml | 3 ++- test/Blocks/sources.jl | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f353c9815..4ac3ebd79 100644 --- a/Project.toml +++ b/Project.toml @@ -36,6 +36,7 @@ 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" @@ -48,4 +49,4 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "LinearAlgebra", "OrdinaryDiffEqDefault", "OrdinaryDiffEq", "Optimization", "SafeTestsets", "Test", "ControlSystemsBase", "DataInterpolations", "SciMLStructures", "SymbolicIndexingInterface", "ForwardDiff"] +test = ["Aqua", "LinearAlgebra", "OrdinaryDiffEqDefault", "OrdinaryDiffEq", "Optimization", "SafeTestsets", "Test", "ControlSystemsBase", "DataFrames", "DataInterpolations", "SciMLStructures", "SymbolicIndexingInterface", "ForwardDiff"] diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index 8437f79df..8598278b6 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -6,6 +6,7 @@ using ModelingToolkitStandardLibrary.Blocks: smooth_sin, smooth_cos, smooth_damp smooth_triangular, triangular, square using OrdinaryDiffEq: ReturnCode.Success using DataInterpolations +using DataFrames using SymbolicIndexingInterface using SciMLStructures: SciMLStructures, Tunable using Optimization From cfee55c3ddd5f64f186ae0f2388c4143682c47b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Thu, 24 Oct 2024 14:17:02 +0300 Subject: [PATCH 44/44] build: add DataFrames compat --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 4ac3ebd79..4bdc82f01 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Aqua = "0.8" ChainRulesCore = "1.24" ControlSystemsBase = "1.4" +DataFrames = "1.7" DataInterpolations = "6" ForwardDiff = "0.10" DiffEqBase = "6.152"