Merge pull request #968 from JuliaControl/pade2

add Padé approximation with different numerator and denominator degree

Author: baggepinnen

lib/ControlSystemsBase/src/delay_systems.jl

@@ -140,16 +140,32 @@ function pade(τ::Real, N::Int)
     return tf(_linscale(Q, -τ), _linscale(Q, τ)) # return Q(-τs)/Q(τs)
 end
 
+# Pade approximation with different degree in numerator and denominator
+"""
+    pade(τ::Real, N_num::Int, N_den::Int)
+
+Compute the Padé approximation of a time-delay of length `τ` with `N_num` and `N_den` degrees in the numerator and denominator, respectively.
+"""
+function pade(τ::Real, m::Int, n::Int)
+    p = [(-1)^i * binomial(m, i) * factorial(m+n-i) / factorial(m+n)
+        for i in 0:m]  |> Polynomials.Polynomial
+
+    q = [binomial(n, i) * factorial(m+n-i) / factorial(m+n)
+        for i in 0:n]  |> Polynomials.Polynomial
+
+    return tf(_linscale(p, τ), _linscale(q, τ))
+end
+
 
 """
     pade(G::DelayLtiSystem, N)
 
 Approximate all time-delays in `G` by Padé approximations of degree `N`.
 """
-function pade(G::DelayLtiSystem, N)
+function pade(G::DelayLtiSystem, N, args...)
     ny, nu = size(G)
     nTau = length(G.Tau)
-    X = append(ss(pade(τ,N)) for τ in G.Tau) # Perhaps append should be renamed blockdiag
+    X = append(ss(pade(τ,N,args...)) for τ in G.Tau) # Perhaps append should be renamed blockdiag
     return lft(G.P.P, X)
 end
 

lib/ControlSystemsBase/test/test_delayed_systems.jl

@@ -229,6 +229,20 @@ P_wb = DemoSystems.woodberry()
 
 @test freqresp(pade(feedback(eye_(2), P_wb), 3), Ω) ≈ freqresp(feedback(eye_(2), P_wb), Ω) atol=1e-4
 
+## Test pade(tau, m, n)
+t = 1.5
+P15 = pade(t, 1, 5)
+dv = denvec(P15)[]
+dv1 = dv[1]/t^5 # Rescale to compare with https://www2.humusoft.cz/www/papers/tcp09/035_hanta.pdf
+@test numvec(P15)[] ./ dv1 ≈ [-120*t, 720]
+@test dv ./ dv1 ≈ [t^5, 10*t^4, 60*t^3, 240*t^2, 600*t, 720]
+
+
+P15 = pade(t, 2, 5)
+# Test vectors computed using RobustPade.robustpade(s->exp(-t*s), 2, 5)
+@test numvec(P15)[] ≈ [0.05357142857142857, -0.42857142857142855, 1.0]
+@test denvec(P15)[] ≈ [0.0030133928571428573, 0.03013392857142857, 0.1607142857142857, 0.5357142857142857, 1.0714285714285714, 1.0]
+
 # test thiran
 for Ts = [1, 1.1]
     z = tf('z', Ts)