add transfer function block

Project.toml

@@ -12,7 +12,14 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 
+[weakdeps]
+ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
+
+[extensions]
+ModelingToolkitSampledDataControlSystemsBaseExt = ["ControlSystemsBase"]
+
 [compat]
+ControlSystemsBase = "1"
 DiffEqCallbacks = "~3.8"
 JuliaSimCompiler = "0.1.19"
 ModelingToolkit = "9"

docs/src/blocks.md

@@ -9,6 +9,7 @@
 - [`Difference`](@ref)
 - [`DiscreteDerivative`](@ref)
 - [`DiscreteIntegrator`](@ref)
+- [`DiscreteTransferFunction`](@ref)
 - [`Sampler`](@ref)
 - [`ZeroOrderHold`](@ref)
 
@@ -24,6 +25,7 @@
 
 ## Discrete-time filters
 - [`ExponentialFilter`](@ref)
+- [`DiscreteTransferFunction`](@ref)
 
 
 ## Docstrings

docs/src/tutorials/noise.md

@@ -88,6 +88,48 @@ plot(figy, figu, plot_title = "DC Motor with Discrete-time Speed Controller")
 You may, e.g.
 - Use [`ExponentialFilter`](@ref) to add exponential filtering using `y(k) ~ (1-a)y(k-1) + a*u(k)`, where `a` is the filter coefficient and `u` is the signal to be filtered.
 - Add moving average filtering using `y(k) ~ 1/N sum(i->u(k-i), i=0:N-1)`, where `N` is the number of samples to average over.
+- Construct a filter transfer function using [`DiscreteTransferFunction`](@ref). If ControlSystemsBase is loaded, `ControlSystemsBase.TransferFunction` objects can be used to construct a [`DiscreteTransferFunction`](@ref). For advanced filter design, see [DSP.jl: Filter design](https://docs.juliadsp.org/stable/filters/#Filter-design), example below.
+
+### Filter design using DSP.jl
+Here, we design a fourth-order Butterworth low-pass filter with a cutoff frequency of 100 Hz and apply it to a noisy signal. The filter is designed using DSP.jl and converted to a transfer function from ControlSystemsBase.jl. The transfer function can be passed to the [`DiscreteTransferFunction`](@ref) block. We also design an [`ExponentialFilter`](@ref) for comparison.
+```@example NOISE
+using ControlSystemsBase: tf
+using ControlSystemsBase.DSP: digitalfilter, Lowpass, Butterworth
+fs = 1000
+cutoff = 100
+z = ShiftIndex(Clock(1/fs))
+F_dsp = digitalfilter(Lowpass(cutoff; fs), Butterworth(4))
+F_tf = tf(F_dsp)
+@mtkmodel FilteredNoise begin
+ @components begin
+ sine = Sine(amplitude=1, frequency=10)
+ noise = NormalNoise(sigma = 0.1, additive = true)
+ filter1 = ExponentialFilter(a = 0.2)
+ filter2 = DiscreteTransferFunction(F_tf; z)
+ end
+ @variables begin
+ x(t) = 0 # Dummy variable to work around a bug for models without continuous-time state
+ end
+ @equations begin
+ connect(sine.output, noise.input)
+ connect(noise.output, filter1.input, filter2.input)
+ D(x) ~ 0 # Dummy equation
+ end
+end
+@named m = FilteredNoise()
+m = complete(m)
+ssys = structural_simplify(IRSystem(m))
+prob = ODEProblem(ssys, [
+ m.filter1.y(z-1) => 0;
+ [m.filter2.y(z-k) => 0 for k = 1:4];
+ [m.filter2.u(z-k) => 0 for k = 1:4];
+ m.noise.y(z-1) => 0;
+], (0.0, 0.1))
+sol = solve(prob, Tsit5())
+plot(sol, idxs=m.noise.y, label="Noisy signal", c=1)
+plot!(sol, idxs=m.filter1.y, label="Exponential filtered signal", c=2)
+plot!(sol, idxs=m.filter2.y, label="Butterworth filtered signal", c=3)
+```
 
 ## Colored noise
 Colored noise can be achieved by filtering white noise through a filter with the desired spectrum. No components are available for this yet.

ext/ModelingToolkitSampledDataControlSystemsBaseExt.jl

@@ -0,0 +1,22 @@
+module ModelingToolkitSampledDataControlSystemsBaseExt
+using ModelingToolkitSampledData
+using ControlSystemsBase: TransferFunction, issiso, numvec, denvec, Discrete
+
+
+"""
+ ModelingToolkitSampledData.DiscreteTransferFunction(G::TransferFunction{<:Discrete}; kwargs...)
+
+Create a DiscreteTransferFunction from a `ControlSystems.TransferFunction`.
+
+The sample time of `G` is used to create a periodic `Clock` object with which the transfer funciton is associated. If this is not desired, pass a custom `ShiftIndex` via the `z` keyword argument. To let the transfer function inherit the sample time of of the surrounding context (clock inference), pass and empty `z = ShiftIndex()`.
+"""
+function ModelingToolkitSampledData.DiscreteTransferFunction(G::TransferFunction{<:Discrete}; z = nothing, kwargs...)
+ issiso(G) || throw(ArgumentError("Only SISO systems are supported"))
+ b,a = numvec(G)[], denvec(G)[]
+ if z === nothing
+ z = ShiftIndex(Clock(G.Ts))
+ end
+ return DiscreteTransferFunction(b, a; z, kwargs...)
+end
+
+end

src/ModelingToolkitSampledData.jl

@@ -9,6 +9,7 @@ export DiscreteIntegrator, DiscreteDerivative, Delay, Difference, ZeroOrderHold,
  DiscretePIDParallel, DiscretePIDStandard, DiscreteStateSpace,
  DiscreteTransferFunction, NormalNoise, UniformNoise, Quantization,
  ExponentialFilter
+export DiscreteTransferFunction
 export DiscreteOnOffController
 include("discrete_blocks.jl")
 

src/discrete_blocks.jl

@@ -718,25 +718,57 @@ With the coeffiencents specified in decreasing orders of ``z``, i.e., ``b = [b_{
 - `output`: Output signal
 
 # Extended help:
-This component supports SISO systems only. To simulate MIMO transfer functions, use [ControlSystemsBase.jl](https://juliacontrol.github.io/ControlSystems.jl/stable/man/creating_systems/) to convert the transfer function to a statespace system, optionally compute a minimal realization using [`minreal`](https://juliacontrol.github.io/ControlSystems.jl/stable/lib/constructors/#ControlSystemsBase.minreal), and then use a [`DiscreteStateSpace`](@ref) component instead.
+This component supports SISO systems only. To simulate MIMO transfer functions, use [ControlSystemsBase.jl](https://juliacontrol.github.io/ControlSystems.jl/stable/man/creating_systems/) to convert the transfer function to a statespace system.
 
-See also [ControlSystemsMTK.jl](https://juliacontrol.github.io/ControlSystemsMTK.jl/stable/) for an interface between [ControlSystems.jl](https://juliacontrol.github.io/ControlSystems.jl/stable/) and ModelingToolkit.jl for advanced manipulation of transfer functions and linear statespace systems. For linearization, see [`linearize`](@ref) and [Linear Analysis](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/).
+See also [ControlSystemsMTK.jl](https://juliacontrol.github.io/ControlSystemsMTK.jl/stable/) for an interface between [ControlSystems.jl](https://juliacontrol.github.io/ControlSystems.jl/stable/) and ModelingToolkit.jl for advanced manipulation of transfer functions and linear statespace systems.
 """
+@component function DiscreteTransferFunction(; a = [1], b = [1],
+ verbose = true,
+ z = ShiftIndex(),
+ name,
+ )
+ na = length(a)
+ nb = length(b)
+ verbose && nb > na && @warn "DiscreteTransferFunction: Numerator degree is larger than denominator degree, this is not a proper transfer function. Simulation of a model including this transfer funciton will require at least $(nb-na) samples additional delay in series. Silence this warning with verbose=false"
+ @parameters begin
+ (b[1:nb] = b), [description = "Numerator coefficients of transfer function (e.g., z - 1 is specified as [1,-1])"]
+ (a[1:na] = a), [description = "Denominator coefficients of transfer function (e.g., z^2 - 0.78z + 0.37 is specified as [1, -0.78, 0.37])"]
+ end
+ a = collect(a)
+ b = collect(b)
+ systems = @named begin
+ input = RealInput()
+ output = RealOutput()
+ end
+ @variables begin
+ (u(t) = 0.0), [description = "Input flowing through connector `input`"]
+ (y(t) = 0.0), [description = "Output flowing through connector `output`"]
+ end
+ equations = [
+ sum(y(z-k+1) * a[k] for k in 1:na) ~ sum(u(z-k+1) * b[k] for k in 1:nb)
+ input.u ~ u
+ output.u ~ y
+ ]
+ return compose(ODESystem(equations, t, [y,u], [b; a]; name), systems)
+end
+
+DiscreteTransferFunction(b, a; kwargs...) = DiscreteTransferFunction(; b, a, kwargs...)
+
 # @mtkmodel DiscreteTransferFunction begin
 # @parameters begin
-# b = [1], [description = "Numerator coefficients of transfer function (e.g., z - 1 is specified as [1,-1])"]
-# a = [1], [description = "Denominator coefficients of transfer function (e.g., z^2 - 0.78z + 0.37 is specified as [1, -0.78, 0.37])"]
+# b[1:1] = [1], [description = "Numerator coefficients of transfer function (e.g., z - 1 is specified as [1,-1])"]
+# a[1:1] = [1], [description = "Denominator coefficients of transfer function (e.g., z^2 - 0.78z + 0.37 is specified as [1, -0.78, 0.37])"]
 # end
 # @structural_parameters begin
 # verbose = true
-# Ts = SampleTime()
-# z = ShiftIndex()
+#  Ts = SampleTime()
+#  z = ShiftIndex()
 # end
 # begin
-# na = length(a)
-# nb = length(b)
+# @show na = length(a)
+# @show nb = length(b)
+# @show a
 # verbose && nb > na && @warn "DiscreteTransferFunction: Numerator degree is larger than denominator degree, this is not a proper transfer function. Simulation of a model including this transfer funciton will require at least $(nb-na) samples additional delay in series. Silence this warning with verbose=false"
-# Ts = SampleTime()
 # end
 # @components begin
 # input = RealInput()
@@ -753,7 +785,7 @@ See also [ControlSystemsMTK.jl](https://juliacontrol.github.io/ControlSystemsMTK
 # end
 # end
 
-# DiscreteTransferFunction(b, a; kwargs...) = DiscreteTransferFunction(; b, a, kwargs...)
+# 
 
 ##
 

test/test_discrete_blocks.jl

@@ -514,4 +514,59 @@ end
  @test sol(0.999, idxs=m.filter.y) == 0
  @test sol(1.1, idxs=m.filter.y) > 0
 
-end
+end
+
+@testset "DiscreteTransferFunction" begin
+ @info "Testing DiscreteTransferFunction"
+ z = ShiftIndex(Clock(0.1))
+ @mtkmodel DiscreteTransferFunctionModel begin
+ @components begin
+ input = Step(start_time=1, smooth=false)
+ filter = DiscreteTransferFunction(; a = [1, -0.9048374180359594], b = [0.09516258196404037], z)
+ end
+ @variables begin
+ x(t) = 0 # Dummy variable to workaround JSCompiler bug
+ end
+ @equations begin
+ connect(input.output, filter.input)
+ D(x) ~ 0
+ end
+ end
+ @named m = DiscreteTransferFunctionModel()
+ m = complete(m)
+ ssys = structural_simplify(IRSystem(m))
+ prob = ODEProblem(ssys, [m.filter.y(z-1) => 0], (0.0, 10.0))
+ sol = solve(prob, Tsit5(), dtmax=0.1)
+ @test sol(10, idxs=m.filter.y) ≈ 1 atol=0.001
+ @test sol(0.999, idxs=m.filter.y) == 0
+ @test sol(1.1, idxs=m.filter.y) > 0
+ @test sol(1:10, idxs=m.filter.y).u ≈ [0.0, 0.6321205588285568, 0.8646647167633857, 0.950212931632134, 0.9816843611112637, 0.9932620530009122, 0.9975212478233313, 0.9990881180344431, 0.999664537372095, 0.9998765901959109]
+
+ import ControlSystemsBase as CS
+ Gc = CS.tf(1, [1, 1])
+ G = CS.c2d(Gc, 0.1)
+
+ @mtkmodel DiscreteTransferFunctionModel begin
+ @components begin
+ input = Step(start_time=1, smooth=false)
+ filter = DiscreteTransferFunction(G)
+ end
+ @variables begin
+ x(t) = 0 # Dummy variable to workaround JSCompiler bug
+ end
+ @equations begin
+ connect(input.output, filter.input)
+ D(x) ~ 0
+ end
+ end
+ @named m = DiscreteTransferFunctionModel()
+ m = complete(m)
+ ssys = structural_simplify(IRSystem(m))
+ prob = ODEProblem(ssys, [m.filter.y(z-1) => 0], (0.0, 10.0))
+ sol = solve(prob, Tsit5(), dtmax=0.1)
+ @test sol(10, idxs=m.filter.y) ≈ 1 atol=0.001
+ @test sol(0.999, idxs=m.filter.y) == 0
+ @test sol(1.1, idxs=m.filter.y) > 0
+ @test sol(1:10, idxs=m.filter.y).u ≈ [0.0, 0.6321205588285568, 0.8646647167633857, 0.950212931632134, 0.9816843611112637, 0.9932620530009122, 0.9975212478233313, 0.9990881180344431, 0.999664537372095, 0.9998765901959109]
+
+end