Repository: python-control/python-control

Issue:
Response Functions as Member Functions?
Hi,
are there some ideological reasons why the response functions are not available as member functions for the respective system classes like they are in sicpy.signal? I feel like this would make the library a bit nicer to use.
I am willing to implement this myself, I just want to ensure that it would actually have a change of getting merged before putting in the effort.
Just to be clear, I mean forwarding the calls or aliasing the methods, not refactoring the code.

Discussion:

Patch:
diff --git a/control/__init__.py b/control/__init__.py
index 40f3a783..1aaaa42e 100644
--- a/control/__init__.py
+++ b/control/__init__.py
@@ -118,5 +118,14 @@ try:
 except ImportError:
     __version__ = "dev"
 
+# patch the LTI class with function aliases for convenience functions
+# this needs to be done after the fact since the modules that contain the functions
+# all heavily depend on the LTI class
+LTI.to_ss = ss
+LTI.to_tf = tf
+LTI.bode_plot = bode_plot
+LTI.nyquist_plot = nyquist_plot
+LTI.nichols_plot = nichols_plot
+
 # Initialize default parameter values
 reset_defaults()
diff --git a/control/bdalg.py b/control/bdalg.py
index 2dbd5c8e..d907cd3c 100644
--- a/control/bdalg.py
+++ b/control/bdalg.py
@@ -600,8 +600,8 @@ def combine_tf(tf_array):
                         f"row {row_index}."
                     )
                 for j_in in range(col.ninputs):
-                    num_row.append(col.num_array[j_out, j_in])
-                    den_row.append(col.den_array[j_out, j_in])
+                    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):
@@ -657,8 +657,8 @@ def split_tf(transfer_function):
         for i_in in range(transfer_function.ninputs):
             row.append(
                 tf.TransferFunction(
-                    transfer_function.num_array[i_out, i_in],
-                    transfer_function.den_array[i_out, i_in],
+                    transfer_function.num[i_out][i_in],
+                    transfer_function.den[i_out][i_in],
                     dt=transfer_function.dt,
                 )
             )
diff --git a/control/flatsys/flatsys.py b/control/flatsys/flatsys.py
index 7d76b9d7..0101d126 100644
--- a/control/flatsys/flatsys.py
+++ b/control/flatsys/flatsys.py
@@ -1,5 +1,42 @@
 # flatsys.py - trajectory generation for differentially flat systems
 # RMM, 10 Nov 2012
+#
+# This file contains routines for computing trajectories for differentially
+# flat nonlinear systems.  It is (very) loosely based on the NTG software
+# package developed by Mark Milam and Kudah Mushambi, but rewritten from
+# scratch in python.
+#
+# Copyright (c) 2012 by California Institute of Technology
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+# 1. Redistributions of source code must retain the above copyright
+#    notice, this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright
+#    notice, this list of conditions and the following disclaimer in the
+#    documentation and/or other materials provided with the distribution.
+#
+# 3. Neither the name of the California Institute of Technology nor
+#    the names of its contributors may be used to endorse or promote
+#    products derived from this software without specific prior
+#    written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+# FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CALTECH
+# OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+# SUCH DAMAGE.
 
 import itertools
 import numpy as np
@@ -17,30 +54,8 @@ class FlatSystem(NonlinearIOSystem):
     """Base class for representing a differentially flat system.
 
     The FlatSystem class is used as a base class to describe differentially
-    flat systems for trajectory generation.  The output of the system does
-    not need to be the differentially flat output.  Flat systems are
-    usually created with the :func:`~control.flatsys.flatsys` factory
-    function.
-
-    Parameters
-    ----------
-    forward : callable
-        A function to compute the flat flag given the states and input.
-    reverse : callable
-        A function to compute the states and input given the flat flag.
-    dt : None, True or float, optional
-        System timebase.
-
-    Attributes
-    ----------
-    ninputs, noutputs, nstates : int
-        Number of input, output and state variables.
-    shape : tuple
-        2-tuple of I/O system dimension, (noutputs, ninputs).
-    input_labels, output_labels, state_labels : list of str
-        Names for the input, output, and state variables.
-    name : string, optional
-        System name.
+    flat systems for trajectory generation.  The output of the system does not
+    need to be the differentially flat output.
 
     Notes
     -----
@@ -219,9 +234,10 @@ def flatsys(*args, updfcn=None, outfcn=None, **kwargs):
         Description of the system states.  Same format as `inputs`.
 
     dt : None, True or float, optional
-        System timebase.  None (default) indicates continuous time, True
-        indicates discrete time with undefined sampling time, positive
-        number is discrete time with specified sampling time.
+        System timebase.  None (default) indicates continuous
+        time, True indicates discrete time with undefined sampling
+        time, positive number is discrete time with specified
+        sampling time.
 
     params : dict, optional
         Parameter values for the systems.  Passed to the evaluation
@@ -236,12 +252,6 @@ def flatsys(*args, updfcn=None, outfcn=None, **kwargs):
     sys: :class:`FlatSystem`
         Flat system.
 
-    Other Parameters
-    ----------------
-    input_prefix, output_prefix, state_prefix : string, optional
-        Set the prefix for input, output, and state signals.  Defaults =
-        'u', 'y', 'x'.
-
     """
     from .linflat import LinearFlatSystem
     from ..statesp import StateSpace
diff --git a/control/frdata.py b/control/frdata.py
index 64a1e822..ac032d3f 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -3,11 +3,11 @@
 # Author: M.M. (Rene) van Paassen (using xferfcn.py as basis)
 # Date: 02 Oct 12
 
-"""Frequency response data representation and functions.
-
-This module contains the FrequencyResponseData (FRD) class and also
-functions that operate on FRD data.
+"""
+Frequency response data representation and functions.
 
+This module contains the FRD class and also functions that operate on
+FRD data.
 """
 
 from collections.abc import Iterable
@@ -21,8 +21,8 @@ from scipy.interpolate import splev, splprep
 
 from . import config
 from .exception import pandas_check
-from .iosys import InputOutputSystem, NamedSignal, _extended_system_name, \
-    _process_iosys_keywords, _process_subsys_index, 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']
@@ -35,29 +35,38 @@ 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, or
+    class constructor, using the :func:`~control.frd` factory function or
     via the :func:`~control.frequency_response` function.
 
     Parameters
     ----------
-    response : 1D or 3D complex array_like
+    d : 1D or 3D complex array_like
         The frequency response at each frequency point.  If 1D, the system is
         assumed to be SISO.  If 3D, the system is MIMO, with the first
         dimension corresponding to the output index of the FRD, the second
         dimension corresponding to the input index, and the 3rd dimension
         corresponding to the frequency points in omega
-    omega : iterable of real frequencies
+    w : iterable of real frequencies
         List of frequency points for which data are available.
+    sysname : str or None
+        Name of the system that generated the data.
     smooth : bool, optional
         If ``True``, create an interpolation function that allows the
         frequency response to be computed at any frequency within the range of
         frequencies give in ``w``.  If ``False`` (default), frequency response
         can only be obtained at the frequencies specified in ``w``.
-    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).
+
+    Attributes
+    ----------
+    ninputs, noutputs : int
+        Number of input and output variables.
+    omega : 1D array
+        Frequency points of the response.
+    fresp : 3D array
+        Frequency response, indexed by output index, input index, and
+        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
@@ -70,46 +79,16 @@ class FrequencyResponseData(LTI):
         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'].
-    sysname : str or None
-        Name of the system that generated the data.
-
-    Attributes
-    ----------
-    fresp : 3D array
-        Frequency response, indexed by output index, input index, and
-        frequency point.
-    frequency : 1D array
-        Array of frequency points for which data are available.
-    ninputs, noutputs : int
-        Number of input and output signals.
-    shape : tuple
-        2-tuple of I/O system dimension, (noutputs, ninputs).
+    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 signals.
-    name : str
-        System name.  For data generated using
-        :func:`~control.frequency_response`, stores the name of the
-        system that created the data.
-    magnitude : array
-        Magnitude of the frequency response, indexed by frequency.
-    phase : array
-        Phase of the frequency response, indexed by frequency.
-
-    Other Parameters
-    ----------------
-    plot_type : str, optional
-        Set the type of plot to generate with ``plot()`` ('bode', 'nichols').
+        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.
-    plot_magnitude, plot_phase : bool, optional
-        If set to `False`, don't plot the magnitude or phase, respectively.
-    return_magphase : bool, optional
-        If True, then a frequency response data object will enumerate as a
-        tuple of the form (mag, phase, omega) where where ``mag`` is the
-        magnitude (absolute value, not dB or log10) of the system
-        frequency response, ``phase`` is the wrapped phase in radians of
-        the system frequency response, and ``omega`` is the (sorted)
-        frequencies at which the response was evaluated.
 
     See Also
     --------
@@ -169,26 +148,22 @@ class FrequencyResponseData(LTI):
     _epsw = 1e-8                #: Bound for exact frequency match
 
     def __init__(self, *args, **kwargs):
-        """FrequencyResponseData(d, w[, dt])
+        """Construct an FRD object.
 
-        Construct a frequency response data (FRD) object.
+        The default constructor is FRD(d, w), where w is an iterable of
+        frequency points, and d is the matching frequency data.
 
-        The default constructor is FrequencyResponseData(d, w), where w is
-        an iterable of frequency points, and d is the matching frequency
-        data.  If d is a single list, 1D array, or tuple, a SISO system
-        description is assumed. d can also be a 2D array, in which case a
-        MIMO response is created.  To call the copy constructor, call
-        FrequencyResponseData(sys), where sys is a FRD object.  The
-        timebase for the frequency response can be provided using an
-        optional third argument or the 'dt' keyword.
+        If d is a single list, 1D array, or tuple, a SISO system description
+        is assumed. d can also be
+
+        To call the copy constructor, call FRD(sys), where sys is a
+        FRD object.
 
-        To construct frequency response data for an existing LTI object,
-        other than an FRD, call FrequencyResponseData(sys, omega).  This
-        functionality can also be obtained using :func:`frequency_response`
-        (which has additional options available).
+        To construct frequency response data for an existing LTI
+        object, other than an FRD, call FRD(sys, omega).
 
-        See :class:`FrequencyResponseData` and :func:`frd` for more
-        information.
+        The timebase for the frequency response can be provided using an
+        optional third argument or the 'dt' keyword.
 
         """
         smooth = kwargs.pop('smooth', False)
@@ -207,12 +182,11 @@ class FrequencyResponseData(LTI):
 
         if len(args) == 2:
             if not isinstance(args[0], FRD) and isinstance(args[0], LTI):
-                # not an FRD, but still an LTI system, second argument
-                # should be the frequency range
+                # not an FRD, but still a system, second argument should be
+                # the frequency range
                 otherlti = args[0]
                 self.omega = sort(np.asarray(args[1], dtype=float))
-
-                # calculate frequency response at specified points
+                # calculate frequency response at my points
                 if otherlti.isctime():
                     s = 1j * self.omega
                     self.fresp = otherlti(s, squeeze=False)
@@ -221,14 +195,6 @@ class FrequencyResponseData(LTI):
                     self.fresp = otherlti(z, squeeze=False)
                 arg_dt = otherlti.dt
 
-                # Copy over signal and system names, if not specified
-                kwargs['inputs'] = kwargs.get('inputs', otherlti.input_labels)
-                kwargs['outputs'] = kwargs.get(
-                    'outputs', otherlti.output_labels)
-                if not otherlti._generic_name_check():
-                    kwargs['name'] = kwargs.get('name', _extended_system_name(
-                        otherlti.name, prefix_suffix_name='sampled'))
-
             else:
                 # The user provided a response and a freq vector
                 self.fresp = array(args[0], dtype=complex, ndmin=1)
@@ -253,10 +219,6 @@ class FrequencyResponseData(LTI):
             self.fresp = args[0].fresp
             arg_dt = args[0].dt
 
-            # Copy over signal and system names, if not specified
-            kwargs['inputs'] = kwargs.get('inputs', args[0].input_labels)
-            kwargs['outputs'] = kwargs.get('outputs', args[0].output_labels)
-
         else:
             raise ValueError(
                 "Needs 1 or 2 arguments; received %i." % len(args))
@@ -287,21 +249,15 @@ class FrequencyResponseData(LTI):
 
         # Process iosys keywords
         defaults = {
-            'inputs': self.fresp.shape[1] if not getattr(
-                self, 'input_index', None) else self.input_labels,
-            'outputs': self.fresp.shape[0] if not getattr(
-                self, 'output_index', None) else self.output_labels,
-            'name': getattr(self, 'name', None)}
+            'inputs': self.fresp.shape[1], 'outputs': self.fresp.shape[0]}
         if arg_dt is not None:
-            if isinstance(args[0], LTI):
-                arg_dt = common_timebase(args[0].dt, arg_dt)
-            kwargs['dt'] = arg_dt
+            defaults['dt'] = arg_dt             # choose compatible timebase
+        name, inputs, outputs, states, dt = _process_iosys_keywords(
+                kwargs, defaults, end=True)
 
         # Process signal names
-        name, inputs, outputs, states, dt = _process_iosys_keywords(
-                kwargs, defaults)
         InputOutputSystem.__init__(
-            self, name=name, inputs=inputs, outputs=outputs, dt=dt, **kwargs)
+            self, name=name, inputs=inputs, outputs=outputs, dt=dt)
 
         # create interpolation functions
         if smooth:
@@ -310,17 +266,17 @@ class FrequencyResponseData(LTI):
                 raise ValueError("can't smooth with only 1 frequency")
             degree = 3 if self.omega.size > 3 else self.omega.size - 1
 
-            self._ifunc = empty((self.fresp.shape[0], self.fresp.shape[1]),
+            self.ifunc = empty((self.fresp.shape[0], self.fresp.shape[1]),
                                dtype=tuple)
             for i in range(self.fresp.shape[0]):
                 for j in range(self.fresp.shape[1]):
-                    self._ifunc[i, j], u = splprep(
+                    self.ifunc[i, j], u = splprep(
                         u=self.omega, x=[real(self.fresp[i, j, :]),
                                          imag(self.fresp[i, j, :])],
                         w=1.0/(absolute(self.fresp[i, j, :]) + 0.001),
                         s=0.0, k=degree)
         else:
-            self._ifunc = None
+            self.ifunc = None
 
     #
     # Frequency response properties
@@ -423,7 +379,7 @@ class FrequencyResponseData(LTI):
         """
         return "FrequencyResponseData({d}, {w}{smooth})".format(
             d=repr(self.fresp), w=repr(self.omega),
-            smooth=(self._ifunc and ", smooth=True") or "")
+            smooth=(self.ifunc and ", smooth=True") or "")
 
     def __neg__(self):
         """Negate a transfer function."""
@@ -482,7 +438,7 @@ class FrequencyResponseData(LTI):
         # Convert the second argument to a transfer function.
         if isinstance(other, (int, float, complex, np.number)):
             return FRD(self.fresp * other, self.omega,
-                       smooth=(self._ifunc is not None))
+                       smooth=(self.ifunc is not None))
         else:
             other = _convert_to_frd(other, omega=self.omega)
 
@@ -500,8 +456,8 @@ class FrequencyResponseData(LTI):
         for i in range(len(self.omega)):
             fresp[:, :, i] = self.fresp[:, :, i] @ other.fresp[:, :, i]
         return FRD(fresp, self.omega,
-                   smooth=(self._ifunc is not None) and
-                          (other._ifunc is not None))
+                   smooth=(self.ifunc is not None) and
+                          (other.ifunc is not None))
 
     def __rmul__(self, other):
         """Right Multiply two LTI objects (serial connection)."""
@@ -509,7 +465,7 @@ class FrequencyResponseData(LTI):
         # Convert the second argument to an frd function.
         if isinstance(other, (int, float, complex, np.number)):
             return FRD(self.fresp * other, self.omega,
-                       smooth=(self._ifunc is not None))
+                       smooth=(self.ifunc is not None))
         else:
             other = _convert_to_frd(other, omega=self.omega)
 
@@ -528,8 +484,8 @@ class FrequencyResponseData(LTI):
         for i in range(len(self.omega)):
             fresp[:, :, i] = other.fresp[:, :, i] @ self.fresp[:, :, i]
         return FRD(fresp, self.omega,
-                   smooth=(self._ifunc is not None) and
-                          (other._ifunc is not None))
+                   smooth=(self.ifunc is not None) and
+                          (other.ifunc is not None))
 
     # TODO: Division of MIMO transfer function objects is not written yet.
     def __truediv__(self, other):
@@ -537,7 +493,7 @@ class FrequencyResponseData(LTI):
 
         if isinstance(other, (int, float, complex, np.number)):
             return FRD(self.fresp * (1/other), self.omega,
-                       smooth=(self._ifunc is not None))
+                       smooth=(self.ifunc is not None))
         else:
             other = _convert_to_frd(other, omega=self.omega)
 
@@ -548,15 +504,15 @@ class FrequencyResponseData(LTI):
                 "systems.")
 
         return FRD(self.fresp/other.fresp, self.omega,
-                   smooth=(self._ifunc is not None) and
-                          (other._ifunc is not None))
+                   smooth=(self.ifunc is not None) and
+                          (other.ifunc is not None))
 
     # TODO: Division of MIMO transfer function objects is not written yet.
     def __rtruediv__(self, other):
         """Right divide two LTI objects."""
         if isinstance(other, (int, float, complex, np.number)):
             return FRD(other / self.fresp, self.omega,
-                       smooth=(self._ifunc is not None))
+                       smooth=(self.ifunc is not None))
         else:
             other = _convert_to_frd(other, omega=self.omega)
 
@@ -573,7 +529,7 @@ class FrequencyResponseData(LTI):
             raise ValueError("Exponent must be an integer")
         if other == 0:
             return FRD(ones(self.fresp.shape), self.omega,
-                       smooth=(self._ifunc is not None))  # unity
+                       smooth=(self.ifunc is not None))  # unity
         if other > 0:
             return self * (self**(other-1))
         if other < 0:
@@ -626,7 +582,7 @@ class FrequencyResponseData(LTI):
         if any(omega_array.imag > 0):
             raise ValueError("FRD.eval can only accept real-valued omega")
 
-        if self._ifunc is None:
+        if self.ifunc is None:
             elements = np.isin(self.omega, omega)  # binary array
             if sum(elements) < len(omega_array):
                 raise ValueError(
@@ -640,7 +596,7 @@ class FrequencyResponseData(LTI):
             for i in range(self.noutputs):
                 for j in range(self.ninputs):
                     for k, w in enumerate(omega_array):
-                        frraw = splev(w, self._ifunc[i, j], der=0)
+                        frraw = splev(w, self.ifunc[i, j], der=0)
                         out[i, j, k] = frraw[0] + 1.0j * frraw[1]
 
         return _process_frequency_response(self, omega, out, squeeze=squeeze)
@@ -795,7 +751,7 @@ class FrequencyResponseData(LTI):
         resfresp = (myfresp @ linalg.inv(I_AB))
         fresp = np.moveaxis(resfresp, 0, 2)
 
-        return FRD(fresp, other.omega, smooth=(self._ifunc is not None))
+        return FRD(fresp, other.omega, smooth=(self.ifunc is not None))
 
     # Plotting interface
     def plot(self, plot_type=None, *args, **kwargs):
@@ -945,8 +901,6 @@ def frd(*args, **kwargs):
 
     Parameters
     ----------
-    sys : LTI (StateSpace or TransferFunction)
-        A linear system that will be evaluated for frequency response data.
     response : array_like or LTI system
         Complex vector with the system response or an LTI system that can
         be used to copmute the frequency response at a list of frequencies.
@@ -963,7 +917,7 @@ def frd(*args, **kwargs):
 
     Returns
     -------
-    sys : FrequencyResponseData
+    sys : :class:`FrequencyResponseData`
         New frequency response data system.
 
     Other Parameters
@@ -972,8 +926,6 @@ def frd(*args, **kwargs):
         List of strings that name the individual signals of the transformed
         system.  If not given, the inputs and outputs are the same as the
         original system.
-    input_prefix, output_prefix : string, optional
-        Set the prefix for input and output signals.  Defaults = 'u', 'y'.
     name : string, optional
         System name. If unspecified, a generic name <sys[id]> is generated
         with a unique integer id.
diff --git a/control/freqplot.py b/control/freqplot.py
index 8d3e6468..54442529 100644
--- a/control/freqplot.py
+++ b/control/freqplot.py
@@ -1305,7 +1305,7 @@ def nyquist_response(
                 "Nyquist plot currently only supports SISO systems.")
 
         # Figure out the frequency range
-        if isinstance(sys, FrequencyResponseData) and sys._ifunc is None \
+        if isinstance(sys, FrequencyResponseData) and sys.ifunc is None \
            and not omega_range_given:
             omega_sys = sys.omega               # use system frequencies
         else:
diff --git a/control/iosys.py b/control/iosys.py
index 2c1f9cea..dd1566eb 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -145,14 +145,20 @@ class InputOutputSystem(object):
     Attributes
     ----------
     ninputs, noutputs, nstates : int
-        Number of input, output, and state variables.
+        Number of input, output and state variables
     input_index, output_index, state_index : dict
-        Dictionary of signal names for the inputs, outputs, and states and
-        the index of the corresponding array.
-    input_labels, output_labels, state_labels : list of str
-        List of signal names for inputs, outputs, and states.
-    shape : tuple
-        2-tuple of I/O system dimension, (noutputs, ninputs).
+        Dictionary of signal names for the inputs, outputs and states and the
+        index of the corresponding array
+    dt : None, True or float
+        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).
+    params : dict, optional
+        Parameter values for the systems.  Passed to the evaluation functions
+        for the system as default values, overriding internal defaults.
+    name : string, optional
+        System name (used for specifying signals)
 
     Other Parameters
     ----------------
@@ -188,7 +194,7 @@ class InputOutputSystem(object):
             raise TypeError("unrecognized keywords: ", str(kwargs))
 
     # Keep track of the keywords that we recognize
-    _kwargs_list = [
+    kwargs_list = [
         'name', 'inputs', 'outputs', 'states', 'input_prefix',
         'output_prefix', 'state_prefix', 'dt']
 
@@ -287,8 +293,13 @@ class InputOutputSystem(object):
         """copy the signal and system name of sys. Name is given as a keyword
         in case a specific name (e.g. append 'linearized') is desired. """
         # Figure out the system name and assign it
-        self.name = _extended_system_name(
-            sys.name, prefix, suffix, prefix_suffix_name)
+        if prefix == "" and prefix_suffix_name is not None:
+            prefix = config.defaults[
+                'iosys.' + prefix_suffix_name + '_system_name_prefix']
+        if suffix == "" and prefix_suffix_name is not None:
+            suffix = config.defaults[
+                'iosys.' + prefix_suffix_name + '_system_name_suffix']
+        self.name = prefix + sys.name + suffix
 
         # Name the inputs, outputs, and states
         self.input_index = sys.input_index.copy()
@@ -421,11 +432,6 @@ class InputOutputSystem(object):
         lambda self: list(self.state_index.keys()),     # getter
         set_states)                                     # setter
 
-    @property
-    def shape(self):
-        """2-tuple of I/O system dimension, (noutputs, ninputs)."""
-        return (self.noutputs, self.ninputs)
-
     # TODO: add dict as a means to selective change names?  [GH #1019]
     def update_names(self, **kwargs):
         """update_names([name, inputs, outputs, states])
@@ -1051,14 +1057,3 @@ def _process_subsys_index(idx, sys_labels, slice_to_list=False):
         idx = range(len(sys_labels))[idx]
 
     return idx, labels
-
-
-# Create an extended system name
-def _extended_system_name(name, prefix="", suffix="", prefix_suffix_name=None):
-    if prefix == "" and prefix_suffix_name is not None:
-        prefix = config.defaults[
-            'iosys.' + prefix_suffix_name + '_system_name_prefix']
-    if suffix == "" and prefix_suffix_name is not None:
-        suffix = config.defaults[
-            'iosys.' + prefix_suffix_name + '_system_name_suffix']
-    return prefix + name + suffix
diff --git a/control/lti.py b/control/lti.py
index 4b354f3f..73240b70 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -11,6 +11,8 @@ from numpy import real, angle, abs
 from warnings import warn
 from . import config
 from .iosys import InputOutputSystem
+import control
+from typing import Callable
 
 __all__ = ['poles', 'zeros', 'damp', 'evalfr', 'frequency_response',
            'freqresp', 'dcgain', 'bandwidth', 'LTI']
@@ -76,9 +78,10 @@ class LTI(InputOutputSystem):
 
         Parameters
         ----------
-        omega : float or 1D array_like
+        omega : float or 1D array_like, optional
             A list, tuple, array, or scalar value of frequencies in
-            radians/sec at which the system will be evaluated.
+            radians/sec at which the system will be evaluated.  If None (default),
+            a set of frequencies is computed based on the system dynamics.
         squeeze : bool, optional
             If squeeze=True, remove single-dimensional entries from the shape
             of the output even if the system is not SISO. If squeeze=False,
@@ -107,6 +110,11 @@ class LTI(InputOutputSystem):
         """
         from .frdata import FrequencyResponseData
 
+        if omega is None:
+            # Use default frequency range
+            from .freqplot import _default_frequency_range
+            omega = _default_frequency_range(self)
+
         omega = np.sort(np.array(omega, ndmin=1))
         if self.isdtime(strict=True):
             # Convert the frequency to discrete time
@@ -202,6 +210,34 @@ class LTI(InputOutputSystem):
         # importing here prevents circular dependancy
         from control.passivity import ispassive
         return ispassive(self)
+    
+    def feedback(self, other=1, sign=-1):
+        raise NotImplementedError(f"feedback not implemented for base {self.__class__.__name__} objects")
+
+    # convenience aliases
+    # most function are only forward declaraed and patched in the __init__.py to avoid circular imports
+
+    # conversions
+    #: Convert to :class:`StateSpace` representation; see :func:`ss`
+    to_ss: Callable
+    #: Convert to :class:`TransferFunction` representation; see :func:`tf`
+    to_tf: Callable
+
+    # freq domain plotting
+    #: Bode plot; see :func:`bode_plot`
+    bode_plot: Callable
+    #: Nyquist plot; see :func:`nyquist_plot`
+    nyquist_plot: Callable
+    #: Nichols plot; see :func:`nichols_plot`
+    nichols_plot: Callable
+
+    # time domain simulation
+    #: Forced response; see :func:`forced_response`
+    forced_response = control.timeresp.forced_response
+    #: Impulse response; see :func:`impulse_response`
+    impulse_response = control.timeresp.impulse_response
+    #: Step response; see :func:`step_response`
+    step_response = control.timeresp.step_response
 
 
 def poles(sys):
@@ -458,7 +494,7 @@ def frequency_response(
     Examples
     --------
     >>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]])
-    >>> mag, phase, omega = ct.frequency_response(G, [0.1, 1., 10.])
+    >>> mag, phase, omega = ct.freqresp(G, [0.1, 1., 10.])
 
     .. todo::
         Add example with MIMO system
@@ -490,8 +526,8 @@ def frequency_response(
 
     responses = []
     for sys_ in syslist:
-        if isinstance(sys_, FrequencyResponseData) and sys_._ifunc is None \
-           and not omega_range_given:
+        if isinstance(sys_, FrequencyResponseData) and sys_.ifunc is None and \
+           not omega_range_given:
             omega_sys = sys_.omega              # use system properties
         else:
             omega_sys = omega_syslist.copy()    # use common omega vector
diff --git a/control/nlsys.py b/control/nlsys.py
index beb2566e..62e4bf78 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -39,11 +39,10 @@ __all__ = ['NonlinearIOSystem', 'InterconnectedSystem', 'nlsys',
 class NonlinearIOSystem(InputOutputSystem):
     """Nonlinear I/O system.
 
-    Creates an :class:`~control.InputOutputSystem` for a nonlinear system
-    by specifying a state update function and an output function.  The new
-    system can be a continuous or discrete time system. Nonlinear I/O
-    systems are usually created with the :func:`~control.nlsys` factory
-    function.
+    Creates an :class:`~control.InputOutputSystem` for a nonlinear system by
+    specifying a state update function and an output function.  The new system
+    can be a continuous or discrete time system (Note: discrete-time systems
+    are not yet supported by most functions.)
 
     Parameters
     ----------
@@ -64,13 +63,20 @@ class NonlinearIOSystem(InputOutputSystem):
 
         where the arguments are the same as for `upfcn`.
 
-    inputs, outputs, states : int, list of str or None, optional
-        Description of the system inputs, outputs, and states.  See
-        :func:`control.nlsys` for more details.
+    inputs : int, list of str or None, optional
+        Description of the system inputs.  This can be given as an integer
+        count or as a list of strings that name the individual signals.
+        If an integer count is specified, the names of the signal will be
+        of the form 's[i]' (where 's' is one of 'u', 'y', or 'x').  If
+        this parameter is not given or given as `None`, the relevant
+        quantity will be determined when possible based on other
+        information provided to functions using the system.
 
-    params : dict, optional
-        Parameter values for the systems.  Passed to the evaluation functions
-        for the system as default values, overriding internal defaults.
+    outputs : int, list of str or None, optional
+        Description of the system outputs.  Same format as `inputs`.
+
+    states : int, list of str, or None, optional
+        Description of the system states.  Same format as `inputs`.
 
     dt : timebase, optional
         The timebase for the system, used to specify whether the system is
@@ -82,16 +88,13 @@ class NonlinearIOSystem(InputOutputSystem):
         * dt = True: discrete time with unspecified sampling period
         * dt = None: no timebase specified
 
-    Attributes
-    ----------
-    ninputs, noutputs, nstates : int
-        Number of input, output and state variables.
-    shape : tuple
-        2-tuple of I/O system dimension, (noutputs, ninputs).
-    input_labels, output_labels, state_labels : list of str
-        Names for the input, output, and state variables.
     name : string, optional
-        System name.
+        System name (used for specifying signals). If unspecified, a
+        generic name <sys[id]> is generated with a unique integer id.
+
+    params : dict, optional
+        Parameter values for the system.  Passed to the evaluation functions
+        for the system as default values, overriding internal defaults.
 
     See Also
     --------
@@ -710,11 +713,6 @@ class InterconnectedSystem(NonlinearIOSystem):
         if outputs is None and outlist is not None:
             outputs = len(outlist)
 
-        if params is None:
-            params = {}
-            for sys in self.syslist:
-                params = params | sys.params
-
         # Create updfcn and outfcn
         def updfcn(t, x, u, params):
             self._update_params(params)
@@ -1222,7 +1220,8 @@ class InterconnectedSystem(NonlinearIOSystem):
         return dropped_inputs, dropped_outputs
 
 
-def nlsys(updfcn, outfcn=None, **kwargs):
+def nlsys(
+        updfcn, outfcn=None, inputs=None, outputs=None, states=None, **kwargs):
     """Create a nonlinear input/output system.
 
     Creates an :class:`~control.InputOutputSystem` for a nonlinear system by
@@ -1231,7 +1230,7 @@ def nlsys(updfcn, outfcn=None, **kwargs):
 
     Parameters
     ----------
-    updfcn : callable (or StateSpace)
+    updfcn : callable
         Function returning the state update function
 
             `updfcn(t, x, u, params) -> array`
@@ -1241,10 +1240,6 @@ def nlsys(updfcn, outfcn=None, **kwargs):
         time, and `params` is a dict containing the values of parameters
         used by the function.
 
-        If a :class:`StateSpace` system is passed as the update function,
-        then a nonlinear I/O system is created that implements the linear
-        dynamics of the state space system.
-
     outfcn : callable
         Function returning the output at the given state
 
@@ -1290,12 +1285,6 @@ def nlsys(updfcn, outfcn=None, **kwargs):
     sys : :class:`NonlinearIOSystem`
         Nonlinear input/output system.
 
-    Other Parameters
-    ----------------
-    input_prefix, output_prefix, state_prefix : string, optional
-        Set the prefix for input, output, and state signals.  Defaults =
-        'u', 'y', 'x'.
-
     See Also
     --------
     ss, tf
@@ -1319,34 +1308,9 @@ def nlsys(updfcn, outfcn=None, **kwargs):
     >>> timepts = np.linspace(0, 10)
     >>> response = ct.input_output_response(
     ...     kincar, timepts, [10, 0.05 * np.sin(timepts)])
-
     """
-    from .statesp import StateSpace
-    from .iosys import _extended_system_name
-
-    if isinstance(updfcn, StateSpace):
-        sys_ss = updfcn
-        kwargs['inputs'] = kwargs.get('inputs', sys_ss.input_labels)
-        kwargs['outputs'] = kwargs.get('outputs', sys_ss.output_labels)
-        kwargs['states'] = kwargs.get('states', sys_ss.state_labels)
-        kwargs['name'] = kwargs.get('name', _extended_system_name(
-            sys_ss.name, prefix_suffix_name='converted'))
-
-        sys_nl = NonlinearIOSystem(
-            lambda t, x, u, params:
-                sys_ss.A @ np.atleast_1d(x) + sys_ss.B @ np.atleast_1d(u),
-            lambda t, x, u, params:
-                sys_ss.C @ np.atleast_1d(x) + sys_ss.D @ np.atleast_1d(u),
-            **kwargs)
-
-        if sys_nl.nstates != sys_ss.nstates or sys_nl.shape != sys_ss.shape:
-            raise ValueError(
-                "new input, output, or state specification "
-                "doesn't match system size")
-
-        return sys_nl
-    else:
-        return NonlinearIOSystem(updfcn, outfcn, **kwargs)
+    return NonlinearIOSystem(
+        updfcn, outfcn, inputs=inputs, outputs=outputs, states=states, **kwargs)
 
 
 def input_output_response(
@@ -2276,8 +2240,7 @@ def interconnect(
 
     params : dict, optional
         Parameter values for the systems.  Passed to the evaluation functions
-        for the system as default values, overriding internal defaults.  If
-        not specified, defaults to parameters from subsystems.
+        for the system as default values, overriding internal defaults.
 
     dt : timebase, optional
         The timebase for the system, used to specify whether the system is
diff --git a/control/optimal.py b/control/optimal.py
index 77cfd370..0eb49c82 100644
--- a/control/optimal.py
+++ b/control/optimal.py
@@ -1203,7 +1203,7 @@ def create_mpc_iosystem(
 
     # Grab the keyword arguments known by this function
     iosys_kwargs = {}
-    for kw in InputOutputSystem._kwargs_list:
+    for kw in InputOutputSystem.kwargs_list:
         if kw in kwargs:
             iosys_kwargs[kw] = kwargs.pop(kw)
 
diff --git a/control/statefbk.py b/control/statefbk.py
index c5c5a803..1dd0be32 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -586,8 +586,7 @@ def create_statefbk_iosystem(
         controller_type=None, xd_labels=None, ud_labels=None, ref_labels=None,
         feedfwd_pattern='trajgen', gainsched_indices=None,
         gainsched_method='linear', control_indices=None, state_indices=None,
-        name=None, inputs=None, outputs=None, states=None, params=None,
-        **kwargs):
+        name=None, inputs=None, outputs=None, states=None, **kwargs):
     r"""Create an I/O system using a (full) state feedback controller.
 
     This function creates an input/output system that implements a
@@ -752,10 +751,6 @@ def create_statefbk_iosystem(
         System name. If unspecified, a generic name <sys[id]> is generated
         with a unique integer id.
 
-    params : dict, optional
-        System parameter values.  By default, these will be copied from
-        `sys` and `ctrl`, but can be overridden with this keyword.
-
     Examples
     --------
     >>> import control as ct
@@ -778,8 +773,7 @@ def create_statefbk_iosystem(
     if not isinstance(sys, NonlinearIOSystem):
         raise ControlArgument("Input system must be I/O system")
 
-    # Process keywords
-    params = sys.params if params is None else params
+    # Process (legacy) keywords
     controller_type = _process_legacy_keyword(
         kwargs, 'type', 'controller_type', controller_type)
     if kwargs:
@@ -976,10 +970,10 @@ def create_statefbk_iosystem(
 
             return u
 
-        ctrl_params = {} if gainsched else {'K': K}
+        params = {} if gainsched else {'K': K}
         ctrl = NonlinearIOSystem(
             _control_update, _control_output, name=name, inputs=inputs,
-            outputs=outputs, states=states, params=ctrl_params)
+            outputs=outputs, states=states, params=params)
 
     elif controller_type == 'iosystem' and feedfwd_pattern == 'trajgen':
         # Use the passed system to compute feedback compensation
@@ -1067,8 +1061,7 @@ def create_statefbk_iosystem(
         [sys, ctrl] if estimator == sys else [sys, ctrl, estimator],
         name=sys.name + "_" + ctrl.name, add_unused=True,
         inplist=inplist, inputs=input_labels,
-        outlist=outlist, outputs=output_labels,
-        params= ctrl.params | params
+        outlist=outlist, outputs=output_labels
     )
     return ctrl, closed
 
diff --git a/control/statesp.py b/control/statesp.py
index 3f53777e..bd312b47 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -1,18 +1,52 @@
-# statesp.py - state space class and related functions
-#
-# Original author: Richard M. Murray
-# Creation date: 24 May 2009
-# Pre-2014 revisions: Kevin K. Chen, Dec 10
-# Use `git shortlog -n -s statesp.py` for full list of contributors
+"""statesp.py
 
-"""State space representation and functions.
+State space representation and functions.
 
-This module contains the StateSpace class, which is used to represent
-linear systems in state space.  This is the primary representation for the
+This file contains the StateSpace class, which is used to represent linear
+systems in state space.  This is the primary representation for the
 python-control library.
 
 """
 
+"""Copyright (c) 2010 by California Institute of Technology
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the California Institute of Technology nor
+   the names of its contributors may be used to endorse or promote
+   products derived from this software without specific prior
+   written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CALTECH
+OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
+
+Author: Richard M. Murray
+Date: 24 May 09
+Revised: Kevin K. Chen, Dec 10
+
+$Id$
+"""
+
 import math
 from collections.abc import Iterable
 from copy import deepcopy
@@ -36,6 +70,7 @@ from .iosys import InputOutputSystem, NamedSignal, _process_dt_keyword, \
     common_timebase, isdtime, issiso
 from .lti import LTI, _process_frequency_response
 from .nlsys import InterconnectedSystem, NonlinearIOSystem
+import control
 
 try:
     from slycot import ab13dd
@@ -67,30 +102,30 @@ class StateSpace(NonlinearIOSystem, LTI):
           dx/dt &= A x + B u \\
               y &= C x + D u
 
-    where `u` is the input, `y` is the output, and `x` is the state.  State
-    space systems are usually created with the :func:`~control.ss` factory
-    function.
+    where `u` is the input, `y` is the output, and `x` is the state.
 
     Parameters
     ----------
-    A, B, C, D : array_like
+    A, B, C, D: array_like
         System matrices of the appropriate dimensions.
     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).
+        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).
 
     Attributes
     ----------
     ninputs, noutputs, nstates : int
         Number of input, output and state variables.
-    shape : tuple
-        2-tuple of I/O system dimension, (noutputs, ninputs).
-    input_labels, output_labels, state_labels : list of str
-        Names for the input, output, and state variables.
-    name : string, optional
-        System name.
+    A, B, C, D : 2D arrays
+        System matrices defining the input/output dynamics.
+    dt : None, True or float
+        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).
 
     Notes
     -----
@@ -130,12 +165,13 @@ class StateSpace(NonlinearIOSystem, LTI):
     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 `control.config.defaults['statesp.latex_num_format']`
-    and `control.config.defaults['statesp.latex_repr_type']`.  The LaTeX
-    output is tailored for MathJax, as used in Jupyter, and may look odd
-    when typeset by non-MathJax LaTeX systems.
+    StateSpace instances have support for IPython LaTeX output,
+    intended for pretty-printing in Jupyter notebooks.  The LaTeX
+    output can be configured using
+    `control.config.defaults['statesp.latex_num_format']` and
+    `control.config.defaults['statesp.latex_repr_type']`.  The LaTeX output is
+    tailored for MathJax, as used in Jupyter, and may look odd when
+    typeset by non-MathJax LaTeX systems.
 
     `control.config.defaults['statesp.latex_num_format']` is a format string
     fragment, specifically the part of the format string after `'{:'`
@@ -159,7 +195,12 @@ class StateSpace(NonlinearIOSystem, LTI):
         True for unspecified sampling time).  To call the copy constructor,
         call StateSpace(sys), where sys is a StateSpace object.
 
-        See :class:`StateSpace` and :func:`ss` for more information.
+        The `remove_useless_states` keyword can be used to scan the A, B, and
+        C matrices for rows or columns of zeros.  If the zeros are such that a
+        particular state has no effect on the input-output dynamics, then that
+        state is removed from the A, B, and C matrices.  If not specified, the
+        value is read from `config.defaults['statesp.remove_useless_states']`
+        (default = False).
 
         """
         #
@@ -1427,6 +1468,9 @@ class StateSpace(NonlinearIOSystem, LTI):
                 raise ValueError("len(u) must be equal to number of inputs")
             return (self.C @ x).reshape((-1,)) \
                 + (self.D @ u).reshape((-1,))  # return as row vector
+        
+    # convenience aliase, import needs to go over the submodule to avoid circular imports
+    initial_response = control.timeresp.initial_response
 
 
 class LinearICSystem(InterconnectedSystem, StateSpace):
@@ -1501,7 +1545,7 @@ def ss(*args, **kwargs):
 
     Create a state space system.
 
-    The function accepts either 1, 4 or 5 positional parameters:
+    The function accepts either 1, 2, 4 or 5 parameters:
 
     ``ss(sys)``
         Convert a linear system into space system form. Always creates a
@@ -1525,11 +1569,11 @@ def ss(*args, **kwargs):
             x[k+1] &= A x[k] + B u[k] \\
               y[k] &= C x[k] + D u[k]
 
-        The matrices can be given as 2D array-like data types.  For SISO
-        systems, `B` and `C` can be given as 1D arrays and D can be given
-        as a scalar.
+        The matrices can be given as *array like* data types or strings.
+        Everything that the constructor of :class:`numpy.matrix` accepts is
+        permissible here too.
 
-    ``ss(*args, inputs=['u1', ..., 'up'], outputs=['y1', ..., 'yq'], states=['x1', ..., 'xn'])``
+    ``ss(args, inputs=['u1', ..., 'up'], outputs=['y1', ..., 'yq'], states=['x1', ..., 'xn'])``
         Create a system with named input, output, and state signals.
 
     Parameters
@@ -1544,33 +1588,23 @@ def ss(*args, **kwargs):
         time, positive number is discrete time with specified
         sampling time, None indicates unspecified timebase (either
         continuous or discrete time).
-    remove_useless_states : bool, optional
-        If `True`, remove states that have no effect on the input/output
-        dynamics.  If not specified, the value is read from
-        `config.defaults['statesp.remove_useless_states']` (default = False).
-    method : str, optional
-        Set the method used for converting a transfer function to a state
-        space system.  Current methods are 'slycot' and 'scipy'.  If set to
-        None (default), try 'slycot' first and then 'scipy' (SISO only).
-
-    Returns
-    -------
-    out : StateSpace
-        Linear input/output system.
-
-    Other Parameters
-    ----------------
     inputs, outputs, states : str, or list of str, optional
         List of strings that name the individual signals.  If this parameter
         is not given or given as `None`, the signal names will be of the
         form `s[i]` (where `s` is one of `u`, `y`, or `x`). See
         :class:`InputOutputSystem` for more information.
-    input_prefix, output_prefix, state_prefix : string, optional
-        Set the prefix for input, output, and state signals.  Defaults =
-        'u', 'y', 'x'.
     name : string, optional
         System name (used for specifying signals). If unspecified, a generic
         name <sys[id]> is generated with a unique integer id.
+    method : str, optional
+        Set the method used for computing the result.  Current methods are
+        'slycot' and 'scipy'.  If set to None (default), try 'slycot' first
+        and then 'scipy' (SISO only).
+
+    Returns
+    -------
+    out: :class:`StateSpace`
+        Linear input/output system.
 
     Raises
     ------
@@ -1579,7 +1613,7 @@ def ss(*args, **kwargs):
 
     See Also
     --------
-    tf, ss2tf, tf2ss, zpk
+    tf, ss2tf, tf2ss
 
     Notes
     -----
@@ -2261,7 +2295,7 @@ def _convert_to_statespace(sys, use_prefix_suffix=False, method=None):
                 D = empty((sys.noutputs, sys.ninputs), dtype=float)
                 for i, j in itertools.product(range(sys.noutputs),
                                               range(sys.ninputs)):
-                    D[i, j] = sys.num_array[i, j][0] / sys.den_array[i, j][0]
+                    D[i, j] = sys.num[i][j][0] / sys.den[i][j][0]
                 newsys = StateSpace([], [], [], D, sys.dt)
             else:
                 if not issiso(sys):
diff --git a/control/timeresp.py b/control/timeresp.py
index 67641d23..072db60d 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -718,10 +718,8 @@ class TimeResponseData:
     def to_pandas(self):
         """Convert response data to pandas data frame.
 
-        Creates a pandas data frame using the input, output, and state labels
-        for the time response.  The column labels are given by the input and
-        output (and state, when present) labels, with time labeled by 'time'
-        and traces (for multi-trace responses) labeled by 'trace'.
+        Creates a pandas data frame using the input, output, and state
+        labels for the time response.
 
         """
         if not pandas_check():
@@ -729,23 +727,16 @@ class TimeResponseData:
         import pandas
 
         # Create a dict for setting up the data frame
-        data = {'time': np.tile(
-            self.time, self.ntraces if self.ntraces > 0 else 1)}
-        if self.ntraces > 0:
-            data['trace'] = np.hstack([
-                np.full(self.time.size, label) for label in self.trace_labels])
+        data = {'time': self.time}
         if self.ninputs > 0:
             data.update(
-                {name: self.u[i].reshape(-1)
-                 for i, name in enumerate(self.input_labels)})
+                {name: self.u[i] for i, name in enumerate(self.input_labels)})
         if self.noutputs > 0:
             data.update(
-                {name: self.y[i].reshape(-1)
-                 for i, name in enumerate(self.output_labels)})
+                {name: self.y[i] for i, name in enumerate(self.output_labels)})
         if self.nstates > 0:
             data.update(
-                {name: self.x[i].reshape(-1)
-                 for i, name in enumerate(self.state_labels)})
+                {name: self.x[i] for i, name in enumerate(self.state_labels)})
 
         return pandas.DataFrame(data)
 
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 5e391d30..969c871f 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -1,29 +1,64 @@
-# xferfcn.py - transfer function class and related functions
-#
-# Original author: Richard M. Murray
-# Creation date: 24 May 2009
-# Pre-2014 revisions: Kevin K. Chen, Dec 2010
-# Use `git shortlog -n -s xferfcn.py` for full list of contributors
+"""xferfcn.py
 
-"""Transfer function representation and functions.
+Transfer function representation and functions.
 
-This module contains the TransferFunction class and also functions that
-operate on transfer functions.  This is the primary representation for the
-python-control library.
+This file contains the TransferFunction class and also functions
+that operate on transfer functions.  This is the primary representation
+for the python-control library.
+"""
+
+"""Copyright (c) 2010 by California Institute of Technology
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the California Institute of Technology nor
+   the names of its contributors may be used to endorse or promote
+   products derived from this software without specific prior
+   written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CALTECH
+OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
+
+Author: Richard M. Murray
+Date: 24 May 09
+Revised: Kevin K. Chen, Dec 10
+
+$Id$
 
 """
 
 from collections.abc import Iterable
 from copy import deepcopy
-from itertools import chain, product
+from itertools import chain
 from re import sub
 from warnings import warn
 
+# External function declarations
 import numpy as np
 import scipy as sp
-from numpy import angle, array, delete, empty, exp, finfo, float64, ndarray, \
-    nonzero, ones, pi, poly, polyadd, polymul, polyval, real, roots, sqrt, \
-    squeeze, where, zeros
+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 scipy.signal import cont2discrete, tf2zpk, zpk2tf
 
@@ -50,68 +85,46 @@ class TransferFunction(LTI):
     A class for representing transfer functions.
 
     The TransferFunction class is used to represent systems in transfer
-    function form.  Transfer functions are usually created with the
-    :func:`~control.tf` factory function.
+    function form.
 
     Parameters
     ----------
-    num : 2D list of coefficient arrays
-        Polynomial coefficients of the numerator.
-    den : 2D list of coefficient arrays
-        Polynomial coefficients of the denominator.
+    num : array_like, or list of list of array_like
+        Polynomial coefficients of the numerator
+    den : array_like, or list of list of array_like
+        Polynomial coefficients of the denominator
     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).
+        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).
+    display_format : None, 'poly' or 'zpk', optional
+        Set the display format used in printing the TransferFunction object.
+        Default behavior is polynomial display and can be changed by
+        changing config.defaults['xferfcn.display_format'].
 
     Attributes
     ----------
-    ninputs, noutputs : int
-        Number of input and output signals.
-    shape : tuple
-        2-tuple of I/O system dimension, (noutputs, ninputs).
-    input_labels, output_labels : list of str
-        Names for the input and output signals.
-    name : string, optional
-        System name.
-    num_array, den_array : 2D array of lists of float
-        Numerator and denominator polynomial coefficients as 2D array
-        of 1D array objects (of varying length).
-    num_list, den_list : 2D list of 1D array
-        Numerator and denominator polynomial coefficients as 2D lists
-        of 1D array objects (of varying length)
-    display_format : None, 'poly' or 'zpk'
-        Display format used in printing the TransferFunction object.
-        Default behavior is polynomial display and can be changed by
-        changing config.defaults['xferfcn.display_format'].
-    s : TransferFunction
-        Represents the continuous time differential operator.
-    z : TransferFunction
-        Represents the discrete time delay operator.
+    ninputs, noutputs, nstates : int
+        Number of input, output and state variables.
+    num, den : 2D list of array
+        Polynomial coefficients of the numerator and denominator.
+    dt : None, True or float
+        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).
 
     Notes
     -----
-    The numerator and denominator polynomials are stored as 2D ndarrays
-    with each element containing a 1D ndarray of coefficients.  These data
-    structures can be retrieved using ``num_array`` and ``den_array``.  For
-    example,
+    The attribues 'num' and 'den' are 2-D lists of arrays containing MIMO
+    numerator and denominator coefficients.  For example,
 
-    >>> sys.num_array[2, 5]         # doctest: +SKIP
+    >>> num[2][5] = numpy.array([1., 4., 8.])                   # doctest: +SKIP
 
-    gives the numerator of the transfer function from the 6th input to the
-    3rd output. (Note: a single 3D ndarray structure cannot be used because
-    the numerators and denominators can have different numbers of
-    coefficients in each entry.)
-
-    The attributes ``num_list`` and ``den_list`` are properties that return
-    2D nested lists containing MIMO numerator and denominator coefficients.
-    For example,
-
-    >>> sys.num_list[2][5]          # doctest: +SKIP
-
-    For legacy purposes, this list-based representation can also be
-    obtained using ``num`` and ``den``.
+    means that the numerator of the transfer function from the 6th input to
+    the 3rd output is set to s^2 + 4s + 8.
 
     A discrete time transfer function is created by specifying a nonzero
     'timebase' dt when the system is constructed:
@@ -159,15 +172,13 @@ class TransferFunction(LTI):
         Construct a transfer function.
 
         The default constructor is TransferFunction(num, den), where num and
-        den are 2D arrays of arrays containing polynomial coefficients.
+        den are lists of lists of arrays containing polynomial coefficients.
         To create a discrete time transfer funtion, use TransferFunction(num,
         den, dt) where 'dt' is the sampling time (or True for unspecified
         sampling time).  To call the copy constructor, call
         TransferFunction(sys), where sys is a TransferFunction object
         (continuous or discrete).
 
-        See :class:`TransferFunction` and :func:`tf` for more information.
-
         """
         #
         # Process positional arguments
@@ -199,8 +210,8 @@ class TransferFunction(LTI):
             raise TypeError("Needs 1, 2 or 3 arguments; received %i."
                              % len(args))
 
-        num = _clean_part(num, "numerator")
-        den = _clean_part(den, "denominator")
+        num = _clean_part(num)
+        den = _clean_part(den)
 
         #
         # Process keyword arguments
@@ -219,16 +230,13 @@ class TransferFunction(LTI):
 
         # Determine if the transfer function is static (needed for dt)
         static = True
-        for arr in [num, den]:
-            for poly in np.nditer(arr, flags=['refs_ok']):
-                if poly.item().size > 1:
+        for col in num + den:
+            for poly in col:
+                if len(poly) > 1:
                     static = False
-                    break
-            if not static:
-                break
 
         defaults = args[0] if len(args) == 1 else \
-            {'inputs': num.shape[1], 'outputs': num.shape[0]}
+            {'inputs': len(num[0]), 'outputs': len(num)}
 
         name, inputs, outputs, states, dt = _process_iosys_keywords(
                 kwargs, defaults, static=static)
@@ -244,17 +252,27 @@ class TransferFunction(LTI):
         # Check to make sure everything is consistent
         #
         # Make sure numerator and denominator matrices have consistent sizes
-        if self.ninputs != den.shape[1]:
+        if self.ninputs != len(den[0]):
             raise ValueError(
                 "The numerator has %i input(s), but the denominator has "
-                "%i input(s)." % (self.ninputs, den.shape[1]))
-        if self.noutputs != den.shape[0]:
+                "%i input(s)." % (self.ninputs, len(den[0])))
+        if self.noutputs != len(den):
             raise ValueError(
                 "The numerator has %i output(s), but the denominator has "
-                "%i output(s)." % (self.noutputs, den.shape[0]))
+                "%i output(s)." % (self.noutputs, len(den)))
 
         # Additional checks/updates on structure of the transfer function
         for i in range(self.noutputs):
+            # Make sure that each row has the same number of columns
+            if len(num[i]) != self.ninputs:
+                raise ValueError(
+                    "Row 0 of the numerator matrix has %i elements, but row "
+                    "%i has %i." % (self.ninputs, i, len(num[i])))
+            if len(den[i]) != self.ninputs:
+                raise ValueError(
+                    "Row 0 of the denominator matrix has %i elements, but row "
+                    "%i has %i." % (self.ninputs, i, len(den[i])))
+
             # Check for zeros in numerator or denominator
             # TODO: Right now these checks are only done during construction.
             # It might be worthwhile to think of a way to perform checks if the
@@ -262,8 +280,8 @@ class TransferFunction(LTI):
             for j in range(self.ninputs):
                 # Check that we don't have any zero denominators.
                 zeroden = True
-                for k in den[i, j]:
-                    if np.any(k):
+                for k in den[i][j]:
+                    if k:
                         zeroden = False
                         break
                 if zeroden:
@@ -273,16 +291,16 @@ class TransferFunction(LTI):
 
                 # If we have zero numerators, set the denominator to 1.
                 zeronum = True
-                for k in num[i, j]:
-                    if np.any(k):
+                for k in num[i][j]:
+                    if k:
                         zeronum = False
                         break
                 if zeronum:
                     den[i][j] = ones(1)
 
         # Store the numerator and denominator
-        self.num_array = num
-        self.den_array = den
+        self.num = num
+        self.den = den
 
         #
         # Final processing
@@ -307,22 +325,29 @@ class TransferFunction(LTI):
     #: :meta hide-value:
     noutputs = 1
 
-    # Numerator and denominator as lists of lists of lists
-    @property
-    def num_list(self):
-        """Numerator polynomial (as 2D nested list of 1D arrays)."""
-        return self.num_array.tolist()
-
-    @property
-    def den_list(self):
-        """Denominator polynomial (as 2D nested lists of 1D arrays)."""
-        return self.den_array.tolist()
+    #: Transfer function numerator polynomial (array)
+    #:
+    #: The numerator of the transfer function is stored as an 2D list of
+    #: arrays containing MIMO numerator coefficients, indexed by outputs and
+    #: inputs.  For example, ``num[2][5]`` is the array of coefficients for
+    #: the numerator of the transfer function from the sixth input to the
+    #: third output.
+    #:
+    #: :meta hide-value:
+    num = [[0]]
 
-    # Legacy versions (TODO: add DeprecationWarning in a later release?)
-    num, den = num_list, den_list
+    #: Transfer function denominator polynomial (array)
+    #:
+    #: The numerator of the transfer function is store as an 2D list of
+    #: arrays containing MIMO numerator coefficients, indexed by outputs and
+    #: inputs.  For example, ``den[2][5]`` is the array of coefficients for
+    #: the denominator of the transfer function from the sixth input to the
+    #: third output.
+    #:
+    #: :meta hide-value:
+    den = [[0]]
 
     def __call__(self, x, squeeze=None, warn_infinite=True):
-
         """Evaluate system's transfer function at complex frequencies.
 
         Returns the complex frequency response `sys(x)` where `x` is `s` for
@@ -402,8 +427,8 @@ class TransferFunction(LTI):
         with np.errstate(all='warn' if warn_infinite else 'ignore'):
             for i in range(self.noutputs):
                 for j in range(self.ninputs):
-                    out[i][j] = (polyval(self.num_array[i, j], x_arr) /
-                                 polyval(self.den_array[i, j], x_arr))
+                    out[i][j] = (polyval(self.num[i][j], x_arr) /
+                                 polyval(self.den[i][j], x_arr))
         return out
 
     def _truncatecoeff(self):
@@ -416,14 +441,14 @@ class TransferFunction(LTI):
         """
 
         # Beware: this is a shallow copy.  This should be okay.
-        data = [self.num_array, self.den_array]
+        data = [self.num, self.den]
         for p in range(len(data)):
             for i in range(self.noutputs):
                 for j in range(self.ninputs):
                     # Find the first nontrivial coefficient.
                     nonzero = None
-                    for k in range(data[p][i, j].size):
-                        if data[p][i, j][k]:
+                    for k in range(data[p][i][j].size):
+                        if data[p][i][j][k]:
                             nonzero = k
                             break
 
@@ -433,7 +458,7 @@ class TransferFunction(LTI):
                     else:
                         # Truncate the trivial coefficients.
                         data[p][i][j] = data[p][i][j][nonzero:]
-        [self.num_array, self.den_array] = data
+        [self.num, self.den] = data
 
     def __str__(self, var=None):
         """String representation of the transfer function.
@@ -453,18 +478,16 @@ class TransferFunction(LTI):
 
                 # Convert the numerator and denominator polynomials to strings.
                 if self.display_format == 'poly':
-                    numstr = _tf_polynomial_to_string(
-                        self.num_array[no, ni], var=var)
-                    denstr = _tf_polynomial_to_string(
-                        self.den_array[no, ni], var=var)
+                    numstr = _tf_polynomial_to_string(self.num[no][ni], var=var)
+                    denstr = _tf_polynomial_to_string(self.den[no][ni], var=var)
                 elif self.display_format == 'zpk':
-                    num = self.num_array[no, ni]
+                    num = self.num[no][ni]
                     if num.size == 1 and num.item() == 0:
                         # Catch a special case that SciPy doesn't handle
-                        z, p, k = tf2zpk([1.], self.den_array[no, ni])
+                        z, p, k = tf2zpk([1.], self.den[no][ni])
                         k = 0
                     else:
-                        z, p, k = tf2zpk(self.num[no][ni], self.den_array[no, ni])
+                        z, p, k = tf2zpk(self.num[no][ni], self.den[no][ni])
                     numstr = _tf_factorized_polynomial_to_string(
                         z, gain=k, var=var)
                     denstr = _tf_factorized_polynomial_to_string(p, var=var)
@@ -489,33 +512,20 @@ class TransferFunction(LTI):
 
     # represent to implement a re-loadable version
     def __repr__(self):
-        """Print transfer function in loadable form."""
+        """Print transfer function in loadable form"""
         if self.issiso():
             return "TransferFunction({num}, {den}{dt})".format(
-                num=self.num_array[0, 0].__repr__(),
-                den=self.den_array[0, 0].__repr__(),
+                num=self.num[0][0].__repr__(), den=self.den[0][0].__repr__(),
                 dt=', {}'.format(self.dt) if isdtime(self, strict=True)
                 else '')
         else:
-            out = "TransferFunction(["
-            for entry in [self.num_array, self.den_array]:
-                for i in range(self.noutputs):
-                    out += "[" if i == 0 else " ["
-                    for j in range(self.ninputs):
-                        out += ", " if j != 0 else ""
-                        numstr = np.array_repr(entry[i, j])
-                        out += numstr
-                    out += "]," if i < self.noutputs - 1 else "]"
-                out += "], [" if entry is self.num_array else "]"
-
-            if config.defaults['control.default_dt'] != self.dt:
-                out += ", {dt}".format(
-                    dt='None' if self.dt is None else self.dt)
-            out += ")"
-            return out
+            return "TransferFunction({num}, {den}{dt})".format(
+                num=self.num.__repr__(), den=self.den.__repr__(),
+                dt=', {}'.format(self.dt) if isdtime(self, strict=True)
+                else '')
 
     def _repr_latex_(self, var=None):
-        """LaTeX representation of transfer function, for Jupyter notebook."""
+        """LaTeX representation of transfer function, for Jupyter notebook"""
 
         mimo = not self.issiso()
 
@@ -531,13 +541,10 @@ class TransferFunction(LTI):
             for ni in range(self.ninputs):
                 # Convert the numerator and denominator polynomials to strings.
                 if self.display_format == 'poly':
-                    numstr = _tf_polynomial_to_string(
-                        self.num_array[no, ni], var=var)
-                    denstr = _tf_polynomial_to_string(
-                        self.den_array[no, ni], var=var)
+                    numstr = _tf_polynomial_to_string(self.num[no][ni], var=var)
+                    denstr = _tf_polynomial_to_string(self.den[no][ni], var=var)
                 elif self.display_format == 'zpk':
-                    z, p, k = tf2zpk(
-                        self.num_array[no, ni], self.den_array[no, ni])
+                    z, p, k = tf2zpk(self.num[no][ni], self.den[no][ni])
                     numstr = _tf_factorized_polynomial_to_string(
                         z, gain=k, var=var)
                     denstr = _tf_factorized_polynomial_to_string(p, var=var)
@@ -566,10 +573,10 @@ class TransferFunction(LTI):
 
     def __neg__(self):
         """Negate a transfer function."""
-        num = deepcopy(self.num_array)
+        num = deepcopy(self.num)
         for i in range(self.noutputs):
             for j in range(self.ninputs):
-                num[i, j] *= -1
+                num[i][j] *= -1
         return TransferFunction(num, self.den, self.dt)
 
     def __add__(self, other):
@@ -599,14 +606,14 @@ class TransferFunction(LTI):
         dt = common_timebase(self.dt, other.dt)
 
         # Preallocate the numerator and denominator of the sum.
-        num = _create_poly_array((self.noutputs, self.ninputs))
-        den = _create_poly_array((self.noutputs, self.ninputs))
+        num = [[[] for j in range(self.ninputs)] for i in range(self.noutputs)]
+        den = [[[] for j in range(self.ninputs)] for i in range(self.noutputs)]
 
         for i in range(self.noutputs):
             for j in range(self.ninputs):
-                num[i, j], den[i, j] = _add_siso(
-                    self.num_array[i, j], self.den_array[i, j],
-                    other.num_array[i, j], other.den_array[i, j])
+                num[i][j], den[i][j] = _add_siso(
+                    self.num[i][j], self.den[i][j],
+                    other.num[i][j], other.den[i][j])
 
         return TransferFunction(num, den, dt)
 
@@ -641,14 +648,14 @@ class TransferFunction(LTI):
                 "C = A * B: A has %i column(s) (input(s)), but B has %i "
                 "row(s)\n(output(s))." % (self.ninputs, other.noutputs))
 
-        ninputs = other.ninputs
-        noutputs = self.noutputs
+        inputs = other.ninputs
+        outputs = self.noutputs
 
         dt = common_timebase(self.dt, other.dt)
 
         # Preallocate the numerator and denominator of the sum.
-        num = _create_poly_array((noutputs, ninputs), [0])
-        den = _create_poly_array((noutputs, ninputs), [1])
+        num = [[[0] for j in range(inputs)] for i in range(outputs)]
+        den = [[[1] for j in range(inputs)] for i in range(outputs)]
 
         # Temporary storage for the summands needed to find the (i, j)th
         # element of the product.
@@ -656,16 +663,17 @@ class TransferFunction(LTI):
         den_summand = [[] for k in range(self.ninputs)]
 
         # Multiply & add.
-        for row in range(noutputs):
-            for col in range(ninputs):
+        for row in range(outputs):
+            for col in range(inputs):
                 for k in range(self.ninputs):
                     num_summand[k] = polymul(
-                        self.num_array[row, k], other.num_array[k, col])
+                        self.num[row][k], other.num[k][col])
                     den_summand[k] = polymul(
-                        self.den_array[row, k], other.den_array[k, col])
-                    num[row, col], den[row, col] = _add_siso(
-                        num[row, col], den[row, col],
+                        self.den[row][k], other.den[k][col])
+                    num[row][col], den[row][col] = _add_siso(
+                        num[row][col], den[row][col],
                         num_summand[k], den_summand[k])
+
         return TransferFunction(num, den, dt)
 
     def __rmul__(self, other):
@@ -684,14 +692,14 @@ class TransferFunction(LTI):
                 "C = A * B: A has %i column(s) (input(s)), but B has %i "
                 "row(s)\n(output(s))." % (other.ninputs, self.noutputs))
 
-        ninputs = self.ninputs
-        noutputs = other.noutputs
+        inputs = self.ninputs
+        outputs = other.noutputs
 
         dt = common_timebase(self.dt, other.dt)
 
         # Preallocate the numerator and denominator of the sum.
-        num = _create_poly_array((noutputs, ninputs), [0])
-        den = _create_poly_array((noutputs, ninputs), [1])
+        num = [[[0] for j in range(inputs)] for i in range(outputs)]
+        den = [[[1] for j in range(inputs)] for i in range(outputs)]
 
         # Temporary storage for the summands needed to find the
         # (i, j)th element
@@ -699,15 +707,13 @@ class TransferFunction(LTI):
         num_summand = [[] for k in range(other.ninputs)]
         den_summand = [[] for k in range(other.ninputs)]
 
-        for i in range(noutputs):  # Iterate through rows of product.
-            for j in range(ninputs):  # Iterate through columns of product.
+        for i in range(outputs):  # Iterate through rows of product.
+            for j in range(inputs):  # Iterate through columns of product.
                 for k in range(other.ninputs):  # Multiply & add.
-                    num_summand[k] = polymul(
-                        other.num_array[i, k], self.num_array[k, j])
-                    den_summand[k] = polymul(
-                        other.den_array[i, k], self.den_array[k, j])
+                    num_summand[k] = polymul(other.num[i][k], self.num[k][j])
+                    den_summand[k] = polymul(other.den[i][k], self.den[k][j])
                     num[i][j], den[i][j] = _add_siso(
-                        num[i, j], den[i, j],
+                        num[i][j], den[i][j],
                         num_summand[k], den_summand[k])
 
         return TransferFunction(num, den, dt)
@@ -730,8 +736,8 @@ class TransferFunction(LTI):
 
         dt = common_timebase(self.dt, other.dt)
 
-        num = polymul(self.num_array[0, 0], other.den_array[0, 0])
-        den = polymul(self.den_array[0, 0], other.num_array[0, 0])
+        num = polymul(self.num[0][0], other.den[0][0])
+        den = polymul(self.den[0][0], other.num[0][0])
 
         return TransferFunction(num, den, dt)
 
@@ -779,14 +785,15 @@ class TransferFunction(LTI):
             indices[1], self.input_labels, slice_to_list=True)
         
         # Construct the transfer function for the subsyste
-        num = _create_poly_array((len(outputs), len(inputs)))
-        den = _create_poly_array(num.shape)
-        for row, i in enumerate(outdx):
-            for col, j in enumerate(inpdx):
-                num[row, col] = self.num_array[i, j]
-                den[row, col] = self.den_array[i, j]
-                col += 1
-            row += 1
+        num, den = [], []
+        for i in outdx:
+            num_i = []
+            den_i = []
+            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)
 
         # Create the system name
         sysname = config.defaults['iosys.indexed_system_name_prefix'] + \
@@ -825,7 +832,7 @@ class TransferFunction(LTI):
                 "for SISO systems.")
         else:
             # for now, just give zeros of a SISO tf
-            return roots(self.num_array[0, 0]).astype(complex)
+            return roots(self.num[0][0]).astype(complex)
 
     def feedback(self, other=1, sign=-1):
         """Feedback interconnection between two LTI objects."""
@@ -839,10 +846,10 @@ class TransferFunction(LTI):
                 "MIMO systems.")
         dt = common_timebase(self.dt, other.dt)
 
-        num1 = self.num_array[0, 0]
-        den1 = self.den_array[0, 0]
-        num2 = other.num_array[0, 0]
-        den2 = other.den_array[0, 0]
+        num1 = self.num[0][0]
+        den1 = self.den[0][0]
+        num2 = other.num[0][0]
+        den2 = other.den[0][0]
 
         num = polymul(num1, den2)
         den = polyadd(polymul(den2, den1), -sign * polymul(num2, num1))
@@ -855,7 +862,7 @@ class TransferFunction(LTI):
         # large.
 
     def minreal(self, tol=None):
-        """Remove cancelling pole/zero pairs from a transfer function."""
+        """Remove cancelling pole/zero pairs from a transfer function"""
         # based on octave minreal
 
         # default accuracy
@@ -863,17 +870,17 @@ class TransferFunction(LTI):
         sqrt_eps = sqrt(float_info.epsilon)
 
         # pre-allocate arrays
-        num = _create_poly_array((self.noutputs, self.ninputs))
-        den = _create_poly_array((self.noutputs, self.ninputs))
+        num = [[[] for j in range(self.ninputs)] for i in range(self.noutputs)]
+        den = [[[] for j in range(self.ninputs)] for i in range(self.noutputs)]
 
         for i in range(self.noutputs):
             for j in range(self.ninputs):
 
                 # split up in zeros, poles and gain
                 newzeros = []
-                zeros = roots(self.num_array[i, j])
-                poles = roots(self.den_array[i, j])
-                gain = self.num_array[i, j][0] / self.den_array[i, j][0]
+                zeros = roots(self.num[i][j])
+                poles = roots(self.den[i][j])
+                gain = self.num[i][j][0] / self.den[i][j][0]
 
                 # check all zeros
                 for z in zeros:
@@ -888,19 +895,19 @@ class TransferFunction(LTI):
                         newzeros.append(z)
 
                 # poly([]) returns a scalar, but we always want a 1d array
-                num[i, j] = np.atleast_1d(gain * real(poly(newzeros)))
-                den[i, j] = np.atleast_1d(real(poly(poles)))
+                num[i][j] = np.atleast_1d(gain * real(poly(newzeros)))
+                den[i][j] = np.atleast_1d(real(poly(poles)))
 
         # end result
         return TransferFunction(num, den, self.dt)
 
     def returnScipySignalLTI(self, strict=True):
-        """Return a 2D array of :class:`scipy.signal.lti` objects.
+        """Return a list of a list of :class:`scipy.signal.lti` objects.
 
         For instance,
 
         >>> out = tfobject.returnScipySignalLTI()               # doctest: +SKIP
-        >>> out[3, 5]                                           # doctest: +SKIP
+        >>> out[3][5]                                           # doctest: +SKIP
 
         is a :class:`scipy.signal.lti` object corresponding to the
         transfer function from the 6th input to the 4th output.
@@ -1108,7 +1115,7 @@ class TransferFunction(LTI):
 
     def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None,
                name=None, copy_names=True, **kwargs):
-        """Convert a continuous-time system to discrete time.
+        """Convert a continuous-time system to discrete time
 
         Creates a discrete-time system from a continuous-time system by
         sampling.  Multiple methods of conversion are supported.
@@ -1255,7 +1262,7 @@ class TransferFunction(LTI):
     # the class attributes are set at the bottom of the file to avoid problems
     # with recursive calls.
 
-    #: Differentation operator (continuous time).
+    #: Differentation operator (continuous time)
     #:
     #: The ``s`` constant can be used to create continuous time transfer
     #: functions using algebraic expressions.
@@ -1268,7 +1275,7 @@ class TransferFunction(LTI):
     #: :meta hide-value:
     s = None
 
-    #: Delay operator (discrete time).
+    #: Delay operator (discrete time)
     #:
     #: The ``z`` constant can be used to create discrete time transfer
     #: functions using algebraic expressions.
@@ -1312,7 +1319,7 @@ def _c2d_matched(sysC, Ts, **kwargs):
 # Utility function to convert a transfer function polynomial to a string
 # Borrowed from poly1d library
 def _tf_polynomial_to_string(coeffs, var='s'):
-    """Convert a transfer function polynomial to a string."""
+    """Convert a transfer function polynomial to a string"""
 
     thestr = "0"
 
@@ -1359,7 +1366,7 @@ def _tf_polynomial_to_string(coeffs, var='s'):
 
 
 def _tf_factorized_polynomial_to_string(roots, gain=1, var='s'):
-    """Convert a factorized polynomial to a string."""
+    """Convert a factorized polynomial to a string"""
 
     if roots.size == 0:
         return _float2str(gain)
@@ -1402,11 +1409,9 @@ def _tf_factorized_polynomial_to_string(roots, gain=1, var='s'):
 
 
 def _tf_string_to_latex(thestr, var='s'):
-    """Superscript all digits in a polynomial string and convert
-    float coefficients in scientific notation to prettier LaTeX
-    representation.
-
-    """
+    """ make sure to superscript all digits in a polynomial string
+        and convert float coefficients in scientific notation
+        to prettier LaTeX representation """
     # TODO: make the multiplication sign configurable
     expmul = r' \\times'
     thestr = sub(var + r'\^(\d{2,})', var + r'^{\1}', thestr)
@@ -1549,20 +1554,15 @@ def tf(*args, **kwargs):
         If `num` and `den` are 1D array_like objects, the function creates a
         SISO system.
 
-        To create a MIMO system, `num` and `den` need to be 2D arrays of
-        of array_like objects (a 3 dimensional data structure in total;
-        for details see note below).  If the denominator for all transfer
-        function is the same, `den` can be specified as a 1D array.
+        To create a MIMO system, `num` and `den` need to be 2D nested lists
+        of array_like objects. (A 3 dimensional data structure in total.)
+        (For details see note below.)
 
     ``tf(num, den, dt)``
         Create a discrete time transfer function system; dt can either be a
         positive number indicating the sampling time or 'True' if no
         specific timebase is given.
 
-    ``tf([[G11, ..., G1m], ..., [Gp1, ..., Gpm]][, dt])``
-        Create a pxm MIMO system from SISO transfer functions Gij.  See
-        :func:`combine_tf` for more details.
-
     ``tf('s')`` or ``tf('z')``
         Create a transfer function representing the differential operator
         ('s') or delay operator ('z').
@@ -1570,27 +1570,20 @@ def tf(*args, **kwargs):
     Parameters
     ----------
     sys : LTI (StateSpace or TransferFunction)
-        A linear system that will be converted to a transfer function.
-    arr : 2D list of TransferFunction
-        2D list of SISO transfer functions to create MIMO transfer function.
+        A linear system
     num : array_like, or list of list of array_like
-        Polynomial coefficients of the numerator.
+        Polynomial coefficients of the numerator
     den : array_like, or list of list of array_like
-        Polynomial coefficients of the denominator.
-    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).
+        Polynomial coefficients of the denominator
     display_format : None, 'poly' or 'zpk'
         Set the display format used in printing the TransferFunction object.
         Default behavior is polynomial display and can be changed by
-        changing config.defaults['xferfcn.display_format'].
+        changing config.defaults['xferfcn.display_format']..
 
     Returns
     -------
-    sys : TransferFunction
-        The new linear system.
+    out: :class:`TransferFunction`
+        The new linear system
 
     Other Parameters
     ----------------
@@ -1598,8 +1591,6 @@ def tf(*args, **kwargs):
         List of strings that name the individual signals of the transformed
         system.  If not given, the inputs and outputs are the same as the
         original system.
-    input_prefix, output_prefix : string, optional
-        Set the prefix for input and output signals.  Defaults = 'u', 'y'.
     name : string, optional
         System name. If unspecified, a generic name <sys[id]> is generated
         with a unique integer id.
@@ -1607,9 +1598,9 @@ def tf(*args, **kwargs):
     Raises
     ------
     ValueError
-        If `num` and `den` have invalid or unequal dimensions.
+        if `num` and `den` have invalid or unequal dimensions
     TypeError
-        If `num` or `den` are of incorrect type.
+        if `num` or `den` are of incorrect type
 
     See Also
     --------
@@ -1620,10 +1611,9 @@ def tf(*args, **kwargs):
 
     Notes
     -----
-    MIMO transfer functions are created by passing a 2D array of coeffients:
     ``num[i][j]`` contains the polynomial coefficients of the numerator
-    for the transfer function from the (j+1)st input to the (i+1)st output,
-    and ``den[i][j]`` works the same way.
+    for the transfer function from the (j+1)st input to the (i+1)st output.
+    ``den[i][j]`` works the same way.
 
     The list ``[2, 3, 4]`` denotes the polynomial :math:`2s^2 + 3s + 4`.
 
@@ -1648,7 +1638,11 @@ def tf(*args, **kwargs):
     >>> sys_tf = ct.tf(sys_ss)
 
     """
-    if len(args) == 1 and isinstance(args[0], str):
+
+    if len(args) == 2 or len(args) == 3:
+        return TransferFunction(*args, **kwargs)
+
+    elif len(args) == 1 and isinstance(args[0], str):
         # Make sure there were no extraneous keywords
         if kwargs:
             raise TypeError("unrecognized keywords: ", str(kwargs))
@@ -1659,53 +1653,19 @@ def tf(*args, **kwargs):
         elif args[0] == 'z':
             return TransferFunction.z
 
-    elif len(args) == 1 and isinstance(args[0], list):
-        # Allow passing an array of SISO transfer functions
-        from .bdalg import combine_tf
-        return combine_tf(*args)
-
     elif len(args) == 1:
         from .statesp import StateSpace
-        if isinstance(sys := args[0], StateSpace):
+        sys = args[0]
+        if isinstance(sys, StateSpace):
             return ss2tf(sys, **kwargs)
         elif isinstance(sys, TransferFunction):
             # Use copy constructor
             return TransferFunction(sys, **kwargs)
-        elif isinstance(data := args[0], np.ndarray) and data.ndim == 2 or \
-             isinstance(data, list) and isinstance(data[0], list):
-            raise NotImplementedError(
-                "arrays of transfer functions not (yet) supported")
         else:
             raise TypeError("tf(sys): sys must be a StateSpace or "
                             "TransferFunction object.   It is %s." % type(sys))
-
-    elif len(args) == 3:
-        if 'dt' in kwargs:
-            warn("received multiple dt arguments, "
-                 f"using positional arg {args[2]}")
-        kwargs['dt'] = args[2]
-        args = args[:2]
-
-    elif len(args) != 2:
-        raise ValueError("Needs 1, 2, or 3 arguments; received %i." % len(args))
-
-    #
-    # Process the numerator and denominator arguments
-    #
-    # If we got through to here, we have two argume nts (num, den) and
-    # the keywords (including dt).  The only thing left to do is look
-    # for some special cases, like having a common denominator.
-    #
-    num, den = args
-
-    num = _clean_part(num, "numerator")
-    den = _clean_part(den, "denominator")
-
-    if den.size == 1 and num.size > 1:
-        # Broadcast denominator to shape of numerator
-        den = np.broadcast_to(den, num.shape).copy()
-
-    return TransferFunction(num, den, **kwargs)
+    else:
+        raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
 
 
 def zpk(zeros, poles, gain, *args, **kwargs):
@@ -1727,7 +1687,7 @@ def zpk(zeros, poles, gain, *args, **kwargs):
     poles : array_like
         Array containing the location of poles.
     gain : float
-        System gain.
+        System gain
     dt : None, True or float, optional
         System timebase. 0 (default) indicates continuous
         time, True indicates discrete time with unspecified sampling
@@ -1749,7 +1709,7 @@ def zpk(zeros, poles, gain, *args, **kwargs):
 
     Returns
     -------
-    out: `TransferFunction`
+    out: :class:`TransferFunction`
         Transfer function with given zeros, poles, and gain.
 
     Examples
@@ -1885,7 +1845,7 @@ def tfdata(sys):
     return tf.num, tf.den
 
 
-def _clean_part(data, name="<unknown>"):
+def _clean_part(data):
     """
     Return a valid, cleaned up numerator or denominator
     for the TransferFunction class.
@@ -1901,37 +1861,27 @@ def _clean_part(data, name="<unknown>"):
     valid_types = (int, float, complex, np.number)
     valid_collection = (list, tuple, ndarray)
 
-    if isinstance(data, np.ndarray) and data.ndim == 2 and \
-       data.dtype == object and isinstance(data[0, 0], np.ndarray):
-        # Data is already in the right format
-        return data
-    elif isinstance(data, ndarray) and data.ndim == 3 and \
-          isinstance(data[0, 0, 0], valid_types):
-        out = np.empty(data.shape[0:2], dtype=np.ndarray)
-        for i, j in product(range(out.shape[0]), range(out.shape[1])):
-            out[i, j] = data[i, j, :]
-    elif (isinstance(data, valid_types) or
+    if (isinstance(data, valid_types) or
             (isinstance(data, ndarray) and data.ndim == 0)):
         # Data is a scalar (including 0d ndarray)
-        out = np.empty((1,1), dtype=np.ndarray)
-        out[0, 0] = array([data])
+        data = [[array([data])]]
+    elif (isinstance(data, ndarray) and data.ndim == 3 and
+          isinstance(data[0, 0, 0], valid_types)):
+        data = [[array(data[i, j])
+                 for j in range(data.shape[1])]
+                for i in range(data.shape[0])]
     elif (isinstance(data, valid_collection) and
             all([isinstance(d, valid_types) for d in data])):
-        out = np.empty((1,1), dtype=np.ndarray)
-        out[0, 0] = array(data)
-    elif isinstance(data, (list, tuple)) and \
-         isinstance(data[0], (list, tuple)) and \
-         (isinstance(data[0][0], valid_collection) and
-          all([isinstance(d, valid_types) for d in data[0][0]]) or \
-          isinstance(data[0][0], valid_types)):
-        out = np.empty((len(data), len(data[0])), dtype=np.ndarray)
-        for i in range(out.shape[0]):
-            if len(data[i]) != out.shape[1]:
-                raise ValueError(
-                    "Row 0 of the %s matrix has %i elements, but row "
-                    "%i has %i." % (name, out.shape[1], i, len(data[i])))
-            for j in range(out.shape[1]):
-                out[i, j] = np.atleast_1d(data[i][j])
+        data = [[array(data)]]
+    elif (isinstance(data, (list, tuple)) and
+          isinstance(data[0], (list, tuple)) and
+          (isinstance(data[0][0], valid_collection) and
+           all([isinstance(d, valid_types) for d in data[0][0]]))):
+        data = list(data)
+        for j in range(len(data)):
+            data[j] = list(data[j])
+            for k in range(len(data[j])):
+                data[j][k] = array(data[j][k])
     else:
         # If the user passed in anything else, then it's unclear what
         # the meaning is.
@@ -1940,36 +1890,20 @@ def _clean_part(data, name="<unknown>"):
             "(for\nSISO), or lists of lists of vectors (for SISO or MIMO).")
 
     # Check for coefficients that are ints and convert to floats
-    for i in range(out.shape[0]):
-        for j in range(out.shape[1]):
-            for k in range(len(out[i, j])):
-                if isinstance(out[i, j][k], (int, np.int32, np.int64)):
-                    out[i, j][k] = float(out[i, j][k])
-    return out
+    for i in range(len(data)):
+        for j in range(len(data[i])):
+            for k in range(len(data[i][j])):
+                if isinstance(data[i][j][k], (int, np.int32, np.int64)):
+                    data[i][j][k] = float(data[i][j][k])
 
+    return data
 
-#
-# Define constants to represent differentiation, unit delay.
-#
-# Set the docstring explicitly to avoid having Sphinx document this as
-# a method instead of a property/attribute.
 
+# Define constants to represent differentiation, unit delay
 TransferFunction.s = TransferFunction([1, 0], [1], 0, name='s')
-TransferFunction.s.__doc__ = "Differentation operator (continuous time)."
-
 TransferFunction.z = TransferFunction([1, 0], [1], True, name='z')
-TransferFunction.z.__doc__ = "Delay operator (discrete time)."
 
 
 def _float2str(value):
     _num_format = config.defaults.get('xferfcn.floating_point_format', ':.4g')
     return f"{value:{_num_format}}"
-
-
-def _create_poly_array(shape, default=None):
-    out = np.empty(shape, dtype=np.ndarray)
-    if default is not None:
-        default = np.array(default)
-        for i, j in product(range(shape[0]), range(shape[1])):
-            out[i, j] = default
-    return out
diff --git a/examples/steering-gainsched.py b/examples/steering-gainsched.py
index 36dafd61..85e8d8bd 100644
--- a/examples/steering-gainsched.py
+++ b/examples/steering-gainsched.py
@@ -46,10 +46,10 @@ def vehicle_output(t, x, u, params):
     return x                            # return x, y, theta (full state)
 
 # Define the vehicle steering dynamics as an input/output system
-vehicle = ct.nlsys(
+vehicle = ct.NonlinearIOSystem(
     vehicle_update, vehicle_output, states=3, name='vehicle',
-    inputs=('v', 'phi'), outputs=('x', 'y', 'theta'),
-    params={'wheelbase': 3, 'maxsteer': 0.5})
+    inputs=('v', 'phi'),
+    outputs=('x', 'y', 'theta'))
 
 #
 # Gain scheduled controller
@@ -89,12 +89,10 @@ def control_output(t, x, u, params):
     return  np.array([v, phi])
 
 # Define the controller as an input/output system
-controller = ct.nlsys(
+controller = ct.NonlinearIOSystem(
     None, control_output, name='controller',        # static system
     inputs=('ex', 'ey', 'etheta', 'vd', 'phid'),    # system inputs
-    outputs=('v', 'phi'),                           # system outputs
-    params={'longpole': -2, 'latpole1': -1/2 + sqrt(-7)/2,
-            'latpole2': -1/2 - sqrt(-7)/2, 'wheelbase': 3}
+    outputs=('v', 'phi')                            # system outputs
 )
 
 #
@@ -115,7 +113,7 @@ def trajgen_output(t, x, u, params):
     return np.array([vref * t, yref, 0, vref, 0])
 
 # Define the trajectory generator as an input/output system
-trajgen = ct.nlsys(
+trajgen = ct.NonlinearIOSystem(
     None, trajgen_output, name='trajgen',
     inputs=('vref', 'yref'),
     outputs=('xd', 'yd', 'thetad', 'vd', 'phid'))
@@ -158,13 +156,10 @@ steering = ct.interconnect(
     inplist=['trajgen.vref', 'trajgen.yref'],
     inputs=['yref', 'vref'],
 
-    # System outputs
+    #  System outputs
     outlist=['vehicle.x', 'vehicle.y', 'vehicle.theta', 'controller.v',
              'controller.phi'],
-    outputs=['x', 'y', 'theta', 'v', 'phi'],
-
-    # Parameters
-    params=trajgen.params | vehicle.params | controller.params,
+    outputs=['x', 'y', 'theta', 'v', 'phi']
 )
 
 # Set up the simulation conditions
@@ -225,8 +220,7 @@ gains = [np.array(ct.lqr(vehicle.linearize(
 # Create the gain scheduled system
 controller, _ = ct.create_statefbk_iosystem(
     vehicle, (gains, points), name='controller', ud_labels=['vd', 'phid'],
-    gainsched_indices=['vd', 'theta'], gainsched_method='linear',
-    params=vehicle.params | controller.params)
+    gainsched_indices=['vd', 'theta'], gainsched_method='linear')
 
 # Connect everything together (note that controller inputs are different)
 steering = ct.interconnect(
@@ -251,13 +245,10 @@ steering = ct.interconnect(
     inplist=['trajgen.vref', 'trajgen.yref'],
     inputs=['yref', 'vref'],
 
-    # System outputs
+    #  System outputs
     outlist=['vehicle.x', 'vehicle.y', 'vehicle.theta', 'controller.v',
              'controller.phi'],
-    outputs=['x', 'y', 'theta', 'v', 'phi'],
-
-    # Parameters
-    params=steering.params
+    outputs=['x', 'y', 'theta', 'v', 'phi']
 )
 
 # Plot the results to compare to the previous case


Test patch:
diff --git a/control/tests/bdalg_test.py b/control/tests/bdalg_test.py
index e0d64f01..8ea67e0f 100644
--- a/control/tests/bdalg_test.py
+++ b/control/tests/bdalg_test.py
@@ -903,15 +903,15 @@ def _tf_close_coeff(tf_a, tf_b, rtol=1e-5, atol=1e-8):
     for i in range(tf_a.noutputs):
         for j in range(tf_a.ninputs):
             if not np.allclose(
-                tf_a.num_array[i, j],
-                tf_b.num_array[i, j],
+                tf_a.num[i][j],
+                tf_b.num[i][j],
                 rtol=rtol,
                 atol=atol,
             ):
                 return False
             if not np.allclose(
-                tf_a.den_array[i, j],
-                tf_b.den_array[i, j],
+                tf_a.den[i][j],
+                tf_b.den[i][j],
                 rtol=rtol,
                 atol=atol,
             ):
diff --git a/control/tests/docstrings_test.py b/control/tests/docstrings_test.py
index 16647895..991ead3e 100644
--- a/control/tests/docstrings_test.py
+++ b/control/tests/docstrings_test.py
@@ -21,9 +21,6 @@ function_skiplist = [
     control.ControlPlot.reshape,                # needed for legacy interface
     control.phase_plot,                         # legacy function
     control.drss,                               # documention in rss
-    control.frd,                                # tested separately below
-    control.ss,                                 # tested separately below
-    control.tf,                                 # tested separately below
 ]
 
 # Checksums to use for checking whether a docstring has changed
@@ -34,10 +31,13 @@ function_docstring_hash = {
     control.dlqr:                       '896cfa651dbbd80e417635904d13c9d6',
     control.lqe:                        '567bf657538935173f2e50700ba87168',
     control.lqr:                        'a3e0a85f781fc9c0f69a4b7da4f0bd22',
+    control.frd:                        '099464bf2d14f25a8769ef951adf658b',
     control.margin:                     'f02b3034f5f1d44ce26f916cc3e51600',
     control.parallel:                   '025c5195a34c57392223374b6244a8c4',
     control.series:                     '9aede1459667738f05cf4fc46603a4f6',
+    control.ss:                         '1b9cfad5dbdf2f474cfdeadf5cb1ad80',
     control.ss2tf:                      '48ff25d22d28e7b396e686dd5eb58831',
+    control.tf:                         '53a13f4a7f75a31c81800e10c88730ef',
     control.tf2ss:                      '086a3692659b7321c2af126f79f4bc11',
     control.markov:                     'a4199c54cb50f07c0163d3790739eafe',
     control.gangof4:                    '0e52eb6cf7ce024f9a41f3ae3ebf04f7',
@@ -260,244 +260,14 @@ def test_deprecated_functions(module, prefix):
                         f"{name} deprecated but w/ non-standard docs/warnings")
                 assert name != 'ss2io'
 
-#
-# Tests for I/O system classes
-#
-# The tests below try to make sure that we document I/O system classes
-# and the factory functions that create them in a uniform way.
-#
-
-ct = control
-fs = control.flatsys
-
-# Dictionary of factory functions associated with primary classes
-class_factory_function = {
-    fs.FlatSystem: fs.flatsys,
-    ct.FrequencyResponseData: ct.frd,
-    ct.InterconnectedSystem: ct.interconnect,
-    ct.LinearICSystem: ct.interconnect,
-    ct.NonlinearIOSystem: ct.nlsys,
-    ct.StateSpace: ct.ss,
-    ct.TransferFunction: ct.tf,
-}
-
-#
-# List of arguments described in class docstrings
-#
-# These are the minimal arguments needed to initialize the class.  Optional
-# arguments should be documented in the factory functions and do not need
-# to be duplicated in the class documentation.
-#
-class_args = {
-    fs.FlatSystem: ['forward', 'reverse'],
-    ct.FrequencyResponseData: ['response', 'omega', 'dt'],
-    ct.NonlinearIOSystem: [
-        'updfcn', 'outfcn', 'inputs', 'outputs', 'states', 'params', 'dt'],
-    ct.StateSpace: ['A', 'B', 'C', 'D', 'dt'],
-    ct.TransferFunction: ['num', 'den', 'dt'],
-}
-
-#
-# List of attributes described in class docstrings
-#
-# This is the list of attributes for the class that are not already listed
-# as parameters used to initialize the class.  These should all be defined
-# in the class docstring.
-#
-# Attributes that are part of all I/O system classes should be listed in
-# `std_class_attributes`.  Attributes that are not commonly needed are
-# defined as part of a parent class can just be documented there, and
-# should be listed in `iosys_parent_attributes` (these will be searched
-# using the MRO).
-
-std_class_attributes = [
-    'ninputs', 'noutputs', 'input_labels', 'output_labels', 'name', 'shape']
-
-# List of attributes defined for specific I/O systems
-class_attributes = {
-    fs.FlatSystem: [],
-    ct.FrequencyResponseData: [],
-    ct.NonlinearIOSystem: ['nstates', 'state_labels'],
-    ct.StateSpace: ['nstates', 'state_labels'],
-    ct.TransferFunction: [],
-}
-
-# List of attributes defined in a parent class (no need to warn)
-iosys_parent_attributes = [
-    'input_index', 'output_index', 'state_index',       # rarely used
-    'states', 'nstates', 'state_labels',                # not need in TF, FRD
-    'params', 'outfcn', 'updfcn'                        # NL I/O, SS overlap
-]
-
-#
-# List of arguments described (only) in factory function docstrings
-#
-# These lists consist of the arguments that should be documented in the
-# factory functions and should not be duplicated in the class
-# documentation, even though in some cases they are actually processed in
-# the class __init__ function.
-#
-std_factory_args = [
-    'inputs', 'outputs', 'name', 'input_prefix', 'output_prefix']
-
-factory_args = {
-    fs.flatsys: ['states', 'state_prefix'],
-    ct.frd: ['sys'],
-    ct.nlsys: ['state_prefix'],
-    ct.ss: ['sys', 'states', 'state_prefix'],
-    ct.tf: ['sys'],
-}
-
-
-@pytest.mark.parametrize(
-    "cls, fcn, args",
-    [(cls, class_factory_function[cls], class_args[cls])
-     for cls in class_args.keys()])
-def test_iosys_primary_classes(cls, fcn, args):
-    docstring = inspect.getdoc(cls)
-
-    # Make sure the typical arguments are there
-    for argname in args + std_class_attributes + class_attributes[cls]:
-        _check_parameter_docs(cls.__name__, argname, docstring)
-
-    # Make sure we reference the factory function
-    if re.search(
-            r"created.*(with|by|using).*the[\s]*"
-            f":func:`~control\\..*{fcn.__name__}`"
-            r"[\s]factory[\s]function", docstring, re.DOTALL) is None:
-        pytest.fail(
-            f"{cls.__name__} does not reference factory function "
-            f"{fcn.__name__}")
-
-    # Make sure we don't reference parameters from the factory function
-    for argname in factory_args[fcn]:
-        if re.search(f"[\\s]+{argname}(, .*)*[\\s]*:", docstring) is not None:
-            pytest.fail(
-                f"{cls.__name__} references factory function parameter "
-                f"'{argname}'")
-
-
-@pytest.mark.parametrize("cls", class_args.keys())
-def test_iosys_attribute_lists(cls, ignore_future_warning):
-    fcn = class_factory_function[cls]
-
-    # Create a system that we can scan for attributes
-    sys = ct.rss(2, 1, 1)
-    ignore_args = []
-    match fcn:
-        case ct.tf:
-            sys = ct.tf(sys)
-            ignore_args = ['state_labels']
-        case ct.frd:
-            sys = ct.frd(sys, [0.1, 1, 10])
-            ignore_args = ['state_labels']
-        case ct.nlsys:
-            sys = ct.nlsys(sys)
-        case fs.flatsys:
-            sys = fs.flatsys(sys)
-            sys = fs.flatsys(sys.forward, sys.reverse)
-
-    docstring = inspect.getdoc(cls)
-    for name, obj in inspect.getmembers(sys):
-        if name.startswith('_') or inspect.ismethod(obj) or name in ignore_args:
-            # Skip hidden variables; class methods are checked elsewhere
-            continue
-
-        # Try to find documentation in primary class
-        if _check_parameter_docs(
-                cls.__name__, name, docstring, fail_if_missing=False):
-            continue
-
-        # Couldn't find in main documentation; look in parent classes
-        for parent in cls.__mro__:
-            if parent == object:
-                pytest.fail(
-                    f"{cls.__name__} attribute '{name}' not documented")
-
-            if _check_parameter_docs(
-                    parent.__name__, name, inspect.getdoc(parent),
-                    fail_if_missing=False):
-                if name not in iosys_parent_attributes + factory_args[fcn]:
-                    warnings.warn(
-                        f"{cls.__name__} attribute '{name}' only documented "
-                        f"in parent class {parent.__name__}")
-                break
-
-
-@pytest.mark.parametrize("cls", [ct.InputOutputSystem, ct.LTI])
-def test_iosys_container_classes(cls):
-    # Create a system that we can scan for attributes
-    sys = cls(states=2, outputs=1, inputs=1)
-
-    docstring = inspect.getdoc(cls)
-    for name, obj in inspect.getmembers(sys):
-        if name.startswith('_') or inspect.ismethod(obj):
-            # Skip hidden variables; class methods are checked elsewhere
-            continue
-
-        # Look through all classes in hierarchy
-        if verbose:
-            print(f"{name=}")
-        for parent in cls.__mro__:
-            if parent == object:
-                pytest.fail(
-                    f"{cls.__name__} attribute '{name}' not documented")
-
-            if verbose:
-                print(f"  {parent=}")
-            if _check_parameter_docs(
-                    parent.__name__, name, inspect.getdoc(parent),
-                    fail_if_missing=False):
-                break
-
-
-@pytest.mark.parametrize("cls", [ct.InterconnectedSystem, ct.LinearICSystem])
-def test_iosys_intermediate_classes(cls):
-    docstring = inspect.getdoc(cls)
-
-    # Make sure there is not a parameters section
-    if re.search(r"\nParameters\n----", docstring) is not None:
-        pytest.fail(f"intermediate {cls} docstring contains Parameters section")
-
-    # Make sure we reference the factory function
-    fcn = class_factory_function[cls]
-    if re.search(f":func:`~control.{fcn.__name__}`", docstring) is None:
-        pytest.fail(
-            f"{cls.__name__} does not reference factory function "
-            f"{fcn.__name__}")
-
-
-@pytest.mark.parametrize("fcn", factory_args.keys())
-def test_iosys_factory_functions(fcn):
-    docstring = inspect.getdoc(fcn)
-    cls = list(class_factory_function.keys())[
-        list(class_factory_function.values()).index(fcn)]
-
-    # Make sure we reference parameters in class and factory function docstring
-    for argname in class_args[cls] + std_factory_args + factory_args[fcn]:
-        _check_parameter_docs(fcn.__name__, argname, docstring)
-
-    # Make sure we don't reference any class attributes
-    for argname in std_class_attributes + class_attributes[cls]:
-        if argname in std_factory_args:
-            continue
-        if re.search(f"[\\s]+{argname}(, .*)*[\\s]*:", docstring) is not None:
-            pytest.fail(
-                f"{fcn.__name__} references class attribute '{argname}'")
-
 
 # Utility function to check for an argument in a docstring
-def _check_parameter_docs(
-        funcname, argname, docstring, prefix="", fail_if_missing=True):
+def _check_parameter_docs(funcname, argname, docstring, prefix=""):
     funcname = prefix + funcname
 
     # Find the "Parameters" section of docstring, where we start searching
-    if not (match := re.search(
-            r"\nParameters\n----", docstring)):
-        if fail_if_missing:
-            pytest.fail(f"{funcname} docstring missing Parameters section")
-        else:
-            return False
+    if not (match := re.search(r"\nParameters\n----", docstring)):
+        pytest.fail(f"{funcname} docstring missing Parameters section")
     else:
         start = match.start()
 
@@ -529,10 +299,7 @@ def _check_parameter_docs(
             docstring)):
         if verbose:
             print(f"      {funcname}: {argname} not documented")
-        if fail_if_missing:
-            pytest.fail(f"{funcname} '{argname}' not documented")
-        else:
-            return False
+        pytest.fail(f"{funcname} '{argname}' not documented")
 
     # Make sure there isn't another instance
     second_match = re.search(
@@ -540,5 +307,3 @@ def _check_parameter_docs(
             docstring[match.end():])
     if second_match:
         pytest.fail(f"{funcname} '{argname}' documented twice")
-
-    return True
diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index b08cd826..c2a29ee2 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -482,7 +482,7 @@ class TestFRD:
         sys1r = eval(repr(sys1))
         np.testing.assert_array_almost_equal(sys1r.fresp, sys1.fresp)
         np.testing.assert_array_almost_equal(sys1r.omega, sys1.omega)
-        assert(sys1._ifunc is not None)
+        assert(sys1.ifunc is not None)
 
         refs = """<FrequencyResponseData>: {sysname}
 Inputs (1): ['u[0]']
diff --git a/control/tests/interconnect_test.py b/control/tests/interconnect_test.py
index d124859f..604488ca 100644
--- a/control/tests/interconnect_test.py
+++ b/control/tests/interconnect_test.py
@@ -666,29 +666,15 @@ def test_interconnect_params():
     # Create a nominally unstable system
     sys1 = ct.nlsys(
         lambda t, x, u, params: params['a'] * x[0] + u[0],
-        states=1, inputs='u', outputs='y', params={'a': 2, 'c':2})
+        states=1, inputs='u', outputs='y', params={'a': 1})
 
     # Simple system for serial interconnection
     sys2 = ct.nlsys(
         None, lambda t, x, u, params: u[0],
-        inputs='r', outputs='u', params={'a': 4, 'b': 3})
+        inputs='r', outputs='u')
 
-    # Make sure default parameters get set as expected
+    # Create a series interconnection
     sys = ct.interconnect([sys1, sys2], inputs='r', outputs='y')
-    assert sys.params == {'a': 4, 'c': 2, 'b': 3}
-    assert sys.dynamics(0, [1], [0]).item() == 4
-
-    # Make sure we can override the parameters
-    sys = ct.interconnect(
-        [sys1, sys2], inputs='r', outputs='y', params={'b': 1})
-    assert sys.params == {'b': 1}
-    assert sys.dynamics(0, [1], [0]).item() == 2
-    assert sys.dynamics(0, [1], [0], params={'a': 5}).item() == 5
-
-    # Create final series interconnection, with proper parameter values
-    sys = ct.interconnect(
-        [sys1, sys2], inputs='r', outputs='y', params={'a': 1})
-    assert sys.params == {'a': 1}
 
     # Make sure we can call the update function
     sys.updfcn(0, [0], [0], {})
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 54d6d56c..0216235b 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -17,8 +17,7 @@ import pytest
 import scipy
 
 import control as ct
-import control.flatsys as fs
-from control.tests.conftest import slycotonly
+
 
 class TestIOSys:
 
@@ -931,8 +930,6 @@ class TestIOSys:
         ios_secord_update = ct.NonlinearIOSystem(
             secord_update, secord_output, inputs=1, outputs=1, states=2,
             params={'omega0':2, 'zeta':0})
-        lin_secord_update = ct.linearize(ios_secord_update, [0, 0], [0])
-        w_update, v_update = np.linalg.eig(lin_secord_update.A)
 
         # Make sure the default parameters haven't changed
         lin_secord_check = ct.linearize(ios_secord_default, [0, 0], [0])
@@ -962,7 +959,7 @@ class TestIOSys:
             ios_series_default_local, [0, 0, 0, 0], [0])
         w, v = np.linalg.eig(lin_series_default_local.A)
         np.testing.assert_array_almost_equal(
-             w, np.concatenate([w_update, w_update]))
+            np.sort(w), np.sort(np.concatenate((w_default, [2j, -2j]))))
 
         # Show that we can change the parameters at linearization
         lin_series_override = ct.linearize(
@@ -2287,55 +2284,3 @@ def test_signal_indexing():
     with pytest.raises(IndexError, match=r"signal name\(s\) not valid"):
         resp.outputs['y[0]', 'u[0]']
 
-@pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.frd, ct.nlsys, fs.flatsys])
-def test_relabeling(fcn):
-    sys = ct.rss(1, 1, 1, name="sys")
-
-    # Rename the inputs, outputs, (states,) system
-    match fcn:
-        case ct.tf:
-            sys = fcn(sys, inputs='u', outputs='y', name='new')
-        case ct.frd:
-            sys = fcn(sys, [0.1, 1, 10], inputs='u', outputs='y', name='new')
-        case _:
-            sys = fcn(sys, inputs='u', outputs='y', states='x', name='new')
-
-    assert sys.input_labels == ['u']
-    assert sys.output_labels == ['y']
-    if sys.nstates:
-        assert sys.state_labels == ['x']
-    assert sys.name == 'new'
-
-
-@pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.frd, ct.nlsys, fs.flatsys])
-def test_signal_prefixing(fcn):
-    sys = ct.rss(2, 1, 1)
-
-    # Recreate the system in different forms, with non-standard prefixes
-    match fcn:
-        case ct.ss:
-            sys = ct.ss(
-                sys.A, sys.B, sys.C, sys.D, state_prefix='xx',
-                input_prefix='uu', output_prefix='yy')
-        case ct.tf:
-            sys = ct.tf(sys)
-            sys = fcn(sys.num, sys.den, input_prefix='uu', output_prefix='yy')
-        case ct.frd:
-            freq = [0.1, 1, 10]
-            data = [sys(w * 1j) for w in freq]
-            sys = fcn(data, freq, input_prefix='uu', output_prefix='yy')
-        case ct.nlsys:
-            sys = ct.nlsys(sys)
-            sys = fcn(
-                sys.updfcn, sys.outfcn, inputs=1, outputs=1, states=2,
-                state_prefix='xx', input_prefix='uu', output_prefix='yy')
-        case fs.flatsys:
-            sys = fs.flatsys(sys)
-            sys = fcn(
-                sys.forward, sys.reverse, inputs=1, outputs=1, states=2,
-                state_prefix='xx', input_prefix='uu', output_prefix='yy')
-
-    assert sys.input_labels == ['uu[0]']
-    assert sys.output_labels == ['yy[0]']
-    if sys.nstates:
-        assert sys.state_labels == ['xx[0]', 'xx[1]']
diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
index b98509d6..95450da0 100644
--- a/control/tests/kwargs_test.py
+++ b/control/tests/kwargs_test.py
@@ -244,6 +244,7 @@ kwarg_unittest = {
     'append': test_unrecognized_kwargs,
     'bode': test_response_plot_kwargs,
     'bode_plot': test_response_plot_kwargs,
+    'LTI.bode_plot': test_response_plot_kwargs, # alias for bode_plot and tested via bode_plot
     'create_estimator_iosystem': stochsys_test.test_estimator_errors,
     'create_statefbk_iosystem': statefbk_test.TestStatefbk.test_statefbk_errors,
     'describing_function_plot': test_matplotlib_kwargs,
@@ -267,11 +268,13 @@ kwarg_unittest = {
     'lqr': test_unrecognized_kwargs,
     'negate': test_unrecognized_kwargs,
     'nichols_plot': test_matplotlib_kwargs,
+    'LTI.nichols_plot': test_matplotlib_kwargs, # alias for nichols_plot and tested via nichols_plot
     'nichols': test_matplotlib_kwargs,
     'nlsys': test_unrecognized_kwargs,
     'nyquist': test_matplotlib_kwargs,
     'nyquist_response': test_response_plot_kwargs,
     'nyquist_plot': test_matplotlib_kwargs,
+    'LTI.nyquist_plot': test_matplotlib_kwargs, # alias for nyquist_plot and tested via nyquist_plot
     'phase_plane_plot': test_matplotlib_kwargs,
     'parallel': test_unrecognized_kwargs,
     'pole_zero_plot': test_unrecognized_kwargs,
diff --git a/control/tests/nlsys_test.py b/control/tests/nlsys_test.py
index 926ca436..7f649e0c 100644
--- a/control/tests/nlsys_test.py
+++ b/control/tests/nlsys_test.py
@@ -7,15 +7,11 @@ The main test functions are contained in iosys_test.py.
 
 """
 
-import math
-import re
-
-import numpy as np
 import pytest
-
+import numpy as np
+import math
 import control as ct
 
-
 # Basic test of nlsys()
 def test_nlsys_basic():
     def kincar_update(t, x, u, params):
@@ -158,44 +154,3 @@ def test_nlsys_empty_io():
 
     resp = ct.forced_response(P, np.linspace(0, 1), 1)
     np.testing.assert_allclose(resp.states[:, -1], 1 - math.exp(-1))
-
-
-def test_ss2io():
-    sys = ct.rss(
-        states=4, inputs=['u1', 'u2'], outputs=['y1', 'y2'], name='sys')
-
-    # Standard conversion
-    nlsys = ct.nlsys(sys)
-    for attr in ['nstates', 'ninputs', 'noutputs']:
-        assert getattr(nlsys, attr) == getattr(sys, attr)
-    assert nlsys.name == 'sys$converted'
-    np.testing.assert_allclose(
-        nlsys.dynamics(0, [1, 2, 3, 4], [0, 0], {}),
-        sys.A @ np.array([1, 2, 3, 4]))
-
-    # Put names back to defaults
-    nlsys = ct.nlsys(
-        sys, inputs=sys.ninputs, outputs=sys.noutputs, states=sys.nstates)
-    for attr, prefix in zip(
-            ['state_labels', 'input_labels', 'output_labels'],
-            ['x', 'u', 'y']):
-        for i in range(len(getattr(nlsys, attr))):
-            assert getattr(nlsys, attr)[i] == f"{prefix}[{i}]"
-    assert re.match(r"sys\$converted", nlsys.name)
-
-    # Override the names with something new
-    nlsys = ct.nlsys(
-        sys, inputs=['U1', 'U2'], outputs=['Y1', 'Y2'],
-        states=['X1', 'X2', 'X3', 'X4'], name='nlsys')
-    for attr, prefix in zip(
-            ['state_labels', 'input_labels', 'output_labels'],
-            ['X', 'U', 'Y']):
-        for i in range(len(getattr(nlsys, attr))):
-            assert getattr(nlsys, attr)[i] == f"{prefix}{i+1}"
-    assert nlsys.name == 'nlsys'
-
-    # Make sure dimension checking works
-    for attr in ['states', 'inputs', 'outputs']:
-        with pytest.raises(ValueError, match=r"new .* doesn't match"):
-            kwargs = {attr: getattr(sys, 'n' + attr) - 1}
-            nlsys = ct.nlsys(sys, **kwargs)
diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py
index 22a946fe..3928fb72 100644
--- a/control/tests/statefbk_test.py
+++ b/control/tests/statefbk_test.py
@@ -930,10 +930,9 @@ def unicycle():
 
     return ct.NonlinearIOSystem(
         unicycle_update, None,
-        inputs=['v', 'phi'],
-        outputs=['x', 'y', 'theta'],
-        states=['x_', 'y_', 'theta_'],
-        params={'a': 1})        # only used for testing params
+        inputs = ['v', 'phi'],
+        outputs = ['x', 'y', 'theta'],
+        states = ['x_', 'y_', 'theta_'])
 
 from math import pi
 
@@ -1195,40 +1194,3 @@ def test_create_statefbk_errors():
 
     with pytest.raises(ControlArgument, match="feedfwd_pattern != 'refgain'"):
         ct.create_statefbk_iosystem(sys, K, Kf, feedfwd_pattern='trajgen')
-
-
-def test_create_statefbk_params(unicycle):
-    # Speeds and angles at which to compute the gains
-    speeds = [1, 5, 10]
-    angles = np.linspace(0, pi/2, 4)
-    points = list(itertools.product(speeds, angles))
-
-    # Gains for each speed (using LQR controller)
-    Q = np.identity(unicycle.nstates)
-    R = np.identity(unicycle.ninputs)
-    gain, _, _ = ct.lqr(unicycle.linearize([0, 0, 0], [5, 0]), Q, R)
-
-    #
-    # Schedule on desired speed and angle
-    #
-
-    # Create a linear controller
-    ctrl, clsys = ct.create_statefbk_iosystem(unicycle, gain)
-    assert [k for k in ctrl.params.keys()] == []
-    assert [k for k in clsys.params.keys()] == ['a']
-    assert clsys.params['a'] == 1
-
-    # Create a nonlinear controller
-    ctrl, clsys = ct.create_statefbk_iosystem(
-        unicycle, gain, controller_type='nonlinear')
-    assert [k for k in ctrl.params.keys()] == ['K']
-    assert [k for k in clsys.params.keys()] == ['K', 'a']
-    assert clsys.params['a'] == 1
-
-    # Override the default parameters
-    ctrl, clsys = ct.create_statefbk_iosystem(
-        unicycle, gain, controller_type='nonlinear', params={'a': 2, 'b': 1})
-    assert [k for k in ctrl.params.keys()] == ['K']
-    assert [k for k in clsys.params.keys()] == ['K', 'a', 'b']
-    assert clsys.params['a'] == 2
-    assert clsys.params['b'] == 1
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 2bb0badc..43dfbbbd 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -127,19 +127,19 @@ class TestStateSpace:
          ((1, 2), TypeError, "1, 4, or 5 arguments"),
          ((np.ones((3, 2)), np.ones((3, 2)),
            np.ones((2, 2)), np.ones((2, 2))), ValueError,
-          r"A is the wrong shape; expected \(3, 3\)"),
+          "A is the wrong shape; expected \(3, 3\)"),
          ((np.ones((3, 3)), np.ones((2, 2)),
            np.ones((2, 3)), np.ones((2, 2))), ValueError,
-          r"B is the wrong shape; expected \(3, 2\)"),
+          "B is the wrong shape; expected \(3, 2\)"),
          ((np.ones((3, 3)), np.ones((3, 2)),
            np.ones((2, 2)), np.ones((2, 2))), ValueError,
-          r"C is the wrong shape; expected \(2, 3\)"),
+          "C is the wrong shape; expected \(2, 3\)"),
          ((np.ones((3, 3)), np.ones((3, 2)),
            np.ones((2, 3)), np.ones((2, 3))), ValueError,
-          r"D is the wrong shape; expected \(2, 2\)"),
+          "D is the wrong shape; expected \(2, 2\)"),
          ((np.ones((3, 3)), np.ones((3, 2)),
            np.ones((2, 3)), np.ones((3, 2))), ValueError,
-          r"D is the wrong shape; expected \(2, 2\)"),
+          "D is the wrong shape; expected \(2, 2\)"),
         ])
     def test_constructor_invalid(self, args, exc, errmsg):
         """Test invalid input to StateSpace() constructor"""
@@ -1281,3 +1281,31 @@ def test_tf2ss_mimo():
     else:
         with pytest.raises(ct.ControlMIMONotImplemented):
             sys_ss = ct.ss(sys_tf)
+
+def test_convenience_aliases():
+    sys = ct.StateSpace(1, 1, 1, 1)
+
+    # test that all the aliases point to the correct function
+    # .__funct__ a bound methods underlying function
+    assert sys.to_ss.__func__ == ct.ss
+    assert sys.to_tf.__func__ == ct.tf
+    assert sys.bode_plot.__func__ == ct.bode_plot
+    assert sys.nyquist_plot.__func__ == ct.nyquist_plot
+    assert sys.nichols_plot.__func__ == ct.nichols_plot
+    assert sys.forced_response.__func__ == ct.forced_response
+    assert sys.impulse_response.__func__ == ct.impulse_response
+    assert sys.step_response.__func__ == ct.step_response
+    assert sys.initial_response.__func__ == ct.initial_response
+
+    # make sure the functions can be used as member function ie they support
+    # an instance of StateSpace as the first argument and that they at least return
+    # the correct type
+    assert isinstance(sys.to_ss(), StateSpace)
+    assert isinstance(sys.to_tf(), TransferFunction)
+    assert isinstance(sys.bode_plot(), ct.ControlPlot)
+    assert isinstance(sys.nyquist_plot(), ct.ControlPlot)
+    assert isinstance(sys.nichols_plot(), ct.ControlPlot)
+    assert isinstance(sys.forced_response([0, 1], [1, 1]), (ct.TimeResponseData, ct.TimeResponseList))
+    assert isinstance(sys.impulse_response(), (ct.TimeResponseData, ct.TimeResponseList))
+    assert isinstance(sys.step_response(), (ct.TimeResponseData, ct.TimeResponseList))
+    assert isinstance(sys.initial_response(X0=1), (ct.TimeResponseData, ct.TimeResponseList))
\ No newline at end of file
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index ae28975e..e5e24b99 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -1247,14 +1247,13 @@ def test_to_pandas():
     np.testing.assert_equal(df['x[1]'], resp.states[1])
 
     # Change the time points
-    sys = ct.rss(2, 1, 2)
+    sys = ct.rss(2, 1, 1)
     T = np.linspace(0, timepts[-1]/2, timepts.size * 2)
-    resp = ct.input_output_response(
-        sys, timepts, [np.sin(timepts), 0], t_eval=T)
+    resp = ct.input_output_response(sys, timepts, np.sin(timepts), t_eval=T)
     df = resp.to_pandas()
     np.testing.assert_equal(df['time'], resp.time)
-    np.testing.assert_equal(df['u[0]'], resp.inputs[0])
-    np.testing.assert_equal(df['y[0]'], resp.outputs[0])
+    np.testing.assert_equal(df['u[0]'], resp.inputs)
+    np.testing.assert_equal(df['y[0]'], resp.outputs)
     np.testing.assert_equal(df['x[0]'], resp.states[0])
     np.testing.assert_equal(df['x[1]'], resp.states[1])
 
@@ -1266,33 +1265,6 @@ def test_to_pandas():
     np.testing.assert_equal(df['u[0]'], resp.inputs)
     np.testing.assert_equal(df['y[0]'], resp.inputs * 5)
 
-    # Multi-trace data
-    # https://github.com/python-control/python-control/issues/1087
-    model = ct.rss(
-        states=['x0', 'x1'], outputs=['y0', 'y1'],
-        inputs=['u0', 'u1'], name='My Model')
-    T = np.linspace(0, 10, 100, endpoint=False)
-    X0 = np.zeros(model.nstates)
-
-    res = ct.step_response(model, T=T, X0=X0, input=0)  # extract single trace
-    df = res.to_pandas()
-    np.testing.assert_equal(
-        df[df['trace'] == 'From u0']['time'], res.time)
-    np.testing.assert_equal(
-        df[df['trace'] == 'From u0']['u0'], res.inputs['u0', 0])
-    np.testing.assert_equal(
-        df[df['trace'] == 'From u0']['y1'], res.outputs['y1', 0])
-
-    res = ct.step_response(model, T=T, X0=X0)           # all traces
-    df = res.to_pandas()
-    for i, label in enumerate(res.trace_labels):
-        np.testing.assert_equal(
-            df[df['trace'] == label]['time'], res.time)
-        np.testing.assert_equal(
-            df[df['trace'] == label]['u1'], res.inputs['u1', i])
-        np.testing.assert_equal(
-            df[df['trace'] == label]['y0'], res.outputs['y0', i])
-
 
 @pytest.mark.skipif(pandas_check(), reason="pandas installed")
 def test_no_pandas():
diff --git a/control/tests/xferfcn_input_test.py b/control/tests/xferfcn_input_test.py
index c375d768..46efbd25 100644
--- a/control/tests/xferfcn_input_test.py
+++ b/control/tests/xferfcn_input_test.py
@@ -64,18 +64,15 @@ def test_clean_part(num, fun, dtype):
     num_ = _clean_part(numa)
     ref_ = np.array(num, dtype=float, ndmin=3)
 
-    assert isinstance(num_, np.ndarray)
-    assert num_.ndim == 2
+    assert isinstance(num_, list)
+    assert np.all([isinstance(part, list) for part in num_])
     for i, numi in enumerate(num_):
         assert len(numi) == ref_.shape[1]
         for j, numj in enumerate(numi):
             np.testing.assert_allclose(numj, ref_[i, j, ...])
 
 
-@pytest.mark.parametrize("badinput", [
-    # [[0., 1.], [2., 3.]],             # OK: treated as static array
-    np.ones((2, 2, 2, 2)),
-    "a"])
+@pytest.mark.parametrize("badinput", [[[0., 1.], [2., 3.]], "a"])
 def test_clean_part_bad_input(badinput):
     """Give the part cleaner invalid input type."""
     with pytest.raises(TypeError):
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index 3f87ef1d..d480cef6 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -30,6 +30,13 @@ class TestXferFcn:
 
     def test_constructor_bad_input_type(self):
         """Give the constructor invalid input types."""
+        # MIMO requires lists of lists of vectors (not lists of vectors)
+        with pytest.raises(TypeError):
+            TransferFunction([[0., 1.], [2., 3.]], [[5., 2.], [3., 0.]])
+        # good input
+        TransferFunction([[[0., 1.], [2., 3.]]],
+                         [[[5., 2.], [3., 0.]]])
+
         # Single argument of the wrong type
         with pytest.raises(TypeError):
             TransferFunction([1])
@@ -183,10 +190,8 @@ class TestXferFcn:
 
         for i in range(sys3.noutputs):
             for j in range(sys3.ninputs):
-                np.testing.assert_allclose(
-                    sys2.num_array[i, j], sys3.num_array[i, j])
-                np.testing.assert_allclose(
-                    sys2.den_array[i, j], sys3.den_array[i, j])
+                np.testing.assert_allclose(sys2.num[i][j], sys3.num[i][j])
+                np.testing.assert_allclose(sys2.den[i][j], sys3.den[i][j])
 
     # Tests for TransferFunction.__add__
 
@@ -231,8 +236,8 @@ class TestXferFcn:
 
         for i in range(sys3.noutputs):
             for j in range(sys3.ninputs):
-                np.testing.assert_allclose(sys3.num_array[i, j], num3[i][j])
-                np.testing.assert_allclose(sys3.den_array[i, j], den3[i][j])
+                np.testing.assert_allclose(sys3.num[i][j], num3[i][j])
+                np.testing.assert_allclose(sys3.den[i][j], den3[i][j])
 
     # Tests for TransferFunction.__sub__
 
@@ -279,8 +284,8 @@ class TestXferFcn:
 
         for i in range(sys3.noutputs):
             for j in range(sys3.ninputs):
-                np.testing.assert_allclose(sys3.num_array[i, j], num3[i][j])
-                np.testing.assert_allclose(sys3.den_array[i, j], den3[i][j])
+                np.testing.assert_allclose(sys3.num[i][j], num3[i][j])
+                np.testing.assert_allclose(sys3.den[i][j], den3[i][j])
 
     # Tests for TransferFunction.__mul__
 
@@ -335,8 +340,8 @@ class TestXferFcn:
 
         for i in range(sys3.noutputs):
             for j in range(sys3.ninputs):
-                np.testing.assert_allclose(sys3.num_array[i, j], num3[i][j])
-                np.testing.assert_allclose(sys3.den_array[i, j], den3[i][j])
+                np.testing.assert_allclose(sys3.num[i][j], num3[i][j])
+                np.testing.assert_allclose(sys3.den[i][j], den3[i][j])
 
     # Tests for TransferFunction.__div__
 
@@ -657,10 +662,10 @@ class TestXferFcn:
 
         for i in range(sys.noutputs):
             for j in range(sys.ninputs):
-                np.testing.assert_array_almost_equal(
-                    tfsys.num_array[i, j], num[i][j])
-                np.testing.assert_array_almost_equal(
-                    tfsys.den_array[i, j], den[i][j])
+                np.testing.assert_array_almost_equal(tfsys.num[i][j],
+                                                     num[i][j])
+                np.testing.assert_array_almost_equal(tfsys.den[i][j],
+                                                     den[i][j])
 
     def test_minreal(self):
         """Try the minreal function, and also test easy entry by creation
@@ -1116,10 +1121,8 @@ class TestXferFcn:
         H2 = eval(H.__repr__())
         for p in range(len(H.num)):
             for m in range(len(H.num[0])):
-                np.testing.assert_array_almost_equal(
-                    H.num_array[p, m], H2.num_array[p, m])
-                np.testing.assert_array_almost_equal(
-                    H.den_array[p, m], H2.den_array[p, m])
+                np.testing.assert_array_almost_equal(H.num[p][m], H2.num[p][m])
+                np.testing.assert_array_almost_equal(H.den[p][m], H2.den[p][m])
             assert H.dt == H2.dt
 
     def test_sample_named_signals(self):
@@ -1177,10 +1180,8 @@ class TestLTIConverter:
         sslti = mimotf.returnScipySignalLTI(strict=False)
         for i in range(2):
             for j in range(3):
-                np.testing.assert_allclose(
-                    sslti[i][j].num, mimotf.num_array[i, j])
-                np.testing.assert_allclose(
-                    sslti[i][j].den, mimotf.den_array[i, j])
+                np.testing.assert_allclose(sslti[i][j].num, mimotf.num[i][j])
+                np.testing.assert_allclose(sslti[i][j].den, mimotf.den[i][j])
                 if mimotf.dt == 0:
                     assert sslti[i][j].dt is None
                 else:
@@ -1285,35 +1286,3 @@ def test_copy_names(create, args, kwargs, convert):
     cpy = convert(sys, inputs='myin', outputs='myout')
     assert cpy.input_labels == ['myin']
     assert cpy.output_labels == ['myout']
-
-s = ct.TransferFunction.s
-@pytest.mark.parametrize("args, num, den", [
-    (('s', ), [[[1, 0]]], [[[1]]]),                     # ctime
-    (('z', ), [[[1, 0]]], [[[1]]]),                     # dtime
-    ((1, 1), [[[1]]], [[[1]]]),                         # scalars as scalars
-    (([[1]], [[1]]), [[[1]]], [[[1]]]),                 # scalars as lists
-    (([[[1, 2]]], [[[3, 4]]]), [[[1, 2]]], [[[3, 4]]]), # SISO as lists
-    (([[np.array([1, 2])]], [[np.array([3, 4])]]),      # SISO as arrays
-     [[[1, 2]]], [[[3, 4]]]),
-    (([[    [1],    [2] ], [[1, 1],      [1, 0] ]],     # MIMO
-      [[ [1, 0], [1, 0] ], [[1, 2],         [1] ]]),
-     [[    [1],    [2] ], [[1, 1],      [1, 0] ]],
-     [[ [1, 0], [1, 0] ], [[1, 2],         [1] ]]),
-    (([[[1, 2], [3, 4]]], [[[5, 6]]]),                  # common denominator
-     [[[1, 2], [3, 4]]], [[[5, 6], [5, 6]]]),
-    (([   [1/s,   2/s],   [(s+1)/(s+2),     s]], ),     # 2x2 from SISO
-     [[    [1],    [2] ], [[1, 1],      [1, 0] ]],      # num
-     [[ [1, 0], [1, 0] ], [[1, 2],         [1] ]]),     # den
-    (([[1, 2], [3, 4]], [[[1, 0], [1, 0]]]), ValueError,
-     r"numerator has 2 output\(s\), but the denominator has 1 output"),
-])
-def test_tf_args(args, num, den):
-    if isinstance(num, type):
-        exception, match = num, den
-        with pytest.raises(exception, match=match):
-            sys = ct.tf(*args)
-    else:
-        sys = ct.tf(*args)
-        chk = ct.tf(num, den)
-        np.testing.assert_equal(sys.num, chk.num)
-        np.testing.assert_equal(sys.den, chk.den)