Merge pull request #905 from JuliaControl/thiran

add Thiran approximation of fractional-delay systems

Author: baggepinnen

docs/make.jl

@@ -22,7 +22,7 @@ makedocs(modules=[ControlSystems, ControlSystemsBase],
     format = Documenter.HTML(prettyurls = haskey(ENV, "CI")),
     sitename="ControlSystems.jl",
     pagesonly = true,
-    draft = false,
+    draft = true,
     pages=[
         "Home" => "index.md",
         "Introductory guide" => Any[
@@ -39,7 +39,7 @@ makedocs(modules=[ControlSystems, ControlSystemsBase],
             "Properties of delay systems" => "examples/delay_systems.md",
             "Automatic differentiation" => "examples/automatic_differentiation.md",
             "Tune a controller using experimental data" => "examples/tuning_from_data.md",
-            "Analysis of hybrid continuous/discrete systems" => "examples/zoh.md",
+            "Analysis of sampled-data (continuous/discrete) systems" => "examples/zoh.md",
         ],
         "Functions" => Any[
             "Constructors" => "lib/constructors.md",

docs/src/examples/delay_systems.md

@@ -32,4 +32,7 @@ For a more advanced example using time delays, see the [Smith predictor](@ref) t
 Time-delay systems are numerically challenging to simulate, if you run into problems, please open an issue with a reproducing example. The [`lsim`](@ref), [`step`](@ref) and [`impulse`](@ref) functions accept keyword arguments that are passed along to the ODE integrator, this can be used to both select integration method and to tweak the integrator options. The documentation for solving delay-differential equations is available [here](https://diffeq.sciml.ai/latest/types/dde_types/) and [here](https://diffeq.sciml.ai/latest/tutorials/dde_example/).
 
 ## Estimation of delay
-See the companion tutorial in ControlSystemIdentification.jl on [Delay estimation](file:///home/fredrikb/.julia/dev/ControlSystemIdentification/docs/build/examples/delayest.html). This tutorial covers the both the detection of the presence of a delay, and estimation of models for systems with delays.
+See the companion tutorial in ControlSystemIdentification.jl on [Delay estimation](file:///home/fredrikb/.julia/dev/ControlSystemIdentification/docs/build/examples/delayest.html). This tutorial covers the both the detection of the presence of a delay, and estimation of models for systems with delays.
+
+## Approximation and discretization of delays
+Delay systems may be approximated as rational functions by means of [Padé approximation](https://en.wikipedia.org/wiki/Pad%C3%A9_approximant) using the function [`pade`](@ref). Pure continuous-time delays can also be discretized using the function [`thiran`](@ref). Continuous-time models with internal delays can be discretized using [`c2d`](@ref), provided that the delay is an integer multiple of the sampling time (fractional delays are not yet supported by [`c2d`](@ref)).

docs/src/examples/zoh.md

@@ -16,6 +16,7 @@ Next, we discretize this system using the standard [`c2d`](@ref) function, which
 ```math
 Z(s) = \dfrac{1 - e^{-T_s s}}{s}
 ```
+[^CCS]: Åström, K. J., & Wittenmark, B. (1997). Computer-Controlled Systems: Theory and Design.
 
 ```@example zoh
 Ts = 1 # Sample interval
@@ -113,7 +114,7 @@ The discrete-time controller ``C(z)`` is a basic PI controller
 Ts = 1
 C = pid(0.01,10; Ts, Tf = 1/100, state_space=true)
 ```
-If we choose option 1. in this case, we get a discretized plant that has a very poor frequency-response fit to the original continuous-time plant, making frequency-domain analysis difficult
+If we choose option 1, we get a discretized plant that has a very poor frequency-response fit to the original continuous-time plant, making frequency-domain analysis difficult
 ```@example zoh
 Pd = c2d(P, Ts)
 Ld = Pd*C
@@ -146,13 +147,11 @@ plot(lsim(PS, disturbance, 0:0.22:3500), lab="Continuous disturbance response")
 plot!(lsim(PSd, disturbance, 3500), lab="Discrete disturbance response")
 hline!([abs(freqresp(PS, ω)[])], l=(:black, :dash), lab="Predicted freq. response amplitude", legend=:bottomleft, fmt=:png)
 ```
-The continuous-time analysis eventually settles at a periodic output that matches the amplitude predicted by the continuous-time frequency response. However, the discrete-time simulation yields, as expected, a very poor result. Noticeable in the simulation is the appearance of a beat frequency, which is expected due to the modulation introduced by sampling.[^CCS]
+The continuous-time analysis eventually settles at a periodic output that matches the amplitude predicted by the continuous-time frequency response. However, the discrete-time simulation yields, as expected, a very poor result. Noticeable in the simulation is the appearance of a beat frequency, which is expected due to the modulation introduced by sampling. [^CCS]
 
 ## Caveats
 - The exact output of the system that was translated from discrete to continuous time is not going to be accurate, so transient properties of the hybrid system cannot be accurately assessed using this kind of analysis. 
-- Interpretation of frequency-responses for sampled-data systems must be done with care. The modulation introduced by sampling implies the creating of additional frequencies in the output. For an input with frequency ``\omega``, the output will contain all frequencies ``\omega ± \omega_s k`` where ``\omega_s`` is the sampling frequency and ``k`` is an integer.[^CCS]
+- Interpretation of frequency-responses for sampled-data systems must be done with care. The modulation introduced by sampling implies the creating of additional frequencies in the output. For an input with frequency ``\omega``, the output will contain all frequencies ``\omega ± \omega_s k`` where ``\omega_s`` is the sampling frequency and ``k`` is an integer. [^CCS]
 
 ## References
 Learn more about sampled-data systems and zero-order hold sampling in chapter 7 of the reference below.
-[^CCS]: Åström, K. J., & Wittenmark, B. (1997). Computer-Controlled Systems: Theory and Design.
-

docs/src/man/creating_systems.md

@@ -123,6 +123,12 @@ Sample Time: 0.1 (seconds)
 Discrete-time transfer function model
 ```
 
+
+## Converting between continuous and discrete time
+A continuous-time system represents differential equations or a transfer function in the [Laplace domain](https://en.wikipedia.org/wiki/Laplace_transform), while a discrete-time system represents difference equations or a transfer function in the [Z-domain](https://en.wikipedia.org/wiki/Z-transform).
+
+The functions [`c2d`](@ref) and [`d2c`](@ref) implement sampling/discretization of continuous-time systems and the inverse mapping from discrete-time to continuous-time systems. 
+
 ## Delay Systems
 The constructor [`delay`](@ref) creates a pure delay, which may be connected to a system by multiplication:
 ```julia
@@ -138,7 +144,7 @@ L = 1.2 # Delay time
 tf(1, [1, 1]) * exp(-L*s)
 ```
 
-Padé approximations of delays can be created using [`pade`](@ref).
+Padé approximations of delays can be created using [`pade`](@ref). Models with delays can be discretized using [`c2d`](@ref), currently, only delays that are integer multiples of the sample time are supported. Pure fractional delays can be approximately discretized using the function [`thiran`](@ref).
 
 A tutorial on delay systems is available here:
 ```@raw html

lib/ControlSystemsBase/src/ControlSystemsBase.jl

@@ -104,6 +104,7 @@ export  LTISystem,
         # delay systems
         delay,
         pade,
+        thiran,
         nonlinearity,
         # demo systems
         ssrand,

lib/ControlSystemsBase/src/delay_systems.jl

@@ -38,15 +38,16 @@ end
 
 """
     delayd_ss(τ, Ts)
-Discrete-time statespace realization of a delay τ sampled with period Ts,
-i.e. of z^-N where N = τ / Ts.
 
-τ must be a multiple of Ts.
+Discrete-time statespace realization of a delay ``τ`` sampled with period ``T_s``,
+i.e. of ``z^{-N}`` where ``N = τ / T_s.``
+
+``τ`` must be a multiple of ``T_s``. See [`thiran`](@ref) for approximate discretization of fractional delays.
 """
 function delayd_ss(τ::Number, Ts::Number)
     n = round(Int, τ / Ts)
     if !(τ - n*Ts ≈ 0)
-        error("The delay τ must be a multiple of the sample time Ts")
+        error("The delay τ must be a multiple of the sample time Ts, use the function `thiran` to approximately discretize fractional delays.")
     end
     ss(diagm(1 => ones(n-1)), [zeros(n-1,1); 1], [1 zeros(1,n-1)], 0, Ts)
 end
@@ -128,6 +129,8 @@ const PADE_Q_COEFFS =  [[2, 1],
     pade(τ::Real, N::Int)
 
 Compute the `N`th order Padé approximation of a time-delay of length `τ`.
+
+See also [`thiran`](@ref) for discretization of delays.
 """
 function pade(τ::Real, N::Int)
     if !(1 <= N <= 10); error("Order of Padé approximation must be between 1 and 10. Got $N."); end
@@ -149,3 +152,25 @@ function pade(G::DelayLtiSystem, N)
     X = append(ss(pade(τ,N)) for τ in G.Tau) # Perhaps append should be renamed blockdiag
     return lft(G.P.P, X)
 end
+
+"""
+    thiran(τ::Real, Ts)
+
+Discretize a potentially fractional delay ``τ`` as a Thiran all-pass filter with sample time `Ts`. 
+
+The Thiran all-pass filter gives an a maximally flat group delay.
+
+If ``τ`` is an integer multiple of ``Ts``, the Thiran all-pass filter reduces to ``z^{-τ/Ts}``.
+
+Ref: T. I. Laakso, V. Valimaki, M. Karjalainen and U. K. Laine, "Splitting the unit delay [FIR/all pass filters design]," in IEEE Signal Processing Magazine, vol. 13, no. 1, 1996.
+"""
+function thiran(τ::Real, Ts)
+    D = τ/Ts
+    N = ceil(Int, D)
+    a = ones(N+1)
+    for k = 1:N
+        P = prod((D-N+n) / (D-N+n+k) for n in 0:N)
+        a[k+1] = (-1)^k * binomial(N, k) * P # Eq 86 in reference
+    end
+    tf(reverse(a), a, Ts)
+end

lib/ControlSystemsBase/test/test_delayed_systems.jl

@@ -229,6 +229,14 @@ P_wb = DemoSystems.woodberry()
 
 @test freqresp(pade(feedback(eye_(2), P_wb), 3), Ω) ≈ freqresp(feedback(eye_(2), P_wb), Ω) atol=1e-4
 
+# test thiran
+for Ts = [1, 1.1]
+    z = tf('z', Ts)
+    @test thiran(2Ts, Ts) == 1/z^2
+end
+
+@test thiran(pi, 1) ≈ tf([-0.00031815668236122736, 0.0042438423976339556, -0.03424682772398137, 0.8290601401044773, 1.0], [1.0, 0.8290601401044773, -0.03424682772398137, 0.0042438423976339556, -0.00031815668236122736], 1)
+
 
 # test automatic frequency selection
 mag, phase, w = bode(DemoSystems.lag()*delay(1))