From a75c48e51d7ce85fccbab41a9a28ae2fbc0a8161 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 7 Feb 2025 07:38:12 +0100
Subject: [PATCH] =?UTF-8?q?add=20Pad=C3=A9=20approximation=20with=20differ?=
 =?UTF-8?q?ent=20numerator=20and=20denomenator=20degree?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 lib/ControlSystemsBase/src/delay_systems.jl   | 20 +++++++++++++++++--
 .../test/test_delayed_systems.jl              | 14 +++++++++++++
 2 files changed, 32 insertions(+), 2 deletions(-)

diff --git a/lib/ControlSystemsBase/src/delay_systems.jl b/lib/ControlSystemsBase/src/delay_systems.jl
index 3bb875bda..d76defbaa 100644
--- a/lib/ControlSystemsBase/src/delay_systems.jl
+++ b/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
 
diff --git a/lib/ControlSystemsBase/test/test_delayed_systems.jl b/lib/ControlSystemsBase/test/test_delayed_systems.jl
index 2213fa8fe..47d836841 100644
--- a/lib/ControlSystemsBase/test/test_delayed_systems.jl
+++ b/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)