Modified _common_den to handle complex roots more accurately.

Kevin Chen · Kevin Chen · commit 2122225174d6 · 2011-02-08T22:18:23.000Z
Minor changes in TestConvert.py and TestBDAlg.py.

Kevin K. Chen &lt;kkchen@princeton.edu&gt;
diff --git a/src/TestBDAlg.py b/src/TestBDAlg.py
@@ -38,7 +38,6 @@ def testScalarSS(self):
         ans1 = feedback(self.x1, self.sys2)
         ans2 = feedback(self.x1, self.sys2, 1.)
         
-        # This one doesn't work yet.  The feedback system has too many states.
         np.testing.assert_array_almost_equal(ans1.A, [[-1.5, 4.], [13., 2.]])
         np.testing.assert_array_almost_equal(ans1.B, [[2.5], [-10.]])
         np.testing.assert_array_almost_equal(ans1.C, [[-2.5, 0.]])
@@ -65,7 +64,6 @@ def testSSScalar(self):
         ans1 = feedback(self.sys2, self.x1)
         ans2 = feedback(self.sys2, self.x1, 1.)
 
-        # This one doesn't work yet.  The feedback system has too many states.
         np.testing.assert_array_almost_equal(ans1.A, [[-1.5, 4.], [13., 2.]])
         np.testing.assert_array_almost_equal(ans1.B, [[1.], [-4.]])
         np.testing.assert_array_almost_equal(ans1.C, [[1., 0.]])
diff --git a/src/TestConvert.py b/src/TestConvert.py
@@ -4,10 +4,13 @@
 
 Test state space and transfer function conversion.
 
-Currently, this unit test script is not complete.  It converts
-several random state spaces back and forth between state space and transfer
-function representations, and asserts that the conversion outputs are correct.
-As it stands, tests pass but there is some rounding error in one of the conversions leading to imaginary numbers. See the warning message. This script may be used to diagnose errors.
+Currently, this unit test script is not complete.  It converts several random
+state spaces back and forth between state space and transfer function
+representations.  Ideally, it should be able to assert that the conversion
+outputs are correct.  This is not yet implemented.
+
+Also, the conversion seems to enter an infinite loop once in a while.  The cause
+of this is unknown.
 
 """
 
@@ -22,11 +25,11 @@ def setUp(self):
         """Set up testing parameters."""
 
         # Number of times to run each of the randomized tests.
-        self.numTests = 1
+        self.numTests = 10
         # Maximum number of states to test + 1
-        self.maxStates = 3
+        self.maxStates = 2
         # Maximum number of inputs and outputs to test + 1
-        self.maxIO = 3
+        self.maxIO = 2
         # Set to True to print systems to the output.
         self.debug = True
 
diff --git a/src/xferfcn.py b/src/xferfcn.py
@@ -75,8 +75,8 @@
 """
 
 # External function declarations
-from numpy import angle, array, empty, finfo, insert, ndarray, ones, polyadd, \
-    polymul, polyval, roots, sort, zeros
+from numpy import angle, any, array, empty, finfo, insert, ndarray, ones, \
+    polyadd, polymul, polyval, roots, sort, sqrt, zeros
 from scipy.signal import lti
 from copy import deepcopy
 from slycot import tb04ad
@@ -491,7 +491,10 @@ def _common_den(self):
         [something] array.
 
         """
-         
+        
+        # Machine precision for floats.
+        eps = finfo(float).eps
+
         # A sorted list to keep track of cumulative poles found as we scan
         # self.den.
         poles = []
@@ -512,8 +515,7 @@ def _common_den(self):
                 # cumulative poles, until one of them reaches the end.  Keep in
                 # mind that both lists are always sorted.
                 while cp_ind < len(currentpoles) and p_ind < len(poles):
-                    if abs(currentpoles[cp_ind] - poles[p_ind]) < (10 *
-                        finfo(float).eps):
+                    if abs(currentpoles[cp_ind] - poles[p_ind]) < (10 * eps):
                         # If the current element of both lists match, then we're
                         # good.  Move to the next pair of elements.
                         cp_ind += 1
@@ -558,8 +560,21 @@ def _common_den(self):
         
         # Construct the common denominator.
         den = 1.
-        for p in poles:
-            den = polymul(den, [1., -p])
+        n = 0
+        while n < len(poles):
+            if abs(poles[n].imag) > 10 * eps:
+                # To prevent buildup of imaginary part error, handle complex
+                # pole pairs together.
+                quad = polymul([1., -poles[n]], [1., -poles[n+1]])
+                assert all(quad.imag < eps), "The quadratic has a nontrivial \
+imaginary part: %g" % quad.imag.max()
+                quad = quad.real
+
+                den = polymul(den, quad)
+                n += 2
+            else:
+                den = polymul(den, [1., -poles[n].real])
+                n += 1
 
         # Modify the numerators so that they each take the common denominator.
         num = deepcopy(self.num)
@@ -585,6 +600,11 @@ def _common_den(self):
                     zeros(largest - len(num[i][j])))
         # Finally, convert the numerator to a 3-D array.
         num = array(num)
+        # Remove trivial imaginary parts.  Check for nontrivial imaginary parts.
+        if any(abs(num.imag) > sqrt(eps)):
+            print ("Warning: The numerator has a nontrivial nontrivial part: %g"
+                % abs(num.imag).max())
+        num = num.real
 
         return num, den
 
