Repository: python-control/python-control

Issue:
Elegant way to generate a transfer function matrix

Hi everybody,
I want to know if there is a more elegant way to generate a transfer-function-matrix for a hinf-controller.

In matlab I just do:
P11 = [0 0; W2*G W2]
P12 = [W1; W2*G]
P21 = [G 1]
P22 = [G]
P_ = [P11 P12; P21 P22]
sys = ss(P_, 'minimal');

In python with control, I habe to do:
import control
from control import *
Gss = ss([[-1.0, -10.0], [1, 0]],  [[10.0], [0]], [0, 1], [0])
G = tf(Gss)
W1 = tf([100], [0.1, 1])
W2 = tf([100], [0.1, 1])
W2G = W2*G
num = [ [ [0], [0], W1.num[0][0] ],  [W2G.num[0][0], W2.num[0][0], G.num[0][0]], [G.num[0][0], [1], G.num[0][0] ] ]   
den = [ [ [1], [1], W1.den[0][0] ],  [W2G.den[0][0], W2.den[0][0], G.den[0][0]], [G.den[0][0], [1], G.den[0][0] ] ] 
P = tf( num, den )
print(P)

The problem is that the return of the numerator and denominator is not directly the list.

And I get a high order system. Should I use just the balanced trunction?
Pss = ss(P)
Pbr = balred(Pss, 4)

________________________
Matlab is using minreal with exist also in "control", but I have to find manually a suitable tolerance...

Discussion:

Solution patch:
diff --git a/control/bdalg.py b/control/bdalg.py
index 024d95fba..d907cd3c5 100644
--- a/control/bdalg.py
+++ b/control/bdalg.py
@@ -10,6 +10,8 @@
 negate
 feedback
 connect
+combine_tf
+split_tf
 
 """
 
@@ -63,7 +65,8 @@
 from . import xferfcn as tf
 from .iosys import InputOutputSystem
 
-__all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect']
+__all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect',
+           'combine_tf', 'split_tf']
 
 
 def series(sys1, *sysn, **kwargs):
@@ -507,3 +510,212 @@ def connect(sys, Q, inputv, outputv):
         Ytrim[i,y-1] = 1.
 
     return Ytrim * sys * Utrim
+
+def combine_tf(tf_array):
+    """Combine array-like of transfer functions into MIMO transfer function.
+
+    Parameters
+    ----------
+    tf_array : list of list of TransferFunction or array_like
+        Transfer matrix represented as a two-dimensional array or list-of-lists
+        containing TransferFunction objects. The TransferFunction objects can
+        have multiple outputs and inputs, as long as the dimensions are
+        compatible.
+
+    Returns
+    -------
+    TransferFunction
+        Transfer matrix represented as a single MIMO TransferFunction object.
+
+    Raises
+    ------
+    ValueError
+        If timesteps of transfer functions do not match.
+    ValueError
+        If ``tf_array`` has incorrect dimensions.
+    ValueError
+        If the transfer functions in a row have mismatched output or input
+        dimensions.
+
+    Examples
+    --------
+    Combine two transfer functions
+
+    >>> s = control.TransferFunction.s
+    >>> control.combine_tf([
+    ...     [1 / (s + 1)],
+    ...     [s / (s + 2)],
+    ... ])
+    TransferFunction([[array([1])], [array([1, 0])]],
+                     [[array([1, 1])], [array([1, 2])]])
+
+    Combine NumPy arrays with transfer functions
+
+    >>> control.combine_tf([
+    ...     [np.eye(2), np.zeros((2, 1))],
+    ...     [np.zeros((1, 2)), control.TransferFunction([1], [1, 0])],
+    ... ])
+    TransferFunction([[array([1.]), array([0.]), array([0.])],
+                      [array([0.]), array([1.]), array([0.])],
+                      [array([0.]), array([0.]), array([1])]],
+                     [[array([1.]), array([1.]), array([1.])],
+                      [array([1.]), array([1.]), array([1.])],
+                      [array([1.]), array([1.]), array([1, 0])]])
+    """
+    # Find common timebase or raise error
+    dt_list = []
+    try:
+        for row in tf_array:
+            for tfn in row:
+                dt_list.append(getattr(tfn, "dt", None))
+    except OSError:
+        raise ValueError("`tf_array` has too few dimensions.")
+    dt_set = set(dt_list)
+    dt_set.discard(None)
+    if len(dt_set) > 1:
+        raise ValueError("Timesteps of transfer functions are "
+                         f"mismatched: {dt_set}")
+    elif len(dt_set) == 0:
+        dt = None
+    else:
+        dt = dt_set.pop()
+    # Convert all entries to transfer function objects
+    ensured_tf_array = []
+    for row in tf_array:
+        ensured_row = []
+        for tfn in row:
+            ensured_row.append(_ensure_tf(tfn, dt))
+        ensured_tf_array.append(ensured_row)
+    # Iterate over
+    num = []
+    den = []
+    for row_index, row in enumerate(ensured_tf_array):
+        for j_out in range(row[0].noutputs):
+            num_row = []
+            den_row = []
+            for col in row:
+                if col.noutputs != row[0].noutputs:
+                    raise ValueError(
+                        "Mismatched number of transfer function outputs in "
+                        f"row {row_index}."
+                    )
+                for j_in in range(col.ninputs):
+                    num_row.append(col.num[j_out][j_in])
+                    den_row.append(col.den[j_out][j_in])
+            num.append(num_row)
+            den.append(den_row)
+    for row_index, row in enumerate(num):
+        if len(row) != len(num[0]):
+            raise ValueError(
+                "Mismatched number transfer function inputs in row "
+                f"{row_index} of numerator."
+            )
+    for row_index, row in enumerate(den):
+        if len(row) != len(den[0]):
+            raise ValueError(
+                "Mismatched number transfer function inputs in row "
+                f"{row_index} of denominator."
+            )
+    return tf.TransferFunction(num, den, dt=dt)
+
+def split_tf(transfer_function):
+    """Split MIMO transfer function into NumPy array of SISO tranfer functions.
+
+    Parameters
+    ----------
+    transfer_function : TransferFunction
+        MIMO transfer function to split.
+
+    Returns
+    -------
+    np.ndarray
+        NumPy array of SISO transfer functions.
+
+    Examples
+    --------
+    Split a MIMO transfer function
+
+    >>> G = control.TransferFunction(
+    ...     [
+    ...         [[87.8], [-86.4]],
+    ...         [[108.2], [-109.6]],
+    ...     ],
+    ...     [
+    ...         [[1, 1], [1, 1]],
+    ...         [[1, 1], [1, 1]],
+    ...     ],
+    ... )
+    >>> control.split_tf(G)
+    array([[TransferFunction(array([87.8]), array([1, 1])),
+            TransferFunction(array([-86.4]), array([1, 1]))],
+           [TransferFunction(array([108.2]), array([1, 1])),
+            TransferFunction(array([-109.6]), array([1, 1]))]], dtype=object)
+    """
+    tf_split_lst = []
+    for i_out in range(transfer_function.noutputs):
+        row = []
+        for i_in in range(transfer_function.ninputs):
+            row.append(
+                tf.TransferFunction(
+                    transfer_function.num[i_out][i_in],
+                    transfer_function.den[i_out][i_in],
+                    dt=transfer_function.dt,
+                )
+            )
+        tf_split_lst.append(row)
+    return np.array(tf_split_lst, dtype=object)
+
+def _ensure_tf(arraylike_or_tf, dt=None):
+    """Convert an array-like to a transfer function.
+
+    Parameters
+    ----------
+    arraylike_or_tf : TransferFunction or array_like
+        Array-like or transfer function.
+    dt : None, True or float, optional
+        System timebase. 0 (default) indicates continuous
+        time, True indicates discrete time with unspecified sampling
+        time, positive number is discrete time with specified
+        sampling time, None indicates unspecified timebase (either
+        continuous or discrete time). If None, timestep is not validated.
+
+    Returns
+    -------
+    TransferFunction
+        Transfer function.
+
+    Raises
+    ------
+    ValueError
+        If input cannot be converted to a transfer function.
+    ValueError
+        If the timesteps do not match.
+    """
+    # If the input is already a transfer function, return it right away
+    if isinstance(arraylike_or_tf, tf.TransferFunction):
+        # If timesteps don't match, raise an exception
+        if (dt is not None) and (arraylike_or_tf.dt != dt):
+            raise ValueError(
+                f"`arraylike_or_tf.dt={arraylike_or_tf.dt}` does not match "
+                f"argument `dt={dt}`."
+            )
+        return arraylike_or_tf
+    if np.ndim(arraylike_or_tf) > 2:
+        raise ValueError(
+            "Array-like must have less than two dimensions to be converted "
+            "into a transfer function."
+        )
+    # If it's not, then convert it to a transfer function
+    arraylike_3d = np.atleast_3d(arraylike_or_tf)
+    try:
+        tfn = tf.TransferFunction(
+            arraylike_3d,
+            np.ones_like(arraylike_3d),
+            dt,
+        )
+    except TypeError:
+        raise ValueError(
+            "`arraylike_or_tf` must only contain array-likes or transfer "
+            "functions."
+        )
+    return tfn
diff --git a/control/exception.py b/control/exception.py
index e4758cc49..add5d01ae 100644
--- a/control/exception.py
+++ b/control/exception.py
@@ -52,6 +52,10 @@ class ControlArgument(TypeError):
     """Raised when arguments to a function are not correct"""
     pass
 
+class ControlIndexError(IndexError):
+    """Raised when arguments to an indexed object are not correct"""
+    pass
+
 class ControlMIMONotImplemented(NotImplementedError):
     """Function is not currently implemented for MIMO systems"""
     pass
diff --git a/control/frdata.py b/control/frdata.py
index 42ecee0d9..1bdf28528 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -10,6 +10,7 @@
 FRD data.
 """
 
+from collections.abc import Iterable
 from copy import copy
 from warnings import warn
 
@@ -20,7 +21,8 @@
 
 from . import config
 from .exception import pandas_check
-from .iosys import InputOutputSystem, _process_iosys_keywords, common_timebase
+from .iosys import InputOutputSystem, NamedSignal, _process_iosys_keywords, \
+    _process_subsys_index, common_timebase
 from .lti import LTI, _process_frequency_response
 
 __all__ = ['FrequencyResponseData', 'FRD', 'frd']
@@ -33,8 +35,8 @@ class FrequencyResponseData(LTI):
 
     The FrequencyResponseData (FRD) class is used to represent systems in
     frequency response data form.  It can be created manually using the
-    class constructor, using the :func:~~control.frd` factory function
-    (preferred), or via the :func:`~control.frequency_response` function.
+    class constructor, using the :func:`~control.frd` factory function or
+    via the :func:`~control.frequency_response` function.
 
     Parameters
     ----------
@@ -65,6 +67,28 @@ class constructor, using the :func:~~control.frd` factory function
         frequency point.
     dt : float, True, or None
         System timebase.
+    squeeze : bool
+        By default, if a system is single-input, single-output (SISO) then
+        the outputs (and inputs) are returned as a 1D array (indexed by
+        frequency) and if a system is multi-input or multi-output, then the
+        outputs are returned as a 2D array (indexed by output and
+        frequency) or a 3D array (indexed by output, trace, and frequency).
+        If ``squeeze=True``, access to the output response will remove
+        single-dimensional entries from the shape of the inputs and outputs
+        even if the system is not SISO. If ``squeeze=False``, the output is
+        returned as a 3D array (indexed by the output, input, and
+        frequency) even if the system is SISO. The default value can be set
+        using config.defaults['control.squeeze_frequency_response'].
+    ninputs, noutputs, nstates : int
+        Number of inputs, outputs, and states of the underlying system.
+    input_labels, output_labels : array of str
+        Names for the input and output variables.
+    sysname : str, optional
+        Name of the system.  For data generated using
+        :func:`~control.frequency_response`, stores the name of the system
+        that created the data.
+    title : str, optional
+        Set the title to use when plotting.
 
     See Also
     --------
@@ -89,6 +113,20 @@ class constructor, using the :func:~~control.frd` factory function
     the imaginary access).  See :meth:`~control.FrequencyResponseData.__call__`
     for a more detailed description.
 
+    A state space system is callable and returns the value of the transfer
+    function evaluated at a point in the complex plane.  See
+    :meth:`~control.StateSpace.__call__` for a more detailed description.
+
+    Subsystem response corresponding to selected input/output pairs can be
+    created by indexing the frequency response data object::
+
+        subsys = sys[output_spec, input_spec]
+
+    The input and output specifications can be single integers, lists of
+    integers, or slices.  In addition, the strings representing the names
+    of the signals can be used and will be replaced with the equivalent
+    signal offsets.
+
     """
     #
     # Class attributes
@@ -243,21 +281,72 @@ def __init__(self, *args, **kwargs):
 
     @property
     def magnitude(self):
-        return np.abs(self.fresp)
+        """Magnitude of the frequency response.
+
+        Magnitude of the frequency response, indexed by either the output
+        and frequency (if only a single input is given) or the output,
+        input, and frequency (for multi-input systems).  See
+        :attr:`FrequencyResponseData.squeeze` for a description of how this
+        can be modified using the `squeeze` keyword.
+
+        Input and output signal names can be used to index the data in
+        place of integer offsets.
+
+        :type: 1D, 2D, or 3D array
+
+        """
+        return NamedSignal(
+            np.abs(self.fresp), self.output_labels, self.input_labels)
 
     @property
     def phase(self):
-        return np.angle(self.fresp)
+        """Phase of the frequency response.
+
+        Phase of the frequency response in radians/sec, indexed by either
+        the output and frequency (if only a single input is given) or the
+        output, input, and frequency (for multi-input systems).  See
+        :attr:`FrequencyResponseData.squeeze` for a description of how this
+        can be modified using the `squeeze` keyword.
+
+        Input and output signal names can be used to index the data in
+        place of integer offsets.
+
+        :type: 1D, 2D, or 3D array
+
+        """
+        return NamedSignal(
+            np.angle(self.fresp), self.output_labels, self.input_labels)
 
     @property
     def frequency(self):
+        """Frequencies at which the response is evaluated.
+
+        :type: 1D array
+
+        """
         return self.omega
 
     @property
     def response(self):
-        return self.fresp
+        """Complex value of the frequency response.
+
+        Value of the frequency response as a complex number, indexed by
+        either the output and frequency (if only a single input is given)
+        or the output, input, and frequency (for multi-input systems).  See
+        :attr:`FrequencyResponseData.squeeze` for a description of how this
+        can be modified using the `squeeze` keyword.
+
+        Input and output signal names can be used to index the data in
+        place of integer offsets.
+
+        :type: 1D, 2D, or 3D array
+
+        """
+        return NamedSignal(
+            self.fresp, self.output_labels, self.input_labels)
 
     def __str__(self):
+
         """String representation of the transfer function."""
 
         mimo = self.ninputs > 1 or self.noutputs > 1
@@ -593,9 +682,25 @@ def __iter__(self):
             return iter((self.omega, fresp))
         return iter((np.abs(fresp), np.angle(fresp), self.omega))
 
-    # Implement (thin) getitem to allow access via legacy indexing
-    def __getitem__(self, index):
-        return list(self.__iter__())[index]
+    def __getitem__(self, key):
+        if not isinstance(key, Iterable) or len(key) != 2:
+            # Implement (thin) getitem to allow access via legacy indexing
+            return list(self.__iter__())[key]
+
+        # Convert signal names to integer offsets (via NamedSignal object)
+        iomap = NamedSignal(
+            self.fresp[:, :, 0], self.output_labels, self.input_labels)
+        indices = iomap._parse_key(key, level=1)  # ignore index checks
+        outdx, outputs = _process_subsys_index(indices[0], self.output_labels)
+        inpdx, inputs = _process_subsys_index(indices[1], self.input_labels)
+
+        # Create the system name
+        sysname = config.defaults['iosys.indexed_system_name_prefix'] + \
+            self.name + config.defaults['iosys.indexed_system_name_suffix']
+
+        return FrequencyResponseData(
+            self.fresp[outdx, :][:, inpdx], self.omega, self.dt,
+            inputs=inputs, outputs=outputs, name=sysname)
 
     # Implement (thin) len to emulate legacy testing interface
     def __len__(self):
diff --git a/control/iosys.py b/control/iosys.py
index 9092b672b..dd1566eb9 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -13,9 +13,10 @@
 import numpy as np
 
 from . import config
+from .exception import ControlIndexError
 
-__all__ = ['InputOutputSystem', 'issiso', 'timebase', 'common_timebase',
-           'isdtime', 'isctime']
+__all__ = ['InputOutputSystem', 'NamedSignal', 'issiso', 'timebase',
+           'common_timebase', 'isdtime', 'isctime']
 
 # Define module default parameter values
 _iosys_defaults = {
@@ -33,6 +34,69 @@
 }
 
 
+# Named signal class
+class NamedSignal(np.ndarray):
+    def __new__(cls, input_array, signal_labels=None, trace_labels=None):
+        # See https://numpy.org/doc/stable/user/basics.subclassing.html
+        obj = np.asarray(input_array).view(cls)     # Cast to our class type
+        obj.signal_labels = signal_labels           # Save signal labels
+        obj.trace_labels = trace_labels             # Save trace labels
+        obj.data_shape = input_array.shape          # Save data shape
+        return obj                                  # Return new object
+
+    def __array_finalize__(self, obj):
+        # See https://numpy.org/doc/stable/user/basics.subclassing.html
+        if obj is None:
+            return
+        self.signal_labels = getattr(obj, 'signal_labels', None)
+        self.trace_labels = getattr(obj, 'trace_labels', None)
+        self.data_shape = getattr(obj, 'data_shape', None)
+
+    def _parse_key(self, key, labels=None, level=0):
+        if labels is None:
+            labels = self.signal_labels
+        try:
+            if isinstance(key, str):
+                key = labels.index(item := key)
+                if level == 0 and len(self.data_shape) < 2:
+                    raise ControlIndexError
+            elif isinstance(key, list):
+                keylist = []
+                for item in key:        # use for loop to save item for error
+                    keylist.append(
+                        self._parse_key(item, labels=labels, level=level+1))
+                if level == 0 and key != keylist and len(self.data_shape) < 2:
+                    raise ControlIndexError
+                key = keylist
+            elif isinstance(key, tuple) and len(key) > 0:
+                keylist = []
+                keylist.append(
+                    self._parse_key(
+                        item := key[0], labels=self.signal_labels,
+                        level=level+1))
+                if len(key) > 1:
+                    keylist.append(
+                        self._parse_key(
+                            item := key[1], labels=self.trace_labels,
+                            level=level+1))
+                if level == 0 and key[:len(keylist)] != tuple(keylist) \
+                   and len(keylist) > len(self.data_shape) - 1:
+                    raise ControlIndexError
+                for i in range(2, len(key)):
+                    keylist.append(key[i])      # pass on remaining elements
+                key = tuple(keylist)
+        except ValueError:
+            raise ValueError(f"unknown signal name '{item}'")
+        except ControlIndexError:
+            raise ControlIndexError(
+                "signal name(s) not valid for squeezed data")
+
+        return key
+
+    def __getitem__(self, key):
+        return super().__getitem__(self._parse_key(key))
+
+
 class InputOutputSystem(object):
     """A class for representing input/output systems.
 
@@ -965,3 +1029,31 @@ def _parse_spec(syslist, spec, signame, dictname=None):
             ValueError(f"signal index '{index}' is out of range")
 
     return system_index, signal_indices, gain
+
+
+#
+# Utility function for processing subsystem indices
+#
+# This function processes an index specification (int, list, or slice) and
+# returns a index specification that can be used to create a subsystem
+#
+def _process_subsys_index(idx, sys_labels, slice_to_list=False):
+    if not isinstance(idx, (slice, list, int)):
+        raise TypeError("system indices must be integers, slices, or lists")
+
+    # Convert singleton lists to integers for proper slicing (below)
+    if isinstance(idx, (list, tuple)) and len(idx) == 1:
+        idx = idx[0]
+
+    # Convert int to slice so that numpy doesn't drop dimension
+    if isinstance(idx, int):
+        idx = slice(idx, idx+1, 1)
+
+    # Get label names (taking care of possibility that we were passed a list)
+    labels = [sys_labels[i] for i in idx] if isinstance(idx, list) \
+        else sys_labels[idx]
+
+    if slice_to_list and isinstance(idx, slice):
+        idx = range(len(sys_labels))[idx]
+
+    return idx, labels
diff --git a/control/nlsys.py b/control/nlsys.py
index 835c16ef6..62e4bf78e 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -193,8 +193,8 @@ def __call__(sys, u, params=None, squeeze=None):
 
         # Evaluate the function on the argument
         out = sys._out(0, np.array((0,)), np.asarray(u))
-        _, out = _process_time_response(
-            None, out, issiso=sys.issiso(), squeeze=squeeze)
+        out = _process_time_response(
+            out, issiso=sys.issiso(), squeeze=squeeze)
         return out
 
     def __mul__(self, other):
diff --git a/control/statesp.py b/control/statesp.py
index aa1c7221b..bfe5f996b 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -48,15 +48,15 @@
 """
 
 import math
+from collections.abc import Iterable
 from copy import deepcopy
 from warnings import warn
-from collections.abc import Iterable
 
 import numpy as np
 import scipy as sp
 import scipy.linalg
-from numpy import (any, asarray, concatenate, cos, delete, empty, exp, eye,
-                   isinf, ones, pad, sin, squeeze, zeros)
+from numpy import any, asarray, concatenate, cos, delete, empty, exp, eye, \
+    isinf, ones, pad, sin, squeeze, zeros
 from numpy.linalg import LinAlgError, eigvals, matrix_rank, solve
 from numpy.random import rand, randn
 from scipy.signal import StateSpace as signalStateSpace
@@ -65,9 +65,9 @@
 from . import config
 from .exception import ControlMIMONotImplemented, ControlSlycot, slycot_check
 from .frdata import FrequencyResponseData
-from .iosys import (InputOutputSystem, _process_dt_keyword,
-                    _process_iosys_keywords, _process_signal_list,
-                    common_timebase, isdtime, issiso)
+from .iosys import InputOutputSystem, NamedSignal, _process_dt_keyword, \
+    _process_iosys_keywords, _process_signal_list, _process_subsys_index, \
+    common_timebase, isdtime, issiso
 from .lti import LTI, _process_frequency_response
 from .nlsys import InterconnectedSystem, NonlinearIOSystem
 
@@ -153,6 +153,17 @@ class StateSpace(NonlinearIOSystem, LTI):
     function evaluated at a point in the complex plane.  See
     :meth:`~control.StateSpace.__call__` for a more detailed description.
 
+    Subsystems corresponding to selected input/output pairs can be
+    created by indexing the state space system::
+
+        subsys = sys[output_spec, input_spec]
+
+    The input and output specifications can be single integers, lists of
+    integers, or slices.  In addition, the strings representing the names
+    of the signals can be used and will be replaced with the equivalent
+    signal offsets.  The subsystem is created by truncating the inputs and
+    outputs, but leaving the full set of system states.
+
     StateSpace instances have support for IPython LaTeX output,
     intended for pretty-printing in Jupyter notebooks.  The LaTeX
     output can be configured using
@@ -1214,25 +1225,25 @@ def append(self, other):
         D[self.noutputs:, self.ninputs:] = other.D
         return StateSpace(A, B, C, D, self.dt)
 
-    def __getitem__(self, indices):
+    def __getitem__(self, key):
         """Array style access"""
-        if not isinstance(indices, Iterable) or len(indices) != 2:
-            raise IOError('must provide indices of length 2 for state space')
-        outdx, inpdx = indices
-        
-        # Convert int to slice to ensure that numpy doesn't drop the dimension
-        if isinstance(outdx, int): outdx = slice(outdx, outdx+1, 1)
-        if isinstance(inpdx, int): inpdx = slice(inpdx, inpdx+1, 1)
+        if not isinstance(key, Iterable) or len(key) != 2:
+            raise IOError("must provide indices of length 2 for state space")
 
-        if not isinstance(outdx, slice) or not isinstance(inpdx, slice):
-            raise TypeError(f"system indices must be integers or slices")
+        # Convert signal names to integer offsets
+        iomap = NamedSignal(self.D, self.output_labels, self.input_labels)
+        indices = iomap._parse_key(key, level=1)  # ignore index checks
+        outdx, output_labels = _process_subsys_index(
+            indices[0], self.output_labels)
+        inpdx, input_labels = _process_subsys_index(
+            indices[1], self.input_labels)
 
         sysname = config.defaults['iosys.indexed_system_name_prefix'] + \
             self.name + config.defaults['iosys.indexed_system_name_suffix']
         return StateSpace(
-            self.A, self.B[:, inpdx], self.C[outdx, :], self.D[outdx, inpdx],
-            self.dt, name=sysname, 
-            inputs=self.input_labels[inpdx], outputs=self.output_labels[outdx])
+            self.A, self.B[:, inpdx], self.C[outdx, :],
+            self.D[outdx, :][:, inpdx], self.dt,
+            name=sysname, inputs=input_labels, outputs=output_labels)
 
     def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None,
                name=None, copy_names=True, **kwargs):
diff --git a/control/timeresp.py b/control/timeresp.py
index 5813c166d..072db60de 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -80,7 +80,7 @@
 
 from . import config
 from .exception import pandas_check
-from .iosys import isctime, isdtime
+from .iosys import NamedSignal, isctime, isdtime
 from .timeplot import time_response_plot
 
 __all__ = ['forced_response', 'step_response', 'step_info',
@@ -204,34 +204,46 @@ class TimeResponseData:
 
     Notes
     -----
-    1. For backward compatibility with earlier versions of python-control,
-       this class has an ``__iter__`` method that allows it to be assigned
-       to a tuple with a variable number of elements.  This allows the
-       following patterns to work:
+    The responses for individual elements of the time response can be
+    accessed using integers, slices, or lists of signal offsets or the
+    names of the appropriate signals::
 
-         t, y = step_response(sys)
-         t, y, x = step_response(sys, return_x=True)
+      sys = ct.rss(4, 2, 1)
+      resp = ct.initial_response(sys, X0=[1, 1, 1, 1])
+      plt.plot(resp.time, resp.outputs['y[0]'])
 
-       When using this (legacy) interface, the state vector is not affected by
-       the `squeeze` parameter.
+    In the case of multi-trace data, the responses should be indexed using
+    the output signal name (or offset) and the input signal name (or
+    offset)::
 
-    2. For backward compatibility with earlier version of python-control,
-       this class has ``__getitem__`` and ``__len__`` methods that allow the
-       return value to be indexed:
+      sys = ct.rss(4, 2, 2, strictly_proper=True)
+      resp = ct.step_response(sys)
+      plt.plot(resp.time, resp.outputs[['y[0]', 'y[1]'], 'u[0]'].T)
 
-         response[0]: returns the time vector
-         response[1]: returns the output vector
-         response[2]: returns the state vector
+    For backward compatibility with earlier versions of python-control,
+    this class has an ``__iter__`` method that allows it to be assigned to
+    a tuple with a variable number of elements.  This allows the following
+    patterns to work::
 
-       When using this (legacy) interface, the state vector is not affected by
-       the `squeeze` parameter.
+       t, y = step_response(sys)
+       t, y, x = step_response(sys, return_x=True)
 
-    3. The default settings for ``return_x``, ``squeeze`` and ``transpose``
-       can be changed by calling the class instance and passing new values:
+    Similarly, the class has ``__getitem__`` and ``__len__`` methods that
+    allow the return value to be indexed:
+
+    * response[0]: returns the time vector
+    * response[1]: returns the output vector
+    * response[2]: returns the state vector
+
+    When using this (legacy) interface, the state vector is not affected
+    by the `squeeze` parameter.
+
+    The default settings for ``return_x``, ``squeeze`` and ``transpose``
+    can be changed by calling the class instance and passing new values::
 
          response(tranpose=True).input
 
-       See :meth:`TimeResponseData.__call__` for more information.
+    See :meth:`TimeResponseData.__call__` for more information.
 
     """
 
@@ -564,13 +576,18 @@ def outputs(self):
         (for multiple traces).  See :attr:`TimeResponseData.squeeze` for a
         description of how this can be modified using the `squeeze` keyword.
 
+        Input and output signal names can be used to index the data in
+        place of integer offsets, with the input signal names being used to
+        access multi-input data.
+
         :type: 1D, 2D, or 3D array
 
         """
-        t, y = _process_time_response(
-            self.t, self.y, issiso=self.issiso,
+        # TODO: move to __init__ to avoid recomputing each time?
+        y = _process_time_response(
+            self.y, issiso=self.issiso,
             transpose=self.transpose, squeeze=self.squeeze)
-        return y
+        return NamedSignal(y, self.output_labels, self.input_labels)
 
     # Getter for states (implements squeeze processing)
     @property
@@ -583,30 +600,25 @@ def states(self):
         for a description of how this can be modified using the `squeeze`
         keyword.
 
+        Input and output signal names can be used to index the data in
+        place of integer offsets, with the input signal names being used to
+        access multi-input data.
+
         :type: 2D or 3D array
 
         """
-        if self.x is None:
-            return None
-
-        elif self.squeeze is True:
-            x = self.x.squeeze()
+        # TODO: move to __init__ to avoid recomputing each time?
+        x = _process_time_response(
+            self.x, transpose=self.transpose,
+            squeeze=self.squeeze, issiso=False)
 
-        elif self.ninputs == 1 and self.noutputs == 1 and \
-             self.ntraces == 1 and self.x.ndim == 3 and \
+        # Special processing for SISO case: always retain state index
+        if self.issiso and self.ntraces == 1 and x.ndim == 3 and \
              self.squeeze is not False:
             # Single-input, single-output system with single trace
-            x = self.x[:, 0, :]
+            x = x[:, 0, :]
 
-        else:
-            # Return the full set of data
-            x = self.x
-
-        # Transpose processing
-        if self.transpose:
-            x = np.transpose(x, np.roll(range(x.ndim), 1))
-
-        return x
+        return NamedSignal(x, self.state_labels, self.input_labels)
 
     # Getter for inputs (implements squeeze processing)
     @property
@@ -621,6 +633,10 @@ def inputs(self):
         the two.  If a 3D vector is passed, then it represents a multi-trace,
         multi-input signal, indexed by input, trace, and time.
 
+        Input and output signal names can be used to index the data in
+        place of integer offsets, with the input signal names being used to
+        access multi-input data.
+
         See :attr:`TimeResponseData.squeeze` for a description of how the
         dimensions of the input vector can be modified using the `squeeze`
         keyword.
@@ -628,15 +644,17 @@ def inputs(self):
         :type: 1D or 2D array
 
         """
+        # TODO: move to __init__ to avoid recomputing each time?
         if self.u is None:
             return None
 
-        t, u = _process_time_response(
-            self.t, self.u, issiso=self.issiso,
+        u = _process_time_response(
+            self.u, issiso=self.issiso,
             transpose=self.transpose, squeeze=self.squeeze)
-        return u
+        return NamedSignal(u, self.input_labels, self.input_labels)
 
     # Getter for legacy state (implements non-standard squeeze processing)
+    # TODO: remove when no longer needed
     @property
     def _legacy_states(self):
         """Time response state vector (legacy version).
@@ -1265,7 +1283,7 @@ def forced_response(sysdata, T=None, U=0., X0=0., transpose=False, params=None,
 
 # Process time responses in a uniform way
 def _process_time_response(
-        tout, yout, issiso=False, transpose=None, squeeze=None):
+        signal, issiso=False, transpose=None, squeeze=None):
     """Process time response signals.
 
     This function processes the outputs (or inputs) of time response
@@ -1273,43 +1291,36 @@ def _process_time_response(
 
     Parameters
     ----------
-    T : 1D array
-        Time values of the output.  Ignored if None.
-
-    yout : ndarray
-        Response of the system.  This can either be a 1D array indexed by time
-        (for SISO systems), a 2D array indexed by output and time (for MIMO
-        systems with no input indexing, such as initial_response or forced
-        response) or a 3D array indexed by output, input, and time.
+    signal : ndarray
+        Data to be processed.  This can either be a 1D array indexed by
+        time (for SISO systems), a 2D array indexed by output and time (for
+        MIMO systems with no input indexing, such as initial_response or
+        forced response) or a 3D array indexed by output, input, and time.
 
     issiso : bool, optional
         If ``True``, process data as single-input, single-output data.
         Default is ``False``.
 
     transpose : bool, optional
-        If True, transpose all input and output arrays (for backward
-        compatibility with MATLAB and :func:`scipy.signal.lsim`).  Default
-        value is False.
+        If True, transpose data (for backward compatibility with MATLAB and
+        :func:`scipy.signal.lsim`).  Default value is False.
 
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then the
-        output response is returned as a 1D array (indexed by time).  If
+        signals are returned as a 1D array (indexed by time).  If
         squeeze=True, remove single-dimensional entries from the shape of the
-        output even if the system is not SISO. If squeeze=False, keep the
-        output as a 3D array (indexed by the output, input, and time) even if
+        signal even if the system is not SISO. If squeeze=False, keep the
+        signal as a 3D array (indexed by the output, input, and time) even if
         the system is SISO. The default value can be set using
         config.defaults['control.squeeze_time_response'].
 
     Returns
     -------
-    T : 1D array
-        Time values of the output.
-
-    yout : ndarray
-        Response of the system.  If the system is SISO and squeeze is not
-        True, the array is 1D (indexed by time).  If the system is not SISO or
-        squeeze is False, the array is either 2D (indexed by output and time)
-        or 3D (indexed by input, output, and time).
+    output : ndarray
+        Processed signal.  If the system is SISO and squeeze is not True,
+        the array is 1D (indexed by time).  If the system is not SISO or
+        squeeze is False, the array is either 2D (indexed by output and
+        time) or 3D (indexed by input, output, and time).
 
     """
     # If squeeze was not specified, figure out the default (might remain None)
@@ -1317,29 +1328,26 @@ def _process_time_response(
         squeeze = config.defaults['control.squeeze_time_response']
 
     # Figure out whether and how to squeeze output data
-    if squeeze is True:         # squeeze all dimensions
-        yout = np.squeeze(yout)
-    elif squeeze is False:      # squeeze no dimensions
+    if squeeze is True:                 # squeeze all dimensions
+        signal = np.squeeze(signal)
+    elif squeeze is False:              # squeeze no dimensions
         pass
-    elif squeeze is None:       # squeeze signals if SISO
+    elif squeeze is None:               # squeeze signals if SISO
         if issiso:
-            if yout.ndim == 3:
-                yout = yout[0][0]       # remove input and output
+            if signal.ndim == 3:
+                signal = signal[0][0]   # remove input and output
             else:
-                yout = yout[0]          # remove input
+                signal = signal[0]      # remove input
     else:
         raise ValueError("Unknown squeeze value")
 
     # See if we need to transpose the data back into MATLAB form
     if transpose:
-        # Transpose time vector in case we are using np.matrix
-        tout = np.transpose(tout)
-
         # For signals, put the last index (time) into the first slot
-        yout = np.transpose(yout, np.roll(range(yout.ndim), 1))
+        signal = np.transpose(signal, np.roll(range(signal.ndim), 1))
 
-    # Return time, output, and (optionally) state
-    return tout, yout
+    # Return output
+    return signal
 
 
 def step_response(
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 499359cbc..56ec7395f 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -48,25 +48,26 @@
 """
 
 from collections.abc import Iterable
+from copy import deepcopy
+from itertools import chain
+from re import sub
+from warnings import warn
 
 # External function declarations
 import numpy as np
-from numpy import angle, array, empty, finfo, ndarray, ones, \
-    polyadd, polymul, polyval, roots, sqrt, zeros, squeeze, exp, pi, \
-    where, delete, real, poly, nonzero
 import scipy as sp
-from scipy.signal import tf2zpk, zpk2tf, cont2discrete
+from numpy import angle, array, delete, empty, exp, finfo, ndarray, nonzero, \
+    ones, pi, poly, polyadd, polymul, polyval, real, roots, sqrt, squeeze, \
+    where, zeros
 from scipy.signal import TransferFunction as signalTransferFunction
-from copy import deepcopy
-from warnings import warn
-from itertools import chain
-from re import sub
-from .lti import LTI, _process_frequency_response
-from .iosys import InputOutputSystem, common_timebase, isdtime, \
-    _process_iosys_keywords
+from scipy.signal import cont2discrete, tf2zpk, zpk2tf
+
+from . import config
 from .exception import ControlMIMONotImplemented
 from .frdata import FrequencyResponseData
-from . import config
+from .iosys import InputOutputSystem, NamedSignal, _process_iosys_keywords, \
+    _process_subsys_index, common_timebase, isdtime
+from .lti import LTI, _process_frequency_response
 
 __all__ = ['TransferFunction', 'tf', 'zpk', 'ss2tf', 'tfdata']
 
@@ -146,6 +147,16 @@ class TransferFunction(LTI):
     function evaluated at a point in the complex plane.  See
     :meth:`~control.TransferFunction.__call__` for a more detailed description.
 
+    Subsystems corresponding to selected input/output pairs can be
+    created by indexing the transfer function::
+
+        subsys = sys[output_spec, input_spec]
+
+    The input and output specifications can be single integers, lists of
+    integers, or slices.  In addition, the strings representing the names
+    of the signals can be used and will be replaced with the equivalent
+    signal offsets.
+
     The TransferFunction class defines two constants ``s`` and ``z`` that
     represent the differentiation and delay operators in continuous and
     discrete time.  These can be used to create variables that allow algebraic
@@ -761,48 +772,30 @@ def __pow__(self, other):
 
     def __getitem__(self, key):
         if not isinstance(key, Iterable) or len(key) != 2:
-            raise IOError('must provide indices of length 2 for transfer functions')
+            raise IOError(
+                "must provide indices of length 2 for transfer functions")
+
+        # Convert signal names to integer offsets (via NamedSignal object)
+        iomap = NamedSignal(
+            np.empty((self.noutputs, self.ninputs)),
+            self.output_labels, self.input_labels)
+        indices = iomap._parse_key(key, level=1)  # ignore index checks
+        outdx, outputs = _process_subsys_index(
+            indices[0], self.output_labels, slice_to_list=True)
+        inpdx, inputs = _process_subsys_index(
+            indices[1], self.input_labels, slice_to_list=True)
         
-        key1, key2 = key
-        if not isinstance(key1, (int, slice)) or not isinstance(key2, (int, slice)):
-            raise TypeError(f"system indices must be integers or slices")
-
-        # pre-process
-        if isinstance(key1, int):
-            key1 = slice(key1, key1 + 1, 1)
-        if isinstance(key2, int):
-            key2 = slice(key2, key2 + 1, 1)
-        # dim1
-        start1, stop1, step1 = key1.start, key1.stop, key1.step
-        if step1 is None:
-            step1 = 1
-        if start1 is None:
-            start1 = 0
-        if stop1 is None:
-            stop1 = len(self.num)
-        # dim1
-        start2, stop2, step2 = key2.start, key2.stop, key2.step
-        if step2 is None:
-            step2 = 1
-        if start2 is None:
-            start2 = 0
-        if stop2 is None:
-            stop2 = len(self.num[0])
-
+        # Construct the transfer function for the subsyste
         num, den = [], []
-        for i in range(start1, stop1, step1):
+        for i in outdx:
             num_i = []
             den_i = []
-            for j in range(start2, stop2, step2):
+            for j in inpdx:
                 num_i.append(self.num[i][j])
                 den_i.append(self.den[i][j])
             num.append(num_i)
             den.append(den_i)
 
-        # Save the label names
-        outputs = [self.output_labels[i] for i in range(start1, stop1, step1)]
-        inputs = [self.input_labels[j] for j in range(start2, stop2, step2)]
-
         # Create the system name
         sysname = config.defaults['iosys.indexed_system_name_prefix'] + \
             self.name + config.defaults['iosys.indexed_system_name_suffix']
diff --git a/doc/conventions.rst b/doc/conventions.rst
index fb1f0715f..d2394e040 100644
--- a/doc/conventions.rst
+++ b/doc/conventions.rst
@@ -98,6 +98,20 @@ argument::
 
     mag, phase, omega = response(squeeze=False)
 
+Frequency response objects are also available as named properties of the
+``response`` object: ``response.magnitude``, ``response.phase``, and
+``response.response`` (for the complex response).  For MIMO systems, these
+elements of the frequency response can be accessed using the names of the
+inputs and outputs::
+
+  response.magnitude['y[0]', 'u[1]']
+
+where the signal names are based on the system that generated the frequency
+response.
+
+Note: The ``fresp`` data member is stored as a NumPy array and cannot be
+accessed with signal names.  Use ``response.response`` to access the
+complex frequency response using signal names.
 
 Discrete time systems
 ---------------------
@@ -132,6 +146,21 @@ constructor for the desired data type using the original system as the sole
 argument or using the explicit conversion functions :func:`ss2tf` and
 :func:`tf2ss`.
 
+Subsystems
+----------
+Subsets of input/output pairs for LTI systems can be obtained by indexing
+the system using either numerical indices (including slices) or signal
+names::
+
+    subsys = sys[[0, 2], 0:2]
+    subsys = sys[['y[0]', 'y[2]'], ['u[0]', 'u[1]']]
+
+Signal names for an indexed subsystem are preserved from the original
+system and the subsystem name is set according to the values of
+``control.config.defaults['iosys.indexed_system_name_prefix']`` and
+``control.config.defaults['iosys.indexed_system_name_suffix']``.  The default
+subsystem name is the original system name with '$indexed' appended.
+
 Simulating LTI systems
 ======================
 
@@ -233,7 +262,7 @@ properties::
 
     sys = ct.rss(4, 1, 1)
     response = ct.step_response(sys)
-    plot(response.time, response.outputs)
+    plt.plot(response.time, response.outputs)
 
 The dimensions of the response properties depend on the function being
 called and whether the system is SISO or MIMO.  In addition, some time
@@ -242,6 +271,17 @@ such as the :func:`step_response` function applied to a MIMO system,
 which will compute the step response for each input/output pair.  See
 :class:`TimeResponseData` for more details.
 
+The input, output, and state elements of the response can be accessed using
+signal names in place of integer offsets::
+
+    plt.plot(response. time, response.states['x[1]']
+
+For multi-trace systems generated by :func:`step_response` and
+:func:`impulse_response`, the input name used to generate the trace can be
+used to access the appropriate input output pair::
+
+    plt.plot(response.time, response.outputs['y[0]', 'u[1]'])
+
 The time response functions can also be assigned to a tuple, which extracts
 the time and output (and optionally the state, if the `return_x` keyword is
 used).  This allows simple commands for plotting::
diff --git a/doc/plotting.rst b/doc/plotting.rst
index 367e2c349..cc292bdbe 100644
--- a/doc/plotting.rst
+++ b/doc/plotting.rst
@@ -79,6 +79,11 @@ the data from the simulation::
       for j in range(2):
           axs[i, j].plot(time, outputs[i, j])
 
+In addition to accessing time response data via integer indices, signal
+names can allow be used::
+
+  plt.plot(response.time, response.outputs['y[0]', 'u[1]'])
+
 A number of options are available in the `plot` method to customize
 the appearance of input output data.  For data produced by the
 :func:`~control.impulse_response` and :func:`~control.step_response`
@@ -278,6 +283,19 @@ maximum frequencies in the (log-spaced) frequency range::
 The number of (log-spaced) points in the frequency can be specified using
 the ``omega_num`` keyword parameter.
 
+Frequency response data can also be accessed directly and plotted manually::
+
+  sys = ct.rss(4, 2, 2, strictly_proper=True)  # 2x2 MIMO system
+  fresp = ct.frequency_response(sys)
+  plt.loglog(fresp.omega, fresp.magnitude['y[1]', 'u[0]'])
+
+Access to frequency response data is available via the attributes
+``omega``, ``magnitude``,` `phase``, and ``response``, where ``response``
+represents the complex value of the frequency response at each frequency.
+The ``magnitude``, ``phase``, and ``response`` arrays can be indexed using
+either input/output indices or signal names, with the first index
+corresponding to the output signal and the second input corresponding to
+the input signal.
 
 Pole/zero data
 ==============


Test patch:
diff --git a/control/tests/bdalg_test.py b/control/tests/bdalg_test.py
index 5629f27f9..8ea67e0f7 100644
--- a/control/tests/bdalg_test.py
+++ b/control/tests/bdalg_test.py
@@ -12,6 +12,7 @@
 from control.statesp import StateSpace
 from control.bdalg import feedback, append, connect
 from control.lti import zeros, poles
+from control.bdalg import _ensure_tf
 
 
 class TestFeedback:
@@ -362,3 +363,557 @@ def test_bdalg_udpate_names_errors():
 
     with pytest.raises(TypeError, match="unrecognized keywords"):
         sys = ctrl.series(sys1, sys2, dt=1)
+
+
+class TestEnsureTf:
+    """Test ``_ensure_tf``."""
+
+    @pytest.mark.parametrize(
+        "arraylike_or_tf, dt, tf",
+        [
+            (
+                ctrl.TransferFunction([1], [1, 2, 3]),
+                None,
+                ctrl.TransferFunction([1], [1, 2, 3]),
+            ),
+            (
+                ctrl.TransferFunction([1], [1, 2, 3]),
+                0,
+                ctrl.TransferFunction([1], [1, 2, 3]),
+            ),
+            (
+                2,
+                None,
+                ctrl.TransferFunction([2], [1]),
+            ),
+            (
+                np.array([2]),
+                None,
+                ctrl.TransferFunction([2], [1]),
+            ),
+            (
+                np.array([[2]]),
+                None,
+                ctrl.TransferFunction([2], [1]),
+            ),
+            (
+                np.array(
+                    [
+                        [2, 0, 3],
+                        [1, 2, 3],
+                    ]
+                ),
+                None,
+                ctrl.TransferFunction(
+                    [
+                        [[2], [0], [3]],
+                        [[1], [2], [3]],
+                    ],
+                    [
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                    ],
+                ),
+            ),
+            (
+                np.array([2, 0, 3]),
+                None,
+                ctrl.TransferFunction(
+                    [
+                        [[2], [0], [3]],
+                    ],
+                    [
+                        [[1], [1], [1]],
+                    ],
+                ),
+            ),
+        ],
+    )
+    def test_ensure(self, arraylike_or_tf, dt, tf):
+        """Test nominal cases"""
+        ensured_tf = _ensure_tf(arraylike_or_tf, dt)
+        assert _tf_close_coeff(tf, ensured_tf)
+
+    @pytest.mark.parametrize(
+        "arraylike_or_tf, dt, exception",
+        [
+            (
+                ctrl.TransferFunction([1], [1, 2, 3]),
+                0.1,
+                ValueError,
+            ),
+            (
+                ctrl.TransferFunction([1], [1, 2, 3], 0.1),
+                0,
+                ValueError,
+            ),
+            (
+                np.ones((1, 1, 1)),
+                None,
+                ValueError,
+            ),
+            (
+                np.ones((1, 1, 1, 1)),
+                None,
+                ValueError,
+            ),
+        ],
+    )
+    def test_error_ensure(self, arraylike_or_tf, dt, exception):
+        """Test error cases"""
+        with pytest.raises(exception):
+            _ensure_tf(arraylike_or_tf, dt)
+
+
+class TestTfCombineSplit:
+    """Test ``combine_tf`` and ``split_tf``."""
+
+    @pytest.mark.parametrize(
+        "tf_array, tf",
+        [
+            # Continuous-time
+            (
+                [
+                    [ctrl.TransferFunction([1], [1, 1])],
+                    [ctrl.TransferFunction([2], [1, 0])],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[1]],
+                        [[2]],
+                    ],
+                    [
+                        [[1, 1]],
+                        [[1, 0]],
+                    ],
+                ),
+            ),
+            # Discrete-time
+            (
+                [
+                    [ctrl.TransferFunction([1], [1, 1], dt=1)],
+                    [ctrl.TransferFunction([2], [1, 0], dt=1)],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[1]],
+                        [[2]],
+                    ],
+                    [
+                        [[1, 1]],
+                        [[1, 0]],
+                    ],
+                    dt=1,
+                ),
+            ),
+            # Scalar
+            (
+                [
+                    [2],
+                    [ctrl.TransferFunction([2], [1, 0])],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[2]],
+                        [[2]],
+                    ],
+                    [
+                        [[1]],
+                        [[1, 0]],
+                    ],
+                ),
+            ),
+            # Matrix
+            (
+                [
+                    [np.eye(3)],
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [0], [3]],
+                                [[1], [2], [3]],
+                            ],
+                            [
+                                [[1], [1], [1]],
+                                [[1], [1], [1]],
+                            ],
+                        )
+                    ],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[1], [0], [0]],
+                        [[0], [1], [0]],
+                        [[0], [0], [1]],
+                        [[2], [0], [3]],
+                        [[1], [2], [3]],
+                    ],
+                    [
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                    ],
+                ),
+            ),
+            # Inhomogeneous
+            (
+                [
+                    [np.eye(3)],
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [0]],
+                                [[1], [2]],
+                            ],
+                            [
+                                [[1], [1]],
+                                [[1], [1]],
+                            ],
+                        ),
+                        ctrl.TransferFunction(
+                            [
+                                [[3]],
+                                [[3]],
+                            ],
+                            [
+                                [[1]],
+                                [[1]],
+                            ],
+                        ),
+                    ],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[1], [0], [0]],
+                        [[0], [1], [0]],
+                        [[0], [0], [1]],
+                        [[2], [0], [3]],
+                        [[1], [2], [3]],
+                    ],
+                    [
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                    ],
+                ),
+            ),
+            # Discrete-time
+            (
+                [
+                    [2],
+                    [ctrl.TransferFunction([2], [1, 0], dt=0.1)],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[2]],
+                        [[2]],
+                    ],
+                    [
+                        [[1]],
+                        [[1, 0]],
+                    ],
+                    dt=0.1,
+                ),
+            ),
+        ],
+    )
+    def test_combine_tf(self, tf_array, tf):
+        """Test combining transfer functions."""
+        tf_combined = ctrl.combine_tf(tf_array)
+        assert _tf_close_coeff(tf_combined, tf)
+
+    @pytest.mark.parametrize(
+        "tf_array, tf",
+        [
+            (
+                np.array(
+                    [
+                        [ctrl.TransferFunction([1], [1, 1])],
+                    ],
+                    dtype=object,
+                ),
+                ctrl.TransferFunction(
+                    [
+                        [[1]],
+                    ],
+                    [
+                        [[1, 1]],
+                    ],
+                ),
+            ),
+            (
+                np.array(
+                    [
+                        [ctrl.TransferFunction([1], [1, 1])],
+                        [ctrl.TransferFunction([2], [1, 0])],
+                    ],
+                    dtype=object,
+                ),
+                ctrl.TransferFunction(
+                    [
+                        [[1]],
+                        [[2]],
+                    ],
+                    [
+                        [[1, 1]],
+                        [[1, 0]],
+                    ],
+                ),
+            ),
+            (
+                np.array(
+                    [
+                        [ctrl.TransferFunction([1], [1, 1], dt=1)],
+                        [ctrl.TransferFunction([2], [1, 0], dt=1)],
+                    ],
+                    dtype=object,
+                ),
+                ctrl.TransferFunction(
+                    [
+                        [[1]],
+                        [[2]],
+                    ],
+                    [
+                        [[1, 1]],
+                        [[1, 0]],
+                    ],
+                    dt=1,
+                ),
+            ),
+            (
+                np.array(
+                    [
+                        [ctrl.TransferFunction([2], [1], dt=0.1)],
+                        [ctrl.TransferFunction([2], [1, 0], dt=0.1)],
+                    ],
+                    dtype=object,
+                ),
+                ctrl.TransferFunction(
+                    [
+                        [[2]],
+                        [[2]],
+                    ],
+                    [
+                        [[1]],
+                        [[1, 0]],
+                    ],
+                    dt=0.1,
+                ),
+            ),
+        ],
+    )
+    def test_split_tf(self, tf_array, tf):
+        """Test splitting transfer functions."""
+        tf_split = ctrl.split_tf(tf)
+        # Test entry-by-entry
+        for i in range(tf_split.shape[0]):
+            for j in range(tf_split.shape[1]):
+                assert _tf_close_coeff(
+                    tf_split[i, j],
+                    tf_array[i, j],
+                )
+        # Test combined
+        assert _tf_close_coeff(
+            ctrl.combine_tf(tf_split),
+            ctrl.combine_tf(tf_array),
+        )
+
+    @pytest.mark.parametrize(
+        "tf_array, exception",
+        [
+            # Wrong timesteps
+            (
+                [
+                    [ctrl.TransferFunction([1], [1, 1], 0.1)],
+                    [ctrl.TransferFunction([2], [1, 0], 0.2)],
+                ],
+                ValueError,
+            ),
+            (
+                [
+                    [ctrl.TransferFunction([1], [1, 1], 0.1)],
+                    [ctrl.TransferFunction([2], [1, 0], 0)],
+                ],
+                ValueError,
+            ),
+            # Too few dimensions
+            (
+                [
+                    ctrl.TransferFunction([1], [1, 1]),
+                    ctrl.TransferFunction([2], [1, 0]),
+                ],
+                ValueError,
+            ),
+            # Too many dimensions
+            (
+                [
+                    [[ctrl.TransferFunction([1], [1, 1], 0.1)]],
+                    [[ctrl.TransferFunction([2], [1, 0], 0)]],
+                ],
+                ValueError,
+            ),
+            # Incompatible dimensions
+            (
+                [
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [
+                                    [1],
+                                ]
+                            ],
+                            [
+                                [
+                                    [1, 1],
+                                ]
+                            ],
+                        ),
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1]],
+                                [[1], [3]],
+                            ],
+                            [
+                                [[1, 0], [1, 0]],
+                                [[1, 0], [1, 0]],
+                            ],
+                        ),
+                    ],
+                ],
+                ValueError,
+            ),
+            (
+                [
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1]],
+                                [[1], [3]],
+                            ],
+                            [
+                                [[1, 0], [1, 0]],
+                                [[1, 0], [1, 0]],
+                            ],
+                        ),
+                        ctrl.TransferFunction(
+                            [
+                                [
+                                    [1],
+                                ]
+                            ],
+                            [
+                                [
+                                    [1, 1],
+                                ]
+                            ],
+                        ),
+                    ],
+                ],
+                ValueError,
+            ),
+            (
+                [
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1]],
+                                [[1], [3]],
+                            ],
+                            [
+                                [[1, 0], [1, 0]],
+                                [[1, 0], [1, 0]],
+                            ],
+                        ),
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1]],
+                                [[1], [3]],
+                            ],
+                            [
+                                [[1, 0], [1, 0]],
+                                [[1, 0], [1, 0]],
+                            ],
+                        ),
+                    ],
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1], [1]],
+                                [[1], [3], [2]],
+                            ],
+                            [
+                                [[1, 0], [1, 0], [1, 0]],
+                                [[1, 0], [1, 0], [1, 0]],
+                            ],
+                        ),
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1]],
+                                [[1], [3]],
+                            ],
+                            [
+                                [[1, 0], [1, 0]],
+                                [[1, 0], [1, 0]],
+                            ],
+                        ),
+                    ],
+                ],
+                ValueError,
+            ),
+        ],
+    )
+    def test_error_combine_tf(self, tf_array, exception):
+        """Test error cases."""
+        with pytest.raises(exception):
+            ctrl.combine_tf(tf_array)
+
+
+def _tf_close_coeff(tf_a, tf_b, rtol=1e-5, atol=1e-8):
+    """Check if two transfer functions have close coefficients.
+
+    Parameters
+    ----------
+    tf_a : TransferFunction
+        First transfer function.
+    tf_b : TransferFunction
+        Second transfer function.
+    rtol : float
+        Relative tolerance for ``np.allclose``.
+    atol : float
+        Absolute tolerance for ``np.allclose``.
+
+    Returns
+    -------
+    bool
+        True if transfer function cofficients are all close.
+    """
+    # Check number of outputs and inputs
+    if tf_a.noutputs != tf_b.noutputs:
+        return False
+    if tf_a.ninputs != tf_b.ninputs:
+        return False
+    # Check timestep
+    if tf_a.dt != tf_b.dt:
+        return False
+    # Check coefficient arrays
+    for i in range(tf_a.noutputs):
+        for j in range(tf_a.ninputs):
+            if not np.allclose(
+                tf_a.num[i][j],
+                tf_b.num[i][j],
+                rtol=rtol,
+                atol=atol,
+            ):
+                return False
+            if not np.allclose(
+                tf_a.den[i][j],
+                tf_b.den[i][j],
+                rtol=rtol,
+                atol=atol,
+            ):
+                return False
+    return True
diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index bae0ec47b..c2a29ee2e 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -609,3 +609,49 @@ def test_frequency_response():
     assert mag_nosq_sq.shape == mag_default.shape
     assert phase_nosq_sq.shape == phase_default.shape
     assert omega_nosq_sq.shape == omega_default.shape
+
+
+def test_signal_labels():
+    # Create a system response for a SISO system
+    sys = ct.rss(4, 1, 1)
+    fresp = ct.frequency_response(sys)
+
+    # Make sure access via strings works
+    np.testing.assert_equal(
+        fresp.magnitude['y[0]'], fresp.magnitude[0])
+    np.testing.assert_equal(
+        fresp.phase['y[0]'], fresp.phase[0])
+
+    # Make sure errors are generated if key is unknown
+    with pytest.raises(ValueError, match="unknown signal name 'bad'"):
+        fresp.magnitude['bad']
+
+    # Create a system response for a MIMO system
+    sys = ct.rss(4, 2, 2)
+    fresp = ct.frequency_response(sys)
+
+    # Make sure access via strings works
+    np.testing.assert_equal(
+        fresp.magnitude['y[0]', 'u[1]'],
+        fresp.magnitude[0, 1])
+    np.testing.assert_equal(
+        fresp.phase['y[0]', 'u[1]'],
+        fresp.phase[0, 1])
+    np.testing.assert_equal(
+        fresp.response['y[0]', 'u[1]'],
+        fresp.response[0, 1])
+
+    # Make sure access via lists of strings works
+    np.testing.assert_equal(
+        fresp.response[['y[1]', 'y[0]'], 'u[0]'],
+        fresp.response[[1, 0], 0])
+
+    # Make sure errors are generated if key is unknown
+    with pytest.raises(ValueError, match="unknown signal name 'bad'"):
+        fresp.magnitude['bad']
+
+    with pytest.raises(ValueError, match="unknown signal name 'bad'"):
+        fresp.response[['y[1]', 'bad']]
+
+    with pytest.raises(ValueError, match=r"unknown signal name 'y\[0\]'"):
+        fresp.response['y[1]', 'y[0]']         # second index = input name
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index baaee03f6..0216235b5 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -2263,3 +2263,24 @@ def test_update_names():
 
     with pytest.raises(TypeError, match=".* takes 1 positional argument"):
         sys.update_names(5)
+
+
+def test_signal_indexing():
+    # Response with two outputs, no traces
+    resp = ct.initial_response(ct.rss(4, 2, 1, strictly_proper=True))
+    assert resp.outputs['y[0]'].shape == resp.outputs.shape[1:]
+    assert resp.outputs[0, 0].item() == 0
+
+    # Implicitly squeezed response
+    resp = ct.step_response(ct.rss(4, 1, 1, strictly_proper=True))
+    for key in ['y[0]', ('y[0]', 'u[0]')]:
+        with pytest.raises(IndexError, match=r"signal name\(s\) not valid"):
+            resp.outputs.__getitem__(key)
+
+    # Explicitly squeezed response
+    resp = ct.step_response(
+        ct.rss(4, 2, 1, strictly_proper=True), squeeze=True)
+    assert resp.outputs['y[0]'].shape == resp.outputs.shape[1:]
+    with pytest.raises(IndexError, match=r"signal name\(s\) not valid"):
+        resp.outputs['y[0]', 'u[0]']
+
diff --git a/control/tests/lti_test.py b/control/tests/lti_test.py
index 734bdb40b..3f001c17b 100644
--- a/control/tests/lti_test.py
+++ b/control/tests/lti_test.py
@@ -303,3 +303,50 @@ def test_squeeze_exceptions(self, fcn):
             sys([[0.1j, 1j], [1j, 10j]])
         with pytest.raises(ValueError, match="must be 1D"):
             evalfr(sys, [[0.1j, 1j], [1j, 10j]])
+
+
+@slycotonly
+@pytest.mark.parametrize(
+    "outdx, inpdx, key",
+    [('y[0]', 'u[1]', (0, 1)),
+     (['y[0]'], ['u[1]'], (0, 1)),
+     (slice(0, 1, 1), slice(1, 2, 1), (0, 1)),
+     (['y[0]', 'y[1]'], ['u[1]', 'u[2]'], ([0, 1], [1, 2])),
+     ([0, 'y[1]'], ['u[1]', 2], ([0, 1], [1, 2])),
+     (slice(0, 2, 1), slice(1, 3, 1), ([0, 1], [1, 2])),
+     (['y[2]', 'y[1]'], ['u[2]', 'u[0]'], ([2, 1], [2, 0])),
+    ])
+@pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.frd])
+def test_subsys_indexing(fcn, outdx, inpdx, key):
+    # Construct the base system and subsystem
+    sys = ct.rss(4, 3, 3)
+    subsys = sys[key]
+
+    # Construct the system to be tested
+    match fcn:
+        case ct.frd:
+            omega = np.logspace(-1, 1)
+            sys = fcn(sys, omega)
+            subsys_chk = fcn(subsys, omega)
+        case _:
+            sys = fcn(sys)
+            subsys_chk = fcn(subsys)
+
+    # Construct the subsystem
+    subsys_fcn = sys[outdx, inpdx]
+
+    # Check to make sure everythng matches up
+    match fcn:
+        case ct.frd:
+            np.testing.assert_almost_equal(
+                subsys_fcn.response, subsys_chk.response)
+        case ct.ss:
+            np.testing.assert_almost_equal(subsys_fcn.A, subsys_chk.A)
+            np.testing.assert_almost_equal(subsys_fcn.B, subsys_chk.B)
+            np.testing.assert_almost_equal(subsys_fcn.C, subsys_chk.C)
+            np.testing.assert_almost_equal(subsys_fcn.D, subsys_chk.D)
+        case ct.tf:
+            omega = np.logspace(-1, 1)
+            np.testing.assert_almost_equal(
+                subsys_fcn.frequency_response(omega).response,
+                subsys_chk.frequency_response(omega).response)
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 2829d6988..cb200c4ab 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -473,18 +473,22 @@ def test_array_access_ss_failure(self):
         with pytest.raises(IOError):
             sys1[0]
 
-    @pytest.mark.parametrize("outdx, inpdx", 
-                             [(0, 1),
-                              (slice(0, 1, 1), 1),
-                              (0, slice(1, 2, 1)),
-                              (slice(0, 1, 1), slice(1, 2, 1)),
-                              (slice(None, None, -1), 1),
-                              (0, slice(None, None, -1)),
-                              (slice(None, 2, None), 1),
-                              (slice(None, None, 1), slice(None, None, 2)),
-                              (0, slice(1, 2, 1)),
-                              (slice(0, 1, 1), slice(1, 2, 1))])
-    def test_array_access_ss(self, outdx, inpdx):
+    @pytest.mark.parametrize(
+        "outdx, inpdx",
+        [(0, 1),
+         (slice(0, 1, 1), 1),
+         (0, slice(1, 2, 1)),
+         (slice(0, 1, 1), slice(1, 2, 1)),
+         (slice(None, None, -1), 1),
+         (0, slice(None, None, -1)),
+         (slice(None, 2, None), 1),
+         (slice(None, None, 1), slice(None, None, 2)),
+         (0, slice(1, 2, 1)),
+         (slice(0, 1, 1), slice(1, 2, 1)),
+         # ([0, 1], [0]),         # lists of indices
+         ])
+    @pytest.mark.parametrize("named", [False, True])
+    def test_array_access_ss(self, outdx, inpdx, named):
         sys1 = StateSpace(
             [[1., 2.], [3., 4.]],
             [[5., 6.], [7., 8.]],
@@ -492,20 +496,22 @@ def test_array_access_ss(self, outdx, inpdx):
             [[13., 14.], [15., 16.]], 1,
             inputs=['u0', 'u1'], outputs=['y0', 'y1'])
 
-        sys1_01 = sys1[outdx, inpdx]
-        
+        if named:
+            # Use names instead of numbers (and re-convert in statesp)
+            outnames = sys1.output_labels[outdx]
+            inpnames = sys1.input_labels[inpdx]
+            sys1_01 = sys1[outnames, inpnames]
+        else:
+            sys1_01 = sys1[outdx, inpdx]
+
         # Convert int to slice to ensure that numpy doesn't drop the dimension
         if isinstance(outdx, int): outdx = slice(outdx, outdx+1, 1)
         if isinstance(inpdx, int): inpdx = slice(inpdx, inpdx+1, 1)
-        
-        np.testing.assert_array_almost_equal(sys1_01.A,
-                                             sys1.A)
-        np.testing.assert_array_almost_equal(sys1_01.B,
-                                             sys1.B[:, inpdx])
-        np.testing.assert_array_almost_equal(sys1_01.C,
-                                             sys1.C[outdx, :])
-        np.testing.assert_array_almost_equal(sys1_01.D,
-                                             sys1.D[outdx, inpdx])
+
+        np.testing.assert_array_almost_equal(sys1_01.A, sys1.A)
+        np.testing.assert_array_almost_equal(sys1_01.B, sys1.B[:, inpdx])
+        np.testing.assert_array_almost_equal(sys1_01.C, sys1.C[outdx, :])
+        np.testing.assert_array_almost_equal(sys1_01.D, sys1.D[outdx, inpdx])
 
         assert sys1.dt == sys1_01.dt
         assert sys1_01.input_labels == sys1.input_labels[inpdx]
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index e2d93be0e..e5e24b990 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -537,7 +537,7 @@ def test_discrete_time_impulse(self, tsystem):
         sysdt = sys.sample(dt, 'impulse')
         np.testing.assert_array_almost_equal(impulse_response(sys, t)[1],
                                              impulse_response(sysdt, t)[1])
-        
+
     def test_discrete_time_impulse_input(self):
         # discrete time impulse input, Only one active input for each trace
         A = [[.5, 0.25],[.0, .5]]
@@ -1318,3 +1318,53 @@ def test_step_info_nonstep():
     assert step_info['Peak'] == 1
     assert step_info['PeakTime'] == 0
     assert isclose(step_info['SteadyStateValue'], 0.96)
+
+
+def test_signal_labels():
+    # Create a system response for a SISO system
+    sys = ct.rss(4, 1, 1)
+    response = ct.step_response(sys)
+
+    # Make sure access via strings works
+    np.testing.assert_equal(response.states['x[2]'], response.states[2])
+
+    # Make sure access via lists of strings works
+    np.testing.assert_equal(
+        response.states[['x[1]', 'x[2]']], response.states[[1, 2]])
+
+    # Make sure errors are generated if key is unknown
+    with pytest.raises(ValueError, match="unknown signal name 'bad'"):
+        response.inputs['bad']
+
+    with pytest.raises(ValueError, match="unknown signal name 'bad'"):
+        response.states[['x[1]', 'bad']]
+
+    # Create a system response for a MIMO system
+    sys = ct.rss(4, 2, 2)
+    response = ct.step_response(sys)
+
+    # Make sure access via strings works
+    np.testing.assert_equal(
+        response.outputs['y[0]', 'u[1]'],
+        response.outputs[0, 1])
+    np.testing.assert_equal(
+        response.states['x[2]', 'u[0]'], response.states[2, 0])
+
+    # Make sure access via lists of strings works
+    np.testing.assert_equal(
+        response.states[['x[1]', 'x[2]'], 'u[0]'],
+        response.states[[1, 2], 0])
+
+    np.testing.assert_equal(
+        response.outputs[['y[1]'], ['u[1]', 'u[0]']],
+        response.outputs[[1], [1, 0]])
+
+    # Make sure errors are generated if key is unknown
+    with pytest.raises(ValueError, match="unknown signal name 'bad'"):
+        response.inputs['bad']
+
+    with pytest.raises(ValueError, match="unknown signal name 'bad'"):
+        response.states[['x[1]', 'bad']]
+
+    with pytest.raises(ValueError, match=r"unknown signal name 'x\[2\]'"):
+        response.states['x[1]', 'x[2]']         # second index = input name
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index 14a11b669..d480cef6e 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -390,19 +390,20 @@ def test_pow(self):
         with pytest.raises(ValueError):
             TransferFunction.__pow__(sys1, 0.5)
 
-    def test_slice(self):
+    @pytest.mark.parametrize("named", [False, True])
+    def test_slice(self, named):
         sys = TransferFunction(
             [ [   [1],    [2],    [3]], [   [3],    [4],    [5]] ],
             [ [[1, 2], [1, 3], [1, 4]], [[1, 4], [1, 5], [1, 6]] ],
             inputs=['u0', 'u1', 'u2'], outputs=['y0', 'y1'], name='sys')
 
-        sys1 = sys[1:, 1:]
+        sys1 = sys[1:, 1:] if not named else sys['y1', ['u1', 'u2']]
         assert (sys1.ninputs, sys1.noutputs) == (2, 1)
         assert sys1.input_labels == ['u1', 'u2']
         assert sys1.output_labels == ['y1']
         assert sys1.name == 'sys$indexed'
 
-        sys2 = sys[:2, :2]
+        sys2 = sys[:2, :2] if not named else sys[['y0', 'y1'], ['u0', 'u1']]
         assert (sys2.ninputs, sys2.noutputs) == (2, 2)
         assert sys2.input_labels == ['u0', 'u1']
         assert sys2.output_labels == ['y0', 'y1']
@@ -411,7 +412,7 @@ def test_slice(self):
         sys = TransferFunction(
             [ [   [1],    [2],    [3]], [   [3],    [4],    [5]] ],
             [ [[1, 2], [1, 3], [1, 4]], [[1, 4], [1, 5], [1, 6]] ], 0.5)
-        sys1 = sys[1:, 1:]
+        sys1 = sys[1:, 1:] if not named else sys[['y[1]'], ['u[1]', 'u[2]']]
         assert (sys1.ninputs, sys1.noutputs) == (2, 1)
         assert sys1.dt == 0.5
         assert sys1.input_labels == ['u[1]', 'u[2]']