GeezerVison is showing github.com — Feature: allow differing numerator and denominator degree in Pade app… · python-control/python-control@2d79f4e
File treeExpand file treeCollapse file tree Open diff view settings Expand file treeCollapse file tree Open diff view settings Collapse file Copy file name to clipboardExpand all lines: control/delay.py +39-12Lines changed: 39 additions & 12 deletions | Original file line number | Diff line number | Diff line change |
|---|
| | 1 | +# -*-coding: utf-8-*- | 1 | 2 | #! TODO: add module docstring | 2 | 3 | # delay.py - functions involving time delays | 3 | 4 | # | | 41 | 42 | # | 42 | 43 | # $Id$ | 43 | 44 |
| 44 | | -# Python 3 compatability (needs to go here) | 45 | | -from __future__ import print_function | | 45 | +from __future__ import division | 46 | 46 |
| 47 | 47 | __all__ = ['pade'] | 48 | 48 |
| 49 | | -def pade(T, n=1): | | 49 | +def pade(T, n=1, numdeg=None): | 50 | 50 | """ | 51 | 51 | Create a linear system that approximates a delay. | 52 | 52 |
| @@ -56,8 +56,12 @@ def pade(T, n=1): | 56 | 56 | ---------- | 57 | 57 | T : number | 58 | 58 | 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 | 61 | 65 |
| 62 | 66 | Returns | 63 | 67 | ------- | @@ -66,22 +70,45 @@ def pade(T, n=1): | 66 | 70 |
| 67 | 71 | Notes | 68 | 72 | ----- | 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 | 71 | 78 | """ | | 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 | + | 72 | 91 | if T == 0: | 73 | 92 | num = [1,] | 74 | 93 | den = [1,] | 75 | 94 | else: | 76 | | - num = [0. for i in range(n+1)] | | 95 | + num = [0. for i in range(numdeg+1)] | 77 | 96 | 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 | + | 78 | 104 | den = [0. for i in range(n+1)] | 79 | 105 | den[-1] = 1. | 80 | | - c = 1. | | 106 | + cd = 1. | 81 | 107 | 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 | + | 85 | 112 | num = [coeff/den[0] for coeff in num] | 86 | 113 | den = [coeff/den[0] for coeff in den] | 87 | 114 | return num, den |
Collapse file Copy file name to clipboard +92Lines changed: 92 additions & 0 deletions | Original file line number | Diff line number | Diff line change |
|---|
| | 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() |
|