From 2d79f4e763a7216e240274c112122e809d8a893b Mon Sep 17 00:00:00 2001 From: Rory Yorke Date: Mon, 4 Jan 2016 08:22:27 +0200 Subject: [PATCH] Feature: allow differing numerator and denominator degree in Pade approx. As suggested by @jason-s in python-control/python-control#73 This is a relatively minimal change, and should do no worse and no better than the previous implementation when numerator and denominator have equal degree. --- control/delay.py | 51 +++++++++++++++----- control/tests/delay_test.py | 92 +++++++++++++++++++++++++++++++++++++ 2 files changed, 131 insertions(+), 12 deletions(-) create mode 100644 control/tests/delay_test.py diff --git a/control/delay.py b/control/delay.py index d99406992..d6350d45b 100644 --- a/control/delay.py +++ b/control/delay.py @@ -1,3 +1,4 @@ +# -*-coding: utf-8-*- #! TODO: add module docstring # delay.py - functions involving time delays # @@ -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. @@ -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 ------- @@ -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 diff --git a/control/tests/delay_test.py b/control/tests/delay_test.py new file mode 100644 index 000000000..17c049d24 --- /dev/null +++ b/control/tests/delay_test.py @@ -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()