Repository: python-control/python-control
Pull request #1024 (2024-07-04 13:59:59)

Implement ERA, change api to TimeResponseData
This PR implements the ERA method. Related to #1002 and #847

- it uses the TimeResponseData Type.


coveralls: 
[![Coverage Status](https://coveralls.io/builds/68469147/badge)](https://coveralls.io/builds/68469147)

coverage: 94.21% (-0.3%) from 94.533%
when pulling **28744f88c691ddf44242a82cbd802a6f5d59dac2 on KybernetikJo:implement_era**
into **4acc78b8aec4f145e7e72c931accfa7e12e5cd4d on python-control:main**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/68472080/badge)](https://coveralls.io/builds/68472080)

coverage: 94.21% (-0.3%) from 94.533%
when pulling **596e5ed7019ae94322fc8e2e1de9e642efe22cd3 on KybernetikJo:implement_era**
into **4acc78b8aec4f145e7e72c931accfa7e12e5cd4d on python-control:main**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/68499241/badge)](https://coveralls.io/builds/68499241)

coverage: 94.374% (-0.2%) from 94.533%
when pulling **55ac0ed10ef05611a52a4830f355c9c8914e5465 on KybernetikJo:implement_era**
into **4acc78b8aec4f145e7e72c931accfa7e12e5cd4d on python-control:main**.


Patch:
diff --git a/control/bdalg.py b/control/bdalg.py
index 7bfd327e..ce800853 100644
--- a/control/bdalg.py
+++ b/control/bdalg.py
@@ -54,19 +54,17 @@ $Id$
 """
 
 from functools import reduce
-from warnings import warn
-
 import numpy as np
-
-from . import frdata as frd
-from . import statesp as ss
+from warnings import warn
 from . import xferfcn as tf
+from . import statesp as ss
+from . import frdata as frd
 from .iosys import InputOutputSystem
 
 __all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect']
 
 
-def series(sys1, *sysn, **kwargs):
+def series(sys1, *sysn):
     r"""series(sys1, sys2, [..., sysn])
 
     Return the series connection (`sysn` \* ...\  \*) `sys2` \* `sys1`.
@@ -81,20 +79,6 @@ def series(sys1, *sysn, **kwargs):
     out : scalar, array, or :class:`InputOutputSystem`
         Series interconnection of the systems.
 
-    Other Parameters
-    ----------------
-    inputs, outputs : str, or list of str, optional
-        List of strings that name the individual signals.  If not given,
-        signal names will be of the form `s[i]` (where `s` is one of `u`,
-        or `y`). See :class:`InputOutputSystem` for more information.
-    states : str, or list of str, optional
-        List of names for system states.  If not given, state names will be
-        of of the form `x[i]` for interconnections of linear systems or
-        '<subsys_name>.<state_name>' for interconnected nonlinear systems.
-    name : string, optional
-        System name (used for specifying signals). If unspecified, a generic
-        name <sys[id]> is generated with a unique integer id.
-
     Raises
     ------
     ValueError
@@ -133,12 +117,10 @@ def series(sys1, *sysn, **kwargs):
     (2, 1, 5)
 
     """
-    sys = reduce(lambda x, y: y * x, sysn, sys1)
-    sys.update_names(**kwargs)
-    return sys
+    return reduce(lambda x, y: y * x, sysn, sys1)
 
 
-def parallel(sys1, *sysn, **kwargs):
+def parallel(sys1, *sysn):
     r"""parallel(sys1, sys2, [..., sysn])
 
     Return the parallel connection `sys1` + `sys2` (+ ...\  + `sysn`).
@@ -153,20 +135,6 @@ def parallel(sys1, *sysn, **kwargs):
     out : scalar, array, or :class:`InputOutputSystem`
         Parallel interconnection of the systems.
 
-    Other Parameters
-    ----------------
-    inputs, outputs : str, or list of str, optional
-        List of strings that name the individual signals.  If not given,
-        signal names will be of the form `s[i]` (where `s` is one of `u`,
-        or `y`). See :class:`InputOutputSystem` for more information.
-    states : str, or list of str, optional
-        List of names for system states.  If not given, state names will be
-        of of the form `x[i]` for interconnections of linear systems or
-        '<subsys_name>.<state_name>' for interconnected nonlinear systems.
-    name : string, optional
-        System name (used for specifying signals). If unspecified, a generic
-        name <sys[id]> is generated with a unique integer id.
-
     Raises
     ------
     ValueError
@@ -203,11 +171,10 @@ def parallel(sys1, *sysn, **kwargs):
     (3, 4, 7)
 
     """
-    sys = reduce(lambda x, y: x + y, sysn, sys1)
-    sys.update_names(**kwargs)
-    return sys
+    return reduce(lambda x, y: x + y, sysn, sys1)
+
 
-def negate(sys, **kwargs):
+def negate(sys):
     """
     Return the negative of a system.
 
@@ -221,29 +188,15 @@ def negate(sys, **kwargs):
     out : scalar, array, or :class:`InputOutputSystem`
         Negated system.
 
-    Other Parameters
-    ----------------
-    inputs, outputs : str, or list of str, optional
-        List of strings that name the individual signals.  If not given,
-        signal names will be of the form `s[i]` (where `s` is one of `u`,
-        or `y`). See :class:`InputOutputSystem` for more information.
-    states : str, or list of str, optional
-        List of names for system states.  If not given, state names will be
-        of of the form `x[i]` for interconnections of linear systems or
-        '<subsys_name>.<state_name>' for interconnected nonlinear systems.
-    name : string, optional
-        System name (used for specifying signals). If unspecified, a generic
-        name <sys[id]> is generated with a unique integer id.
-
-    See Also
-    --------
-    append, feedback, interconnect, parallel, series
-
     Notes
     -----
     This function is a wrapper for the __neg__ function in the StateSpace and
     TransferFunction classes.  The output type is the same as the input type.
 
+    See Also
+    --------
+    append, feedback, interconnect, parallel, series
+
     Examples
     --------
     >>> G = ct.tf([2], [1, 1])
@@ -255,12 +208,11 @@ def negate(sys, **kwargs):
     np.float64(-2.0)
 
     """
-    sys = -sys
-    sys.update_names(**kwargs)
-    return sys
+    return -sys
 
 #! TODO: expand to allow sys2 default to work in MIMO case?
-def feedback(sys1, sys2=1, sign=-1, **kwargs):
+#! TODO: allow renaming of signals (for all bdalg operations)
+def feedback(sys1, sys2=1, sign=-1):
     """Feedback interconnection between two I/O systems.
 
     Parameters
@@ -277,20 +229,6 @@ def feedback(sys1, sys2=1, sign=-1, **kwargs):
     out : scalar, array, or :class:`InputOutputSystem`
         Feedback interconnection of the systems.
 
-    Other Parameters
-    ----------------
-    inputs, outputs : str, or list of str, optional
-        List of strings that name the individual signals.  If not given,
-        signal names will be of the form `s[i]` (where `s` is one of `u`,
-        or `y`). See :class:`InputOutputSystem` for more information.
-    states : str, or list of str, optional
-        List of names for system states.  If not given, state names will be
-        of of the form `x[i]` for interconnections of linear systems or
-        '<subsys_name>.<state_name>' for interconnected nonlinear systems.
-    name : string, optional
-        System name (used for specifying signals). If unspecified, a generic
-        name <sys[id]> is generated with a unique integer id.
-
     Raises
     ------
     ValueError
@@ -323,7 +261,7 @@ def feedback(sys1, sys2=1, sign=-1, **kwargs):
     # Allow anything with a feedback function to call that function
     # TODO: rewrite to allow __rfeedback__
     try:
-        return sys1.feedback(sys2, sign, **kwargs)
+        return sys1.feedback(sys2, sign)
     except (AttributeError, TypeError):
         pass
 
@@ -346,11 +284,9 @@ def feedback(sys1, sys2=1, sign=-1, **kwargs):
         else:
             sys1 = ss._convert_to_statespace(sys1)
 
-    sys = sys1.feedback(sys2, sign)
-    sys.update_names(**kwargs)
-    return sys
+    return sys1.feedback(sys2, sign)
 
-def append(*sys, **kwargs):
+def append(*sys):
     """append(sys1, sys2, [..., sysn])
 
     Group LTI state space models by appending their inputs and outputs.
@@ -363,20 +299,6 @@ def append(*sys, **kwargs):
     sys1, sys2, ..., sysn: scalar, array, or :class:`StateSpace`
         I/O systems to combine.
 
-    Other Parameters
-    ----------------
-    inputs, outputs : str, or list of str, optional
-        List of strings that name the individual signals.  If not given,
-        signal names will be of the form `s[i]` (where `s` is one of `u`,
-        or `y`). See :class:`InputOutputSystem` for more information.
-    states : str, or list of str, optional
-        List of names for system states.  If not given, state names will be
-        of of the form `x[i]` for interconnections of linear systems or
-        '<subsys_name>.<state_name>' for interconnected nonlinear systems.
-    name : string, optional
-        System name (used for specifying signals). If unspecified, a generic
-        name <sys[id]> is generated with a unique integer id.
-
     Returns
     -------
     out: :class:`StateSpace`
@@ -405,7 +327,6 @@ def append(*sys, **kwargs):
     s1 = ss._convert_to_statespace(sys[0])
     for s in sys[1:]:
         s1 = s1.append(s)
-    s1.update_names(**kwargs)
     return s1
 
 def connect(sys, Q, inputv, outputv):
@@ -449,12 +370,6 @@ def connect(sys, Q, inputv, outputv):
     --------
     append, feedback, interconnect, negate, parallel, series
 
-    Notes
-    -----
-    The :func:`~control.interconnect` function in the :ref:`input/output
-    systems <iosys-module>` module allows the use of named signals and
-    provides an alternative method for interconnecting multiple systems.
-
     Examples
     --------
     >>> G = ct.rss(7, inputs=2, outputs=2)
@@ -463,6 +378,12 @@ def connect(sys, Q, inputv, outputv):
     >>> T.ninputs, T.noutputs, T.nstates
     (1, 2, 7)
 
+    Notes
+    -----
+    The :func:`~control.interconnect` function in the :ref:`input/output
+    systems <iosys-module>` module allows the use of named signals and
+    provides an alternative method for interconnecting multiple systems.
+
     """
     # TODO: maintain `connect` for use in MATLAB submodule (?)
     warn("`connect` is deprecated; use `interconnect`", DeprecationWarning)
diff --git a/control/iosys.py b/control/iosys.py
index d00dade6..fbd5c1db 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -6,12 +6,10 @@
 # FrequencyResponseData, InterconnectedSystem and other similar classes
 # that allow naming of signals.
 
-import re
+import numpy as np
 from copy import deepcopy
 from warnings import warn
-
-import numpy as np
-
+import re
 from . import config
 
 __all__ = ['InputOutputSystem', 'issiso', 'timebase', 'common_timebase',
@@ -368,51 +366,6 @@ class InputOutputSystem(object):
         lambda self: list(self.state_index.keys()),     # getter
         set_states)                                     # setter
 
-    # TODO: add dict as a means to selective change names?  [GH #1019]
-    def update_names(self, **kwargs):
-        """update_names([name, inputs, outputs, states])
-
-        Update signal and system names for an I/O system.
-
-        Parameters
-        ----------
-        name : str, optional
-            New system name.
-        inputs : list of str, int, or None, optional
-            List of strings that name the individual input signals.  If
-            given as an integer or None, signal names default to the form
-            `u[i]`.  See :class:`InputOutputSystem` for more information.
-        outputs : list of str, int, or None, optional
-            Description of output signals; defaults to `y[i]`.
-        states : int, list of str, int, or None, optional
-            Description of system states; defaults to `x[i]`.
-
-        """
-        self.name = kwargs.pop('name', self.name)
-        if 'inputs' in kwargs:
-            ninputs, input_index = _process_signal_list(
-                kwargs.pop('inputs'), prefix=kwargs.pop('input_prefix', 'u'))
-            if self.ninputs and self.ninputs != ninputs:
-                raise ValueError("number of inputs does not match system size")
-            self.input_index = input_index
-        if 'outputs' in kwargs:
-            noutputs, output_index = _process_signal_list(
-                kwargs.pop('outputs'), prefix=kwargs.pop('output_prefix', 'y'))
-            if self.noutputs and self.noutputs != noutputs:
-                raise ValueError("number of outputs does not match system size")
-            self.output_index = output_index
-        if 'states' in kwargs:
-            nstates, state_index = _process_signal_list(
-                kwargs.pop('states'), prefix=kwargs.pop('state_prefix', 'x'))
-            if self.nstates != nstates:
-                raise ValueError("number of states does not match system size")
-            self.state_index = state_index
-
-        # Make sure we processed all of the arguments
-        if kwargs:
-            raise TypeError("unrecognized keywords: ", str(kwargs))
-
-
     def isctime(self, strict=False):
         """
         Check to see if a system is a continuous-time system.
@@ -870,6 +823,7 @@ def _process_labels(labels, name, default):
 # This function returns the subsystem index, a list of indices for the
 # system signals, and the gain to use for that set of signals.
 #
+import re
 
 def _parse_spec(syslist, spec, signame, dictname=None):
     """Parse a signal specification, returning system and signal index."""
diff --git a/control/modelsimp.py b/control/modelsimp.py
index 06c3d350..4b72e2bb 100644
--- a/control/modelsimp.py
+++ b/control/modelsimp.py
@@ -48,8 +48,9 @@ from .exception import ControlSlycot, ControlMIMONotImplemented, \
 from .iosys import isdtime, isctime
 from .statesp import StateSpace
 from .statefbk import gram
+from .timeresp import TimeResponseData
 
-__all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal']
+__all__ = ['hsvd', 'balred', 'modred', 'eigensys_realization', 'markov', 'minreal', 'era']
 
 
 # Hankel Singular Value Decomposition
@@ -368,38 +369,129 @@ def minreal(sys, tol=None, verbose=True):
     return sysr
 
 
-def era(YY, m, n, nin, nout, r):
-    """Calculate an ERA model of order `r` based on the impulse-response data
+def _block_hankel(Y, m, n):
+    """Create a block Hankel matrix from impulse response"""
+    q, p, _ = Y.shape
+    YY = Y.transpose(0,2,1) # transpose for reshape
+    
+    H = np.zeros((q*m,p*n))
+    
+    for r in range(m):
+        # shift and add row to Hankel matrix
+        new_row = YY[:,r:r+n,:]
+        H[q*r:q*(r+1),:] = new_row.reshape((q,p*n))
+            
+    return H
+
+
+def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
+    r"""eigensys_realization(YY, r)
+
+    Calculate an ERA model of order `r` based on the impulse-response data
     `YY`.
 
-    .. note:: This function is not implemented yet.
+    This function computes a discrete time system
+
+    .. math::
+
+        x[k+1] &= A x[k] + B u[k] \\\\
+        y[k] &= C x[k] + D u[k]
+
+    for a given impulse-response data (see [1]_).
+
+    The function can be called with 2 arguments:
+
+    * ``sysd, S = eigensys_realization(data, r)``
+    * ``sysd, S = eigensys_realization(YY, r)``
+
+    where `data` is a `TimeResponseData` object, `YY` is a 1D or 3D
+    array, and r is an integer.
 
     Parameters
     ----------
-    YY: array
-        `nout` x `nin` dimensional impulse-response data
-    m: integer
-        Number of rows in Hankel matrix
-    n: integer
-        Number of columns in Hankel matrix
-    nin: integer
-        Number of input variables
-    nout: integer
-        Number of output variables
-    r: integer
-        Order of model
+    YY : array_like
+        Impulse response from which the StateSpace model is estimated, 1D
+        or 3D array.
+    data : TimeResponseData
+        Impulse response from which the StateSpace model is estimated.
+    r : integer
+        Order of model.
+    m : integer, optional
+        Number of rows in Hankel matrix. Default is 2*r.
+    n : integer, optional
+        Number of columns in Hankel matrix. Default is 2*r.
+    dt : True or float, optional
+        True indicates discrete time with unspecified sampling time and a
+        positive float is discrete time with the specified sampling time.
+        It can be used to scale the StateSpace model in order to match the
+        unit-area impulse response of python-control. Default is True.
+    transpose : bool, optional
+        Assume that input data is transposed relative to the standard
+        :ref:`time-series-convention`. For TimeResponseData this parameter 
+        is ignored. Default is False.
 
     Returns
     -------
-    sys: StateSpace
-        A reduced order model sys=ss(Ar,Br,Cr,Dr)
+    sys : StateSpace
+        A reduced order model sys=StateSpace(Ar,Br,Cr,Dr,dt).
+    S : array
+        Singular values of Hankel matrix. Can be used to choose a good r
+        value.
+
+    References
+    ----------
+    .. [1] Samet Oymak and Necmiye Ozay, Non-asymptotic Identification of
+       LTI Systems from a Single Trajectory.
+       https://arxiv.org/abs/1806.05722
 
     Examples
     --------
-    >>> rsys = era(YY, m, n, nin, nout, r)                      # doctest: +SKIP
+    >>> T = np.linspace(0, 10, 100)
+    >>> _, YY = ct.impulse_response(ct.tf([1], [1, 0.5], True), T)
+    >>> sysd, _ = ct.eigensys_realization(YY, r=1)
 
+    >>> T = np.linspace(0, 10, 100)
+    >>> response = ct.impulse_response(ct.tf([1], [1, 0.5], True), T)
+    >>> sysd, _ = ct.eigensys_realization(response, r=1)
     """
-    raise NotImplementedError('This function is not implemented yet.')
+    if isinstance(arg, TimeResponseData):
+        YY = np.array(arg.outputs, ndmin=3)
+        if arg.transpose:
+            YY = np.transpose(YY)
+    else:
+        YY = np.array(arg, ndmin=3)
+        if transpose:
+            YY = np.transpose(YY)
+    
+    q, p, l = YY.shape
+
+    if m is None:
+        m = 2*r
+    if n is None:
+        n = 2*r
+
+    if m*q < r or n*p < r:
+        raise ValueError("Hankel parameters are to small")
+    
+    if (l-1) < m+n:
+        raise ValueError("not enough data for requested number of parameters")
+    
+    H = _block_hankel(YY[:,:,1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
+    Hf = H[:,:-p] # first p*n columns of H
+    Hl = H[:,p:] # last p*n columns of H
+    
+    U,S,Vh = np.linalg.svd(Hf, True)
+    Ur =U[:,0:r]
+    Vhr =Vh[0:r,:]
+
+    # balanced realizations
+    Sigma_inv = np.diag(1./np.sqrt(S[0:r]))
+    Ar = Sigma_inv @ Ur.T @ Hl @ Vhr.T @ Sigma_inv
+    Br = Sigma_inv @ Ur.T @ Hf[:,0:p]*dt # dt scaling for unit-area impulse
+    Cr = Hf[0:q,:] @ Vhr.T @ Sigma_inv
+    Dr = YY[:,:,0]
+
+    return StateSpace(Ar,Br,Cr,Dr,dt), S
 
 
 def markov(Y, U, m=None, transpose=False):
@@ -556,3 +648,6 @@ def markov(Y, U, m=None, transpose=False):
 
     # Return the first m Markov parameters
     return H if transpose else np.transpose(H)
+
+# Function aliases
+era = eigensys_realization
diff --git a/control/nlsys.py b/control/nlsys.py
index d6d3b1b7..a82d3328 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -143,9 +143,6 @@ class NonlinearIOSystem(InputOutputSystem):
                 self.outfcn = lambda t, x, u, params: np.zeros(0)
             elif self.noutputs is None and self.nstates is not None:
                 self.noutputs = self.nstates
-                if len(self.output_index) == 0:
-                    # Use state names for outputs
-                    self.output_index = self.state_index
             elif self.noutputs is not None and self.noutputs == self.nstates:
                 # Number of outputs = number of states => all is OK
                 pass
@@ -517,7 +514,7 @@ class NonlinearIOSystem(InputOutputSystem):
         return newsys
 
     def linearize(self, x0, u0, t=0, params=None, eps=1e-6,
-                  copy_names=False, **kwargs):
+                  name=None, copy_names=False, **kwargs):
         """Linearize an input/output system at a given state and input.
 
         Return the linearization of an input/output system at a given state
@@ -532,16 +529,26 @@ class NonlinearIOSystem(InputOutputSystem):
         # numerical linearization use the `_rhs()` and `_out()` member
         # functions.
         #
-        # Process nominal states and inputs
-        x0, nstates = _process_vector_argument(x0, "x0", self.nstates)
-        u0, ninputs = _process_vector_argument(u0, "u0", self.ninputs)
+        # If x0 and u0 are specified as lists, concatenate the elements
+        x0 = _concatenate_list_elements(x0, 'x0')
+        u0 = _concatenate_list_elements(u0, 'u0')
 
-        # Update the current parameters (prior to calling _out())
-        self._update_params(params)
+        # Figure out dimensions if they were not specified.
+        nstates = _find_size(self.nstates, x0, "states")
+        ninputs = _find_size(self.ninputs, u0, "inputs")
+
+        # Convert x0, u0 to arrays, if needed
+        if np.isscalar(x0):
+            x0 = np.ones((nstates,)) * x0
+        if np.isscalar(u0):
+            u0 = np.ones((ninputs,)) * u0
 
         # Compute number of outputs by evaluating the output function
         noutputs = _find_size(self.noutputs, self._out(t, x0, u0), "outputs")
 
+        # Update the current parameters
+        self._update_params(params)
+
         # Compute the nominal value of the update law and output
         F0 = self._rhs(t, x0, u0)
         H0 = self._out(t, x0, u0)
@@ -572,6 +579,8 @@ class NonlinearIOSystem(InputOutputSystem):
         # Set the system name, inputs, outputs, and states
         if copy_names:
             linsys._copy_names(self, prefix_suffix_name='linearized')
+            if name is not None:
+                linsys.name = name
 
         # re-init to include desired signal names if names were provided
         return StateSpace(linsys, **kwargs)
@@ -1468,16 +1477,7 @@ def input_output_response(
         # Use the input time points as the output time points
         t_eval = T
 
-    #
-    # Process input argument
-    #
-    # The input argument is interpreted very flexibly, allowing the
-    # use of lists and/or tuples of mixed scalar and vector elements.
-    #
-    # Much of the processing here is similar to the processing in
-    # _process_vector_argument, but applied to a time series.
-
-    # If we were passed a list of inputs, concatenate them (w/ broadcast)
+    # If we were passed a list of input, concatenate them (w/ broadcast)
     if isinstance(U, (tuple, list)) and len(U) != ntimepts:
         U_elements = []
         for i, u in enumerate(U):
@@ -1501,44 +1501,62 @@ def input_output_response(
         # Save the newly created input vector
         U = np.vstack(U_elements)
 
-    # Figure out the number of inputs
-    if sys.ninputs is None:
-        if isinstance(U, np.ndarray):
-            ninputs = U.shape[0] if U.size > 1 else U.size
-        else:
-            ninputs = 1
-    else:
-        ninputs = sys.ninputs
-
     # Make sure the input has the right shape
-    if ninputs is None or ninputs == 1:
+    if sys.ninputs is None or sys.ninputs == 1:
         legal_shapes = [(ntimepts,), (1, ntimepts)]
     else:
-        legal_shapes = [(ninputs, ntimepts)]
+        legal_shapes = [(sys.ninputs, ntimepts)]
 
-    U = _check_convert_array(
-        U, legal_shapes, 'Parameter ``U``: ', squeeze=False)
+    U = _check_convert_array(U, legal_shapes,
+                             'Parameter ``U``: ', squeeze=False)
 
     # Always store the input as a 2D array
     U = U.reshape(-1, ntimepts)
     ninputs = U.shape[0]
 
-    # Process initial states
-    X0, nstates = _process_vector_argument(X0, "X0", sys.nstates)
+    # If we were passed a list of initial states, concatenate them
+    X0 = _concatenate_list_elements(X0, 'X0')
 
-    # Update the parameter values (prior to evaluating outfcn)
-    sys._update_params(params)
+    # If the initial state is too short, make it longer (NB: sys.nstates
+    # could be None if nstates comes from size of initial condition)
+    if sys.nstates and isinstance(X0, np.ndarray) and X0.size < sys.nstates:
+        if X0[-1] != 0:
+            warn("initial state too short; padding with zeros")
+        X0 = np.hstack([X0, np.zeros(sys.nstates - X0.size)])
+
+    # If we were passed a list of initial states, concatenate them
+    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)
+    if sys.nstates and isinstance(X0, np.ndarray) and X0.size < sys.nstates:
+        if X0[-1] != 0:
+            warn("initial state too short; padding with zeros")
+        X0 = np.hstack([X0, np.zeros(sys.nstates - X0.size)])
+
+    # Compute the number of states
+    nstates = _find_size(sys.nstates, X0, "states")
+
+    # create X0 if not given, test if X0 has correct shape
+    X0 = _check_convert_array(X0, [(nstates,), (nstates, 1)],
+                              'Parameter ``X0``: ', squeeze=True)
 
     # Figure out the number of outputs
-    if sys.outfcn is None:
-        noutputs = nstates if sys.noutputs is None else sys.noutputs
-    else:
+    if sys.noutputs is None:
+        # Evaluate the output function to find number of outputs
         noutputs = np.shape(sys._out(T[0], X0, U[:, 0]))[0]
+    else:
+        noutputs = sys.noutputs
 
-    if sys.noutputs is not None and sys.noutputs != noutputs:
-        raise RuntimeError(
-            f"inconsistent size of outputs; system specified {sys.noutputs}, "
-            f"output function returned {noutputs}")
+    # Update the parameter values
+    sys._update_params(params)
 
     #
     # Define a function to evaluate the input at an arbitrary time
@@ -1745,9 +1763,17 @@ def find_eqpt(sys, x0, u0=None, y0=None, t=0, params=None,
     from scipy.optimize import root
 
     # Figure out the number of states, inputs, and outputs
-    x0, nstates = _process_vector_argument(x0, "x0", sys.nstates)
-    u0, ninputs = _process_vector_argument(u0, "u0", sys.ninputs)
-    y0, noutputs = _process_vector_argument(y0, "y0", sys.noutputs)
+    nstates = _find_size(sys.nstates, x0, "states")
+    ninputs = _find_size(sys.ninputs, u0, "inputs")
+    noutputs = _find_size(sys.noutputs, y0, "outputs")
+
+    # Convert x0, u0, y0 to arrays, if needed
+    if np.isscalar(x0):
+        x0 = np.ones((nstates,)) * x0
+    if np.isscalar(u0):
+        u0 = np.ones((ninputs,)) * u0
+    if np.isscalar(y0):
+        y0 = np.ones((ninputs,)) * y0
 
     # Make sure the input arguments match the sizes of the system
     if len(x0) != nstates or \
@@ -1982,7 +2008,7 @@ def linearize(sys, xeq, ueq=None, t=0, params=None, **kw):
     return sys.linearize(xeq, ueq, t=t, params=params, **kw)
 
 
-def _find_size(sysval, vecval, name="system component"):
+def _find_size(sysval, vecval, label):
     """Utility function to find the size of a system parameter
 
     If both parameters are not None, they must be consistent.
@@ -1990,8 +2016,7 @@ def _find_size(sysval, vecval, name="system component"):
     if hasattr(vecval, '__len__'):
         if sysval is not None and sysval != len(vecval):
             raise ValueError(
-                f"inconsistent information to determine size of {name}; "
-                f"expected {sysval} values, received {len(vecval)}")
+                f"inconsistent information for number of {label}")
         return len(vecval)
     # None or 0, which is a valid value for "a (sysval, ) vector of zeros".
     if not vecval:
@@ -1999,7 +2024,7 @@ def _find_size(sysval, vecval, name="system component"):
     elif sysval == 1:
         # (1, scalar) is also a valid combination from legacy code
         return 1
-    raise ValueError(f"can't determine size of {name}")
+    raise ValueError(f"can't determine number of {label}")
 
 
 # Function to create an interconnected system
@@ -2571,59 +2596,18 @@ def interconnect(
     return newsys
 
 
-def _process_vector_argument(arg, name, size):
-    """Utility function to process vector elements (states, inputs)
-
-    Process state and input arguments to turn them into lists of the
-    appropriate length.
-
-    Parameters
-    ----------
-    arg : array_like
-        Value of the parameter passed to the function.  Can be a list,
-        tuple, ndarray, scalar, or None.
-    name : string
-        Name of the argument being processed.  Used in errors/warnings.
-    size : int or None
-        Size of the element.  If None, size is determined by arg.
-
-    Returns
-    -------
-    val : array or None
-        Value of the element, zero-padded to proper length.
-    nelem : int or None
-        Number of elements in the returned value.
-
-    Warns
-    -----
-    UserWarning : "{name} too short; padding with zeros"
-        If argument is too short and last value in arg is not 0.
-
-    """
-    # Allow and expand list
-    if isinstance(arg, (tuple, list)):
-        val_list = []
-        for i, v in enumerate(arg):
-            v = np.array(v).reshape(-1)             # convert to 1D array
-            val_list += v.tolist()                  # add elements to list
-        val = np.array(val_list)
-    elif np.isscalar(arg) and size is not None:     # extend scalars
-        val = np.ones((size, )) * arg
-    elif np.isscalar(arg) and size is None:         # single scalar
-        val = np.array([arg])
-    elif isinstance(arg, np.ndarray):
-        val = arg.reshape(-1)                       # convert to 1D array
-    else:
-        val = arg                                   # return what we were given
-
-    if size is not None and isinstance(val, np.ndarray) and val.size < size:
-        # If needed, extend the size of the vector to match desired size
-        if val[-1] != 0:
-            warn(f"{name} too short; padding with zeros")
-        val = np.hstack([val, np.zeros(size - val.size)])
+# Utility function to allow lists of 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)
 
-    nelem = _find_size(size, val, name)                 # determine size
-    return val, nelem
+    # Otherwise, do nothing
+    return X
 
 
 # Utility function to create an I/O system from a static gain
diff --git a/doc/control.rst b/doc/control.rst
index efd643d8..30ae1a03 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -137,7 +137,7 @@ Model simplification tools
     balred
     hsvd
     modred
-    era
+    eigensys_realization
     markov
 
 Nonlinear system support
diff --git a/doc/era_msd.py b/doc/era_msd.py
new file mode 120000
index 00000000..0cf6a528
--- /dev/null
+++ b/doc/era_msd.py
@@ -0,0 +1,1 @@
+../examples/era_msd.py
\ No newline at end of file
diff --git a/doc/era_msd.rst b/doc/era_msd.rst
new file mode 100644
index 00000000..de702406
--- /dev/null
+++ b/doc/era_msd.rst
@@ -0,0 +1,15 @@
+ERA example, mass spring damper system
+--------------------------------------
+
+Code
+....
+.. literalinclude:: era_msd.py
+   :language: python
+   :linenos:
+
+
+Notes
+.....
+
+1. The environment variable `PYCONTROL_TEST_EXAMPLES` is used for
+testing to turn off plotting of the outputs.0
\ No newline at end of file
diff --git a/doc/examples.rst b/doc/examples.rst
index 21364157..25e16115 100644
--- a/doc/examples.rst
+++ b/doc/examples.rst
@@ -35,6 +35,7 @@ other sources.
    kincar-flatsys
    mrac_siso_mit
    mrac_siso_lyapunov
+   era_msd
 
 Jupyter notebooks
 =================
diff --git a/doc/iosys.rst b/doc/iosys.rst
index eb4311e0..c0c2cca3 100644
--- a/doc/iosys.rst
+++ b/doc/iosys.rst
@@ -191,86 +191,7 @@ Additional features
 ===================
 
 The I/O systems module has a number of other features that can be used to
-simplify the creation and use of interconnected input/output systems.
-
-Vector elements processing
---------------------------
-
-Several I/O system commands perform processing of vector elements
-(such as initial states or input vectors) and broadcast these to the
-proper shape.
-
-For static elements, such as the initial state in a simulation or the
-nominal state and input for a linearization), the following processing
-is done:
-
-* Scalars are automatically converted to a vector of the appropriate
-  size consisting of the scalar value. This is commonly used when
-  specifying the origin ('0') or a step input ('1').
-
-* Lists of values are concatenated into a single vector.  This is
-  often used when you have an interconnected system and you need to
-  specify the initial condition or input value for each subsystem
-  (e.g., [X1eq, X2eq, ...]).
-
-* Vector elements are zero padded to the required length. If you
-  specify only a portion of the values for states or inputs, the
-  remaining values are taken as zero.  (If the final element in the
-  given vector is non-zero, a warning is issued.)
-
-Similar processing is done for input time series, used for the
-:func:`~control.input_output_response` and
-:func:`~control.forced_response` commands, with the following
-additional feature:
-
-* Time series elements are broadcast to match the number of time points
-  specified.  If a list of time series and static elements are given (as a
-  list), static elements are broadcast to the proper number of time points,
-  and the overall list of elements concatenated to provide the full input
-  vector.
-
-As an example, suppose we have an interconnected system consisting of three
-subsystems, a controlled process, an estimator, and a (static) controller::
-
-  proc = ct.nlsys(...,
-      states=2, inputs=['u1', 'u2', 'd'], outputs='y')
-  estim = ct.nlsys(...,
-      states=2, inputs='y', outputs=['xhat[0]', 'xhat[1]')
-  ctrl = ct.nlsys(...,
-      states=0, inputs=['r', 'xhat[0]', 'xhat[1]'], outputs=['u1', 'u2'])
-
-  clsys = ct.interconnect(
-      [proc, estim, ctrl], inputs=['r', 'd'], outputs=['y', 'u1', 'u2'])
-
-To linearize the system around the origin, we can utilize the scalar
-processing feature of vector elements::
-
-  P = proc.linearize(0, 0)
-
-In this command, the states and the inputs are broadcast to the size of the
-state and input vectors, respectively.
-
-If we want to linearize the closed loop system around a process state
-``x0`` (with two elements) and an estimator state ``0`` (for both states),
-we can use the list processing feature::
-
-  H = clsys.linearize([x0, 0], 0)
-
-Note that this also utilizes the zero-padding functionality, since the
-second argument in the list ``[x0, 0]`` is a scalar and so the vector
-``[x0, 0]`` only has three elements instead of the required four.
-
-To run an input/output simulation with a sinusoidal signal for the first
-input, a constant for the second input, and no external disturbance, we can
-use the list processing feature combined with time series broadcasting::
-
-  timepts = np.linspace(0, 10)
-  u1 = np.sin(timepts)
-  u2 = 1
-  resp = ct.input_output_response(clsys, timepts, [u1, u2, 0])
-
-In this command, the second and third arguments will be broadcast to match
-the number of time points.
+simplify the creation of interconnected input/output systems.
 
 Summing junction
 ----------------
diff --git a/examples/era_msd.py b/examples/era_msd.py
new file mode 100644
index 00000000..10193343
--- /dev/null
+++ b/examples/era_msd.py
@@ -0,0 +1,63 @@
+# era_msd.py
+# Johannes Kaisinger, 4 July 2024
+#
+# Demonstrate estimation of State Space model from impulse response.
+# SISO, SIMO, MISO, MIMO case
+
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+
+import control as ct
+
+# set up a mass spring damper system (2dof, MIMO case)
+# Mechanical Vibrations: Theory and Application, SI Edition, 1st ed.
+# Figure 6.5 / Example 6.7
+# m q_dd + c q_d + k q = f
+m1, k1, c1 = 1., 4., 1.
+m2, k2, c2 = 2., 2., 1.
+k3, c3 = 6., 2.
+
+A = np.array([
+    [0., 0., 1., 0.],
+    [0., 0., 0., 1.],
+    [-(k1+k2)/m1, (k2)/m1, -(c1+c2)/m1, c2/m1],
+    [(k2)/m2, -(k2+k3)/m2, c2/m2, -(c2+c3)/m2]
+])
+B = np.array([[0.,0.],[0.,0.],[1/m1,0.],[0.,1/m2]])
+C = np.array([[1.0, 0.0, 0.0, 0.0],[0.0, 1.0, 0.0, 0.0]])
+D = np.zeros((2,2))
+
+xixo_list = ["SISO","SIMO","MISO","MIMO"]
+xixo = xixo_list[3] # choose a system for estimation
+match xixo:
+    case "SISO":
+        sys = ct.StateSpace(A, B[:,0], C[0,:], D[0,0])
+    case "SIMO":
+        sys = ct.StateSpace(A, B[:,:1], C, D[:,:1])
+    case "MISO":
+        sys = ct.StateSpace(A, B, C[:1,:], D[:1,:])
+    case "MIMO":
+        sys = ct.StateSpace(A, B, C, D)
+
+
+dt = 0.1
+sysd = sys.sample(dt, method='zoh')
+response = ct.impulse_response(sysd)
+response.plot()
+plt.show()
+
+sysd_est, _ = ct.eigensys_realization(response,r=4,dt=dt)
+
+step_true = ct.step_response(sysd)
+step_true.sysname="H_true"
+step_est = ct.step_response(sysd_est)
+step_est.sysname="H_est"
+
+step_true.plot(title=xixo)
+step_est.plot(color='orange',linestyle='dashed')
+
+plt.show()
+
+if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
+    plt.show()
\ No newline at end of file


Test patch:
diff --git a/control/tests/bdalg_test.py b/control/tests/bdalg_test.py
index b9e26e8c..2ed793ef 100644
--- a/control/tests/bdalg_test.py
+++ b/control/tests/bdalg_test.py
@@ -316,49 +316,3 @@ class TestFeedback:
                 connect(sys, Q, [2], [1, 0])
             with pytest.raises(IndexError):
                 connect(sys, Q, [2], [1, -1])
-
-
-@pytest.mark.parametrize(
-    "op, nsys, ninputs, noutputs, nstates", [
-        (ctrl.series, 2, 1, 1, 4),
-        (ctrl.parallel, 2, 1, 1, 4),
-        (ctrl.feedback, 2, 1, 1, 4),
-        (ctrl.append, 2, 2, 2, 4),
-        (ctrl.negate, 1, 1, 1, 2),
-    ])
-def test_bdalg_update_names(op, nsys, ninputs, noutputs, nstates):
-    syslist = [ctrl.rss(2, 1, 1), ctrl.rss(2, 1, 1)]
-    inputs = ['in1', 'in2']
-    outputs = ['out1', 'out2']
-    states = ['x1', 'x2', 'x3', 'x4']
-
-    newsys = op(
-        *syslist[:nsys], name='newsys', inputs=inputs[:ninputs],
-        outputs=outputs[:noutputs], states=states[:nstates])
-    assert newsys.name == 'newsys'
-    assert newsys.ninputs == ninputs
-    assert newsys.input_labels == inputs[:ninputs]
-    assert newsys.noutputs == noutputs
-    assert newsys.output_labels == outputs[:noutputs]
-    assert newsys.nstates == nstates
-    assert newsys.state_labels == states[:nstates]
-
-
-def test_bdalg_udpate_names_errors():
-    sys1 = ctrl.rss(2, 1, 1)
-    sys2 = ctrl.rss(2, 1, 1)
-
-    with pytest.raises(ValueError, match="number of inputs does not match"):
-        sys = ctrl.series(sys1, sys2, inputs=2)
-
-    with pytest.raises(ValueError, match="number of outputs does not match"):
-        sys = ctrl.series(sys1, sys2, outputs=2)
-
-    with pytest.raises(ValueError, match="number of states does not match"):
-        sys = ctrl.series(sys1, sys2, states=2)
-
-    with pytest.raises(ValueError, match="number of states does not match"):
-        sys = ctrl.series(ctrl.tf(sys1), ctrl.tf(sys2), states=2)
-
-    with pytest.raises(TypeError, match="unrecognized keywords"):
-        sys = ctrl.series(sys1, sys2, dt=1)
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index cf4e3dd4..f3693cf0 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -231,20 +231,11 @@ class TestIOSys:
             linearized.C, [[1, 0, 0], [0, 1, 0]])
         np.testing.assert_array_almost_equal(linearized.D, np.zeros((2,2)))
 
-        # Pass fewer than the required elements
-        padded = iosys.linearize([0, 0], np.array([0]))
-        assert padded.nstates == linearized.nstates
-        assert padded.ninputs == linearized.ninputs
-
-        # Check for warning if last element before padding is nonzero
-        with pytest.warns(UserWarning, match="x0 too short; padding"):
-            padded = iosys.linearize([0, 1], np.array([0]))
-
     @pytest.mark.usefixtures("editsdefaults")
     def test_linearize_named_signals(self, kincar):
         # Full form of the call
-        linearized = kincar.linearize(
-            [0, 0, 0], [0, 0], copy_names=True, name='linearized')
+        linearized = kincar.linearize([0, 0, 0], [0, 0], copy_names=True,
+                                      name='linearized')
         assert linearized.name == 'linearized'
         assert linearized.find_input('v') == 0
         assert linearized.find_input('phi') == 1
@@ -265,8 +256,8 @@ class TestIOSys:
         assert lin_nocopy.find_state('x') is None
 
         # if signal names are provided, they should override those of kincar
-        linearized_newnames = kincar.linearize(
-            [0, 0, 0], [0, 0], name='linearized',
+        linearized_newnames = kincar.linearize([0, 0, 0], [0, 0],
+            name='linearized',
             copy_names=True, inputs=['v2', 'phi2'], outputs=['x2','y2'])
         assert linearized_newnames.name == 'linearized'
         assert linearized_newnames.find_input('v2') == 0
@@ -278,11 +269,6 @@ class TestIOSys:
         assert linearized_newnames.find_output('x') is None
         assert linearized_newnames.find_output('y') is None
 
-        # if system name is provided but copy_names is false, override name
-        linearized_newsysname = kincar.linearize(
-            [0, 0, 0], [0, 0], name='newname', copy_names=False)
-        assert linearized_newsysname.name == 'newname'
-
         # Test legacy version as well
         with pytest.warns(UserWarning, match="NumPy matrix class no longer"):
             ct.use_legacy_defaults('0.8.4')
@@ -1431,7 +1417,7 @@ class TestIOSys:
     def test_neg_badsize(self):
         # Create a system of unspecified size
         sys = ct.NonlinearIOSystem(lambda t, x, u, params: -x)
-        with pytest.raises(ValueError, match="Can't determine number"):
+        with pytest.raises(ValueError, match="Can't determine"):
             -sys
 
     def test_bad_signal_list(self):
@@ -1895,7 +1881,7 @@ def test_input_output_broadcasting():
     np.testing.assert_equal(resp_cov0.states, resp_init.states)
 
     # Specify only some of the initial conditions
-    with pytest.warns(UserWarning, match="X0 too short; padding"):
+    with pytest.warns(UserWarning, match="initial state too short; padding"):
         resp_short = ct.input_output_response(sys, T, [U[0], [0, 1]], [X0, 1])
 
     # Make sure that inconsistent settings don't work
@@ -2086,7 +2072,6 @@ def test_find_eqpt(x0, ix, u0, iu, y0, iy, dx0, idx, dt, x_expect, u_expect):
     np.testing.assert_allclose(np.array(xeq), x_expect, atol=1e-6)
     np.testing.assert_allclose(np.array(ueq), u_expect, atol=1e-6)
 
-
 def test_iosys_sample():
     csys = ct.rss(2, 1, 1)
     dsys = csys.sample(0.1)
@@ -2097,74 +2082,3 @@ def test_iosys_sample():
     dsys = ct.sample_system(csys, 0.1)
     assert isinstance(dsys, ct.StateSpace)
     assert dsys.dt == 0.1
-
-
-# Make sure that we can determine system sizes automatically
-def test_find_size():
-    # Create a nonlinear system with no size information
-    sys = ct.nlsys(
-        lambda t, x, u, params: -x + u,
-        lambda t, x, u, params: x[:1])
-
-    # Run a simulation with size set by parameters
-    timepts = np.linspace(0, 1)
-    resp = ct.input_output_response(sys, timepts, [0, 1], X0=[0, 0])
-    assert resp.states.shape[0] == 2
-    assert resp.inputs.shape[0] == 2
-    assert resp.outputs.shape[0] == 1
-
-    #
-    # Make sure we get warnings if things are inconsistent
-    #
-
-    # Define a system of fixed size
-    sys = ct.nlsys(
-        lambda t, x, u, params: -x + u,
-        lambda t, x, u, params: x[:1],
-        inputs=2, states=2)
-
-    with pytest.raises(ValueError, match="inconsistent .* size of X0"):
-        resp = ct.input_output_response(sys, timepts, [0, 1], X0=[0, 0, 1])
-
-    with pytest.raises(ValueError, match=".*U.* Wrong shape"):
-        resp = ct.input_output_response(sys, timepts, [0, 1, 2], X0=[0, 0])
-
-    with pytest.raises(RuntimeError, match="inconsistent size of outputs"):
-        sys = ct.nlsys(
-            lambda t, x, u, params: -x + u,
-            lambda t, x, u, params: x[:1],
-            inputs=2, states=2, outputs=2)
-        resp = ct.input_output_response(sys, timepts, [0, 1], X0=[0, 0])
-
-
-def test_update_names():
-    sys = ct.rss(['x1', 'x2'], 2, 2)
-    sys.update_names(
-        name='new', states=2, inputs=['u1', 'u2'],
-        outputs=2, output_prefix='yy')
-    assert sys.name == 'new'
-    assert sys.ninputs == 2
-    assert sys.input_labels == ['u1', 'u2']
-    assert sys.ninputs == 2
-    assert sys.output_labels == ['yy[0]', 'yy[1]']
-    assert sys.state_labels == ['x[0]', 'x[1]']
-
-    # Generate some error conditions
-    with pytest.raises(ValueError, match="number of inputs does not match"):
-        sys.update_names(inputs=3)
-
-    with pytest.raises(ValueError, match="number of outputs does not match"):
-        sys.update_names(outputs=3)
-
-    with pytest.raises(ValueError, match="number of states does not match"):
-        sys.update_names(states=3)
-
-    with pytest.raises(ValueError, match="number of states does not match"):
-        siso = ct.tf([1], [1, 2, 1])
-        ct.tf(siso).update_names(states=2)
-
-    with pytest.raises(TypeError, match="unrecognized keywords"):
-        sys.update_names(dt=1)
-
-    with pytest.raises(TypeError, match=".* takes 1 positional argument"):
-        sys.update_names(5)
diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
index 4d252ab1..0e0abc0d 100644
--- a/control/tests/kwargs_test.py
+++ b/control/tests/kwargs_test.py
@@ -92,25 +92,20 @@ def test_kwarg_search(module, prefix):
 
 @pytest.mark.parametrize(
     "function, nsssys, ntfsys, moreargs, kwargs",
-    [(control.append, 2, 0, (), {}),
-     (control.dlqe, 1, 0, ([[1]], [[1]]), {}),
+    [(control.dlqe, 1, 0, ([[1]], [[1]]), {}),
      (control.dlqr, 1, 0, ([[1, 0], [0, 1]], [[1]]), {}),
      (control.drss, 0, 0, (2, 1, 1), {}),
-     (control.feedback, 2, 0, (), {}),
      (control.flatsys.flatsys, 1, 0, (), {}),
      (control.input_output_response, 1, 0, ([0, 1, 2], [1, 1, 1]), {}),
      (control.lqe, 1, 0, ([[1]], [[1]]), {}),
      (control.lqr, 1, 0, ([[1, 0], [0, 1]], [[1]]), {}),
      (control.linearize, 1, 0, (0, 0), {}),
-     (control.negate, 1, 0, (), {}),
      (control.nlsys, 0, 0, (lambda t, x, u, params: np.array([0]),), {}),
-     (control.parallel, 2, 0, (), {}),
      (control.pzmap, 1, 0, (), {}),
      (control.rlocus, 0, 1, (), {}),
      (control.root_locus, 0, 1, (), {}),
      (control.root_locus_plot, 0, 1, (), {}),
      (control.rss, 0, 0, (2, 1, 1), {}),
-     (control.series, 2, 0, (), {}),
      (control.set_defaults, 0, 0, ('control',), {'default_dt': True}),
      (control.ss, 0, 0, (0, 0, 0, 0), {'dt': 1}),
      (control.ss2io, 1, 0,  (), {}),
@@ -127,7 +122,6 @@ def test_kwarg_search(module, prefix):
      (control.LTI, 0, 0, (),
       {'inputs': 1, 'outputs': 1, 'states': 1}),
      (control.flatsys.LinearFlatSystem, 1, 0, (), {}),
-     (control.InputOutputSystem.update_names, 1, 0, (), {}),
      (control.NonlinearIOSystem.linearize, 1, 0, (0, 0), {}),
      (control.StateSpace.sample, 1, 0, (0.1,), {}),
      (control.StateSpace, 0, 0,
@@ -238,7 +232,6 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
 #
 
 kwarg_unittest = {
-    'append': test_unrecognized_kwargs,
     'bode': test_response_plot_kwargs,
     'bode_plot': test_response_plot_kwargs,
     'create_estimator_iosystem': stochsys_test.test_estimator_errors,
@@ -249,7 +242,6 @@ kwarg_unittest = {
     'dlqe': test_unrecognized_kwargs,
     'dlqr': test_unrecognized_kwargs,
     'drss': test_unrecognized_kwargs,
-    'feedback': test_unrecognized_kwargs,
     'flatsys.flatsys': test_unrecognized_kwargs,
     'frd': frd_test.TestFRD.test_unrecognized_keyword,
     'gangof4': test_matplotlib_kwargs,
@@ -260,7 +252,6 @@ kwarg_unittest = {
     'linearize': test_unrecognized_kwargs,
     'lqe': test_unrecognized_kwargs,
     'lqr': test_unrecognized_kwargs,
-    'negate': test_unrecognized_kwargs,
     'nichols_plot': test_matplotlib_kwargs,
     'nichols': test_matplotlib_kwargs,
     'nlsys': test_unrecognized_kwargs,
@@ -268,14 +259,12 @@ kwarg_unittest = {
     'nyquist_response': test_response_plot_kwargs,
     'nyquist_plot': test_matplotlib_kwargs,
     'phase_plane_plot': test_matplotlib_kwargs,
-    'parallel': test_unrecognized_kwargs,
     'pole_zero_plot': test_unrecognized_kwargs,
     'pzmap': test_unrecognized_kwargs,
     'rlocus': test_unrecognized_kwargs,
     'root_locus': test_unrecognized_kwargs,
     'root_locus_plot': test_unrecognized_kwargs,
     'rss': test_unrecognized_kwargs,
-    'series': test_unrecognized_kwargs,
     'set_defaults': test_unrecognized_kwargs,
     'singular_values_plot': test_matplotlib_kwargs,
     'ss': test_unrecognized_kwargs,
@@ -303,7 +292,6 @@ kwarg_unittest = {
     'DescribingFunctionResponse.plot':
         descfcn_test.test_describing_function_exceptions,
     'InputOutputSystem.__init__': test_unrecognized_kwargs,
-    'InputOutputSystem.update_names': test_unrecognized_kwargs,
     'LTI.__init__': test_unrecognized_kwargs,
     'flatsys.LinearFlatSystem.__init__': test_unrecognized_kwargs,
     'NonlinearIOSystem.linearize': test_unrecognized_kwargs,
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
index 49c2afd5..dc50ce96 100644
--- a/control/tests/modelsimp_test.py
+++ b/control/tests/modelsimp_test.py
@@ -7,10 +7,10 @@ import numpy as np
 import pytest
 
 
-from control import StateSpace, forced_response, tf, rss, c2d
+from control import StateSpace, impulse_response, step_response, forced_response, tf, rss, c2d
 from control.exception import ControlMIMONotImplemented
 from control.tests.conftest import slycotonly
-from control.modelsimp import balred, hsvd, markov, modred
+from control.modelsimp import balred, hsvd, markov, modred, eigensys_realization
 
 
 class TestModelsimp:
@@ -111,6 +111,85 @@ class TestModelsimp:
         # for k=5, m=n=10: 0.015 %
         np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
 
+    def testERASignature(self):
+
+        # test siso
+        # Katayama, Subspace Methods for System Identification
+        # Example 6.1, Fibonacci sequence
+        H_true = np.array([0.,1.,1.,2.,3.,5.,8.,13.,21.,34.])
+
+        # A realization of fibonacci impulse response
+        A = np.array([[0., 1.],[1., 1.,]])
+        B = np.array([[1.],[1.,]])
+        C = np.array([[1., 0.,]])
+        D = np.array([[0.,]])
+        
+        T = np.arange(0,10,1)
+        sysd_true = StateSpace(A,B,C,D,True)
+        ir_true = impulse_response(sysd_true,T=T)
+
+        # test TimeResponseData
+        sysd_est, _  = eigensys_realization(ir_true,r=2)
+        ir_est = impulse_response(sysd_est, T=T)
+        _, H_est = ir_est
+    
+        np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
+
+        # test ndarray
+        _, YY_true = ir_true
+        sysd_est, _  = eigensys_realization(YY_true,r=2)
+        ir_est = impulse_response(sysd_est, T=T)
+        _, H_est = ir_est
+    
+        np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
+
+        # test mimo
+        # Mechanical Vibrations: Theory and Application, SI Edition, 1st ed.
+        # Figure 6.5 / Example 6.7
+        # m q_dd + c q_d + k q = f
+        m1, k1, c1 = 1., 4., 1.
+        m2, k2, c2 = 2., 2., 1.
+        k3, c3 = 6., 2.
+
+        A = np.array([
+            [0., 0., 1., 0.],
+            [0., 0., 0., 1.],
+            [-(k1+k2)/m1, (k2)/m1, -(c1+c2)/m1, c2/m1],
+            [(k2)/m2, -(k2+k3)/m2, c2/m2, -(c2+c3)/m2]
+        ])
+        B = np.array([[0.,0.],[0.,0.],[1/m1,0.],[0.,1/m2]])
+        C = np.array([[1.0, 0.0, 0.0, 0.0],[0.0, 1.0, 0.0, 0.0]])
+        D = np.zeros((2,2))
+
+        sys = StateSpace(A, B, C, D)
+
+        dt = 0.1
+        T = np.arange(0,10,dt)
+        sysd_true = sys.sample(dt, method='zoh')
+        ir_true = impulse_response(sysd_true, T=T)
+
+        # test TimeResponseData
+        sysd_est, _ = eigensys_realization(ir_true,r=4,dt=dt)
+
+        step_true = step_response(sysd_true)
+        step_est = step_response(sysd_est)
+
+        np.testing.assert_allclose(step_true.outputs, 
+                                   step_est.outputs,
+                                   rtol=1e-6, atol=1e-8)
+        
+        # test ndarray
+        _, YY_true = ir_true
+        sysd_est, _  = eigensys_realization(YY_true,r=4,dt=dt)
+
+        step_true = step_response(sysd_true, T=T)
+        step_est = step_response(sysd_est, T=T)
+
+        np.testing.assert_allclose(step_true.outputs, 
+                                   step_est.outputs,
+                                   rtol=1e-6, atol=1e-8)
+
+
     def testModredMatchDC(self):
         #balanced realization computed in matlab for the transfer function:
         # num = [1 11 45 32], den = [1 15 60 200 60]
diff --git a/control/tests/nlsys_test.py b/control/tests/nlsys_test.py
index 7f649e0c..80baa646 100644
--- a/control/tests/nlsys_test.py
+++ b/control/tests/nlsys_test.py
@@ -46,7 +46,7 @@ def test_nlsys_basic():
     ])
 def test_lti_nlsys_response(nin, nout, input, output):
     sys_ss = ct.rss(4, nin, nout, strictly_proper=True)
-    sys_ss.A = np.diag([-1, -2, -3, -4])        # avoid random numerical errors
+    sys_ss.A = np.diag([-1, -2, -3, -4])        # avoid random noise errors
     sys_nl = ct.nlsys(
         lambda t, x, u, params: sys_ss.A @ x + sys_ss.B @ u,
         lambda t, x, u, params: sys_ss.C @ x + sys_ss.D @ u,