add transfer function block

Author: baggepinnen

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

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