python-control · bnavigator · Aug 17, 2020 · Jul 22, 2020 · Jul 23, 2020 · Jul 23, 2020
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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])