Repository: python-control/python-control
Pull request #913 (2023-06-10 17:18:14)

Remove NumPy matrix class
This PR removes the use of the NumPy matrix class from python-control.  The main changes are:

* The `use_numpy_matrix` function is no longer available.
* The use of MATLAB-style matrix strings ("[[1 2; 3 4]]") is no longer supported.
* The `use_legacy_defaults` function returns a warning if the version is pre-0.9 (since those versions used the `np.matrix` type for state space systems).
* Unit tests and CI tests were updated to remove testing against `np.matrix`.
* Docstrings and documentation were updated to be compatible with the above changes.


Patch:
diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml
index cea5e542..04b46a46 100644
--- a/.github/workflows/python-package-conda.yml
+++ b/.github/workflows/python-package-conda.yml
@@ -9,7 +9,6 @@ jobs:
       ${{ matrix.slycot || 'no' }} Slycot;
       ${{ matrix.pandas || 'no' }} Pandas;
       ${{ matrix.cvxopt || 'no' }} CVXOPT
-      ${{ matrix.array-and-matrix == 1 && '; array and matrix' || '' }}
       ${{ matrix.mplbackend && format('; {0}', matrix.mplbackend) }}
     runs-on: ubuntu-latest
 
@@ -22,14 +21,12 @@ jobs:
         pandas: [""]
         cvxopt: ["", "conda"]
         mplbackend: [""]
-        array-and-matrix: [0]
         include:
           - python-version: '3.11'
             slycot: conda
             pandas: conda
             cvxopt: conda
             mplbackend: QtAgg
-            array-and-matrix: 1
 
     steps:
     - uses: actions/checkout@v3
@@ -63,7 +60,6 @@ jobs:
     - name: Test with pytest
       shell: bash -l {0}
       env:
-        PYTHON_CONTROL_ARRAY_AND_MATRIX: ${{ matrix.array-and-matrix }}
         MPLBACKEND: ${{ matrix.mplbackend }}
       run: pytest -v --cov=control --cov-config=.coveragerc control/tests
 
diff --git a/control/config.py b/control/config.py
index f75bd52d..9009a30f 100644
--- a/control/config.py
+++ b/control/config.py
@@ -14,7 +14,7 @@ from .exception import ControlArgument
 
 __all__ = ['defaults', 'set_defaults', 'reset_defaults',
            'use_matlab_defaults', 'use_fbs_defaults',
-           'use_legacy_defaults', 'use_numpy_matrix']
+           'use_legacy_defaults']
 
 # Package level default values
 _control_defaults = {
@@ -202,7 +202,6 @@ def use_matlab_defaults():
     The following conventions are used:
         * Bode plots plot gain in dB, phase in degrees, frequency in
           rad/sec, with grids
-        * State space class and functions use Numpy matrix objects
 
     Examples
     --------
@@ -211,7 +210,6 @@ def use_matlab_defaults():
 
     """
     set_defaults('freqplot', dB=True, deg=True, Hz=False, grid=True)
-    set_defaults('statesp', use_numpy_matrix=True)
 
 
 # Set defaults to match FBS (Astrom and Murray)
@@ -233,41 +231,6 @@ def use_fbs_defaults():
     set_defaults('nyquist', mirror_style='--')
 
 
-# Decide whether to use numpy.matrix for state space operations
-def use_numpy_matrix(flag=True, warn=True):
-    """Turn on/off use of Numpy `matrix` class for state space operations.
-
-    Parameters
-    ----------
-    flag : bool
-        If flag is `True` (default), use the deprecated Numpy
-        `matrix` class to represent matrices in the `~control.StateSpace`
-        class and functions.  If flat is `False`, then matrices are
-        represented by a 2D `ndarray` object.
-
-    warn : bool
-        If flag is `True` (default), issue a warning when turning on the use
-        of the Numpy `matrix` class.  Set `warn` to false to omit display of
-        the warning message.
-
-    Notes
-    -----
-    Prior to release 0.9.x, the default type for 2D arrays is the Numpy
-    `matrix` class.  Starting in release 0.9.0, the default type for state
-    space operations is a 2D array.
-
-    Examples
-    --------
-    >>> ct.use_numpy_matrix(True, False)
-    >>> # do some legacy calculations using np.matrix
-
-    """
-    if flag and warn:
-        warnings.warn("Return type numpy.matrix is deprecated.",
-                      stacklevel=2, category=DeprecationWarning)
-    set_defaults('statesp', use_numpy_matrix=flag)
-
-
 def use_legacy_defaults(version):
     """ Sets the defaults to whatever they were in a given release.
 
@@ -331,7 +294,7 @@ def use_legacy_defaults(version):
     # Version 0.9.0:
     if major == 0 and minor < 9:
         # switched to 'array' as default for state space objects
-        set_defaults('statesp', use_numpy_matrix=True)
+        warnings.warn("NumPy matrix class no longer supported")
 
         # switched to 0 (=continuous) as default timestep
         set_defaults('control', default_dt=None)
diff --git a/control/iosys.py b/control/iosys.py
index 6c6458f3..00ed288f 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -2227,8 +2227,6 @@ def ss(*args, **kwargs):
               y[k] &= C x[k] + D u[k]
 
         The matrices can be given as *array like* data types or strings.
-        Everything that the constructor of :class:`numpy.matrix` accepts is
-        permissible here too.
 
     ``ss(args, inputs=['u1', ..., 'up'], outputs=['y1', ..., 'yq'], states=['x1', ..., 'xn'])``
         Create a system with named input, output, and state signals.
diff --git a/control/mateqn.py b/control/mateqn.py
index 1cf2e65d..339f1a88 100644
--- a/control/mateqn.py
+++ b/control/mateqn.py
@@ -126,14 +126,9 @@ def lyap(A, Q, C=None, E=None, method=None):
 
     Returns
     -------
-    X : 2D array (or matrix)
+    X : 2D array
         Solution to the Lyapunov or Sylvester equation
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     """
     # Decide what method to use
     method = _slycot_or_scipy(method)
@@ -260,11 +255,6 @@ def dlyap(A, Q, C=None, E=None, method=None):
     X : 2D array (or matrix)
         Solution to the Lyapunov or Sylvester equation
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     """
     # Decide what method to use
     method = _slycot_or_scipy(method)
@@ -395,11 +385,6 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True, method=None,
     G : 2D array (or matrix)
         Gain matrix
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     """
     # Decide what method to use
     method = _slycot_or_scipy(method)
@@ -554,11 +539,6 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True, method=None,
     G : 2D array (or matrix)
         Gain matrix
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     """
     # Decide what method to use
     method = _slycot_or_scipy(method)
diff --git a/control/matlab/wrappers.py b/control/matlab/wrappers.py
index e7d75724..d98dcabf 100644
--- a/control/matlab/wrappers.py
+++ b/control/matlab/wrappers.py
@@ -48,7 +48,7 @@ def bode(*args, **kwargs):
     --------
     >>> from control.matlab import ss, bode
 
-    >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
+    >>> sys = ss([[1, -2], [3, -4]], [[5], [7]], [[6, 8]], 9)
     >>> mag, phase, omega = bode(sys)
 
     .. todo::
diff --git a/control/statefbk.py b/control/statefbk.py
index f9897419..50bad75c 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -110,9 +110,6 @@ def place(A, B, p):
         The algorithm will not place poles at the same location more
         than rank(B) times.
 
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     References
     ----------
     .. [1] A.L. Tits and Y. Yang, "Globally convergent algorithms for robust
@@ -193,11 +190,6 @@ def place_varga(A, B, p, dtime=False, alpha=None):
     [1] Varga A. "A Schur method for pole assignment."  IEEE Trans. Automatic
         Control, Vol. AC-26, pp. 517-519, 1981.
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     Examples
     --------
     >>> A = [[-1, -1], [0, 1]]
@@ -279,10 +271,6 @@ def acker(A, B, poles):
     K : 2D array (or matrix)
         Gains such that A - B K has given eigenvalues
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
     """
     # Convert the inputs to matrices
     a = _ssmatrix(A)
@@ -366,13 +354,10 @@ def lqr(*args, **kwargs):
 
     Notes
     -----
-    1. If the first argument is an LTI object, then this object will be used
-       to define the dynamics and input matrices.  Furthermore, if the LTI
-       object corresponds to a discrete time system, the ``dlqr()`` function
-       will be called.
-
-    2. The return type for 2D arrays depends on the default class set for
-       state space operations.  See :func:`~control.use_numpy_matrix`.
+    If the first argument is an LTI object, then this object will be used
+    to define the dynamics and input matrices.  Furthermore, if the LTI
+    object corresponds to a discrete time system, the ``dlqr()`` function
+    will be called.
 
     Examples
     --------
@@ -514,11 +499,6 @@ def dlqr(*args, **kwargs):
     --------
     lqr, lqe, dlqe
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     Examples
     --------
     >>> K, S, E = dlqr(dsys, Q, R, [N])                         # doctest: +SKIP
@@ -971,11 +951,6 @@ def ctrb(A, B):
     C : 2D array (or matrix)
         Controllability matrix
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     Examples
     --------
     >>> G = ct.tf2ss([1], [1, 2, 3])
@@ -1010,11 +985,6 @@ def obsv(A, C):
     O : 2D array (or matrix)
         Observability matrix
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     Examples
     --------
     >>> G = ct.tf2ss([1], [1, 2, 3])
@@ -1063,11 +1033,6 @@ def gram(sys, type):
         if slycot routine sb03md cannot be found
         if slycot routine sb03od cannot be found
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     Examples
     --------
     >>> G = ct.rss(4)
diff --git a/control/statesp.py b/control/statesp.py
index b2a3934a..ae2c32c5 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -77,7 +77,6 @@ __all__ = ['StateSpace', 'tf2ss', 'ssdata', 'linfnorm']
 
 # Define module default parameter values
 _statesp_defaults = {
-    'statesp.use_numpy_matrix': False,  # False is default in 0.9.0 and above
     'statesp.remove_useless_states': False,
     'statesp.latex_num_format': '.3g',
     'statesp.latex_repr_type': 'partitioned',
@@ -104,14 +103,8 @@ def _ssmatrix(data, axis=1):
     arr : 2D array, with shape (0, 0) if a is empty
 
     """
-    # 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']:
-        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)
+    # Convert the data into an array
+    arr = np.array(data, dtype=float)
     ndim = arr.ndim
     shape = arr.shape
 
@@ -205,12 +198,7 @@ class StateSpace(LTI):
     -----
     The main data members in the ``StateSpace`` class are the A, B, C, and D
     matrices.  The class also keeps track of the number of states (i.e.,
-    the size of A).  The data format used to store state space matrices is
-    set using the value of `config.defaults['use_numpy_matrix']`.  If True
-    (default), the state space elements are stored as `numpy.matrix` objects;
-    otherwise they are `numpy.ndarray` objects.  The
-    :func:`~control.use_numpy_matrix` function can be used to set the storage
-    type.
+    the size of A).
 
     A discrete time system is created by specifying a nonzero 'timebase', dt
     when the system is constructed:
@@ -358,10 +346,8 @@ class StateSpace(LTI):
         elif kwargs:
             raise TypeError("unrecognized keyword(s): ", str(kwargs))
 
-        # Reset shapes (may not be needed once np.matrix support is removed)
+        # Reset shape if system is static
         if self._isstatic():
-            # static gain
-            # matrix's default "empty" shape is 1x0
             A.shape = (0, 0)
             B.shape = (0, self.ninputs)
             C.shape = (self.noutputs, 0)
@@ -467,10 +453,6 @@ class StateSpace(LTI):
         """
 
         # Search for useless states and get indices of these states.
-        #
-        # Note: shape from np.where depends on whether we are storing state
-        # space objects as np.matrix or np.array.  Code below will work
-        # correctly in either case.
         ax1_A = np.where(~self.A.any(axis=1))[0]
         ax1_B = np.where(~self.B.any(axis=1))[0]
         ax0_A = np.where(~self.A.any(axis=0))[-1]
@@ -502,12 +484,11 @@ class StateSpace(LTI):
         return string
 
     # represent to implement a re-loadable version
-    # TODO: remove the conversion to array when matrix is no longer used
     def __repr__(self):
         """Print state-space system in loadable form."""
         return "StateSpace({A}, {B}, {C}, {D}{dt})".format(
-            A=asarray(self.A).__repr__(), B=asarray(self.B).__repr__(),
-            C=asarray(self.C).__repr__(), D=asarray(self.D).__repr__(),
+            A=self.A.__repr__(), B=self.B.__repr__(),
+            C=self.C.__repr__(), D=self.D.__repr__(),
             dt=(isdtime(self, strict=True) and ", {}".format(self.dt)) or '')
 
     def _latex_partitioned_stateless(self):
@@ -930,18 +911,17 @@ class StateSpace(LTI):
         x_arr = np.atleast_1d(x).astype(complex, copy=False)
 
         # return fast on systems with 0 or 1 state
-        if not config.defaults['statesp.use_numpy_matrix']:
-            if self.nstates == 0:
-                return self.D[:, :, np.newaxis] \
-                    * np.ones_like(x_arr, dtype=complex)
-            if self.nstates == 1:
-                with np.errstate(divide='ignore', invalid='ignore'):
-                    out = self.C[:, :, np.newaxis] \
-                          / (x_arr - self.A[0, 0]) \
-                          * self.B[:, :, np.newaxis] \
-                          + self.D[:, :, np.newaxis]
-                out[np.isnan(out)] = complex(np.inf, np.nan)
-                return out
+        if self.nstates == 0:
+            return self.D[:, :, np.newaxis] \
+                * np.ones_like(x_arr, dtype=complex)
+        elif self.nstates == 1:
+            with np.errstate(divide='ignore', invalid='ignore'):
+                out = self.C[:, :, np.newaxis] \
+                    / (x_arr - self.A[0, 0]) \
+                    * self.B[:, :, np.newaxis] \
+                    + self.D[:, :, np.newaxis]
+            out[np.isnan(out)] = complex(np.inf, np.nan)
+            return out
 
         try:
             out = self.slycot_laub(x_arr)
diff --git a/control/stochsys.py b/control/stochsys.py
index 663b09ec..85e10833 100644
--- a/control/stochsys.py
+++ b/control/stochsys.py
@@ -87,9 +87,9 @@ def lqe(*args, **kwargs):
 
     Returns
     -------
-    L : 2D array (or matrix)
+    L : 2D array
         Kalman estimator gain
-    P : 2D array (or matrix)
+    P : 2D array
         Solution to Riccati equation
 
         .. math::
@@ -101,13 +101,10 @@ def lqe(*args, **kwargs):
 
     Notes
     -----
-    1. If the first argument is an LTI object, then this object will be used
-       to define the dynamics, noise and output matrices.  Furthermore, if
-       the LTI object corresponds to a discrete time system, the ``dlqe()``
-       function will be called.
-
-    2. The return type for 2D arrays depends on the default class set for
-       state space operations.  See :func:`~control.use_numpy_matrix`.
+    If the first argument is an LTI object, then this object will be used
+    to define the dynamics, noise and output matrices.  Furthermore, if the
+    LTI object corresponds to a discrete time system, the ``dlqe()``
+    function will be called.
 
     Examples
     --------
@@ -224,9 +221,9 @@ def dlqe(*args, **kwargs):
 
     Returns
     -------
-    L : 2D array (or matrix)
+    L : 2D array
         Kalman estimator gain
-    P : 2D array (or matrix)
+    P : 2D array
         Solution to Riccati equation
 
         .. math::
@@ -236,11 +233,6 @@ def dlqe(*args, **kwargs):
     E : 1D array
         Eigenvalues of estimator poles eig(A - L C)
 
-    Notes
-    -----
-    The return type for 2D arrays depends on the default class set for
-    state space operations.  See :func:`~control.use_numpy_matrix`.
-
     Examples
     --------
     >>> L, P, E = dlqe(A, G, C, QN, RN)                         # doctest: +SKIP
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 56006d69..7303d5e7 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -1640,8 +1640,8 @@ def tf(*args, **kwargs):
     >>> G  = (s + 1)/(s**2 + 2*s + 1)
 
     >>> # Convert a StateSpace to a TransferFunction object.
-    >>> sys_ss = ct.ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
-    >>> sys2 = ct.tf(sys1)
+    >>> sys_ss = ct.ss([[1, -2], [3, -4]], [[5], [7]], [[6, 8]], 9)
+    >>> sys_tf = ct.tf(sys_ss)
 
     """
 
@@ -1801,7 +1801,7 @@ def ss2tf(*args, **kwargs):
     >>> sys1 = ct.ss2tf(A, B, C, D)
 
     >>> sys_ss = ct.ss(A, B, C, D)
-    >>> sys2 = ct.ss2tf(sys_ss)
+    >>> sys_tf = ct.ss2tf(sys_ss)
 
     """
 
diff --git a/doc/control.rst b/doc/control.rst
index 8dc8a09a..ca46043d 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -197,6 +197,5 @@ Utility functions and conversions
     unwrap
     use_fbs_defaults
     use_matlab_defaults
-    use_numpy_matrix
 
 
diff --git a/doc/conventions.rst b/doc/conventions.rst
index 7c9c1ec6..b5073b8e 100644
--- a/doc/conventions.rst
+++ b/doc/conventions.rst
@@ -303,9 +303,6 @@ Selected variables that can be configured, along with their default values:
   * freqplot.feature_periphery_decade (1.0): How many decades to include in
     the frequency range on both sides of features (poles, zeros).
 
-  * statesp.use_numpy_matrix (True): set the return type for state space
-    matrices to `numpy.matrix` (verus numpy.ndarray)
-
   * statesp.default_dt and xferfcn.default_dt (None): set the default value
     of dt when constructing new LTI systems
 
@@ -322,5 +319,4 @@ Functions that can be used to set standard configurations:
     reset_defaults
     use_fbs_defaults
     use_matlab_defaults
-    use_numpy_matrix
     use_legacy_defaults


Test patch:
diff --git a/.github/workflows/doctest.yml b/.github/workflows/doctest.yml
index 62638d10..49455a5c 100644
--- a/.github/workflows/doctest.yml
+++ b/.github/workflows/doctest.yml
@@ -33,7 +33,6 @@ jobs:
     - name: Run doctest
       shell: bash -l {0}
       env:
-        PYTHON_CONTROL_ARRAY_AND_MATRIX: ${{ matrix.array-and-matrix }}
         MPLBACKEND: ${{ matrix.mplbackend }}
       working-directory: doc
       run: |
diff --git a/control/tests/config_test.py b/control/tests/config_test.py
index 15229139..1547b7e2 100644
--- a/control/tests/config_test.py
+++ b/control/tests/config_test.py
@@ -242,15 +242,13 @@ class TestConfig:
         assert ct.config.defaults['freqplot.feature_periphery_decades'] == 1.0
 
     def test_legacy_defaults(self):
-        with pytest.deprecated_call():
+        with pytest.warns(UserWarning, match="NumPy matrix class no longer"):
             ct.use_legacy_defaults('0.8.3')
-            assert(isinstance(ct.ss(0, 0, 0, 1).D, np.matrix))
-        ct.reset_defaults()
-        assert isinstance(ct.ss(0, 0, 0, 1).D, np.ndarray)
-        assert not isinstance(ct.ss(0, 0, 0, 1).D, np.matrix)
+            ct.reset_defaults()
 
-        ct.use_legacy_defaults('0.8.4')
-        assert ct.config.defaults['forced_response.return_x'] is True
+        with pytest.warns(UserWarning, match="NumPy matrix class no longer"):
+            ct.use_legacy_defaults('0.8.4')
+            assert ct.config.defaults['forced_response.return_x'] is True
 
         ct.use_legacy_defaults('0.9.0')
         assert isinstance(ct.ss(0, 0, 0, 1).D, np.ndarray)
diff --git a/control/tests/conftest.py b/control/tests/conftest.py
index b63db3e1..c5ab6cb8 100644
--- a/control/tests/conftest.py
+++ b/control/tests/conftest.py
@@ -13,24 +13,21 @@ TEST_MATRIX_AND_ARRAY = os.getenv("PYTHON_CONTROL_ARRAY_AND_MATRIX") == "1"
 
 # some common pytest marks. These can be used as test decorators or in
 # pytest.param(marks=)
-slycotonly = pytest.mark.skipif(not control.exception.slycot_check(),
-                                reason="slycot not installed")
-cvxoptonly = pytest.mark.skipif(not control.exception.cvxopt_check(),
-                                reason="cvxopt not installed")
-matrixfilter = pytest.mark.filterwarnings("ignore:.*matrix subclass:"
-                                          "PendingDeprecationWarning")
-matrixerrorfilter = pytest.mark.filterwarnings("error:.*matrix subclass:"
-                                               "PendingDeprecationWarning")
+slycotonly = pytest.mark.skipif(
+    not control.exception.slycot_check(), reason="slycot not installed")
+cvxoptonly = pytest.mark.skipif(
+    not control.exception.cvxopt_check(), reason="cvxopt not installed")
 
 
 @pytest.fixture(scope="session", autouse=True)
 def control_defaults():
     """Make sure the testing session always starts with the defaults.
 
-    This should be the first fixture initialized,
-    so that all other fixtures see the general defaults (unless they set them
-    themselves) even before importing control/__init__. Enforce this by adding
-    it as an argument to all other session scoped fixtures.
+    This should be the first fixture initialized, so that all other
+    fixtures see the general defaults (unless they set them themselves)
+    even before importing control/__init__. Enforce this by adding it as an
+    argument to all other session scoped fixtures.
+
     """
     control.reset_defaults()
     the_defaults = control.config.defaults.copy()
@@ -39,63 +36,6 @@ def control_defaults():
     assert control.config.defaults == the_defaults
 
 
-@pytest.fixture(scope="function", autouse=TEST_MATRIX_AND_ARRAY,
-                params=[pytest.param("arrayout", marks=matrixerrorfilter),
-                        pytest.param("matrixout", marks=matrixfilter)])
-def matarrayout(request):
-    """Switch the config to use np.ndarray and np.matrix as returns."""
-    restore = control.config.defaults['statesp.use_numpy_matrix']
-    control.use_numpy_matrix(request.param == "matrixout", warn=False)
-    yield
-    control.use_numpy_matrix(restore, warn=False)
-
-
-def ismatarrayout(obj):
-    """Test if the returned object has the correct type as configured.
-
-    note that isinstance(np.matrix(obj), np.ndarray) is True
-    """
-    use_matrix = control.config.defaults['statesp.use_numpy_matrix']
-    return (isinstance(obj, np.ndarray)
-            and isinstance(obj, np.matrix) == use_matrix)
-
-
-def asmatarrayout(obj):
-    """Return a object according to the configured default."""
-    use_matrix = control.config.defaults['statesp.use_numpy_matrix']
-    matarray = np.asmatrix if use_matrix else np.asarray
-    return matarray(obj)
-
-
-@contextmanager
-def check_deprecated_matrix():
-    """Check that a call produces a deprecation warning because of np.matrix."""
-    use_matrix = control.config.defaults['statesp.use_numpy_matrix']
-    if use_matrix:
-        with pytest.deprecated_call():
-            try:
-                yield
-            finally:
-                pass
-    else:
-        yield
-
-
-@pytest.fixture(scope="function",
-                params=[p for p, usebydefault in
-                        [(pytest.param(np.array,
-                                       id="arrayin"),
-                          True),
-                         (pytest.param(np.matrix,
-                                       id="matrixin",
-                                       marks=matrixfilter),
-                          False)]
-                        if usebydefault or TEST_MATRIX_AND_ARRAY])
-def matarrayin(request):
-    """Use array and matrix to construct input data in tests."""
-    return request.param
-
-
 @pytest.fixture(scope="function")
 def editsdefaults():
     """Make sure any changes to the defaults only last during a test."""
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 4012770b..1fe57d57 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -17,7 +17,6 @@ from math import sqrt
 
 import control as ct
 from control import iosys as ios
-from control.tests.conftest import matrixfilter
 
 
 class TestIOSys:
@@ -253,10 +252,10 @@ class TestIOSys:
         assert linearized_newnames.find_output('y') is None
 
         # Test legacy version as well
-        ct.use_legacy_defaults('0.8.4')
-        ct.config.use_numpy_matrix(False)       # np.matrix deprecated
-        linearized = kincar.linearize([0, 0, 0], [0, 0], copy_names=True)
-        assert linearized.name == kincar.name + '_linearized'
+        with pytest.warns(UserWarning, match="NumPy matrix class no longer"):
+            ct.use_legacy_defaults('0.8.4')
+            linearized = kincar.linearize([0, 0, 0], [0, 0], copy_names=True)
+            assert linearized.name == kincar.name + '_linearized'
 
     def test_connect(self, tsys):
         # Define a couple of (linear) systems to interconnection
@@ -1060,8 +1059,8 @@ class TestIOSys:
         """Enforce generic system names 'sys[i]' to be present when systems are
         created without explicit names."""
 
-        ct.config.use_legacy_defaults('0.8.4')  # changed delims in 0.9.0
-        ct.config.use_numpy_matrix(False)       # np.matrix deprecated
+        with pytest.warns(UserWarning, match="NumPy matrix class no longer"):
+            ct.config.use_legacy_defaults('0.8.4')  # changed delims in 0.9.0
 
         # Create a system with a known ID
         ct.namedio.NamedIOSystem._idCounter = 0
@@ -1128,8 +1127,8 @@ class TestIOSys:
         output: 'y[i]'
         """
 
-        ct.config.use_legacy_defaults('0.8.4')  # changed delims in 0.9.0
-        ct.config.use_numpy_matrix(False)       # np.matrix deprecated
+        with pytest.warns(UserWarning, match="NumPy matrix class no longer"):
+            ct.config.use_legacy_defaults('0.8.4')  # changed delims in 0.9.0
 
         # Create a system with a known ID
         ct.namedio.NamedIOSystem._idCounter = 0
@@ -1434,8 +1433,8 @@ class TestIOSys:
             ios_series = nlios * nlios
 
         # Nonduplicate objects
-        ct.config.use_legacy_defaults('0.8.4')  # changed delims in 0.9.0
-        ct.config.use_numpy_matrix(False)       # np.matrix deprecated
+        with pytest.warns(UserWarning, match="NumPy matrix class no longer"):
+            ct.config.use_legacy_defaults('0.8.4')  # changed delims in 0.9.0
         nlios1 = nlios.copy()
         nlios2 = nlios.copy()
         with pytest.warns(UserWarning, match="duplicate name"):
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
index 0746e3fe..49c2afd5 100644
--- a/control/tests/modelsimp_test.py
+++ b/control/tests/modelsimp_test.py
@@ -9,7 +9,7 @@ import pytest
 
 from control import StateSpace, forced_response, tf, rss, c2d
 from control.exception import ControlMIMONotImplemented
-from control.tests.conftest import slycotonly, matarrayin
+from control.tests.conftest import slycotonly
 from control.modelsimp import balred, hsvd, markov, modred
 
 
@@ -17,11 +17,11 @@ class TestModelsimp:
     """Test model reduction functions"""
 
     @slycotonly
-    def testHSVD(self, matarrayout, matarrayin):
-        A = matarrayin([[1., -2.], [3., -4.]])
-        B = matarrayin([[5.], [7.]])
-        C = matarrayin([[6., 8.]])
-        D = matarrayin([[9.]])
+    def testHSVD(self):
+        A = np.array([[1., -2.], [3., -4.]])
+        B = np.array([[5.], [7.]])
+        C = np.array([[6., 8.]])
+        D = np.array([[9.]])
         sys = StateSpace(A, B, C, D)
         hsv = hsvd(sys)
         hsvtrue = np.array([24.42686, 0.5731395])  # from MATLAB
@@ -32,8 +32,8 @@ class TestModelsimp:
         assert isinstance(hsv, np.ndarray)
         assert not isinstance(hsv, np.matrix)
 
-    def testMarkovSignature(self, matarrayout, matarrayin):
-        U = matarrayin([[1., 1., 1., 1., 1.]])
+    def testMarkovSignature(self):
+        U = np.array([[1., 1., 1., 1., 1.]])
         Y = U
         m = 3
         H = markov(Y, U, m, transpose=False)
@@ -111,17 +111,17 @@ class TestModelsimp:
         # for k=5, m=n=10: 0.015 %
         np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
 
-    def testModredMatchDC(self, matarrayin):
+    def testModredMatchDC(self):
         #balanced realization computed in matlab for the transfer function:
         # num = [1 11 45 32], den = [1 15 60 200 60]
-        A = matarrayin(
+        A = np.array(
             [[-1.958, -1.194, 1.824, -1.464],
              [-1.194, -0.8344, 2.563, -1.351],
              [-1.824, -2.563, -1.124, 2.704],
              [-1.464, -1.351, -2.704, -11.08]])
-        B = matarrayin([[-0.9057], [-0.4068], [-0.3263], [-0.3474]])
-        C = matarrayin([[-0.9057, -0.4068, 0.3263, -0.3474]])
-        D = matarrayin([[0.]])
+        B = np.array([[-0.9057], [-0.4068], [-0.3263], [-0.3474]])
+        C = np.array([[-0.9057, -0.4068, 0.3263, -0.3474]])
+        D = np.array([[0.]])
         sys = StateSpace(A, B, C, D)
         rsys = modred(sys,[2, 3],'matchdc')
         Artrue = np.array([[-4.431, -4.552], [-4.552, -5.361]])
@@ -133,30 +133,30 @@ class TestModelsimp:
         np.testing.assert_array_almost_equal(rsys.C, Crtrue, decimal=3)
         np.testing.assert_array_almost_equal(rsys.D, Drtrue, decimal=2)
 
-    def testModredUnstable(self, matarrayin):
+    def testModredUnstable(self):
         """Check if an error is thrown when an unstable system is given"""
-        A = matarrayin(
+        A = np.array(
             [[4.5418, 3.3999, 5.0342, 4.3808],
              [0.3890, 0.3599, 0.4195, 0.1760],
              [-4.2117, -3.2395, -4.6760, -4.2180],
              [0.0052, 0.0429, 0.0155, 0.2743]])
-        B = matarrayin([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0], [4.0, 4.0]])
-        C = matarrayin([[1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 4.0]])
-        D = matarrayin([[0.0, 0.0], [0.0, 0.0]])
+        B = np.array([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0], [4.0, 4.0]])
+        C = np.array([[1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 4.0]])
+        D = np.array([[0.0, 0.0], [0.0, 0.0]])
         sys = StateSpace(A, B, C, D)
         np.testing.assert_raises(ValueError, modred, sys, [2, 3])
 
-    def testModredTruncate(self, matarrayin):
+    def testModredTruncate(self):
         #balanced realization computed in matlab for the transfer function:
         # num = [1 11 45 32], den = [1 15 60 200 60]
-        A = matarrayin(
+        A = np.array(
             [[-1.958, -1.194, 1.824, -1.464],
              [-1.194, -0.8344, 2.563, -1.351],
              [-1.824, -2.563, -1.124, 2.704],
              [-1.464, -1.351, -2.704, -11.08]])
-        B = matarrayin([[-0.9057], [-0.4068], [-0.3263], [-0.3474]])
-        C = matarrayin([[-0.9057, -0.4068, 0.3263, -0.3474]])
-        D = matarrayin([[0.]])
+        B = np.array([[-0.9057], [-0.4068], [-0.3263], [-0.3474]])
+        C = np.array([[-0.9057, -0.4068, 0.3263, -0.3474]])
+        D = np.array([[0.]])
         sys = StateSpace(A, B, C, D)
         rsys = modred(sys,[2, 3],'truncate')
         Artrue = np.array([[-1.958, -1.194], [-1.194, -0.8344]])
@@ -170,18 +170,18 @@ class TestModelsimp:
 
 
     @slycotonly
-    def testBalredTruncate(self, matarrayin):
+    def testBalredTruncate(self):
         # controlable canonical realization computed in matlab for the transfer
         # function:
         # num = [1 11 45 32], den = [1 15 60 200 60]
-        A = matarrayin(
+        A = np.array(
             [[-15., -7.5, -6.25, -1.875],
              [8., 0., 0., 0.],
              [0., 4., 0., 0.],
              [0., 0., 1., 0.]])
-        B = matarrayin([[2.], [0.], [0.], [0.]])
-        C = matarrayin([[0.5, 0.6875, 0.7031, 0.5]])
-        D = matarrayin([[0.]])
+        B = np.array([[2.], [0.], [0.], [0.]])
+        C = np.array([[0.5, 0.6875, 0.7031, 0.5]])
+        D = np.array([[0.]])
         
         sys = StateSpace(A, B, C, D)
         orders = 2
@@ -211,18 +211,18 @@ class TestModelsimp:
         np.testing.assert_array_almost_equal(Dr, Drtrue, decimal=4)
 
     @slycotonly
-    def testBalredMatchDC(self, matarrayin):
+    def testBalredMatchDC(self):
         # controlable canonical realization computed in matlab for the transfer
         # function:
         # num = [1 11 45 32], den = [1 15 60 200 60]
-        A = matarrayin(
+        A = np.array(
             [[-15., -7.5, -6.25, -1.875],
              [8., 0., 0., 0.],
              [0., 4., 0., 0.],
              [0., 0., 1., 0.]])
-        B = matarrayin([[2.], [0.], [0.], [0.]])
-        C = matarrayin([[0.5, 0.6875, 0.7031, 0.5]])
-        D = matarrayin([[0.]])
+        B = np.array([[2.], [0.], [0.], [0.]])
+        C = np.array([[0.5, 0.6875, 0.7031, 0.5]])
+        D = np.array([[0.]])
         
         sys = StateSpace(A, B, C, D)
         orders = 2
diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py
index 951c817f..dc72c072 100644
--- a/control/tests/statefbk_test.py
+++ b/control/tests/statefbk_test.py
@@ -16,8 +16,7 @@ from control.exception import ControlDimension, ControlSlycot, \
 from control.mateqn import care, dare
 from control.statefbk import (ctrb, obsv, place, place_varga, lqr, dlqr,
                               gram, acker)
-from control.tests.conftest import (slycotonly, check_deprecated_matrix,
-                                    ismatarrayout, asmatarrayout)
+from control.tests.conftest import slycotonly
 
 
 @pytest.fixture
@@ -36,48 +35,37 @@ class TestStatefbk:
     # Set to True to print systems to the output.
     debug = False
 
-    def testCtrbSISO(self, matarrayin, matarrayout):
-        A = matarrayin([[1., 2.], [3., 4.]])
-        B = matarrayin([[5.], [7.]])
+    def testCtrbSISO(self):
+        A = np.array([[1., 2.], [3., 4.]])
+        B = np.array([[5.], [7.]])
         Wctrue = np.array([[5., 19.], [7., 43.]])
-
-        with check_deprecated_matrix():
-            Wc = ctrb(A, B)
-        assert ismatarrayout(Wc)
-
+        Wc = ctrb(A, B)
         np.testing.assert_array_almost_equal(Wc, Wctrue)
 
-    def testCtrbMIMO(self, matarrayin):
-        A = matarrayin([[1., 2.], [3., 4.]])
-        B = matarrayin([[5., 6.], [7., 8.]])
+    def testCtrbMIMO(self):
+        A = np.array([[1., 2.], [3., 4.]])
+        B = np.array([[5., 6.], [7., 8.]])
         Wctrue = np.array([[5., 6., 19., 22.], [7., 8., 43., 50.]])
         Wc = ctrb(A, B)
         np.testing.assert_array_almost_equal(Wc, Wctrue)
 
-        # Make sure default type values are correct
-        assert ismatarrayout(Wc)
-
-    def testObsvSISO(self, matarrayin):
-        A = matarrayin([[1., 2.], [3., 4.]])
-        C = matarrayin([[5., 7.]])
+    def testObsvSISO(self):
+        A = np.array([[1., 2.], [3., 4.]])
+        C = np.array([[5., 7.]])
         Wotrue = np.array([[5., 7.], [26., 38.]])
         Wo = obsv(A, C)
         np.testing.assert_array_almost_equal(Wo, Wotrue)
 
-        # Make sure default type values are correct
-        assert ismatarrayout(Wo)
-
-
-    def testObsvMIMO(self, matarrayin):
-        A = matarrayin([[1., 2.], [3., 4.]])
-        C = matarrayin([[5., 6.], [7., 8.]])
+    def testObsvMIMO(self):
+        A = np.array([[1., 2.], [3., 4.]])
+        C = np.array([[5., 6.], [7., 8.]])
         Wotrue = np.array([[5., 6.], [7., 8.], [23., 34.], [31., 46.]])
         Wo = obsv(A, C)
         np.testing.assert_array_almost_equal(Wo, Wotrue)
 
-    def testCtrbObsvDuality(self, matarrayin):
-        A = matarrayin([[1.2, -2.3], [3.4, -4.5]])
-        B = matarrayin([[5.8, 6.9], [8., 9.1]])
+    def testCtrbObsvDuality(self):
+        A = np.array([[1.2, -2.3], [3.4, -4.5]])
+        B = np.array([[5.8, 6.9], [8., 9.1]])
         Wc = ctrb(A, B)
         A = np.transpose(A)
         C = np.transpose(B)
@@ -85,59 +73,55 @@ class TestStatefbk:
         np.testing.assert_array_almost_equal(Wc,Wo)
 
     @slycotonly
-    def testGramWc(self, matarrayin, matarrayout):
-        A = matarrayin([[1., -2.], [3., -4.]])
-        B = matarrayin([[5., 6.], [7., 8.]])
-        C = matarrayin([[4., 5.], [6., 7.]])
-        D = matarrayin([[13., 14.], [15., 16.]])
+    def testGramWc(self):
+        A = np.array([[1., -2.], [3., -4.]])
+        B = np.array([[5., 6.], [7., 8.]])
+        C = np.array([[4., 5.], [6., 7.]])
+        D = np.array([[13., 14.], [15., 16.]])
         sys = ss(A, B, C, D)
         Wctrue = np.array([[18.5, 24.5], [24.5, 32.5]])
-
-        with check_deprecated_matrix():
-            Wc = gram(sys, 'c')
-
-        assert ismatarrayout(Wc)
+        Wc = gram(sys, 'c')
         np.testing.assert_array_almost_equal(Wc, Wctrue)
 
     @slycotonly
-    def testGramRc(self, matarrayin):
-        A = matarrayin([[1., -2.], [3., -4.]])
-        B = matarrayin([[5., 6.], [7., 8.]])
-        C = matarrayin([[4., 5.], [6., 7.]])
-        D = matarrayin([[13., 14.], [15., 16.]])
+    def testGramRc(self):
+        A = np.array([[1., -2.], [3., -4.]])
+        B = np.array([[5., 6.], [7., 8.]])
+        C = np.array([[4., 5.], [6., 7.]])
+        D = np.array([[13., 14.], [15., 16.]])
         sys = ss(A, B, C, D)
         Rctrue = np.array([[4.30116263, 5.6961343], [0., 0.23249528]])
         Rc = gram(sys, 'cf')
         np.testing.assert_array_almost_equal(Rc, Rctrue)
 
     @slycotonly
-    def testGramWo(self, matarrayin):
-        A = matarrayin([[1., -2.], [3., -4.]])
-        B = matarrayin([[5., 6.], [7., 8.]])
-        C = matarrayin([[4., 5.], [6., 7.]])
-        D = matarrayin([[13., 14.], [15., 16.]])
+    def testGramWo(self):
+        A = np.array([[1., -2.], [3., -4.]])
+        B = np.array([[5., 6.], [7., 8.]])
+        C = np.array([[4., 5.], [6., 7.]])
+        D = np.array([[13., 14.], [15., 16.]])
         sys = ss(A, B, C, D)
         Wotrue = np.array([[257.5, -94.5], [-94.5, 56.5]])
         Wo = gram(sys, 'o')
         np.testing.assert_array_almost_equal(Wo, Wotrue)
 
     @slycotonly
-    def testGramWo2(self, matarrayin):
-        A = matarrayin([[1., -2.], [3., -4.]])
-        B = matarrayin([[5.], [7.]])
-        C = matarrayin([[6., 8.]])
-        D = matarrayin([[9.]])
+    def testGramWo2(self):
+        A = np.array([[1., -2.], [3., -4.]])
+        B = np.array([[5.], [7.]])
+        C = np.array([[6., 8.]])
+        D = np.array([[9.]])
         sys = ss(A,B,C,D)
         Wotrue = np.array([[198., -72.], [-72., 44.]])
         Wo = gram(sys, 'o')
         np.testing.assert_array_almost_equal(Wo, Wotrue)
 
     @slycotonly
-    def testGramRo(self, matarrayin):
-        A = matarrayin([[1., -2.], [3., -4.]])
-        B = matarrayin([[5., 6.], [7., 8.]])
-        C = matarrayin([[4., 5.], [6., 7.]])
-        D = matarrayin([[13., 14.], [15., 16.]])
+    def testGramRo(self):
+        A = np.array([[1., -2.], [3., -4.]])
+        B = np.array([[5., 6.], [7., 8.]])
+        C = np.array([[4., 5.], [6., 7.]])
+        D = np.array([[13., 14.], [15., 16.]])
         sys = ss(A, B, C, D)
         Rotrue = np.array([[16.04680654, -5.8890222], [0., 4.67112593]])
         Ro = gram(sys, 'of')
@@ -195,19 +179,18 @@ class TestStatefbk:
         P_placed.sort()
         np.testing.assert_array_almost_equal(P_expected, P_placed)
 
-    def testPlace(self, matarrayin):
+    def testPlace(self):
         # Matrices shamelessly stolen from scipy example code.
-        A = matarrayin([[1.380, -0.2077, 6.715, -5.676],
+        A = np.array([[1.380, -0.2077, 6.715, -5.676],
                         [-0.5814, -4.290, 0, 0.6750],
                         [1.067, 4.273, -6.654, 5.893],
                         [0.0480, 4.273, 1.343, -2.104]])
-        B = matarrayin([[0, 5.679],
+        B = np.array([[0, 5.679],
                         [1.136, 1.136],
                         [0, 0],
                         [-3.146, 0]])
-        P = matarrayin([-0.5 + 1j, -0.5 - 1j, -5.0566, -8.6659])
+        P = np.array([-0.5 + 1j, -0.5 - 1j, -5.0566, -8.6659])
         K = place(A, B, P)
-        assert ismatarrayout(K)
         P_placed = np.linalg.eigvals(A - B @ K)
         self.checkPlaced(P, P_placed)
 
@@ -219,17 +202,17 @@ class TestStatefbk:
 
         # Check that we get an error if we ask for too many poles in the same
         # location. Here, rank(B) = 2, so lets place three at the same spot.
-        P_repeated = matarrayin([-0.5, -0.5, -0.5, -8.6659])
+        P_repeated = np.array([-0.5, -0.5, -0.5, -8.6659])
         with pytest.raises(ValueError):
             place(A, B, P_repeated)
 
     @slycotonly
-    def testPlace_varga_continuous(self, matarrayin):
+    def testPlace_varga_continuous(self):
         """
         Check that we can place eigenvalues for dtime=False
         """
-        A = matarrayin([[1., -2.], [3., -4.]])
-        B = matarrayin([[5.], [7.]])
+        A = np.array([[1., -2.], [3., -4.]])
+        B = np.array([[5.], [7.]])
 
         P = [-2., -2.]
         K = place_varga(A, B, P)
@@ -242,26 +225,26 @@ class TestStatefbk:
 
         # Regression test against bug #177
         # https://github.com/python-control/python-control/issues/177
-        A = matarrayin([[0, 1], [100, 0]])
-        B = matarrayin([[0], [1]])
-        P = matarrayin([-20 + 10*1j, -20 - 10*1j])
+        A = np.array([[0, 1], [100, 0]])
+        B = np.array([[0], [1]])
+        P = np.array([-20 + 10*1j, -20 - 10*1j])
         K = place_varga(A, B, P)
         P_placed = np.linalg.eigvals(A - B @ K)
         self.checkPlaced(P, P_placed)
 
 
     @slycotonly
-    def testPlace_varga_continuous_partial_eigs(self, matarrayin):
+    def testPlace_varga_continuous_partial_eigs(self):
         """
         Check that we are able to use the alpha parameter to only place
         a subset of the eigenvalues, for the continous time case.
         """
         # A matrix has eigenvalues at s=-1, and s=-2. Choose alpha = -1.5
         # and check that eigenvalue at s=-2 stays put.
-        A = matarrayin([[1., -2.], [3., -4.]])
-        B = matarrayin([[5.], [7.]])
+        A = np.array([[1., -2.], [3., -4.]])
+        B = np.array([[5.], [7.]])
 
-        P = matarrayin([-3.])
+        P = np.array([-3.])
         P_expected = np.array([-2.0, -3.0])
         alpha = -1.5
         K = place_varga(A, B, P, alpha=alpha)
@@ -271,30 +254,30 @@ class TestStatefbk:
         self.checkPlaced(P_expected, P_placed)
 
     @slycotonly
-    def testPlace_varga_discrete(self, matarrayin):
+    def testPlace_varga_discrete(self):
         """
         Check that we can place poles using dtime=True (discrete time)
         """
-        A = matarrayin([[1., 0], [0, 0.5]])
-        B = matarrayin([[5.], [7.]])
+        A = np.array([[1., 0], [0, 0.5]])
+        B = np.array([[5.], [7.]])
 
-        P = matarrayin([0.5, 0.5])
+        P = np.array([0.5, 0.5])
         K = place_varga(A, B, P, dtime=True)
         P_placed = np.linalg.eigvals(A - B @ K)
         # No guarantee of the ordering, so sort them
         self.checkPlaced(P, P_placed)
 
     @slycotonly
-    def testPlace_varga_discrete_partial_eigs(self, matarrayin):
+    def testPlace_varga_discrete_partial_eigs(self):
         """"
         Check that we can only assign a single eigenvalue in the discrete
         time case.
         """
         # A matrix has eigenvalues at 1.0 and 0.5. Set alpha = 0.51, and
         # check that the eigenvalue at 0.5 is not moved.
-        A = matarrayin([[1., 0], [0, 0.5]])
-        B = matarrayin([[5.], [7.]])
-        P = matarrayin([0.2, 0.6])
+        A = np.array([[1., 0], [0, 0.5]])
+        B = np.array([[5.], [7.]])
+        P = np.array([0.2, 0.6])
         P_expected = np.array([0.5, 0.6])
         alpha = 0.51
         K = place_varga(A, B, P, dtime=True, alpha=alpha)
@@ -302,49 +285,49 @@ class TestStatefbk:
         self.checkPlaced(P_expected, P_placed)
 
     def check_LQR(self, K, S, poles, Q, R):
-        S_expected = asmatarrayout(np.sqrt(Q @ R))
-        K_expected = asmatarrayout(S_expected / R)
+        S_expected = np.sqrt(Q @ R)
+        K_expected = S_expected / R
         poles_expected = -np.squeeze(np.asarray(K_expected))
         np.testing.assert_array_almost_equal(S, S_expected)
         np.testing.assert_array_almost_equal(K, K_expected)
         np.testing.assert_array_almost_equal(poles, poles_expected)
 
     def check_DLQR(self, K, S, poles, Q, R):
-        S_expected = asmatarrayout(Q)
-        K_expected = asmatarrayout(0)
+        S_expected = Q
+        K_expected = 0
         poles_expected = -np.squeeze(np.asarray(K_expected))
         np.testing.assert_array_almost_equal(S, S_expected)
         np.testing.assert_array_almost_equal(K, K_expected)
         np.testing.assert_array_almost_equal(poles, poles_expected)
 
     @pytest.mark.parametrize("method", [None, 'slycot', 'scipy'])
-    def test_LQR_integrator(self, matarrayin, matarrayout, method):
+    def test_LQR_integrator(self, method):
         if method == 'slycot' and not slycot_check():
             return
-        A, B, Q, R = (matarrayin([[X]]) for X in [0., 1., 10., 2.])
+        A, B, Q, R = (np.array([[X]]) for X in [0., 1., 10., 2.])
         K, S, poles = lqr(A, B, Q, R, method=method)
         self.check_LQR(K, S, poles, Q, R)
 
     @pytest.mark.parametrize("method", [None, 'slycot', 'scipy'])
-    def test_LQR_3args(self, matarrayin, matarrayout, method):
+    def test_LQR_3args(self, method):
         if method == 'slycot' and not slycot_check():
             return
         sys = ss(0., 1., 1., 0.)
-        Q, R = (matarrayin([[X]]) for X in [10., 2.])
+        Q, R = (np.array([[X]]) for X in [10., 2.])
         K, S, poles = lqr(sys, Q, R, method=method)
         self.check_LQR(K, S, poles, Q, R)
 
     @pytest.mark.parametrize("method", [None, 'slycot', 'scipy'])
-    def test_DLQR_3args(self, matarrayin, matarrayout, method):
+    def test_DLQR_3args(self, method):
         if method == 'slycot' and not slycot_check():
             return
         dsys = ss(0., 1., 1., 0., .1)
-        Q, R = (matarrayin([[X]]) for X in [10., 2.])
+        Q, R = (np.array([[X]]) for X in [10., 2.])
         K, S, poles = dlqr(dsys, Q, R, method=method)
         self.check_DLQR(K, S, poles, Q, R)
 
-    def test_DLQR_4args(self, matarrayin, matarrayout):
-        A, B, Q, R = (matarrayin([[X]]) for X in [0., 1., 10., 2.])
+    def test_DLQR_4args(self):
+        A, B, Q, R = (np.array([[X]]) for X in [0., 1., 10., 2.])
         K, S, poles = dlqr(A, B, Q, R)
         self.check_DLQR(K, S, poles, Q, R)
 
@@ -443,14 +426,14 @@ class TestStatefbk:
         with pytest.warns(UserWarning):
             (K, S, E) = dlqr(A, B, Q, R, N)
 
-    def test_care(self, matarrayin):
+    def test_care(self):
         """Test stabilizing and anti-stabilizing feedback, continuous"""
-        A = matarrayin(np.diag([1, -1]))
-        B = matarrayin(np.identity(2))
-        Q = matarrayin(np.identity(2))
-        R = matarrayin(np.identity(2))
-        S = matarrayin(np.zeros((2, 2)))
-        E = matarrayin(np.identity(2))
+        A = np.diag([1, -1])
+        B = np.identity(2)
+        Q = np.identity(2)
+        R = np.identity(2)
+        S = np.zeros((2, 2))
+        E = np.identity(2)
 
         X, L, G = care(A, B, Q, R, S, E, stabilizing=True)
         assert np.all(np.real(L) < 0)
@@ -465,14 +448,14 @@ class TestStatefbk:
     @pytest.mark.parametrize(
         "stabilizing",
         [True, pytest.param(False, marks=slycotonly)])
-    def test_dare(self, matarrayin, stabilizing):
+    def test_dare(self, stabilizing):
         """Test stabilizing and anti-stabilizing feedback, discrete"""
-        A = matarrayin(np.diag([0.5, 2]))
-        B = matarrayin(np.identity(2))
-        Q = matarrayin(np.identity(2))
-        R = matarrayin(np.identity(2))
-        S = matarrayin(np.zeros((2, 2)))
-        E = matarrayin(np.identity(2))
+        A = np.diag([0.5, 2])
+        B = np.identity(2)
+        Q = np.identity(2)
+        R = np.identity(2)
+        S = np.zeros((2, 2))
+        E = np.identity(2)
 
         X, L, G = dare(A, B, Q, R, S, E, stabilizing=stabilizing)
         sgn = {True: -1, False: 1}[stabilizing]
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 1182674c..83dc58b4 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -21,7 +21,7 @@ from control.lti import evalfr
 from control.statesp import StateSpace, _convert_to_statespace, tf2ss, \
     _statesp_defaults, _rss_generate, linfnorm
 from control.iosys import ss, rss, drss
-from control.tests.conftest import ismatarrayout, slycotonly
+from control.tests.conftest import slycotonly
 from control.xferfcn import TransferFunction, ss2tf
 
 
@@ -196,17 +196,6 @@ class TestStateSpace:
         sys = StateSpace(sysin)
         assert sys.dt is None
 
-    def test_matlab_style_constructor(self):
-        """Use (deprecated) matrix-style construction string"""
-        with pytest.deprecated_call():
-            sys = StateSpace("-1 1; 0 2", "0; 1", "1, 0", "0")
-        assert sys.A.shape == (2, 2)
-        assert sys.B.shape == (2, 1)
-        assert sys.C.shape == (1, 2)
-        assert sys.D.shape == (1, 1)
-        for X in [sys.A, sys.B, sys.C, sys.D]:
-            assert ismatarrayout(X)
-
     def test_D_broadcast(self, sys623):
         """Test broadcast of D=0 to the right shape"""
         # Giving D as a scalar 0 should broadcast to the right shape
@@ -1017,13 +1006,7 @@ class TestLTIConverter:
 
 class TestStateSpaceConfig:
     """Test the configuration of the StateSpace module"""
-
-    @pytest.fixture
-    def matarrayout(self):
-        """Override autoused global fixture within this class"""
-        pass
-
-    def test_statespace_defaults(self, matarrayout):
+    def test_statespace_defaults(self):
         """Make sure the tests are run with the configured defaults"""
         for k, v in _statesp_defaults.items():
             assert defaults[k] == v, \
diff --git a/control/tests/stochsys_test.py b/control/tests/stochsys_test.py
index b2d90e2a..91f4a1a0 100644
--- a/control/tests/stochsys_test.py
+++ b/control/tests/stochsys_test.py
@@ -3,7 +3,6 @@
 
 import numpy as np
 import pytest
-from control.tests.conftest import asmatarrayout
 
 import control as ct
 import control.optimal as opt
@@ -12,8 +11,8 @@ from math import log, pi
 
 # Utility function to check LQE answer
 def check_LQE(L, P, poles, G, QN, RN):
-    P_expected = asmatarrayout(np.sqrt(G @ QN @ G @ RN))
-    L_expected = asmatarrayout(P_expected / RN)
+    P_expected = np.sqrt(G @ QN @ G @ RN)
+    L_expected = P_expected / RN
     poles_expected = -np.squeeze(np.asarray(L_expected))
     np.testing.assert_almost_equal(P, P_expected)
     np.testing.assert_almost_equal(L, L_expected)
@@ -21,19 +20,19 @@ def check_LQE(L, P, poles, G, QN, RN):
 
 # Utility function to check discrete LQE solutions
 def check_DLQE(L, P, poles, G, QN, RN):
-    P_expected = asmatarrayout(G.dot(QN).dot(G))
-    L_expected = asmatarrayout(0)
+    P_expected = G.dot(QN).dot(G)
+    L_expected = 0
     poles_expected = -np.squeeze(np.asarray(L_expected))
     np.testing.assert_almost_equal(P, P_expected)
     np.testing.assert_almost_equal(L, L_expected)
     np.testing.assert_almost_equal(poles, poles_expected)
 
 @pytest.mark.parametrize("method", [None, 'slycot', 'scipy'])
-def test_LQE(matarrayin, method):
+def test_LQE(method):
     if method == 'slycot' and not slycot_check():
         return
 
-    A, G, C, QN, RN = (matarrayin([[X]]) for X in [0., .1, 1., 10., 2.])
+    A, G, C, QN, RN = (np.array([[X]]) for X in [0., .1, 1., 10., 2.])
     L, P, poles = lqe(A, G, C, QN, RN, method=method)
     check_LQE(L, P, poles, G, QN, RN)
 
@@ -80,11 +79,11 @@ def test_lqe_call_format(cdlqe):
         L, P, E = cdlqe(sys_tf, Q, R)
 
 @pytest.mark.parametrize("method", [None, 'slycot', 'scipy'])
-def test_DLQE(matarrayin, method):
+def test_DLQE(method):
     if method == 'slycot' and not slycot_check():
         return
 
-    A, G, C, QN, RN = (matarrayin([[X]]) for X in [0., .1, 1., 10., 2.])
+    A, G, C, QN, RN = (np.array([[X]]) for X in [0., .1, 1., 10., 2.])
     L, P, poles = dlqe(A, G, C, QN, RN, method=method)
     check_DLQE(L, P, poles, G, QN, RN)
 
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 124e16c1..ccd80816 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -857,7 +857,7 @@ class TestTimeresp:
                                      initial_response,
                                      forced_response])
     @pytest.mark.parametrize("squeeze", [None, True, False])
-    def test_time_vector(self, tsystem, fun, squeeze, matarrayout):
+    def test_time_vector(self, tsystem, fun, squeeze):
         """Test time vector handling and correct output convention
 
         gh-239, gh-295
@@ -1130,8 +1130,8 @@ class TestTimeresp:
     ])
     def test_squeeze_0_8_4(self, nstate, nout, ninp, squeeze, shape):
         # Set defaults to match release 0.8.4
-        ct.config.use_legacy_defaults('0.8.4')
-        ct.config.use_numpy_matrix(False)
+        with pytest.warns(UserWarning, match="NumPy matrix class no longer"):
+            ct.config.use_legacy_defaults('0.8.4')
 
         # Generate system, time, and input vectors
         sys = ct.rss(nstate, nout, ninp, strictly_proper=True)
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index eec30580..fd1076db 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -14,7 +14,7 @@ from control import isctime, isdtime, sample_system
 from control import defaults, reset_defaults, set_defaults
 from control.statesp import _convert_to_statespace
 from control.xferfcn import _convert_to_transfer_function
-from control.tests.conftest import slycotonly, matrixfilter
+from control.tests.conftest import slycotonly
 
 
 class TestXferFcn:
@@ -749,20 +749,16 @@ class TestXferFcn:
         np.testing.assert_array_almost_equal(sys.num[1][1], tm.num[1][2])
         np.testing.assert_array_almost_equal(sys.den[1][1], tm.den[1][2])
 
-    @pytest.mark.parametrize(
-        "matarrayin",
-        [pytest.param(np.array,
-                      id="arrayin",
-                      marks=[pytest.mark.skip(".__matmul__ not implemented")]),
-         pytest.param(np.matrix,
-                      id="matrixin",
-                      marks=matrixfilter)],
-        indirect=True)
-    @pytest.mark.parametrize("X_, ij",
-                             [([[2., 0., ]], 0),
-                              ([[0., 2., ]], 1)])
-    def test_matrix_array_multiply(self, matarrayin, X_, ij):
-        """Test mulitplication of MIMO TF with matrix and matmul with array"""
+    @pytest.mark.parametrize("op", [
+        pytest.param('mul'),
+        pytest.param(
+            'matmul', marks=pytest.mark.skip(".__matmul__ not implemented")),
+    ])
+    @pytest.mark.parametrize("X, ij",
+                             [(np.array([[2., 0., ]]), 0),
+                              (np.array([[0., 2., ]]), 1)])
+    def test_matrix_array_multiply(self, op, X, ij):
+        """Test mulitplication of MIMO TF with matrix"""
         # 2 inputs, 2 outputs with prime zeros so they do not cancel
         n = 2
         p = [3, 5, 7, 11, 13, 17, 19, 23]
@@ -771,13 +767,12 @@ class TestXferFcn:
              for i in range(n)],
             [[[1, -1]] * n] * n)
 
-        X = matarrayin(X_)
-
-        if matarrayin is np.matrix:
+        if op == 'matmul':
+            XH = X @ H
+        elif op == 'mul':
             XH = X * H
         else:
-            # XH = X @ H
-            XH = np.matmul(X, H)
+            assert NotImplemented(f"unknown operator '{op}'")
         XH = XH.minreal()
         assert XH.ninputs == n
         assert XH.noutputs == X.shape[0]
@@ -790,11 +785,12 @@ class TestXferFcn:
         np.testing.assert_allclose(2. * H.num[ij][1], XH.num[0][1], rtol=1e-4)
         np.testing.assert_allclose(     H.den[ij][1], XH.den[0][1], rtol=1e-4)
 
-        if matarrayin is np.matrix:
+        if op == 'matmul':
+            HXt = H @ X.T
+        elif op == 'mul':
             HXt = H * X.T
         else:
-            # HXt = H @ X.T
-            HXt = np.matmul(H, X.T)
+            assert NotImplemented(f"unknown operator '{op}'")
         HXt = HXt.minreal()
         assert HXt.ninputs == X.T.shape[1]
         assert HXt.noutputs == n