Merge pull request #242 from murrayrm/fix_common_den

fix issue #240: _common_den was not padding properly

control/tests/convert_test.py

@@ -245,6 +245,29 @@ def testTf2SsDuplicatePoles(self):
         except ImportError:
             print("Slycot not present, skipping")
 
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def test_tf2ss_robustness(self):
+        """Unit test to make sure that tf2ss is working correctly.
+         Source: https://github.com/python-control/python-control/issues/240
+        """
+        import control
+
+        num =  [ [[0], [1]],           [[1],   [0]] ]
+        den1 = [ [[1], [1,1]],         [[1,4], [1]] ]
+        sys1tf = control.tf(num, den1)
+        sys1ss = control.tf2ss(sys1tf)
+
+        # slight perturbation
+        den2 = [ [[1], [1e-10, 1, 1]], [[1,4], [1]] ]
+        sys2tf = control.tf(num, den2)
+        sys2ss = control.tf2ss(sys2tf)
+
+        # Make sure that the poles match for StateSpace and TransferFunction
+        np.testing.assert_array_almost_equal(np.sort(sys1tf.pole()),
+                                             np.sort(sys1ss.pole()))
+        np.testing.assert_array_almost_equal(np.sort(sys2tf.pole()),
+                                             np.sort(sys2ss.pole()))
+
 def suite():
    return unittest.TestLoader().loadTestsFromTestCase(TestConvert)
 

control/tests/timeresp_test.py

@@ -15,6 +15,7 @@
 from control.statesp import *
 from control.xferfcn import TransferFunction, _convertToTransferFunction
 from control.dtime import c2d
+from control.exception import slycot_check
 
 class TestTimeresp(unittest.TestCase):
     def setUp(self):
@@ -233,6 +234,23 @@ def test_discrete_initial(self):
         t, yout = impulse_response(h1, np.arange(4))
         np.testing.assert_array_equal(yout[0], [0., 1., 0., 0.])
 
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def test_step_robustness(self):
+        "Unit test: https://github.com/python-control/python-control/issues/240"
+        # Create 2 input, 2 output system
+        num =  [ [[0], [1]],           [[1],   [0]] ]
+
+        den1 = [ [[1], [1,1]],         [[1,4], [1]] ]
+        sys1 = TransferFunction(num, den1)
+
+        den2 = [ [[1], [1e-10, 1, 1]], [[1,4], [1]] ]   # slight perturbation
+        sys2 = TransferFunction(num, den2)
+
+        # Compute step response from input 1 to output 1, 2
+        t1, y1 = step_response(sys1, input=0)
+        t2, y2 = step_response(sys2, input=0)
+        np.testing.assert_array_almost_equal(y1, y2)
+
 def suite():
     return unittest.TestLoader().loadTestsFromTestCase(TestTimeresp)
 

control/xferfcn.py

@@ -840,9 +840,13 @@ def _common_den(self, imag_tol=None):
                     num[i,j,0] = poleset[i][j][2]
             else:
                 # create the denominator matching this input
+                # polyfromroots gives coeffs in opposite order from what we use
+                # coefficients should be padded on right, ending at np
                 np = len(poles[j])
                 den[j,np::-1] = polyfromroots(poles[j]).real
                 denorder[j] = np
+
+                # now create the numerator, also padded on the right
                 for i in range(self.outputs):
                     # start with the current set of zeros for this output
                     nwzeros = list(poleset[i][j][0])
@@ -851,14 +855,15 @@ def _common_den(self, imag_tol=None):
                     for ip in chain(poleset[i][j][3],
                                     range(poleset[i][j][4],np)):
                         nwzeros.append(poles[j][ip])
-                    
+
                     numpoly = poleset[i][j][2] * polyfromroots(nwzeros).real 
-                    m = npmax - len(numpoly)
-                    #print(j,i,m,len(numpoly),len(poles[j]))
-                    if m < 0:
-                        num[i,j,::-1] = numpoly
-                    else:
-                        num[i,j,:m:-1] = numpoly   
+                    # print(numpoly, den[j])
+                    # polyfromroots gives coeffs in opposite order => invert
+                    # numerator polynomial should be padded on left and right
+                    #   ending at np to line up with what td04ad expects...
+                    num[i, j, np+1-len(numpoly):np+1] = numpoly[::-1]
+                    # print(num[i, j])
+
         if (abs(den.imag) > epsnm).any():
             print("Warning: The denominator has a nontrivial imaginary part: %f"
                       % abs(den.imag).max())