Skip to content

Commit 0b26e06

Browse files
slivingstonslivingston
committed
Merge pull request #76 from roryyorke/rory-pade-spec-num
Feature: allow differing numerator and denominator degree in Padé approximations
2 parents 49045ca + 2d79f4e commit 0b26e06

File tree

2 files changed

+131
-12
lines changed

2 files changed

+131
-12
lines changed

‎control/delay.py

Lines changed: 39 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
# -*-coding: utf-8-*-
12
#! TODO: add module docstring
23
# delay.py - functions involving time delays
34
#
@@ -41,12 +42,11 @@
4142
#
4243
# $Id$
4344

44-
# Python 3 compatability (needs to go here)
45-
from __future__ import print_function
45+
from __future__ import division
4646

4747
__all__ = ['pade']
4848

49-
def pade(T, n=1):
49+
def pade(T, n=1, numdeg=None):
5050
"""
5151
Create a linear system that approximates a delay.
5252
@@ -56,8 +56,12 @@ def pade(T, n=1):
5656
----------
5757
T : number
5858
time delay
59-
n : integer
60-
order of approximation
59+
n : positive integer
60+
degree of denominator of approximation
61+
numdeg: integer, or None (the default)
62+
If None, numerator degree equals denominator degree
63+
If >= 0, specifies degree of numerator
64+
If < 0, numerator degree is n+numdeg
6165
6266
Returns
6367
-------
@@ -66,22 +70,45 @@ def pade(T, n=1):
6670
6771
Notes
6872
-----
69-
Based on an algorithm in Golub and van Loan, "Matrix Computation" 3rd.
70-
Ed. pp. 572-574.
73+
Based on:
74+
1. Algorithm 11.3.1 in Golub and van Loan, "Matrix Computation" 3rd.
75+
Ed. pp. 572-574
76+
2. M. Vajta, "Some remarks on Padé-approximations",
77+
3rd TEMPUS-INTCOM Symposium
7178
"""
79+
if numdeg is None:
80+
numdeg = n
81+
elif numdeg < 0:
82+
numdeg += n
83+
84+
if not T >= 0:
85+
raise ValueError("require T >= 0")
86+
if not n >= 0:
87+
raise ValueError("require n >= 0")
88+
if not (0 <= numdeg <= n):
89+
raise ValueError("require 0 <= numdeg <= n")
90+
7291
if T == 0:
7392
num = [1,]
7493
den = [1,]
7594
else:
76-
num = [0. for i in range(n+1)]
95+
num = [0. for i in range(numdeg+1)]
7796
num[-1] = 1.
97+
cn = 1.
98+
for k in range(1, numdeg+1):
99+
# derived from Gloub and van Loan eq. for Dpq(z) on p. 572
100+
# this accumulative style follows Alg 11.3.1
101+
cn *= -T * (numdeg - k + 1)/(numdeg + n - k + 1)/k
102+
num[numdeg-k] = cn
103+
78104
den = [0. for i in range(n+1)]
79105
den[-1] = 1.
80-
c = 1.
106+
cd = 1.
81107
for k in range(1, n+1):
82-
c = T * c * (n - k + 1)/(2 * n - k + 1)/k
83-
num[n - k] = c * (-1)**k
84-
den[n - k] = c
108+
# see cn above
109+
cd *= T * (n - k + 1)/(numdeg + n - k + 1)/k
110+
den[n-k] = cd
111+
85112
num = [coeff/den[0] for coeff in num]
86113
den = [coeff/den[0] for coeff in den]
87114
return num, den

‎control/tests/delay_test.py

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
#!/usr/bin/env python -*-coding: utf-8-*-
2+
#
3+
# Test Pade approx
4+
#
5+
# Primitive; ideally test to numerical limits
6+
7+
from __future__ import division
8+
9+
import unittest
10+
11+
import numpy as np
12+
13+
from control.delay import pade
14+
15+
16+
class TestPade(unittest.TestCase):
17+
18+
# Reference data from Miklos Vajta's paper "Some remarks on
19+
# Padé-approximations", Table 1, with corrections. The
20+
# corrections are to highest power coeff in numerator for
21+
# (ddeg,ndeg)=(4,3) and (5,4); use Eq (12) in the paper to verify
22+
23+
# all for T = 1
24+
ref = [
25+
# dendeg numdeg den num
26+
( 1, 1, [1,2], [-1,2]),
27+
( 1, 0, [1,1], [1]),
28+
( 2, 2, [1,6,12], [1,-6,12]),
29+
( 2, 1, [1,4,6], [-2,6]),
30+
( 3, 3, [1,12,60,120], [-1,12,-60,120]),
31+
( 3, 2, [1,9,36,60], [3,-24,60]),
32+
( 4, 4, [1,20,180,840,1680], [1,-20,180,-840,1680]),
33+
( 4, 3, [1,16,120,480,840], [-4,60,-360,840]),
34+
( 5, 5, [1,30,420,3360,15120,30240], [-1,30,-420,3360,-15120,30240]),
35+
( 5, 4, [1,25,300,2100,8400,15120,], [5,-120,1260,-6720,15120]),
36+
]
37+
38+
def testRefs(self):
39+
"test reference cases for T=1"
40+
T = 1
41+
for dendeg, numdeg, refden, refnum in self.ref:
42+
num, den = pade(T, dendeg, numdeg)
43+
np.testing.assert_array_almost_equal_nulp(np.array(refden), den, nulp=2)
44+
np.testing.assert_array_almost_equal_nulp(np.array(refnum), num, nulp=2)
45+
46+
def testTvalues(self):
47+
"test reference cases for T!=1"
48+
Ts = [1/53, 21.95]
49+
for dendeg, numdeg, baseden, basenum in self.ref:
50+
for T in Ts:
51+
refden = T**np.arange(dendeg, -1, -1)*baseden
52+
refnum = T**np.arange(numdeg, -1, -1)*basenum
53+
refnum /= refden[0]
54+
refden /= refden[0]
55+
num, den = pade(T, dendeg, numdeg)
56+
np.testing.assert_array_almost_equal_nulp(refden, den, nulp=2)
57+
np.testing.assert_array_almost_equal_nulp(refnum, num, nulp=2)
58+
59+
def testErrors(self):
60+
"ValueError raised for invalid arguments"
61+
self.assertRaises(ValueError,pade,-1,1) # T<0
62+
self.assertRaises(ValueError,pade,1,-1) # dendeg < 0
63+
self.assertRaises(ValueError,pade,1,2,-3) # numdeg < 0
64+
self.assertRaises(ValueError,pade,1,2,3) # numdeg > dendeg
65+
66+
def testNumdeg(self):
67+
"numdeg argument follows docs"
68+
# trivialish - interface check, not math check
69+
T = 1
70+
dendeg = 5
71+
ref = [pade(T,dendeg,numdeg)
72+
for numdeg in range(0,dendeg+1)]
73+
testneg = [pade(T,dendeg,numdeg)
74+
for numdeg in range(-dendeg,0)]
75+
self.assertEqual(ref[:-1],testneg)
76+
self.assertEqual(ref[-1], pade(T,dendeg,dendeg))
77+
self.assertEqual(ref[-1], pade(T,dendeg,None))
78+
self.assertEqual(ref[-1], pade(T,dendeg))
79+
80+
def testT0(self):
81+
"T=0 always returns [1],[1]"
82+
T = 0
83+
refnum = [1.0]
84+
refden = [1.0]
85+
for dendeg in range(1, 6):
86+
for numdeg in range(0, dendeg+1):
87+
num, den = pade(T, dendeg, numdeg)
88+
np.testing.assert_array_almost_equal_nulp(np.array(refnum), np.array(num))
89+
np.testing.assert_array_almost_equal_nulp(np.array(refden), np.array(den))
90+
91+
if __name__ == '__main__':
92+
unittest.main()

0 commit comments

Comments
 (0)