diff --git a/.travis.yml b/.travis.yml index 022e48c6f..ec615501d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -49,6 +49,11 @@ jobs: services: xvfb python: "3.8" env: SCIPY=scipy SLYCOT=source + - name: "use numpy matrix" + dist: xenial + services: xvfb + python: "3.8" + env: SCIPY=scipy SLYCOT=source PYTHON_CONTROL_STATESPACE_ARRAY=1 # Exclude combinations that are very unlikely (and don't work) exclude: diff --git a/control/canonical.py b/control/canonical.py index b578418bd..bd9ee4a94 100644 --- a/control/canonical.py +++ b/control/canonical.py @@ -6,7 +6,8 @@ from .statesp import StateSpace from .statefbk import ctrb, obsv -from numpy import zeros, shape, poly, iscomplex, hstack, dot, transpose +from numpy import zeros, zeros_like, shape, poly, iscomplex, vstack, hstack, dot, \ + transpose, empty from numpy.linalg import solve, matrix_rank, eig __all__ = ['canonical_form', 'reachable_form', 'observable_form', 'modal_form', @@ -70,9 +71,9 @@ def reachable_form(xsys): zsys = StateSpace(xsys) # Generate the system matrices for the desired canonical form - zsys.B = zeros(shape(xsys.B)) + zsys.B = zeros_like(xsys.B) zsys.B[0, 0] = 1.0 - zsys.A = zeros(shape(xsys.A)) + zsys.A = zeros_like(xsys.A) Apoly = poly(xsys.A) # characteristic polynomial for i in range(0, xsys.states): zsys.A[0, i] = -Apoly[i+1] / Apoly[0] @@ -124,9 +125,9 @@ def observable_form(xsys): zsys = StateSpace(xsys) # Generate the system matrices for the desired canonical form - zsys.C = zeros(shape(xsys.C)) + zsys.C = zeros_like(xsys.C) zsys.C[0, 0] = 1 - zsys.A = zeros(shape(xsys.A)) + zsys.A = zeros_like(xsys.A) Apoly = poly(xsys.A) # characteristic polynomial for i in range(0, xsys.states): zsys.A[i, 0] = -Apoly[i+1] / Apoly[0] @@ -144,7 +145,7 @@ def observable_form(xsys): raise ValueError("Transformation matrix singular to working precision.") # Finally, compute the output matrix - zsys.B = Tzx * xsys.B + zsys.B = Tzx.dot(xsys.B) return zsys, Tzx @@ -174,9 +175,9 @@ def modal_form(xsys): # Calculate eigenvalues and matrix of eigenvectors Tzx, eigval, eigvec = eig(xsys.A) - # Eigenvalues and according eigenvectors are not sorted, + # Eigenvalues and corresponding eigenvectors are not sorted, # thus modal transformation is ambiguous - # Sorting eigenvalues and respective vectors by largest to smallest eigenvalue + # Sort eigenvalues and vectors from largest to smallest eigenvalue idx = eigval.argsort()[::-1] eigval = eigval[idx] eigvec = eigvec[:,idx] @@ -189,23 +190,18 @@ def modal_form(xsys): # Keep track of complex conjugates (need only one) lst_conjugates = [] - Tzx = None + Tzx = empty((0, xsys.A.shape[0])) # empty zero-height row matrix for val, vec in zip(eigval, eigvec.T): if iscomplex(val): if val not in lst_conjugates: lst_conjugates.append(val.conjugate()) - if Tzx is not None: - Tzx = hstack((Tzx, hstack((vec.real.T, vec.imag.T)))) - else: - Tzx = hstack((vec.real.T, vec.imag.T)) + Tzx = vstack((Tzx, vec.real, vec.imag)) else: # if conjugate has already been seen, skip this eigenvalue lst_conjugates.remove(val) else: - if Tzx is not None: - Tzx = hstack((Tzx, vec.real.T)) - else: - Tzx = vec.real.T + Tzx = vstack((Tzx, vec.real)) + Tzx = Tzx.T # Generate the system matrices for the desired canonical form zsys.A = solve(Tzx, xsys.A).dot(Tzx) diff --git a/control/statesp.py b/control/statesp.py index 522d187a9..526fc129b 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -77,7 +77,10 @@ def _ssmatrix(data, axis=1): - """Convert argument to a (possibly empty) state space matrix. + """Convert argument to a (possibly empty) 2D state space matrix. + + The axis keyword argument makes it convenient to specify that if the input + is a vector, it is a row (axis=1) or column (axis=0) vector. Parameters ---------- @@ -94,8 +97,10 @@ 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): + if config.defaults['statesp.use_numpy_matrix']: arr = np.matrix(data, dtype=float) + elif isinstance(data, str): + arr = np.array(np.matrix(data, dtype=float)) else: arr = np.array(data, dtype=float) ndim = arr.ndim @@ -195,12 +200,20 @@ def __init__(self, *args, **kw): raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) # Process keyword arguments - remove_useless = kw.get('remove_useless', config.defaults['statesp.remove_useless_states']) + remove_useless = kw.get('remove_useless', + config.defaults['statesp.remove_useless_states']) # Convert all matrices to standard form A = _ssmatrix(A) - B = _ssmatrix(B, axis=0) - C = _ssmatrix(C, axis=1) + # if B is a 1D array, turn it into a column vector if it fits + if np.asarray(B).ndim == 1 and len(B) == A.shape[0]: + B = _ssmatrix(B, axis=0) + else: + B = _ssmatrix(B) + if np.asarray(C).ndim == 1 and len(C) == A.shape[0]: + C = _ssmatrix(C, axis=1) + else: + C = _ssmatrix(C, axis=0) #if this doesn't work, error below if np.isscalar(D) and D == 0 and B.shape[1] > 0 and C.shape[0] > 0: # If D is a scalar zero, broadcast it to the proper size D = np.zeros((C.shape[0], B.shape[1])) @@ -1240,8 +1253,8 @@ def _mimo2simo(sys, input, warn_conversion=False): "Only input {i} is used." .format(i=input)) # $X = A*X + B*U # Y = C*X + D*U - new_B = sys.B[:, input] - new_D = sys.D[:, input] + new_B = sys.B[:, input:input+1] + new_D = sys.D[:, input:input+1] sys = StateSpace(sys.A, new_B, sys.C, new_D, sys.dt) return sys diff --git a/control/tests/canonical_test.py b/control/tests/canonical_test.py index 3172f13b7..7d4ae4e27 100644 --- a/control/tests/canonical_test.py +++ b/control/tests/canonical_test.py @@ -22,13 +22,13 @@ def test_reachable_form(self): D_true = 42.0 # Perform a coordinate transform with a random invertible matrix - T_true = np.matrix([[-0.27144004, -0.39933167, 0.75634684, 0.44135471], + T_true = np.array([[-0.27144004, -0.39933167, 0.75634684, 0.44135471], [-0.74855725, -0.39136285, -0.18142339, -0.50356997], [-0.40688007, 0.81416369, 0.38002113, -0.16483334], [-0.44769516, 0.15654653, -0.50060858, 0.72419146]]) - A = np.linalg.solve(T_true, A_true)*T_true + A = np.linalg.solve(T_true, A_true).dot(T_true) B = np.linalg.solve(T_true, B_true) - C = C_true*T_true + C = C_true.dot(T_true) D = D_true # Create a state space system and convert it to the reachable canonical form @@ -69,11 +69,11 @@ def test_modal_form(self): D_true = 42.0 # Perform a coordinate transform with a random invertible matrix - T_true = np.matrix([[-0.27144004, -0.39933167, 0.75634684, 0.44135471], + T_true = np.array([[-0.27144004, -0.39933167, 0.75634684, 0.44135471], [-0.74855725, -0.39136285, -0.18142339, -0.50356997], [-0.40688007, 0.81416369, 0.38002113, -0.16483334], [-0.44769516, 0.15654653, -0.50060858, 0.72419146]]) - A = np.linalg.solve(T_true, A_true)*T_true + A = np.linalg.solve(T_true, A_true).dot(T_true) B = np.linalg.solve(T_true, B_true) C = C_true*T_true D = D_true @@ -98,9 +98,9 @@ def test_modal_form(self): C_true = np.array([[1, 0, 0, 1]]) D_true = np.array([[0]]) - A = np.linalg.solve(T_true, A_true) * T_true + A = np.linalg.solve(T_true, A_true).dot(T_true) B = np.linalg.solve(T_true, B_true) - C = C_true * T_true + C = C_true.dot(T_true) D = D_true # Create state space system and convert to modal canonical form @@ -132,9 +132,9 @@ def test_modal_form(self): C_true = np.array([[0, 1, 0, 1]]) D_true = np.array([[0]]) - A = np.linalg.solve(T_true, A_true) * T_true + A = np.linalg.solve(T_true, A_true).dot(T_true) B = np.linalg.solve(T_true, B_true) - C = C_true * T_true + C = C_true.dot(T_true) D = D_true # Create state space system and convert to modal canonical form @@ -173,13 +173,13 @@ def test_observable_form(self): D_true = 42.0 # Perform a coordinate transform with a random invertible matrix - T_true = np.matrix([[-0.27144004, -0.39933167, 0.75634684, 0.44135471], + T_true = np.array([[-0.27144004, -0.39933167, 0.75634684, 0.44135471], [-0.74855725, -0.39136285, -0.18142339, -0.50356997], [-0.40688007, 0.81416369, 0.38002113, -0.16483334], [-0.44769516, 0.15654653, -0.50060858, 0.72419146]]) - A = np.linalg.solve(T_true, A_true)*T_true + A = np.linalg.solve(T_true, A_true).dot(T_true) B = np.linalg.solve(T_true, B_true) - C = C_true*T_true + C = C_true.dot(T_true) D = D_true # Create a state space system and convert it to the observable canonical form diff --git a/control/tests/conftest.py b/control/tests/conftest.py new file mode 100755 index 000000000..e98bbe1d7 --- /dev/null +++ b/control/tests/conftest.py @@ -0,0 +1,13 @@ +# contest.py - pytest local plugins and fixtures + +import control +import os + +import pytest + + +@pytest.fixture(scope="session", autouse=True) +def use_numpy_ndarray(): + """Switch the config to use ndarray instead of matrix""" + if os.getenv("PYTHON_CONTROL_STATESPACE_ARRAY") == "1": + control.config.defaults['statesp.use_numpy_matrix'] = False diff --git a/control/tests/discrete_test.py b/control/tests/discrete_test.py index 6598e3a81..9c1928dab 100644 --- a/control/tests/discrete_test.py +++ b/control/tests/discrete_test.py @@ -353,7 +353,7 @@ def test_sample_ss(self): for sys in (sys1, sys2): for h in (0.1, 0.5, 1, 2): Ad = I + h * sys.A - Bd = h * sys.B + 0.5 * h**2 * (sys.A * sys.B) + Bd = h * sys.B + 0.5 * h**2 * np.dot(sys.A, sys.B) sysd = sample_system(sys, h, method='zoh') np.testing.assert_array_almost_equal(sysd.A, Ad) np.testing.assert_array_almost_equal(sysd.B, Bd) diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py index 0738e8b18..22f8307d2 100644 --- a/control/tests/iosys_test.py +++ b/control/tests/iosys_test.py @@ -53,7 +53,7 @@ def test_linear_iosys(self): for x, u in (([0, 0], 0), ([1, 0], 0), ([0, 1], 0), ([0, 0], 1)): np.testing.assert_array_almost_equal( np.reshape(iosys._rhs(0, x, u), (-1,1)), - linsys.A * np.reshape(x, (-1, 1)) + linsys.B * u) + np.dot(linsys.A, np.reshape(x, (-1, 1))) + np.dot(linsys.B, u)) # Make sure that simulations also line up T, U, X0 = self.T, self.U, self.X0 @@ -151,9 +151,9 @@ def test_nonlinear_iosys(self): # Create a nonlinear system with the same dynamics nlupd = lambda t, x, u, params: \ - np.reshape(linsys.A * np.reshape(x, (-1, 1)) + linsys.B * u, (-1,)) + np.reshape(np.dot(linsys.A, np.reshape(x, (-1, 1))) + np.dot(linsys.B, u), (-1,)) nlout = lambda t, x, u, params: \ - np.reshape(linsys.C * np.reshape(x, (-1, 1)) + linsys.D * u, (-1,)) + np.reshape(np.dot(linsys.C, np.reshape(x, (-1, 1))) + np.dot(linsys.D, u), (-1,)) nlsys = ios.NonlinearIOSystem(nlupd, nlout) # Make sure that simulations also line up @@ -747,8 +747,8 @@ def test_named_signals(self): + np.dot(self.mimo_linsys1.B, np.reshape(u, (-1, 1))) ).reshape(-1,), outfcn = lambda t, x, u, params: np.array( - self.mimo_linsys1.C * np.reshape(x, (-1, 1)) \ - + self.mimo_linsys1.D * np.reshape(u, (-1, 1)) + np.dot(self.mimo_linsys1.C, np.reshape(x, (-1, 1))) \ + + np.dot(self.mimo_linsys1.D, np.reshape(u, (-1, 1))) ).reshape(-1,), inputs = ('u[0]', 'u[1]'), outputs = ('y[0]', 'y[1]'), diff --git a/control/tests/statesp_array_test.py b/control/tests/statesp_array_test.py index a2d034075..f0574cf24 100644 --- a/control/tests/statesp_array_test.py +++ b/control/tests/statesp_array_test.py @@ -13,6 +13,7 @@ from control.lti import evalfr from control.exception import slycot_check from control.config import use_numpy_matrix, reset_defaults +from control.config import defaults class TestStateSpace(unittest.TestCase): """Tests for the StateSpace class.""" @@ -74,8 +75,12 @@ def test_matlab_style_constructor(self): self.assertEqual(sys.B.shape, (2, 1)) self.assertEqual(sys.C.shape, (1, 2)) self.assertEqual(sys.D.shape, (1, 1)) - for X in [sys.A, sys.B, sys.C, sys.D]: - self.assertTrue(isinstance(X, np.matrix)) + if defaults['statesp.use_numpy_matrix']: + for X in [sys.A, sys.B, sys.C, sys.D]: + self.assertTrue(isinstance(X, np.matrix)) + else: + for X in [sys.A, sys.B, sys.C, sys.D]: + self.assertTrue(isinstance(X, np.ndarray)) def test_pole(self): """Evaluate the poles of a MIMO system.""" diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py index 96404d79f..34a17f992 100644 --- a/control/tests/statesp_test.py +++ b/control/tests/statesp_test.py @@ -323,9 +323,9 @@ def test_array_access_ss(self): np.testing.assert_array_almost_equal(sys1_11.A, sys1.A) np.testing.assert_array_almost_equal(sys1_11.B, - sys1.B[:, 1]) + sys1.B[:, 1:2]) np.testing.assert_array_almost_equal(sys1_11.C, - sys1.C[0, :]) + sys1.C[0:1, :]) np.testing.assert_array_almost_equal(sys1_11.D, sys1.D[0, 1]) diff --git a/control/tests/test_control_matlab.py b/control/tests/test_control_matlab.py index e45b52523..aa8633e7c 100644 --- a/control/tests/test_control_matlab.py +++ b/control/tests/test_control_matlab.py @@ -11,7 +11,7 @@ from numpy.testing import assert_array_almost_equal from numpy import array, asarray, matrix, asmatrix, zeros, ones, linspace,\ all, hstack, vstack, c_, r_ -from matplotlib.pylab import show, figure, plot, legend, subplot2grid +from matplotlib.pyplot import show, figure, plot, legend, subplot2grid from control.matlab import ss, step, impulse, initial, lsim, dcgain, \ ss2tf from control.statesp import _mimo2siso @@ -24,29 +24,13 @@ class TestControlMatlab(unittest.TestCase): def setUp(self): pass - def plot_matrix(self): - #Test: can matplotlib correctly plot matrices? - #Yes, but slightly inconvenient - figure() - t = matrix([[ 1.], - [ 2.], - [ 3.], - [ 4.]]) - y = matrix([[ 1., 4.], - [ 4., 5.], - [ 9., 6.], - [16., 7.]]) - plot(t, y) - #plot(asarray(t)[0], asarray(y)[0]) - - def make_SISO_mats(self): """Return matrices for a SISO system""" - A = matrix([[-81.82, -45.45], + A = array([[-81.82, -45.45], [ 10., -1. ]]) - B = matrix([[9.09], + B = array([[9.09], [0. ]]) - C = matrix([[0, 0.159]]) + C = array([[0, 0.159]]) D = zeros((1, 1)) return A, B, C, D @@ -181,7 +165,7 @@ def test_impulse(self): #Test MIMO system A, B, C, D = self.make_MIMO_mats() - sys = ss(A, B, C, D) + sys = ss(A, B, C, D) t, y = impulse(sys) plot(t, y, label='MIMO System') @@ -202,7 +186,7 @@ def test_initial(self): #X0=[1,1] : produces a spike subplot2grid(plot_shape, (0, 1)) - t, y = initial(sys, X0=matrix("1; 1")) + t, y = initial(sys, X0=array(matrix("1; 1"))) plot(t, y) #Test MIMO system @@ -318,21 +302,11 @@ def test_lsim(self): plot(t, y, label='y') legend(loc='best') - #Test with matrices - subplot2grid(plot_shape, (1, 0)) - t = matrix(linspace(0, 1, 100)) - u = matrix(r_[1:1:50j, 0:0:50j]) - x0 = matrix("0.; 0") - y, t_out, _x = lsim(sys, u, t, x0) - plot(t_out, y, label='y') - plot(t_out, asarray(u/10)[0], label='u/10') - legend(loc='best') - #Test with MIMO system subplot2grid(plot_shape, (1, 1)) A, B, C, D = self.make_MIMO_mats() sys = ss(A, B, C, D) - t = matrix(linspace(0, 1, 100)) + t = array(linspace(0, 1, 100)) u = array([r_[1:1:50j, 0:0:50j], r_[0:1:50j, 0:0:50j]]) x0 = [0, 0, 0, 0] @@ -404,12 +378,12 @@ def test_convert_MIMO_to_SISO(self): #Test with additional systems -------------------------------------------- #They have crossed inputs and direct feedthrough #SISO system - As = matrix([[-81.82, -45.45], + As = array([[-81.82, -45.45], [ 10., -1. ]]) - Bs = matrix([[9.09], + Bs = array([[9.09], [0. ]]) - Cs = matrix([[0, 0.159]]) - Ds = matrix([[0.02]]) + Cs = array([[0, 0.159]]) + Ds = array([[0.02]]) sys_siso = ss(As, Bs, Cs, Ds) # t, y = step(sys_siso) # plot(t, y, label='sys_siso d=0.02') @@ -428,7 +402,7 @@ def test_convert_MIMO_to_SISO(self): [0 , 0 ]]) Cm = array([[0, 0, 0, 0.159], [0, 0.159, 0, 0 ]]) - Dm = matrix([[0, 0.02], + Dm = array([[0, 0.02], [0.02, 0 ]]) sys_mimo = ss(Am, Bm, Cm, Dm)