diff --git a/control/bdalg.py b/control/bdalg.py
index ce8008537..7bfd327eb 100644
--- a/control/bdalg.py
+++ b/control/bdalg.py
@@ -54,17 +54,19 @@
 """
 
 from functools import reduce
-import numpy as np
 from warnings import warn
-from . import xferfcn as tf
-from . import statesp as ss
+
+import numpy as np
+
 from . import frdata as frd
+from . import statesp as ss
+from . import xferfcn as tf
 from .iosys import InputOutputSystem
 
 __all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect']
 
 
-def series(sys1, *sysn):
+def series(sys1, *sysn, **kwargs):
     r"""series(sys1, sys2, [..., sysn])
 
     Return the series connection (`sysn` \* ...\  \*) `sys2` \* `sys1`.
@@ -79,6 +81,20 @@ def series(sys1, *sysn):
     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
@@ -117,10 +133,12 @@ def series(sys1, *sysn):
     (2, 1, 5)
 
     """
-    return reduce(lambda x, y: y * x, sysn, sys1)
+    sys = reduce(lambda x, y: y * x, sysn, sys1)
+    sys.update_names(**kwargs)
+    return sys
 
 
-def parallel(sys1, *sysn):
+def parallel(sys1, *sysn, **kwargs):
     r"""parallel(sys1, sys2, [..., sysn])
 
     Return the parallel connection `sys1` + `sys2` (+ ...\  + `sysn`).
@@ -135,6 +153,20 @@ def parallel(sys1, *sysn):
     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
@@ -171,10 +203,11 @@ def parallel(sys1, *sysn):
     (3, 4, 7)
 
     """
-    return reduce(lambda x, y: x + y, sysn, sys1)
-
+    sys = reduce(lambda x, y: x + y, sysn, sys1)
+    sys.update_names(**kwargs)
+    return sys
 
-def negate(sys):
+def negate(sys, **kwargs):
     """
     Return the negative of a system.
 
@@ -188,15 +221,29 @@ def negate(sys):
     out : scalar, array, or :class:`InputOutputSystem`
         Negated system.
 
-    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.
+    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.
+
     Examples
     --------
     >>> G = ct.tf([2], [1, 1])
@@ -208,11 +255,12 @@ def negate(sys):
     np.float64(-2.0)
 
     """
-    return -sys
+    sys = -sys
+    sys.update_names(**kwargs)
+    return sys
 
 #! TODO: expand to allow sys2 default to work in MIMO case?
-#! TODO: allow renaming of signals (for all bdalg operations)
-def feedback(sys1, sys2=1, sign=-1):
+def feedback(sys1, sys2=1, sign=-1, **kwargs):
     """Feedback interconnection between two I/O systems.
 
     Parameters
@@ -229,6 +277,20 @@ def feedback(sys1, sys2=1, sign=-1):
     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
@@ -261,7 +323,7 @@ def feedback(sys1, sys2=1, sign=-1):
     # Allow anything with a feedback function to call that function
     # TODO: rewrite to allow __rfeedback__
     try:
-        return sys1.feedback(sys2, sign)
+        return sys1.feedback(sys2, sign, **kwargs)
     except (AttributeError, TypeError):
         pass
 
@@ -284,9 +346,11 @@ def feedback(sys1, sys2=1, sign=-1):
         else:
             sys1 = ss._convert_to_statespace(sys1)
 
-    return sys1.feedback(sys2, sign)
+    sys = sys1.feedback(sys2, sign)
+    sys.update_names(**kwargs)
+    return sys
 
-def append(*sys):
+def append(*sys, **kwargs):
     """append(sys1, sys2, [..., sysn])
 
     Group LTI state space models by appending their inputs and outputs.
@@ -299,6 +363,20 @@ def append(*sys):
     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`
@@ -327,6 +405,7 @@ def append(*sys):
     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):
@@ -370,6 +449,12 @@ 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)
@@ -378,12 +463,6 @@ 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 fbd5c1dba..d00dade65 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -6,10 +6,12 @@
 # FrequencyResponseData, InterconnectedSystem and other similar classes
 # that allow naming of signals.
 
-import numpy as np
+import re
 from copy import deepcopy
 from warnings import warn
-import re
+
+import numpy as np
+
 from . import config
 
 __all__ = ['InputOutputSystem', 'issiso', 'timebase', 'common_timebase',
@@ -366,6 +368,51 @@ def find_states(self, name_list):
         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.
@@ -823,7 +870,6 @@ 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/nlsys.py b/control/nlsys.py
index 358c4b125..139853fe6 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -143,6 +143,9 @@ def __init__(self, updfcn, outfcn=None, params=None, **kwargs):
                 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
@@ -514,7 +517,7 @@ def feedback(self, other=1, sign=-1, params=None):
         return newsys
 
     def linearize(self, x0, u0, t=0, params=None, eps=1e-6,
-                  name=None, copy_names=False, **kwargs):
+                  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
@@ -529,26 +532,16 @@ def linearize(self, x0, u0, t=0, params=None, eps=1e-6,
         # numerical linearization use the `_rhs()` and `_out()` member
         # functions.
         #
-        # If x0 and u0 are specified as lists, concatenate the elements
-        x0 = _concatenate_list_elements(x0, 'x0')
-        u0 = _concatenate_list_elements(u0, 'u0')
+        # Process nominal states and inputs
+        x0, nstates = _process_vector_argument(x0, "x0", self.nstates)
+        u0, ninputs = _process_vector_argument(u0, "u0", self.ninputs)
 
-        # 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
+        # Update the current parameters (prior to calling _out())
+        self._update_params(params)
 
         # 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)
@@ -579,8 +572,6 @@ def linearize(self, x0, u0, t=0, params=None, eps=1e-6,
         # 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)
@@ -1467,7 +1458,16 @@ def input_output_response(
         # Use the input time points as the output time points
         t_eval = T
 
-    # If we were passed a list of input, concatenate them (w/ broadcast)
+    #
+    # 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 isinstance(U, (tuple, list)) and len(U) != ntimepts:
         U_elements = []
         for i, u in enumerate(U):
@@ -1491,62 +1491,44 @@ 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 sys.ninputs is None or sys.ninputs == 1:
+    if ninputs is None or ninputs == 1:
         legal_shapes = [(ntimepts,), (1, ntimepts)]
     else:
-        legal_shapes = [(sys.ninputs, ntimepts)]
+        legal_shapes = [(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]
 
-    # If we were passed a list of initial states, concatenate them
-    X0 = _concatenate_list_elements(X0, 'X0')
-
-    # 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)])
+    # Process initial states
+    X0, nstates = _process_vector_argument(X0, "X0", sys.nstates)
 
-    # 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)
+    # Update the parameter values (prior to evaluating outfcn)
+    sys._update_params(params)
 
     # Figure out the number of outputs
-    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]
+    if sys.outfcn is None:
+        noutputs = nstates if sys.noutputs is None else sys.noutputs
     else:
-        noutputs = sys.noutputs
+        noutputs = np.shape(sys._out(T[0], X0, U[:, 0]))[0]
 
-    # Update the parameter values
-    sys._update_params(params)
+    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}")
 
     #
     # Define a function to evaluate the input at an arbitrary time
@@ -1753,17 +1735,9 @@ 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
-    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
+    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)
 
     # Make sure the input arguments match the sizes of the system
     if len(x0) != nstates or \
@@ -1998,7 +1972,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, label):
+def _find_size(sysval, vecval, name="system component"):
     """Utility function to find the size of a system parameter
 
     If both parameters are not None, they must be consistent.
@@ -2006,7 +1980,8 @@ def _find_size(sysval, vecval, label):
     if hasattr(vecval, '__len__'):
         if sysval is not None and sysval != len(vecval):
             raise ValueError(
-                f"inconsistent information for number of {label}")
+                f"inconsistent information to determine size of {name}; "
+                f"expected {sysval} values, received {len(vecval)}")
         return len(vecval)
     # None or 0, which is a valid value for "a (sysval, ) vector of zeros".
     if not vecval:
@@ -2014,7 +1989,7 @@ def _find_size(sysval, vecval, label):
     elif sysval == 1:
         # (1, scalar) is also a valid combination from legacy code
         return 1
-    raise ValueError(f"can't determine number of {label}")
+    raise ValueError(f"can't determine size of {name}")
 
 
 # Function to create an interconnected system
@@ -2572,18 +2547,59 @@ def interconnect(
     return newsys
 
 
-# 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)
+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)])
 
-    # Otherwise, do nothing
-    return X
+    nelem = _find_size(size, val, name)                 # determine size
+    return val, nelem
 
 
 # Utility function to create an I/O system from a static gain
diff --git a/control/tests/bdalg_test.py b/control/tests/bdalg_test.py
index 2ed793ef2..b9e26e8c0 100644
--- a/control/tests/bdalg_test.py
+++ b/control/tests/bdalg_test.py
@@ -316,3 +316,49 @@ def testConnect(self, tsys):
                 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 f3693cf00..cf4e3dd43 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -231,11 +231,20 @@ def test_linearize(self, tsys, kincar):
             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
@@ -256,8 +265,8 @@ def test_linearize_named_signals(self, kincar):
         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
@@ -269,6 +278,11 @@ def test_linearize_named_signals(self, kincar):
         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')
@@ -1417,7 +1431,7 @@ def test_operand_badtype(self, C, op):
     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"):
+        with pytest.raises(ValueError, match="Can't determine number"):
             -sys
 
     def test_bad_signal_list(self):
@@ -1881,7 +1895,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="initial state too short; padding"):
+    with pytest.warns(UserWarning, match="X0 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
@@ -2072,6 +2086,7 @@ 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)
@@ -2082,3 +2097,74 @@ 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 d6bd06487..73dbe3229 100644
--- a/control/tests/kwargs_test.py
+++ b/control/tests/kwargs_test.py
@@ -92,20 +92,25 @@ def test_kwarg_search(module, prefix):
 
 @pytest.mark.parametrize(
     "function, nsssys, ntfsys, moreargs, kwargs",
-    [(control.dlqe, 1, 0, ([[1]], [[1]]), {}),
+    [(control.append, 2, 0, (), {}),
+     (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,  (), {}),
@@ -122,6 +127,7 @@ 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,
@@ -232,6 +238,7 @@ 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,
@@ -242,6 +249,7 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
     '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,
@@ -252,6 +260,7 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
     '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,
@@ -259,12 +268,14 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
     '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,
@@ -292,6 +303,7 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
     '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/nlsys_test.py b/control/tests/nlsys_test.py
index 80baa646f..7f649e0cc 100644
--- a/control/tests/nlsys_test.py
+++ b/control/tests/nlsys_test.py
@@ -46,7 +46,7 @@ def kincar_output(t, x, u, params):
     ])
 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 noise errors
+    sys_ss.A = np.diag([-1, -2, -3, -4])        # avoid random numerical 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,
diff --git a/doc/iosys.rst b/doc/iosys.rst
index c0c2cca31..eb4311e05 100644
--- a/doc/iosys.rst
+++ b/doc/iosys.rst
@@ -191,7 +191,86 @@ Additional features
 ===================
 
 The I/O systems module has a number of other features that can be used to
-simplify the creation of interconnected input/output systems.
+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.
 
 Summing junction
 ----------------