diff --git a/control/statesp.py b/control/statesp.py index 1779dfbfd..405cdc62f 100644 --- 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 index 9273877af..8c94e7d76 100644 --- 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."""