python-control · lytex · Feb 9, 2020 · Feb 29, 2020 · Feb 29, 2020 · Feb 29, 2020
diff --git a/control/statesp.py b/control/statesp.py
@@ -93,9 +93,15 @@ def _ssmatrix(data, axis=1):
     # Convert the data into an array or matrix, as configured
     # If data is passed as a string, use (deprecated?) matrix constructor
     if config.defaults['statesp.use_numpy_matrix'] or isinstance(data, str):
-        arr = np.matrix(data, dtype=float)
+        if np.isrealobj(data):
+            arr = np.matrix(data, dtype=float)
+        elif np.iscomplexobj(data):
+            arr = np.matrix(data, dtype=complex)
     else:
-        arr = np.array(data, dtype=float)
+        if np.isrealobj(data):
+            arr = np.array(data, dtype=float)
+        elif np.iscomplexobj(data):
+            arr = np.array(data, dtype=complex)
     ndim = arr.ndim
     shape = arr.shape
 
@@ -572,10 +578,15 @@ def zero(self):
         # Use AB08ND from Slycot if it's available, otherwise use
         # scipy.lingalg.eigvals().
         try:
-            from slycot import ab08nd
-
-            out = ab08nd(self.A.shape[0], self.B.shape[1], self.C.shape[0],
-                         self.A, self.B, self.C, self.D)
+            if np.isrealobj(self.A) and np.isrealobj(self.B) and\
+            np.isrealobj(self.C) and np.isrealobj(self.D):
+                from slycot import ab08nd
+                out = ab08nd(self.A.shape[0], self.B.shape[1], self.C.shape[0],
+                self.A, self.B, self.C, self.D)
+            else:
+                from slycot import ab08nz
+                out = ab08nz(self.A.shape[0], self.B.shape[1], self.C.shape[0],
+                            self.A, self.B, self.C, self.D)
             nu = out[0]
             if nu == 0:
                 return np.array([])

diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
@@ -7,13 +7,12 @@
 import numpy as np
 from numpy.linalg import solve
 from scipy.linalg import eigvals, block_diag
-from control import matlab
+from control import matlab, series
 from control.statesp import StateSpace, _convertToStateSpace, tf2ss
 from control.xferfcn import TransferFunction, ss2tf
 from control.lti import evalfr
 from control.exception import slycot_check
 
-
 class TestStateSpace(unittest.TestCase):
     """Tests for the StateSpace class."""
 
@@ -63,6 +62,50 @@ def setUp(self):
         D623 = np.zeros((3, 2))
         self.sys623 = StateSpace(A623, B623, C623, D623)
 
+        # Systems with complex matrices, sysC322, sysC222, sysC623
+        # These systems have the same shape as the previous ones
+
+        A = np.array([[6.0 + 9.0j, 8.0 + 9.0j,  -4.0 - 7.0j],
+                      [8.0 - 7.0j, 3.0,         1.0 - 1.0j],
+                      [-7.0 + 9.0j, -8.0 + 6.0j, 9.0 + 8.0j]])
+        B = np.array([[6.0 + 3.0j, -9.0 - 2.0j],
+                      [9.0 + 5.0j,  7.0 + 3.0j],
+                      [3.0 + 5.0j, 8.0 - 6.0j]])
+        C = np.array([[4.0 + 4.0j, -4.0 + 9.0j, -8.0 - 1.0j],
+        [-9.0 - 3.0j, -9.0 - 9.0j, 6.0 - 2.0j]])
+        D = np.array([[5.0 - 1.0j, -6.0 + 4.0j],
+        [6.0 + 3.0j, 5.0j]])
+        self.sysC322 = StateSpace(A, B, C, D)
+
+        A = np.array([[-4.0 - 7.0j, 3.0 + 9.0j],
+        [3.0, -6.0 - 3.0j]])
+        B = np.array([[2.0, 5.0 + 7.0j],
+        [-5.0 + 4.0j, -5.0 + 9.0j]])
+        C = np.array([[1.0 + 6.0j, -7.0 + 6.0j],
+        [-7.0 - 5.0j, -5.0 - 5.0j]])
+        D = np.array([[8.0 + 2.0j, -6.0 - 3.0j],
+        [-3.0 - 1.0j, -5.0 + 6.0j]])
+        self.sysC222 = StateSpace(A, B, C, D)
+
+        A = np.array(
+        [[2.0 - 8.0j, -2.0 + 6.0j, 8.0 - 1.0j, -6.0 + 7.0j, -5.0 - 3.0j, -5.0 - 6.0j],
+        [1.0 - 1.0j, 1.0 + 7.0j, -7.0 + 8.0j, 6.0 + 2.0j, 3.0, 8.0 - 5.0j],
+        [8.0 - 7.0j, -8.0 - 8.0j, 1.0 - 6.0j, -4.0 + 1.0j, 4.0 - 2.0j, -7.0 - 2.0j],
+        [-4.0 + 9.0j, -8.0 - 2.0j, -1.0 - 4.0j, 1.0 - 7.0j, 5.0 - 8.0j, 6.0 - 9.0j],
+        [5.0 - 9.0j, 1.0 - 5.0j, -9.0 - 7.0j, -6.0 + 7.0j, -1.0 - 5.0j, 1.0 + 8.0j],
+        [5.0 + 5.0j, 5.0 + 6.0j, -3.0 - 7.0j, 2.0 + 2.0j, -8.0 - 7.0j, 9.0 + 8.0j]])
+        B = np.array([[6.0j, 5.0 + 3.0j, 8.0 + 4.0j],
+        [-9.0j, -2.0 - 1.0j, 9.0 - 6.0j],
+        [-3.0 - 9.0j, -5.0 + 1.0j, 1.0 - 2.0j],
+        [8.0 - 6.0j, -2.0 - 4.0j, -8.0 + 2.0j],
+        [-2.0 + 3.0j, -8.0 + 5.0j, -5.0 + 5.0j],
+        [-7.0 + 4.0j, -7.0 - 6.0j, -3.0 - 8.0j]])
+        C = np.array([[8.0 + 6.0j, -3.0j, -1.0 + 7.0j, 2.0j, 6.0 - 6.0j, 3.0 - 1.0j],
+        [5.0 + 1.0j, -1.0 + 8.0j, -4.0 + 1.0j, 2.0j, 6.0 - 4.0j, -2.0 - 5.0j]])
+        D = np.array([[7.0 - 4.0j, -5.0 - 1.0j, -5.0 + 8.0j],
+        [-6.0 + 8.0j, -6.0 - 6.0j, -1.0 + 9.0j]])
+        self.sysC623 = StateSpace(A, B, C, D)
+
     def test_D_broadcast(self):
         """Test broadcast of D=0 to the right shape"""
         # Giving D as a scalar 0 should broadcast to the right shape
@@ -519,6 +562,179 @@ def test_lft(self):
         np.testing.assert_allclose(np.array(pk.C).reshape(-1), Cmatlab)
         np.testing.assert_allclose(np.array(pk.D).reshape(-1), Dmatlab)
 
+    def test_pole_complex(self):
+        """Evaluate the poles of a complex MIMO system."""
+
+        p = np.sort(self.sysC322.pole())
+        true_p = np.sort([21.7521962567499187 +7.8381752267389437j
+        -6.4620734598457208 +2.8701578198630169j,
+        2.7098772030958163 +6.2916669533980274j])
+
+        np.testing.assert_array_almost_equal(p, true_p)
+
+
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def test_zero_siso_complex(self):
+        """Evaluate the zeros of a complex SISO system."""
+        # extract only first input / first output system of sysC222. This system is denoted sysC111
+        #  or tfC111
+        tfC111 = ss2tf(self.sysC222)
+        sysC111 = tf2ss(tfC111[0, 0])
+
+        # compute zeros as root of the characteristic polynomial at the numerator of tfC111
+        # this method is simple and assumed as valid in this test
+        true_z = np.sort(tfC111[0, 0].zero())
+        # Compute the zeros through ab08nd, which is tested here
+        z = np.sort(sysC111.zero())
+
+        np.testing.assert_almost_equal(true_z, z)
+
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def test_zero_mimo_sysC322_square_complex(self):
+        """Evaluate the zeros of a square complex MIMO system."""
+
+        z = np.sort(self.sysC322.zero())
+        true_z = np.sort([36.493759 + 8.073864j,
+                -4.7860079 + 29.3266582j,
+                -7.4509692 +5.5262915j])
+        np.testing.assert_array_almost_equal(z, true_z)
+
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def test_zero_mimo_sysC222_square_complex(self):
+        """Evaluate the zeros of a square complex MIMO system."""
+
+        z = np.sort(self.sysC222.zero())
+        true_z = np.sort([-10.56850055,
+                        3.3685005])
+        np.testing.assert_array_almost_equal(z, true_z)
+
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def test_zero_mimo_sysC623_non_square_complex(self):
+        """Evaluate the zeros of a non square complex MIMO system."""
+
+        z = np.sort(self.sysC623.zero())
+        true_z = np.sort([]) # System has no transmission zeros, not sure what slycot returns
+        np.testing.assert_array_almost_equal(z, true_z)
+
+    def test_add_ss_complex(self):
+        """Add two complex MIMO systems."""
+
+        A = np.array([[6.0 +9.0j, 8.0 +9.0j, -4.0 -7.0j, 0.0, 0.0],
+                    [8.0 -7.0j, 3.0, 1.0 -1.0j, 0.0, 0.0],
+                    [-7.0 +9.0j, -8.0 +6.0j, 9.0 +8.0j, 0.0, 0.0],
+                    [0.0, 0.0, 0.0, -4.0 -7.0j, 3.0 +9.0j],
+                    [0.0, 0.0, 0.0, 3.0, -6.0 -3.0j]])
+        B = np.array([[6.0 +3.0j, -9.0 -2.0j],
+                    [9.0 +5.0j, 7.0 +3.0j],
+                    [3.0 +5.0j, 8.0 -6.0j],
+                    [2.0, 5.0 +7.0j],
+                    [-5.0 +4.0j, -5.0 +9.0j]])
+        C = np.array([[4.0 +4.0j, -4.0 +9.0j, -8.0 -1.0j, 1.0 +6.0j, -7.0 +6.0j],
+                    [-9.0 -3.0j, -9.0 -9.0j, 6.0 -2.0j, -7.0 -5.0j, -5.0 -5.0j]])
+        D = np.array([[13.0 +1.0j, -12.0 +1.0j],
+                    [3.0 +2.0j, -5.0 +11.0j]])
+
+        sys = self.sysC322 + self.sysC222
+
+        np.testing.assert_array_almost_equal(sys.A, A)
+        np.testing.assert_array_almost_equal(sys.B, B)
+        np.testing.assert_array_almost_equal(sys.C, C)
+        np.testing.assert_array_almost_equal(sys.D, D)
+
+    def test_subtract_ss_complex(self):
+        """Subtract two complex MIMO systems."""
+
+        A = np.array([[6.0 +9.0j, 8.0 +9.0j, -4.0 -7.0j, 0.0, 0.0],
+                    [8.0 -7.0j, 3.0, 1.0 -1.0j, 0.0, 0.0],
+                    [-7.0 +9.0j, -8.0 +6.0j, 9.0 +8.0j, 0.0, 0.0],
+                    [0.0, 0.0, 0.0, -4.0 -7.0j, 3.0 +9.0j],
+                    [0.0, 0.0, 0.0, 3.0, -6.0 -3.0j]])
+        B = np.array([[6.0 +3.0j, -9.0 -2.0j],
+                    [9.0 +5.0j, 7.0 +3.0j],
+                    [3.0 +5.0j, 8.0 -6.0j],
+                    [2.0, 5.0 +7.0j],
+                    [-5.0 +4.0j, -5.0 +9.0j]])
+        C = np.array([[4.0 +4.0j, -4.0 +9.0j, -8.0 -1.0j, -1.0 -6.0j, 7.0 -6.0j],
+                    [-9.0 -3.0j, -9.0 -9.0j, 6.0 -2.0j, 7.0 +5.0j, 5.0 +5.0j]])
+        D = np.array([[-3.0 -3.0j, 7.0j],
+                    [9.0 +4.0j, 5.0 -1.0j]])
+
+        sys = self.sysC322 - self.sysC222
+
+        np.testing.assert_array_almost_equal(sys.A, A)
+        np.testing.assert_array_almost_equal(sys.B, B)
+        np.testing.assert_array_almost_equal(sys.C, C)
+        np.testing.assert_array_almost_equal(sys.D, D)
+
+    def test_multiply_ss_complex(self):
+        """Multiply two complex MIMO systems."""
+
+        A = np.array([[6.0 +9.0j, 8.0 +9.0j, -4.0 -7.0j, 41.0 +98.0j, -25.0 +70.0j],
+                    [8.0 -7.0j, 3.0, 1.0 -1.0j, -55.0 +3.0j, -113.0 -31.0j],
+                    [-7.0 +9.0j, -8.0 +6.0j, 9.0 +8.0j, -113.0 +25.0j, -121.0 -27.0j],
+                    [0.0, 0.0, 0.0, -4.0 -7.0j, 3.0 +9.0j],
+                    [0.0, 0.0, 0.0, 3.0, -6.0 -3.0j]])
+        B = np.array([[67.0 +51.0j, 30.0 -80.0j],
+                    [44.0 +42.0j, -92.0 -30.0j],
+                    [-16.0 +56.0j, -7.0 +39.0j],
+                    [2.0, 5.0 +7.0j],
+                    [-5.0 +4.0j, -5.0 +9.0j]])
+        C = np.array([[4.0 +4.0j, -4.0 +9.0j, -8.0 -1.0j, 73.0 +31.0j, 21.0 +47.0j],
+                    [-9.0 -3.0j, -9.0 -9.0j, 6.0 -2.0j, 13.0 +4.0j, -35.0 -10.0j]])
+        D = np.array([[64.0 -4.0j, -27.0 -65.0j],
+                    [47.0 +21.0j, -57.0 -61.0j]])
+
+        sys_aux = StateSpace(A, B, C, D)
+
+        sys = series(self.sysC322, self.sysC222)
+
+        np.testing.assert_array_almost_equal(sys.A, sys_aux.A)
+        np.testing.assert_array_almost_equal(sys.B, sys_aux.B)
+        np.testing.assert_array_almost_equal(sys.C, sys_aux.C)
+        np.testing.assert_array_almost_equal(sys.D, sys_aux.D)
+
+    def test_evalfr_complex(self):
+        """Evaluate the frequency response at one frequency."""
+
+        resp = [[4.6799374 - 34.9854626j,
+                -10.8392352 - 10.3778031j],
+                [28.8313352 + 17.1145433j,
+                5.6628990 + 8.8759694j]]
+
+        # Correct versions of the call
+        np.testing.assert_almost_equal(evalfr(self.sysC322, 1j), resp)
+        np.testing.assert_almost_equal(self.sysC322._evalfr(1.), resp)
+
+        # Deprecated version of the call (should generate warning)
+        import warnings
+        with warnings.catch_warnings(record=True) as w:
+            # Set up warnings filter to only show warnings in control module
+            warnings.filterwarnings("ignore")
+            warnings.filterwarnings("always", module="control")
+
+            # Make sure that we get a pending deprecation warning
+            self.sysC322.evalfr(1.)
+            assert len(w) == 1
+            assert issubclass(w[-1].category, PendingDeprecationWarning)
+
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def test_freq_resp_complex(self):
+        """Evaluate the frequency response at multiple frequencies."""
+
+        true_mag = [[15.2721178, 3.9176691, 20.5865790],
+        [24.5384389, 2.8374330, 18.2268344]]
+
+        true_phase = [[1.03455334, 2.2133291, 2.6715324],
+        [1.8217663, -2.8266936, 2.2694910]]
+        true_omega = [0.1, 10, 0.01j, -1j];
+
+        mag, phase, omega = self.sysC623.freqresp(true_omega)
+
+        np.testing.assert_almost_equal(mag, true_mag)
+        np.testing.assert_almost_equal(phase, true_phase)
+        np.testing.assert_equal(omega, true_omega)
+
+
 class TestRss(unittest.TestCase):
     """These are tests for the proper functionality of statesp.rss."""