Skip to content

Feature: allow differing numerator and denominator degree in Pade app… #76

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Mar 8, 2016
Merged

Feature: allow differing numerator and denominator degree in Pade app… #76

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 39 additions & 12 deletions control/delay.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
# -*-coding: utf-8-*-
#! TODO: add module docstring
# delay.py - functions involving time delays
#
Expand Down Expand Up @@ -41,12 +42,11 @@
#
# $Id$

# Python 3 compatability (needs to go here)
from __future__ import print_function
from __future__ import division

__all__ = ['pade']

def pade(T, n=1):
def pade(T, n=1, numdeg=None):
"""
Create a linear system that approximates a delay.

Expand All @@ -56,8 +56,12 @@ def pade(T, n=1):
----------
T : number
time delay
n : integer
order of approximation
n : positive integer
degree of denominator of approximation
numdeg: integer, or None (the default)
If None, numerator degree equals denominator degree
If >= 0, specifies degree of numerator
If < 0, numerator degree is n+numdeg

Returns
-------
Expand All @@ -66,22 +70,45 @@ def pade(T, n=1):

Notes
-----
Based on an algorithm in Golub and van Loan, "Matrix Computation" 3rd.
Ed. pp. 572-574.
Based on:
1. Algorithm 11.3.1 in Golub and van Loan, "Matrix Computation" 3rd.
Ed. pp. 572-574
2. M. Vajta, "Some remarks on Padé-approximations",
3rd TEMPUS-INTCOM Symposium
"""
if numdeg is None:
numdeg = n
elif numdeg < 0:
numdeg += n

if not T >= 0:
raise ValueError("require T >= 0")
if not n >= 0:
raise ValueError("require n >= 0")
if not (0 <= numdeg <= n):
raise ValueError("require 0 <= numdeg <= n")

if T == 0:
num = [1,]
den = [1,]
else:
num = [0. for i in range(n+1)]
num = [0. for i in range(numdeg+1)]
num[-1] = 1.
cn = 1.
for k in range(1, numdeg+1):
# derived from Gloub and van Loan eq. for Dpq(z) on p. 572
# this accumulative style follows Alg 11.3.1
cn *= -T * (numdeg - k + 1)/(numdeg + n - k + 1)/k
num[numdeg-k] = cn

den = [0. for i in range(n+1)]
den[-1] = 1.
c = 1.
cd = 1.
for k in range(1, n+1):
c = T * c * (n - k + 1)/(2 * n - k + 1)/k
num[n - k] = c * (-1)**k
den[n - k] = c
# see cn above
cd *= T * (n - k + 1)/(numdeg + n - k + 1)/k
den[n-k] = cd

num = [coeff/den[0] for coeff in num]
den = [coeff/den[0] for coeff in den]
return num, den
92 changes: 92 additions & 0 deletions control/tests/delay_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,92 @@
#!/usr/bin/env python -*-coding: utf-8-*-
#
# Test Pade approx
#
# Primitive; ideally test to numerical limits

from __future__ import division

import unittest

import numpy as np

from control.delay import pade


class TestPade(unittest.TestCase):

# Reference data from Miklos Vajta's paper "Some remarks on
# Padé-approximations", Table 1, with corrections. The
# corrections are to highest power coeff in numerator for
# (ddeg,ndeg)=(4,3) and (5,4); use Eq (12) in the paper to verify

# all for T = 1
ref = [
# dendeg numdeg den num
( 1, 1, [1,2], [-1,2]),
( 1, 0, [1,1], [1]),
( 2, 2, [1,6,12], [1,-6,12]),
( 2, 1, [1,4,6], [-2,6]),
( 3, 3, [1,12,60,120], [-1,12,-60,120]),
( 3, 2, [1,9,36,60], [3,-24,60]),
( 4, 4, [1,20,180,840,1680], [1,-20,180,-840,1680]),
( 4, 3, [1,16,120,480,840], [-4,60,-360,840]),
( 5, 5, [1,30,420,3360,15120,30240], [-1,30,-420,3360,-15120,30240]),
( 5, 4, [1,25,300,2100,8400,15120,], [5,-120,1260,-6720,15120]),
]

def testRefs(self):
"test reference cases for T=1"
T = 1
for dendeg, numdeg, refden, refnum in self.ref:
num, den = pade(T, dendeg, numdeg)
np.testing.assert_array_almost_equal_nulp(np.array(refden), den, nulp=2)
np.testing.assert_array_almost_equal_nulp(np.array(refnum), num, nulp=2)

def testTvalues(self):
"test reference cases for T!=1"
Ts = [1/53, 21.95]
for dendeg, numdeg, baseden, basenum in self.ref:
for T in Ts:
refden = T**np.arange(dendeg, -1, -1)*baseden
refnum = T**np.arange(numdeg, -1, -1)*basenum
refnum /= refden[0]
refden /= refden[0]
num, den = pade(T, dendeg, numdeg)
np.testing.assert_array_almost_equal_nulp(refden, den, nulp=2)
np.testing.assert_array_almost_equal_nulp(refnum, num, nulp=2)

def testErrors(self):
"ValueError raised for invalid arguments"
self.assertRaises(ValueError,pade,-1,1) # T<0
self.assertRaises(ValueError,pade,1,-1) # dendeg < 0
self.assertRaises(ValueError,pade,1,2,-3) # numdeg < 0
self.assertRaises(ValueError,pade,1,2,3) # numdeg > dendeg

def testNumdeg(self):
"numdeg argument follows docs"
# trivialish - interface check, not math check
T = 1
dendeg = 5
ref = [pade(T,dendeg,numdeg)
for numdeg in range(0,dendeg+1)]
testneg = [pade(T,dendeg,numdeg)
for numdeg in range(-dendeg,0)]
self.assertEqual(ref[:-1],testneg)
self.assertEqual(ref[-1], pade(T,dendeg,dendeg))
self.assertEqual(ref[-1], pade(T,dendeg,None))
self.assertEqual(ref[-1], pade(T,dendeg))

def testT0(self):
"T=0 always returns [1],[1]"
T = 0
refnum = [1.0]
refden = [1.0]
for dendeg in range(1, 6):
for numdeg in range(0, dendeg+1):
num, den = pade(T, dendeg, numdeg)
np.testing.assert_array_almost_equal_nulp(np.array(refnum), np.array(num))
np.testing.assert_array_almost_equal_nulp(np.array(refden), np.array(den))

if __name__ == '__main__':
unittest.main()