Repository: python-control/python-control
Pull request #829 (2022-12-30 20:37:55)

continuous time system support for create_estimator_iosystem
This PR adds support for continuous time systems in the `create_estimator_iosystem` function + unit tests.

(I also had to update some of the xfail exception conditions since it looks like NumPy is now at v1.24.1, but there are still some failures under a small number of conditions.)

This should be the last update before putting out control-0.9.3, which I hope to do in the next day or two.


Patch:
diff --git a/control/iosys.py b/control/iosys.py
index e16a15ed..a9dc4514 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -536,10 +536,6 @@ class InputOutputSystem(NamedIOSystem):
         # functions.
         #
 
-        # If x0 and u0 are specified as lists, concatenate the elements
-        x0 = _concatenate_list_elements(x0, 'x0')
-        u0 = _concatenate_list_elements(u0, 'u0')
-
         # Figure out dimensions if they were not specified.
         nstates = _find_size(self.nstates, x0)
         ninputs = _find_size(self.ninputs, u0)
@@ -1766,7 +1762,14 @@ def input_output_response(
     ninputs = U.shape[0]
 
     # If we were passed a list of initial states, concatenate them
-    X0 = _concatenate_list_elements(X0, 'X0')
+    if isinstance(X0, (tuple, list)):
+        X0_list = []
+        for i, x0 in enumerate(X0):
+            x0 = np.array(x0).reshape(-1)       # convert everyting to 1D array
+            X0_list += x0.tolist()              # add elements to initial state
+
+        # Save the newly created input vector
+        X0 = np.array(X0_list)
 
     # If the initial state is too short, make it longer (NB: sys.nstates
     # could be None if nstates comes from size of initial condition)
@@ -2390,19 +2393,6 @@ def ss(*args, **kwargs):
     return sys
 
 
-# Utility function to allow lists states, inputs
-def _concatenate_list_elements(X, name='X'):
-    # If we were passed a list, concatenate the elements together
-    if isinstance(X, (tuple, list)):
-        X_list = []
-        for i, x in enumerate(X):
-            x = np.array(x).reshape(-1)         # convert everyting to 1D array
-            X_list += x.tolist()                # add elements to initial state
-        return np.array(X_list)
-
-    # Otherwise, do nothing
-    return X
-
 def rss(states=1, outputs=1, inputs=1, strictly_proper=False, **kwargs):
     """Create a stable random state space object.
 
diff --git a/control/statefbk.py b/control/statefbk.py
index 1f61134a..c88b98c5 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -41,7 +41,6 @@
 
 # External packages and modules
 import numpy as np
-import scipy as sp
 
 from . import statesp
 from .mateqn import care, dare, _check_shape
@@ -72,7 +71,7 @@ except ImportError:
     sb03od = None
 
 
-__all__ = ['ctrb', 'obsv', 'gram', 'place', 'place_varga', 'lqr',
+__all__ = ['ctrb', 'obsv', 'gram', 'place', 'place_varga', 'lqr', 
            'dlqr', 'acker', 'create_statefbk_iosystem']
 
 
@@ -597,10 +596,8 @@ def dlqr(*args, **kwargs):
 
 # Function to create an I/O sytems representing a state feedback controller
 def create_statefbk_iosystem(
-        sys, gain, integral_action=None, estimator=None, type=None,
-        xd_labels='xd[{i}]', ud_labels='ud[{i}]', gainsched_indices=None,
-        gainsched_method='linear', name=None, inputs=None, outputs=None,
-        states=None):
+        sys, K, integral_action=None, xd_labels='xd[{i}]', ud_labels='ud[{i}]',
+        estimator=None, type='linear'):
     """Create an I/O system using a (full) state feedback controller
 
     This function creates an input/output system that implements a
@@ -616,47 +613,31 @@ def create_statefbk_iosystem(
     feedback gain (eg, from LQR).  The function returns the controller
     ``ctrl`` and the closed loop systems ``clsys``, both as I/O systems.
 
-    A gain scheduled controller can also be created, by passing a list of
-    gains and a corresponding list of values of a set of scheduling
-    variables.  In this case, the controller has the form
-
-        u = ud - K_p(mu) (x - xd) - K_i(mu) integral(C x - C x_d)
-
-    where mu represents the scheduling variable.
-
     Parameters
     ----------
     sys : InputOutputSystem
         The I/O system that represents the process dynamics.  If no estimator
         is given, the output of this system should represent the full state.
 
-    gain : ndarray or tuple
-        If a array is give, it represents the state feedback gain (K).
-        This matrix defines the gains to be applied to the system.  If
-        ``integral_action`` is None, then the dimensions of this array
-        should be (sys.ninputs, sys.nstates).  If `integral action` is
-        set to a matrix or a function, then additional columns
-        represent the gains of the integral states of the controller.
-
-        If a tuple is given, then it specifies a gain schedule.  The tuple
-        should be of the form ``(gains, points)`` where gains is a list of
-        gains :math:`K_j` and points is a list of values :math:`\\mu_j` at
-        which the gains are computed.  The `gainsched_indices` parameter
-        should be used to specify the scheduling variables.
+    K : ndarray
+        The state feedback gain.  This matrix defines the gains to be
+        applied to the system.  If ``integral_action`` is None, then the
+        dimensions of this array should be (sys.ninputs, sys.nstates).  If
+        `integral action` is set to a matrix or a function, then additional
+        columns represent the gains of the integral states of the
+        controller.
 
     xd_labels, ud_labels : str or list of str, optional
         Set the name of the signals to use for the desired state and inputs.
         If a single string is specified, it should be a format string using
-        the variable ``i`` as an index.  Otherwise, a list of strings
-        matching the size of xd and ud, respectively, should be used.
-        Default is ``'xd[{i}]'`` for xd_labels and ``'ud[{i}]'`` for
-        ud_labels.  These settings can also be overriden using the `inputs`
-        keyword.
+        the variable ``i`` as an index.  Otherwise, a list of strings matching
+        the size of xd and ud, respectively, should be used.  Default is
+        ``'xd[{i}]'`` for xd_labels and ``'xd[{i}]'`` for ud_labels.
 
     integral_action : None, ndarray, or func, optional
         If this keyword is specified, the controller can include integral
-        action in addition to state feedback.  If ``integral_action`` is a
-        matrix, it will be multiplied by the current and desired state to
+        action in addition to state feedback.  If ``integral_action`` is an
+        ndarray, it will be multiplied by the current and desired state to
         generate the error for the internal integrator states of the control
         law.  If ``integral_action`` is a function ``h``, that function will
         be called with the signature h(t, x, u, params) to obtain the
@@ -665,63 +646,30 @@ def create_statefbk_iosystem(
         ``K`` matrix.
 
     estimator : InputOutputSystem, optional
-        If an estimator is provided, use the states of the estimator as
+        If an estimator is provided, using the states of the estimator as
         the system inputs for the controller.
 
-    gainsched_indices : list of int or str, optional
-        If a gain scheduled controller is specified, specify the indices of
-        the controller input to use for scheduling the gain. The input to
-        the controller is the desired state xd, the desired input ud, and
-        the system state x (or state estimate xhat, if an estimator is
-        given). The indices can either be specified as integer offsets into
-        the input vector or as strings matching the signal names of the
-        input vector. The default is to use the desire state xd.
-
-    gainsched_method : str, optional
-        The method to use for gain scheduling.  Possible values are 'linear'
-        (default), 'nearest', and 'cubic'.  More information is available in
-        :func:`scipy.interpolate.griddata`. For points outside of the convex
-        hull of the scheduling points, the gain at the nearest point is
-        used.
-
-    type : 'linear' or 'nonlinear', optional
-        Set the type of controller to create. The default for a linear gain
-        is a linear controller implementing the LQR regulator. If the type
-        is 'nonlinear', a :class:NonlinearIOSystem is created instead, with
-        the gain ``K`` as a parameter (allowing modifications of the gain at
-        runtime). If the gain parameter is a tuple, then a nonlinear,
-        gain-scheduled controller is created.
+    type : 'nonlinear' or 'linear', optional
+        Set the type of controller to create. The default is a linear
+        controller implementing the LQR regulator. If the type is 'nonlinear',
+        a :class:NonlinearIOSystem is created instead, with the gain ``K`` as
+        a parameter (allowing modifications of the gain at runtime).
 
     Returns
     -------
     ctrl : InputOutputSystem
         Input/output system representing the controller.  This system takes
-        as inputs the desired state ``xd``, the desired input ``ud``, and
-        either the system state ``x`` or the estimated state ``xhat``.  It
-        outputs the controller action u according to the formula :math:`u =
-        u_d - K(x - x_d)`.  If the keyword ``integral_action`` is specified,
-        then an additional set of integrators is included in the control
-        system (with the gain matrix ``K`` having the integral gains
-        appended after the state gains).  If a gain scheduled controller is
-        specified, the gain (proportional and integral) are evaluated using
-        the scheduling variables specified by ``gainsched_indices``.
+        as inputs the desired state xd, the desired input ud, and the system
+        state x.  It outputs the controller action u according to the
+        formula u = ud - K(x - xd).  If the keyword `integral_action` is
+        specified, then an additional set of integrators is included in the
+        control system (with the gain matrix K having the integral gains
+        appended after the state gains).
 
     clsys : InputOutputSystem
         Input/output system representing the closed loop system.  This
-        systems takes as inputs the desired trajectory ``(xd, ud)`` and
-        outputs the system state ``x`` and the applied input ``u``
-        (vertically stacked).
-
-    Other Parameters
-    ----------------
-    inputs, outputs : str, or list of str, optional
-        List of strings that name the individual signals of the transformed
-        system.  If not given, the inputs and outputs are the same as the
-        original system.
-    
-    name : string, optional
-        System name. If unspecified, a generic name <sys[id]> is generated
-        with a unique integer id.
+        systems takes as inputs the desired trajectory (xd, ud) and outputs
+        the system state x and the applied input u (vertically stacked).
 
     """
     # Make sure that we were passed an I/O system as an input
@@ -757,127 +705,53 @@ def create_statefbk_iosystem(
         C = np.zeros((0, sys.nstates))
 
     # Check to make sure that state feedback has the right shape
-    if isinstance(gain, np.ndarray):
-        K = gain
-        if K.shape != (sys.ninputs, estimator.noutputs + nintegrators):
-            raise ControlArgument(
-                f'Control gain must be an array of size {sys.ninputs}'
-                f'x {sys.nstates}' +
-                (f'+{nintegrators}' if nintegrators > 0 else ''))
-        gainsched = False
-
-    elif isinstance(gain, tuple):
-        # Check for gain scheduled controller
-        if len(gain) != 2:
-            raise ControlArgument("gain must be a 2-tuple for gain scheduling")
-        gains, points = gain[0:2]
-
-        # Stack gains and points if past as a list
-        gains = np.stack(gains)
-        points = np.stack(points)
-        gainsched=True
-
-    else:
-        raise ControlArgument("gain must be an array or a tuple")
-
-    # Decide on the type of system to create
-    if gainsched and type == 'linear':
+    if not isinstance(K, np.ndarray) or \
+       K.shape != (sys.ninputs, estimator.noutputs + nintegrators):
         raise ControlArgument(
-            "type 'linear' not allowed for gain scheduled controller")
-    elif type is None:
-        type = 'nonlinear' if gainsched else 'linear'
-    elif type not in {'linear', 'nonlinear'}:
-        raise ControlArgument(f"unknown type '{type}'")
+            f'Control gain must be an array of size {sys.ninputs}'
+            f'x {sys.nstates}' +
+            (f'+{nintegrators}' if nintegrators > 0 else ''))
 
     # Figure out the labels to use
     if isinstance(xd_labels, str):
-        # Generate the list of labels using the argument as a format string
+        # Gnerate the list of labels using the argument as a format string
         xd_labels = [xd_labels.format(i=i) for i in range(sys.nstates)]
 
     if isinstance(ud_labels, str):
-        # Generate the list of labels using the argument as a format string
+        # Gnerate the list of labels using the argument as a format string
         ud_labels = [ud_labels.format(i=i) for i in range(sys.ninputs)]
 
-    # Create the signal and system names
-    if inputs is None:
-        inputs = xd_labels + ud_labels + estimator.output_labels
-    if outputs is None:
-        outputs = list(sys.input_index.keys())
-    if states is None:
-        states = nintegrators
-
-    # Process gainscheduling variables, if present
-    if gainsched:
-        # Create a copy of the scheduling variable indices (default = xd)
-        gainsched_indices = range(sys.nstates) if gainsched_indices is None \
-            else list(gainsched_indices)
-
-        # Make sure the scheduling variable indices are the right length
-        if len(gainsched_indices) != points.shape[1]:
-            raise ControlArgument(
-                "length of gainsched_indices must match dimension of"
-                " scheduling variables")
-
-        # Process scheduling variables
-        for i, idx in enumerate(gainsched_indices):
-            if isinstance(idx, str):
-                gainsched_indices[i] = inputs.index(gainsched_indices[i])
-
-        # Create interpolating function
-        if gainsched_method == 'nearest':
-            _interp = sp.interpolate.NearestNDInterpolator(points, gains)
-            def _nearest(mu):
-                raise SystemError(f"could not find nearest gain at mu = {mu}")
-        elif gainsched_method == 'linear':
-            _interp = sp.interpolate.LinearNDInterpolator(points, gains)
-            _nearest = sp.interpolate.NearestNDInterpolator(points, gains)
-        elif gainsched_method == 'cubic':
-            _interp = sp.interpolate.CloughTocher2DInterpolator(points, gains)
-            _nearest = sp.interpolate.NearestNDInterpolator(points, gains)
-        else:
-            raise ControlArgument(
-                f"unknown gain scheduling method '{gainsched_method}'")
-
-        def _compute_gain(mu):
-            K = _interp(mu)
-            if np.isnan(K).any():
-                K = _nearest(mu)
-            return K
-
     # Define the controller system
     if type == 'nonlinear':
         # Create an I/O system for the state feedback gains
-        def _control_update(t, states, inputs, params):
+        def _control_update(t, x, inputs, params):
             # Split input into desired state, nominal input, and current state
             xd_vec = inputs[0:sys.nstates]
             x_vec = inputs[-estimator.nstates:]
 
             # Compute the integral error in the xy coordinates
-            return C @ (x_vec - xd_vec)
+            return C @ x_vec - C @ xd_vec
 
-        def _control_output(t, states, inputs, params):
-            if gainsched:
-                mu = inputs[gainsched_indices]
-                K_ = _compute_gain(mu)
-            else:
-                K_ = params.get('K')
+        def _control_output(t, e, z, params):
+            K = params.get('K')
 
             # Split input into desired state, nominal input, and current state
-            xd_vec = inputs[0:sys.nstates]
-            ud_vec = inputs[sys.nstates:sys.nstates + sys.ninputs]
-            x_vec = inputs[-sys.nstates:]
+            xd_vec = z[0:sys.nstates]
+            ud_vec = z[sys.nstates:sys.nstates + sys.ninputs]
+            x_vec = z[-sys.nstates:]
 
             # Compute the control law
-            u = ud_vec - K_[:, 0:sys.nstates] @ (x_vec - xd_vec)
+            u = ud_vec - K[:, 0:sys.nstates] @ (x_vec - xd_vec)
             if nintegrators > 0:
-                u -= K_[:, sys.nstates:] @ states
+                u -= K[:, sys.nstates:] @ e
 
             return u
 
-        params = {} if gainsched else {'K': K}
         ctrl = NonlinearIOSystem(
-            _control_update, _control_output, name=name, inputs=inputs,
-            outputs=outputs, states=states, params=params)
+            _control_update, _control_output, name='control',
+            inputs=xd_labels + ud_labels + estimator.output_labels,
+            outputs=list(sys.input_index.keys()), params={'K': K},
+            states=nintegrators)
 
     elif type == 'linear' or type is None:
         # Create the matrices implementing the controller
@@ -894,8 +768,9 @@ def create_statefbk_iosystem(
         ])
 
         ctrl = ss(
-            A_lqr, B_lqr, C_lqr, D_lqr, dt=sys.dt, name=name,
-            inputs=inputs, outputs=outputs, states=states)
+            A_lqr, B_lqr, C_lqr, D_lqr, dt=sys.dt, name='control',
+            inputs=xd_labels + ud_labels + estimator.output_labels,
+            outputs=list(sys.input_index.keys()), states=nintegrators)
 
     else:
         raise ControlArgument(f"unknown type '{type}'")
@@ -904,7 +779,7 @@ def create_statefbk_iosystem(
     closed = interconnect(
         [sys, ctrl] if estimator == sys else [sys, ctrl, estimator],
         name=sys.name + "_" + ctrl.name,
-        inplist=inputs[:-sys.nstates], inputs=inputs[:-sys.nstates],
+        inplist=xd_labels + ud_labels, inputs=xd_labels + ud_labels,
         outlist=sys.output_labels + sys.input_labels,
         outputs=sys.output_labels + sys.input_labels
     )
diff --git a/control/stochsys.py b/control/stochsys.py
index 2b823307..3b297800 100644
--- a/control/stochsys.py
+++ b/control/stochsys.py
@@ -20,7 +20,7 @@ import numpy as np
 import scipy as sp
 from math import sqrt
 
-from .iosys import InputOutputSystem, NonlinearIOSystem
+from .iosys import InputOutputSystem, LinearIOSystem, NonlinearIOSystem
 from .lti import LTI
 from .namedio import isctime, isdtime
 from .mateqn import care, dare, _check_shape
@@ -314,7 +314,13 @@ def create_estimator_iosystem(
     """Create an I/O system implementing a linqear quadratic estimator
 
     This function creates an input/output system that implements a
-    state estimator of the form
+    continuous time state estimator of the form
+
+        \dot xhat = A x + B u - L (C xhat - y)
+        \dot P = A P + P A^T + F QN F^T - P C^T RN^{-1} C P
+        L = P C^T RN^{-1}
+
+    or a discrete time state estimator of the form
 
         xhat[k + 1] = A x[k] + B u[k] - L (C xhat[k] - y[k])
         P[k + 1] = A P A^T + F QN F^T - A P C^T Reps^{-1} C P A
@@ -359,8 +365,9 @@ def create_estimator_iosystem(
     Returns
     -------
     estim : InputOutputSystem
-        Input/output system representing the estimator.  This system takes the
-        system input and output and generates the estimated state.
+        Input/output system representing the estimator.  This system takes
+        the system output y and input u and generates the estimated state
+        xhat.
 
     Notes
     -----
@@ -384,8 +391,8 @@ def create_estimator_iosystem(
     """
 
     # Make sure that we were passed an I/O system as an input
-    if not isinstance(sys, InputOutputSystem):
-        raise ControlArgument("Input system must be I/O system")
+    if not isinstance(sys, LinearIOSystem):
+        raise ControlArgument("Input system must be a linear I/O system")
 
     # Extract the matrices that we need for easy reference
     A, B = sys.A, sys.B
@@ -409,7 +416,7 @@ def create_estimator_iosystem(
     # Initialize the covariance matrix
     if P0 is None:
         # Initalize P0 to the steady state value
-        _, P0, _ = lqe(A, G, C, QN, RN)
+        L0, P0, _ = lqe(A, G, C, QN, RN)
 
     # Figure out the labels to use
     if isinstance(state_labels, str):
@@ -432,16 +439,46 @@ def create_estimator_iosystem(
         sensor_labels = [sensor_labels.format(i=i) for i in range(C.shape[0])]
 
     if isctime(sys):
-        raise NotImplementedError("continuous time not yet implemented")
-
-    else:
         # Create an I/O system for the state feedback gains
         # Note: reshape vectors into column vectors for legacy np.matrix
+
+        R_inv = np.linalg.inv(RN)
+        Reps_inv = C.T @ R_inv @ C
+
+        def _estim_update(t, x, u, params):
+            # See if we are estimating or predicting
+            correct = params.get('correct', True)
+
+            # Get the state of the estimator
+            xhat = x[0:sys.nstates].reshape(-1, 1)
+            P = x[sys.nstates:].reshape(sys.nstates, sys.nstates)
+
+            # Extract the inputs to the estimator
+            y = u[0:C.shape[0]].reshape(-1, 1)
+            u = u[C.shape[0]:].reshape(-1, 1)
+
+            # Compute the optimal gain
+            L = P @ C.T @ R_inv
+
+            # Update the state estimate
+            dxhat = A @ xhat + B @ u            # prediction
+            if correct:
+                dxhat -= L @ (C @ xhat - y)     # correction
+
+            # Update the covariance
+            dP = A @ P + P @ A.T + G @ QN @ G.T
+            if correct:
+                dP -= P @ Reps_inv @ P
+
+            # Return the update
+            return np.hstack([dxhat.reshape(-1), dP.reshape(-1)])
+
+    else:
         def _estim_update(t, x, u, params):
             # See if we are estimating or predicting
             correct = params.get('correct', True)
 
-            # Get the state of the estimator 
+            # Get the state of the estimator
             xhat = x[0:sys.nstates].reshape(-1, 1)
             P = x[sys.nstates:].reshape(sys.nstates, sys.nstates)
 
@@ -449,7 +486,7 @@ def create_estimator_iosystem(
             y = u[0:C.shape[0]].reshape(-1, 1)
             u = u[C.shape[0]:].reshape(-1, 1)
 
-            # Compute the optimal again
+            # Compute the optimal gain
             Reps_inv = np.linalg.inv(RN + C @ P @ C.T)
             L = A @ P @ C.T @ Reps_inv
 
@@ -466,8 +503,8 @@ def create_estimator_iosystem(
             # Return the update
             return np.hstack([dxhat.reshape(-1), dP.reshape(-1)])
 
-        def _estim_output(t, x, u, params):
-            return x[0:sys.nstates]
+    def _estim_output(t, x, u, params):
+        return x[0:sys.nstates]
 
     # Define the estimator system
     return NonlinearIOSystem(
diff --git a/doc/iosys.rst b/doc/iosys.rst
index 1f5f21e6..1da7f588 100644
--- a/doc/iosys.rst
+++ b/doc/iosys.rst
@@ -73,18 +73,18 @@ We begin by defining the dynamics of the system
       d = params.get('d', 0.56)
       k = params.get('k', 125)
       r = params.get('r', 1.6)
-
+      
       # Map the states into local variable names
       H = x[0]
       L = x[1]
 
       # Compute the control action (only allow addition of food)
       u_0 = u if u > 0 else 0
-
+  
       # Compute the discrete updates
       dH = (r + u_0) * H * (1 - H/k) - (a * H * L)/(c + H)
       dL = b * (a * H *  L)/(c + H) - d * L
-
+  
       return [dH, dL]
 
 We now create an input/output system using these dynamics:
@@ -105,10 +105,10 @@ of the system:
 
   X0 = [25, 20]                 # Initial H, L
   T = np.linspace(0, 70, 500)   # Simulation 70 years of time
-
+  
   # Simulate the system
   t, y = ct.input_output_response(io_predprey, T, 0, X0)
-
+  
   # Plot the response
   plt.figure(1)
   plt.plot(t, y[0])
@@ -168,14 +168,14 @@ function:
     inplist=['control.Ld'],
     outlist=['predprey.H', 'predprey.L', 'control.y[0]']
   )
-
+       
 Finally, we simulate the closed loop system:
 
 .. code-block:: python
 
   # Simulate the system
   t, y = ct.input_output_response(io_closed, T, 30, [15, 20])
-
+  
   # Plot the response
   plt.figure(2)
   plt.subplot(2, 1, 1)
@@ -198,7 +198,7 @@ Summing junction
 
 The :func:`~control.summing_junction` function can be used to create an
 input/output system that takes the sum of an arbitrary number of inputs.  For
-example, to create an input/output system that takes the sum of three inputs,
+ezample, to create an input/output system that takes the sum of three inputs,
 use the command
 
 .. code-block:: python
@@ -256,90 +256,6 @@ of the interconnected system) is not found, but inputs and outputs of
 individual systems that are not connected to other systems are left
 unconnected (so be careful!).
 
-Automated creation of state feedback systems
---------------------------------------------
-
-The :func:`~control.create_statefbk_iosystem` function can be used to
-create an I/O system consisting of a state feedback gain (with
-optional integral action and gain scheduling) and an estimator.  A
-basic state feedback controller of the form
-
-.. math::
-
-  u = u_\text{d} - K (x - x_\text{d})
-
-can be created with the command::
-
-  ctrl, clsys = ct.create_statefbk_iosystem(sys, K)
-
-where `sys` is the process dynamics and `K` is the state feedback gain
-(e.g., from LQR).  The function returns the controller `ctrl` and the
-closed loop systems `clsys`, both as I/O systems.  The input to the
-controller is the vector of desired states :math:`x_\text{d}`, desired
-inputs :math:`u_\text{d}`, and system states :math:`x`.
-
-If the full system state is not available, the output of a state
-estimator can be used to construct the controller using the command::
-
-  ctrl, clsys = ct.create_statefbk_iosystem(sys, K, estimator=estim)
-
-where `estim` is the state estimator I/O system.  The controller will
-have the same form as above, but with the system state :math:`x`
-replaced by the estimated state :math:`\hat x` (output of `estim`).
-The closed loop controller will include both the state feedback and
-the estimator.
-
-Integral action can be included using the `integral_action` keyword.
-The value of this keyword can either be an matrix (ndarray) or a
-function.  If a matrix :math:`C` is specified, the difference between
-the desired state and system state will be multiplied by this matrix
-and integrated.  The controller gain should then consist of a set of
-proportional gains :math:`K_\text{p}` and integral gains
-:math:`K_\text{i}` with
-
-.. math::
-
-   K = \begin{bmatrix} K_\text{p} \\ K_\text{i} \end{bmatrix}
-
-and the control action will be given by
-
-.. math::
-
-  u = u_\text{d} - K\text{p} (x - x_\text{d}) -
-      K_\text{i} \int C (x - x_\text{d}) dt.
-
-If `integral_action` is a function `h`, that function will be called
-with the signature `h(t, x, u, params)` to obtain the outputs that
-should be integrated.  The number of outputs that are to be integrated
-must match the number of additional columns in the `K` matrix.  If an
-estimator is specified, :math:`\hat x` will be used in place of
-:math:`x`.
-
-Finally, gain scheduling on the desired state, desired input, or
-system state can be implemented by setting the gain to a 2-tuple
-consisting of a list of gains and a list of points at which the gains
-were computed, as well as a description of the scheduling variables::
-
-  ctrl, clsys = ct.create_statefbk_iosystem(
-      sys, ([g1, ..., gN], [p1, ..., pN]), gainsched_indices=[s1, ..., sq])
-
-The list of indices can either be integers indicating the offset into
-the controller input vector :math:`(x_\text{d}, u_\text{d}, x)` or a
-list of strings matching the names of the input signals.  The
-controller implemented in this case has the form
-
-.. math::
-
-  u = u_\text{d} - K(\mu) (x - x_\text{d})
-
-where :math:`\mu` represents the scheduling variables.  See
-:ref:`steering-gainsched.py` for an example implementation of a gain
-scheduled controller (in the alternative formulation section at the
-bottom of the file).
-
-Integral action and state estimation can also be used with gain
-scheduled controllers.
-
 
 Module classes and functions
 ============================
diff --git a/doc/steering-gainsched.rst b/doc/steering-gainsched.rst
index a5ec2e0c..511f76b8 100644
--- a/doc/steering-gainsched.rst
+++ b/doc/steering-gainsched.rst
@@ -1,5 +1,3 @@
-.. _steering-gainsched.py:
-
 Gain scheduled control for vehicle steeering (I/O system)
 ---------------------------------------------------------
 
diff --git a/examples/steering-gainsched.py b/examples/steering-gainsched.py
index 88eed9a9..7ddc6b5b 100644
--- a/examples/steering-gainsched.py
+++ b/examples/steering-gainsched.py
@@ -70,7 +70,7 @@ def control_output(t, x, u, params):
     latpole1 = params.get('latpole1', -1/2 + sqrt(-7)/2)
     latpole2 = params.get('latpole2', -1/2 - sqrt(-7)/2)
     l = params.get('wheelbase', 3)
-
+    
     # Extract the system inputs
     ex, ey, etheta, vd, phid = u
 
@@ -85,7 +85,7 @@ def control_output(t, x, u, params):
     else:
         # We aren't moving, so don't turn the steering wheel
         phi = phid
-
+    
     return  np.array([v, phi])
 
 # Define the controller as an input/output system
@@ -135,7 +135,7 @@ trajgen = ct.NonlinearIOSystem(
 #                          +----------- [-1] -----------+
 #
 # We construct the system using the InterconnectedSystem constructor and using
-# signal labels to keep track of everything.
+# signal labels to keep track of everything.  
 
 steering = ct.interconnect(
     # List of subsystems
@@ -186,93 +186,8 @@ for vref in [8, 10, 12]:
     plt.plot([0, 5], [vref, vref], 'k--')
 
     # Plot the system output
-    y_line, = plt.plot(tout, yout[y_index], 'r')  # lateral position
-    v_line, = plt.plot(tout, yout[v_index], 'b')  # vehicle velocity
-
-# Add axis labels
-plt.xlabel('Time (s)')
-plt.ylabel('x vel (m/s), y pos (m)')
-plt.legend((v_line, y_line), ('v', 'y'), loc='center right', frameon=False)
-
-#
-# Alternative formulation, using create_statefbk_iosystem()
-#
-# A different way to implement gain scheduling is to use the gain scheduling
-# functionality built into the create_statefbk_iosystem() function, where we
-# pass a table of gains instead of a single gain.  To generate a more
-# interesting plot, we scale the feedforward input to generate some error.
-#
-
-import itertools
-from math import pi
-
-# Define the points for the scheduling variables
-speeds = [1, 10, 20]
-angles = np.linspace(-pi, pi, 4)
-points = list(itertools.product(speeds, angles))
-
-# Create controllers at each scheduling point
-Q = np.diag([1, 1, 1])
-R = np.diag([0.1, 0.1])
-gains = [np.array(ct.lqr(vehicle.linearize(
-    [0, 0, angle], [speed, 0]), Q, R)[0]) for speed, angle in points]
-
-# Create the gain scheduled system
-controller, _ = ct.create_statefbk_iosystem(
-    vehicle, (gains, points), name='controller', ud_labels=['vd', 'phid'],
-    gainsched_indices=['vd', 'theta'], gainsched_method='linear')
-
-# Connect everything together (note that controller inputs are different
-steering = ct.interconnect(
-    # List of subsystems
-    (trajgen, controller, vehicle), name='steering',
-
-    # Interconnections between  subsystems
-    connections=(
-        ['controller.xd[0]', 'trajgen.xd'],
-        ['controller.xd[1]', 'trajgen.yd'],
-        ['controller.xd[2]', 'trajgen.thetad'],
-        ['controller.x', 'vehicle.x'],
-        ['controller.y', 'vehicle.y'],
-        ['controller.theta', 'vehicle.theta'],
-        ['controller.vd', ('trajgen', 'vd', 0.2)],      # create error
-        ['controller.phid', 'trajgen.phid'],
-        ['vehicle.v', 'controller.v'],
-        ['vehicle.phi', 'controller.phi']
-    ),
-
-    # System inputs
-    inplist=['trajgen.vref', 'trajgen.yref'],
-    inputs=['yref', 'vref'],
-
-    #  System outputs
-    outlist=['vehicle.x', 'vehicle.y', 'vehicle.theta', 'controller.v',
-             'controller.phi'],
-    outputs=['x', 'y', 'theta', 'v', 'phi']
-)
-
-# Plot the results to compare to the previous case
-plt.figure();
-
-# Plot the reference trajectory for the y position
-plt.plot([0, 5], [yref, yref], 'k--')
-
-# Find the signals we want to plot
-y_index = steering.find_output('y')
-v_index = steering.find_output('v')
-
-# Do an iteration through different speeds
-for vref in [8, 10, 12]:
-    # Simulate the closed loop controller response
-    tout, yout = ct.input_output_response(
-        steering, T, [vref * np.ones(len(T)), yref * np.ones(len(T))])
-
-    # Plot the reference speed
-    plt.plot([0, 5], [vref, vref], 'k--')
-
-    # Plot the system output
-    y_line, = plt.plot(tout, yout[y_index], 'r')  # lateral position
-    v_line, = plt.plot(tout, yout[v_index], 'b')  # vehicle velocity
+    y_line, = plt.plot(tout, yout[y_index, :], 'r')  # lateral position
+    v_line, = plt.plot(tout, yout[v_index, :], 'b')  # vehicle velocity
 
 # Add axis labels
 plt.xlabel('Time (s)')


Test patch:
diff --git a/control/tests/flatsys_test.py b/control/tests/flatsys_test.py
index 196cddbd..95fb8cf7 100644
--- a/control/tests/flatsys_test.py
+++ b/control/tests/flatsys_test.py
@@ -212,15 +212,17 @@ class TestFlatSys:
             elif re.match("Iteration limit.*", traj_ocp.message) and \
                  re.match(
                      "conda ubuntu-3.* Generic", os.getenv('JOBNAME', '')) and \
-                 np.__version__ == '1.24.0':
+                 re.match("1.24.[01]", np.__version__):
                 pytest.xfail("gh820: iteration limit exceeded")
 
             else:
                 # Dump out information to allow creation of an exception
-                print("Platform: ", platform.platform())
-                print("Python: ", platform.python_version())
+                print("Message:", traj_ocp.message)
+                print("Platform:", platform.platform())
+                print("Python:", platform.python_version())
+                print("NumPy version:", np.__version__)
                 np.show_config()
-                print("JOBNAME: ", os.getenv('JOBNAME'))
+                print("JOBNAME:", os.getenv('JOBNAME'))
 
                 pytest.fail(
                     "unknown failure; view output to identify configuration")
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index fa26ded4..527f597b 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -1188,7 +1188,7 @@ class TestIOSys:
             assert "copy of namedsys.x0" in same_name_series.state_index
 
     def test_named_signals_linearize_inconsistent(self, tsys):
-        """Make sure that providing inputs or outputs not consistent with
+        """Mare sure that providing inputs or outputs not consistent with
            updfcn or outfcn fail
         """
 
@@ -1233,17 +1233,6 @@ class TestIOSys:
             with pytest.raises(ValueError):
                 sys2.linearize(x0, u0)
 
-    def test_linearize_concatenation(self, kincar):
-        # Create a simple nonlinear system to check (kinematic car)
-        iosys = kincar
-        linearized = iosys.linearize([0, np.array([0, 0])], [0, 0])
-        np.testing.assert_array_almost_equal(linearized.A, np.zeros((3,3)))
-        np.testing.assert_array_almost_equal(
-            linearized.B, [[1, 0], [0, 0], [0, 1]])
-        np.testing.assert_array_almost_equal(
-            linearized.C, [[1, 0, 0], [0, 1, 0]])
-        np.testing.assert_array_almost_equal(linearized.D, np.zeros((2,2)))
-
     def test_lineariosys_statespace(self, tsys):
         """Make sure that a LinearIOSystem is also a StateSpace object"""
         iosys_siso = ct.LinearIOSystem(tsys.siso_linsys, name='siso')
diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py
index 9e8feb4c..13f164e1 100644
--- a/control/tests/statefbk_test.py
+++ b/control/tests/statefbk_test.py
@@ -5,8 +5,6 @@ RMM, 30 Mar 2011 (based on TestStatefbk from v0.4a)
 
 import numpy as np
 import pytest
-import itertools
-from math import pi, atan
 
 import control as ct
 from control import lqe, dlqe, poles, rss, ss, tf
@@ -779,201 +777,3 @@ class TestStatefbk:
         with pytest.raises(ControlArgument, match="must be an array of size"):
             ctrl, clsys = ct.create_statefbk_iosystem(
                 sys, K, integral_action=C_int)
-
-
-# Kinematic car example for testing gain scheduling
-@pytest.fixture
-def unicycle():
-    # Create a simple nonlinear system to check (kinematic car)
-    def unicycle_update(t, x, u, params):
-        return np.array([np.cos(x[2]) * u[0], np.sin(x[2]) * u[0], u[1]])
-
-    def unicycle_output(t, x, u, params):
-        return x
-
-    return ct.NonlinearIOSystem(
-        unicycle_update, unicycle_output,
-        inputs = ['v', 'phi'],
-        outputs = ['x', 'y', 'theta'],
-        states = ['x_', 'y_', 'theta_'])
-
-from math import pi
-
-@pytest.mark.parametrize("method", ['nearest', 'linear', 'cubic'])
-def test_gainsched_unicycle(unicycle, method):
-    # Speeds and angles at which to compute the gains
-    speeds = [1, 5, 10]
-    angles = np.linspace(0, pi/2, 4)
-    points = list(itertools.product(speeds, angles))
-
-    # Gains for each speed (using LQR controller)
-    Q = np.identity(unicycle.nstates)
-    R = np.identity(unicycle.ninputs)
-    gains = [np.array(ct.lqr(unicycle.linearize(
-        [0, 0, angle], [speed, 0]), Q, R)[0]) for speed, angle in points]
-
-    #
-    # Schedule on desired speed and angle
-    #
-
-    # Create gain scheduled controller
-    ctrl, clsys = ct.create_statefbk_iosystem(
-        unicycle, (gains, points),
-        gainsched_indices=[3, 2], gainsched_method=method)
-
-    # Check the gain at the selected points
-    for speed, angle in points:
-        # Figure out the desired state and input
-        xe, ue = np.array([0, 0, angle]), np.array([speed, 0])
-        xd, ud = np.array([0, 0, angle]), np.array([speed, 0])
-
-        # Check the control system at the scheduling points
-        ctrl_lin = ctrl.linearize([], [xd, ud, xe*0])
-        K, S, E = ct.lqr(unicycle.linearize(xd, ud), Q, R)
-        np.testing.assert_allclose(
-            ctrl_lin.D[-xe.size:, -xe.size:], -K, rtol=1e-2)
-
-        # Check the closed loop system at the scheduling points
-        clsys_lin = clsys.linearize(xe, [xd, ud])
-        np.testing.assert_allclose(
-            np.sort(clsys_lin.poles()), np.sort(E), rtol=1e-2)
-
-    # Check the gain at an intermediate point and confirm stability
-    speed, angle = 2, pi/3
-    xe, ue = np.array([0, 0, angle]), np.array([speed, 0])
-    xd, ud = np.array([0, 0, angle]), np.array([speed, 0])
-    clsys_lin = clsys.linearize(xe, [xd, ud])
-    assert np.all(np.real(clsys_lin.poles()) < 0)
-
-    # Make sure that gains are different from 'nearest'
-    if method is not None and method != 'nearest':
-        ctrl_nearest, clsys_nearest = ct.create_statefbk_iosystem(
-            unicycle, (gains, points),
-            gainsched_indices=['ud[0]', 2], gainsched_method='nearest')
-        nearest_lin = clsys_nearest.linearize(xe, [xd, ud])
-        assert not np.allclose(
-            np.sort(clsys_lin.poles()), np.sort(nearest_lin.poles()), rtol=1e-2)
-
-    # Run a simulation following a curved path
-    T = 10                      # length of the trajectory [sec]
-    r = 10                      # radius of the circle [m]
-    timepts = np.linspace(0, T, 50)
-    Xd = np.vstack([
-        r * np.cos(timepts/T * pi/2 + 3*pi/2),
-        r * np.sin(timepts/T * pi/2 + 3*pi/2) + r,
-        timepts/T * pi/2
-    ])
-    Ud = np.vstack([
-        np.ones_like(timepts) * (r * pi/2) / T,
-        np.ones_like(timepts) * (pi / 2) / T
-    ])
-    X0 = Xd[:, 0] + np.array([-0.1, -0.1, -0.1])
-
-    resp = ct.input_output_response(clsys, timepts, [Xd, Ud], X0)
-    # plt.plot(resp.states[0], resp.states[1])
-    np.testing.assert_allclose(
-        resp.states[:, -1], Xd[:, -1], atol=1e-2, rtol=1e-2)
-
-    #
-    # Schedule on actual speed
-    #
-
-    # Create gain scheduled controller
-    ctrl, clsys = ct.create_statefbk_iosystem(
-        unicycle, (gains, points),
-        ud_labels=['vd', 'phid'], gainsched_indices=['vd', 'theta'])
-
-    # Check the gain at the selected points
-    for speed, angle in points:
-        # Figure out the desired state and input
-        xe, ue = np.array([0, 0, angle]), np.array([speed, 0])
-        xd, ud = np.array([0, 0, angle]), np.array([speed, 0])
-
-        # Check the control system at the scheduling points
-        ctrl_lin = ctrl.linearize([], [xd*0, ud, xe])
-        K, S, E = ct.lqr(unicycle.linearize(xe, ue), Q, R)
-        np.testing.assert_allclose(
-            ctrl_lin.D[-xe.size:, -xe.size:], -K, rtol=1e-2)
-
-        # Check the closed loop system at the scheduling points
-        clsys_lin = clsys.linearize(xe, [xd, ud])
-        np.testing.assert_allclose(np.sort(
-            clsys_lin.poles()), np.sort(E), rtol=1e-2)
-
-    # Run a simulation following a curved path
-    resp = ct.input_output_response(clsys, timepts, [Xd, Ud], X0)
-    np.testing.assert_allclose(
-        resp.states[:, -1], Xd[:, -1], atol=1e-2, rtol=1e-2)
-
-
-def test_gainsched_default_indices():
-    # Define a linear system to test
-    sys = ct.ss([[-1, 0.1], [0, -2]], [[0], [1]], np.eye(2), 0)
-
-    # Define gains at origin + corners of unit cube
-    points = [[0, 0]] + list(itertools.product([-1, 1], [-1, 1]))
-
-    # Define gain to be constant
-    K, _, _ = ct.lqr(sys, np.eye(sys.nstates), np.eye(sys.ninputs))
-    gains = [K for p in points]
-
-    # Define the paramters for the simulations
-    timepts = np.linspace(0, 10, 100)
-    X0 = np.ones(sys.nstates) * 0.9
-
-    # Create a controller and simulate the initial response
-    gs_ctrl, gs_clsys = ct.create_statefbk_iosystem(sys, (gains, points))
-    gs_resp = ct.input_output_response(gs_clsys, timepts, 0, X0)
-
-    # Verify that we get the same result as a constant gain
-    ck_clsys = ct.ss(sys.A - sys.B @ K, sys.B, sys.C, 0)
-    ck_resp = ct.input_output_response(ck_clsys, timepts, 0, X0)
-
-    np.testing.assert_allclose(gs_resp.states, ck_resp.states)
-
-
-def test_gainsched_errors(unicycle):
-    # Set up gain schedule (same as previous test)
-    speeds = [1, 5, 10]
-    angles = np.linspace(0, pi/2, 4)
-    points = list(itertools.product(speeds, angles))
-
-    Q = np.identity(unicycle.nstates)
-    R = np.identity(unicycle.ninputs)
-    gains = [np.array(ct.lqr(unicycle.linearize(
-        [0, 0, angle], [speed, 0]), Q, R)[0]) for speed, angle in points]
-
-    # Make sure the generic case works OK
-    ctrl, clsys = ct.create_statefbk_iosystem(
-        unicycle, (gains, points), gainsched_indices=[3, 2])
-    xd, ud = np.array([0, 0, angles[0]]), np.array([speeds[0], 0])
-    ctrl_lin = ctrl.linearize([], [xd, ud, xd*0])
-    K, S, E = ct.lqr(unicycle.linearize(xd, ud), Q, R)
-    np.testing.assert_allclose(
-        ctrl_lin.D[-xd.size:, -xd.size:], -K, rtol=1e-2)
-
-    # Wrong type of gain schedule argument
-    with pytest.raises(ControlArgument, match="gain must be an array"):
-        ctrl, clsys = ct.create_statefbk_iosystem(
-            unicycle, [gains, points], gainsched_indices=[3, 2])
-
-    # Wrong number of gain schedule argument
-    with pytest.raises(ControlArgument, match="gain must be a 2-tuple"):
-        ctrl, clsys = ct.create_statefbk_iosystem(
-            unicycle, (gains, speeds, angles), gainsched_indices=[3, 2])
-
-    # Mismatched dimensions for gains and points
-    with pytest.raises(ControlArgument, match="length of gainsched_indices"):
-        ctrl, clsys = ct.create_statefbk_iosystem(
-            unicycle, (gains, [speeds]), gainsched_indices=[3, 2])
-
-    # Unknown gain scheduling variable label
-    with pytest.raises(ValueError, match=".* not in list"):
-        ctrl, clsys = ct.create_statefbk_iosystem(
-            unicycle, (gains, points), gainsched_indices=['stuff', 2])
-
-    # Unknown gain scheduling method
-    with pytest.raises(ControlArgument, match="unknown gain scheduling method"):
-        ctrl, clsys = ct.create_statefbk_iosystem(
-            unicycle, (gains, points),
-            gainsched_indices=[3, 2], gainsched_method='unknown')
diff --git a/control/tests/stochsys_test.py b/control/tests/stochsys_test.py
index 11084d9d..75e5a510 100644
--- a/control/tests/stochsys_test.py
+++ b/control/tests/stochsys_test.py
@@ -7,6 +7,7 @@ from control.tests.conftest import asmatarrayout
 
 import control as ct
 from control import lqe, dlqe, rss, drss, tf, ss, ControlArgument, slycot_check
+from math import log, pi
 
 # Utility function to check LQE answer
 def check_LQE(L, P, poles, G, QN, RN):
@@ -48,7 +49,7 @@ def test_lqe_call_format(cdlqe):
 
     # Standard calling format
     Lref, Pref, Eref = cdlqe(sys.A, sys.B, sys.C, Q, R)
-    
+
     # Call with system instead of matricees
     L, P, E = cdlqe(sys, Q, R)
     np.testing.assert_almost_equal(Lref, L)
@@ -58,15 +59,15 @@ def test_lqe_call_format(cdlqe):
     # Make sure we get an error if we specify N
     with pytest.raises(ct.ControlNotImplemented):
         L, P, E = cdlqe(sys, Q, R, N)
-        
+
     # Inconsistent system dimensions
     with pytest.raises(ct.ControlDimension, match="Incompatible"):
         L, P, E = cdlqe(sys.A, sys.C, sys.B, Q, R)
-        
+
     # Incorrect covariance matrix dimensions
     with pytest.raises(ct.ControlDimension, match="Incompatible"):
         L, P, E = cdlqe(sys.A, sys.B, sys.C, R, Q)
-        
+
     # Too few input arguments
     with pytest.raises(ct.ControlArgument, match="not enough input"):
         L, P, E = cdlqe(sys.A, sys.C)
@@ -99,14 +100,14 @@ def test_lqe_discrete():
     np.testing.assert_almost_equal(K_csys, K_expl)
     np.testing.assert_almost_equal(S_csys, S_expl)
     np.testing.assert_almost_equal(E_csys, E_expl)
-    
+
     # Calling lqe() with a discrete time system should call dlqe()
     K_lqe, S_lqe, E_lqe = ct.lqe(dsys, Q, R)
     K_dlqe, S_dlqe, E_dlqe = ct.dlqe(dsys, Q, R)
     np.testing.assert_almost_equal(K_lqe, K_dlqe)
     np.testing.assert_almost_equal(S_lqe, S_dlqe)
     np.testing.assert_almost_equal(E_lqe, E_dlqe)
-    
+
     # Calling lqe() with no timebase should call lqe()
     asys = ct.ss(csys.A, csys.B, csys.C, csys.D, dt=None)
     K_asys, S_asys, E_asys = ct.lqe(asys, Q, R)
@@ -114,11 +115,11 @@ def test_lqe_discrete():
     np.testing.assert_almost_equal(K_asys, K_expl)
     np.testing.assert_almost_equal(S_asys, S_expl)
     np.testing.assert_almost_equal(E_asys, E_expl)
-    
+
     # Calling dlqe() with a continuous time system should raise an error
     with pytest.raises(ControlArgument, match="called with a continuous"):
         K, S, E = ct.dlqe(csys, Q, R)
-        
+
 def test_estimator_iosys():
     sys = ct.drss(4, 2, 2, strictly_proper=True)
 
@@ -129,7 +130,7 @@ def test_estimator_iosys():
     QN = np.eye(sys.ninputs)
     RN = np.eye(sys.noutputs)
     estim = ct.create_estimator_iosystem(sys, QN, RN, P0)
-        
+
     ctrl, clsys = ct.create_statefbk_iosystem(sys, K, estimator=estim)
 
     # Extract the elements of the estimator
@@ -162,20 +163,76 @@ def test_estimator_iosys():
     np.testing.assert_almost_equal(cls.D, D_clchk)
 
 
+@pytest.mark.parametrize("sys_args", [
+    ([[-1]], [[1]], [[1]], 0),                          # scalar system
+    ([[-1, 0.1], [0, -2]], [[0], [1]], [[1, 0]], 0),    # SISO, 2 state
+    ([[-1, 0.1], [0, -2]], [[1, 0], [0, 1]], [[1, 0]], 0),    # 2i, 1o, 2s
+    ([[-1, 0.1, 0.1], [0, -2, 0], [0.1, 0, -3]],        # 2i, 2o, 3s
+     [[1, 0], [0, 0.1], [0, 1]],
+     [[1, 0, 0.1], [0, 1, 0.1]], 0),
+])
+def test_estimator_iosys_ctime(sys_args):
+    # Define the system we want to test
+    sys = ct.ss(*sys_args)
+    T = 10 * log(1e-2) / np.max(sys.poles().real)
+    assert T > 0
+
+    # Create nonlinear version of the system to match integration methods
+    nl_sys = ct.NonlinearIOSystem(
+        lambda t, x, u, params : sys.A @ x + sys.B @ u,
+        lambda t, x, u, params : sys.C @ x + sys.D @ u,
+        inputs=sys.ninputs, outputs=sys.noutputs, states=sys.nstates)
+
+    # Define an initial condition, inputs (small, to avoid integration errors)
+    timepts = np.linspace(0, T, 500)
+    U = 2e-2 * np.array([np.sin(timepts + i*pi/3) for i in range(sys.ninputs)])
+    X0 = np.ones(sys.nstates)
+
+    # Set up the parameters for the filter
+    P0 = np.eye(sys.nstates)
+    QN = np.eye(sys.ninputs)
+    RN = np.eye(sys.noutputs)
+
+    # Construct the estimator
+    estim = ct.create_estimator_iosystem(sys, QN, RN)
+
+    # Compute the system response and the optimal covariance
+    sys_resp = ct.input_output_response(nl_sys, timepts, U, X0)
+    _, Pf, _ = ct.lqe(sys, QN, RN)
+    Pf = np.array(Pf)           # convert from matrix, if needed
+
+    # Make sure that we converge to the optimal estimate
+    estim_resp = ct.input_output_response(
+        estim, timepts, [sys_resp.outputs, U], [0*X0, P0])
+    np.testing.assert_allclose(
+        estim_resp.states[0:sys.nstates, -1], sys_resp.states[:, -1],
+        atol=1e-6, rtol=1e-3)
+    np.testing.assert_allclose(
+        estim_resp.states[sys.nstates:, -1], Pf.reshape(-1),
+        atol=1e-6, rtol=1e-3)
+
+    # Make sure that optimal estimate is an eq pt
+    ss_resp = ct.input_output_response(
+        estim, timepts, [sys_resp.outputs, U], [X0, Pf])
+    np.testing.assert_allclose(
+        ss_resp.states[sys.nstates:],
+        np.outer(Pf.reshape(-1), np.ones_like(timepts)),
+        atol=1e-4, rtol=1e-2)
+    np.testing.assert_allclose(
+        ss_resp.states[0:sys.nstates], sys_resp.states,
+        atol=1e-4, rtol=1e-2)
+
+
 def test_estimator_errors():
     sys = ct.drss(4, 2, 2, strictly_proper=True)
     P0 = np.eye(sys.nstates)
     QN = np.eye(sys.ninputs)
     RN = np.eye(sys.noutputs)
 
-    with pytest.raises(ct.ControlArgument, match="Input system must be I/O"):
+    with pytest.raises(ct.ControlArgument, match=".* system must be a linear"):
         sys_tf = ct.tf([1], [1, 1], dt=True)
         estim = ct.create_estimator_iosystem(sys_tf, QN, RN)
-            
-    with pytest.raises(NotImplementedError, match="continuous time not"):
-        sys_ct = ct.rss(4, 2, 2, strictly_proper=True)
-        estim = ct.create_estimator_iosystem(sys_ct, QN, RN)
-            
+
     with pytest.raises(ValueError, match="output must be full state"):
         C = np.eye(2, 4)
         estim = ct.create_estimator_iosystem(sys, QN, RN, C=C)
@@ -246,7 +303,7 @@ def test_correlation():
     # Try passing a second argument
     tau, Rneg = ct.correlation(T, V, -V)
     np.testing.assert_equal(Rtau, -Rneg)
-    
+
     # Test error conditions
     with pytest.raises(ValueError, match="Time vector T must be 1D"):
         tau, Rtau = ct.correlation(V, V)