Repository: python-control/python-control
Pull request #1082 (2024-12-20 21:07:27)

I/O system enhancements
This PR provides a number of enhancements to I/O systems classes, most of which are internal but a few of which are visible to the user:

* Updated factory functions to provide consistent functionality in terms of copying signal/system names, overriding system/signal names, and converting between classes.
* Updated docstrings for all I/O classes and factory functions to provide consistent information and added unit tests to ensure consistency.  The main changes were to make sure that all factory function parameters are documented while having the class docstrings focus on the primary parameters + attributes (with a reference to the factory function for more detailed documentation).
* Updated the `tf()` factory function to allow a 2D list of SISO transfer functions to be given as a means of creating a MIMO transfer function (use the new `combine_tf` function).  This is redundant with `combine_tf`, but seemed useful to include in the factory function for the special case of creating a MIMO transfer function via a set of SISO expressions.  E.g.:
  ```
  s = ct.tf('s')
  sys = ct.tf([[1/s, 0], [(s+1)/(s**2 + 2s + 1), 1/(s+1)]])
  ```
* [8f3615b, 22 Dec '24] Added processing of the `params` keyword to `create_statefbk_iosystem`.
* [e4d373c, 22 Dec '24] Updated `interconnect()` to combine parameter lists of subsystems if `params` is not specified.
* Changed some methods and attributes to be hidden so that they don't get flagged in docstring testing.  These are functions that are not needed by the user (eg, `FrequencyResponseData._ifunc`, which is the internally generated interpolation function if the `smooth` keyword is used).
* Changed the internal representation of transfer function numerators and denominators to a 2D array of arrays instead of a 2D nested list of arrays.  This allows user of NumPy routines when accessing input/output pairs (eg, `sys.num_array[i, j]` instead of `sys.num_list[i][j]`).  The nested list representation is still available via the `{num,den}_list` properties and the `num` and `den` attributes retrieve the nested list structure for backward compatibility.  (We may want to change that in the future, but it could break user code where someone is accessing the numerator and denominator polynomials directly.)
* Add the ability for `nlsys` to create a `NonlinearIOSystem` representation of a `StateSpace` system.  It's not clear that this is so useful for a user (since `StateSpace` is a subclass of `NonlinearIOSystem`), but it is handy for testing since you can quickly create a nonlinear I/O system from a linear system but have it treated as a "pure" `NonlinearIOSystem` instead of a `StateSpace` system.
* Got rid of redundant licensing information at the top of files that were modified (the information is contained in the `LICENSE` file in the main directory).

Note: this PR is likely to conflict with #1081, so the plan is to review/merge #1081 first, then rebase this PR on top of that one.


Patch:
diff --git a/control/bdalg.py b/control/bdalg.py
index d907cd3c..2dbd5c8e 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[j_out][j_in])
-                    den_row.append(col.den[j_out][j_in])
+                    num_row.append(col.num_array[j_out, j_in])
+                    den_row.append(col.den_array[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[i_out][i_in],
-                    transfer_function.den[i_out][i_in],
+                    transfer_function.num_array[i_out, i_in],
+                    transfer_function.den_array[i_out, i_in],
                     dt=transfer_function.dt,
                 )
             )
diff --git a/control/flatsys/flatsys.py b/control/flatsys/flatsys.py
index 0101d126..7d76b9d7 100644
--- a/control/flatsys/flatsys.py
+++ b/control/flatsys/flatsys.py
@@ -1,42 +1,5 @@
 # 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
@@ -54,8 +17,30 @@ 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 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.
 
     Notes
     -----
@@ -234,10 +219,9 @@ 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
@@ -252,6 +236,12 @@ 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 ac032d3f..64a1e822 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.
+"""Frequency response data representation and functions.
+
+This module contains the FrequencyResponseData (FRD) class and also
+functions that operate on FRD data.
 
-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, _process_iosys_keywords, \
-    _process_subsys_index, common_timebase
+from .iosys import InputOutputSystem, NamedSignal, _extended_system_name, \
+    _process_iosys_keywords, _process_subsys_index, common_timebase
 from .lti import LTI, _process_frequency_response
 
 __all__ = ['FrequencyResponseData', 'FRD', 'frd']
@@ -35,38 +35,29 @@ 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
     ----------
-    d : 1D or 3D complex array_like
+    response : 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
-    w : iterable of real frequencies
+    omega : 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``.
-
-    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.
+    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).
     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
@@ -79,16 +70,46 @@ 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'].
-    ninputs, noutputs, nstates : int
-        Number of inputs, outputs, and states of the underlying system.
+    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).
     input_labels, output_labels : array of str
-        Names for the input and output variables.
-    sysname : str, optional
-        Name of the system.  For data generated using
-        :func:`~control.frequency_response`, stores the name of the system
-        that created the data.
+        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').
     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
     --------
@@ -148,22 +169,26 @@ class FrequencyResponseData(LTI):
     _epsw = 1e-8                #: Bound for exact frequency match
 
     def __init__(self, *args, **kwargs):
-        """Construct an FRD object.
+        """FrequencyResponseData(d, w[, dt])
 
-        The default constructor is FRD(d, w), where w is an iterable of
-        frequency points, and d is the matching frequency data.
+        Construct a frequency response data (FRD) object.
 
-        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.
+        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.
 
-        To construct frequency response data for an existing LTI
-        object, other than an FRD, call FRD(sys, omega).
+        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).
 
-        The timebase for the frequency response can be provided using an
-        optional third argument or the 'dt' keyword.
+        See :class:`FrequencyResponseData` and :func:`frd` for more
+        information.
 
         """
         smooth = kwargs.pop('smooth', False)
@@ -182,11 +207,12 @@ class FrequencyResponseData(LTI):
 
         if len(args) == 2:
             if not isinstance(args[0], FRD) and isinstance(args[0], LTI):
-                # not an FRD, but still a system, second argument should be
-                # the frequency range
+                # not an FRD, but still an LTI system, second argument
+                # should be the frequency range
                 otherlti = args[0]
                 self.omega = sort(np.asarray(args[1], dtype=float))
-                # calculate frequency response at my points
+
+                # calculate frequency response at specified points
                 if otherlti.isctime():
                     s = 1j * self.omega
                     self.fresp = otherlti(s, squeeze=False)
@@ -195,6 +221,14 @@ 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)
@@ -219,6 +253,10 @@ 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))
@@ -249,15 +287,21 @@ class FrequencyResponseData(LTI):
 
         # Process iosys keywords
         defaults = {
-            'inputs': self.fresp.shape[1], 'outputs': self.fresp.shape[0]}
+            '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)}
         if arg_dt is not None:
-            defaults['dt'] = arg_dt             # choose compatible timebase
-        name, inputs, outputs, states, dt = _process_iosys_keywords(
-                kwargs, defaults, end=True)
+            if isinstance(args[0], LTI):
+                arg_dt = common_timebase(args[0].dt, arg_dt)
+            kwargs['dt'] = arg_dt
 
         # Process signal names
+        name, inputs, outputs, states, dt = _process_iosys_keywords(
+                kwargs, defaults)
         InputOutputSystem.__init__(
-            self, name=name, inputs=inputs, outputs=outputs, dt=dt)
+            self, name=name, inputs=inputs, outputs=outputs, dt=dt, **kwargs)
 
         # create interpolation functions
         if smooth:
@@ -266,17 +310,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
@@ -379,7 +423,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."""
@@ -438,7 +482,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)
 
@@ -456,8 +500,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)."""
@@ -465,7 +509,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)
 
@@ -484,8 +528,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):
@@ -493,7 +537,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)
 
@@ -504,15 +548,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)
 
@@ -529,7 +573,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:
@@ -582,7 +626,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(
@@ -596,7 +640,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)
@@ -751,7 +795,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):
@@ -901,6 +945,8 @@ 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.
@@ -917,7 +963,7 @@ def frd(*args, **kwargs):
 
     Returns
     -------
-    sys : :class:`FrequencyResponseData`
+    sys : FrequencyResponseData
         New frequency response data system.
 
     Other Parameters
@@ -926,6 +972,8 @@ 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 54442529..8d3e6468 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 dd1566eb..2c1f9cea 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -145,20 +145,14 @@ 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
-    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)
+        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).
 
     Other Parameters
     ----------------
@@ -194,7 +188,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']
 
@@ -293,13 +287,8 @@ 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
-        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
+        self.name = _extended_system_name(
+            sys.name, prefix, suffix, prefix_suffix_name)
 
         # Name the inputs, outputs, and states
         self.input_index = sys.input_index.copy()
@@ -432,6 +421,11 @@ 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])
@@ -1057,3 +1051,14 @@ 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 e9455aed..4b354f3f 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -458,7 +458,7 @@ def frequency_response(
     Examples
     --------
     >>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]])
-    >>> mag, phase, omega = ct.freqresp(G, [0.1, 1., 10.])
+    >>> mag, phase, omega = ct.frequency_response(G, [0.1, 1., 10.])
 
     .. todo::
         Add example with MIMO system
@@ -490,8 +490,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 62e4bf78..beb2566e 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -39,10 +39,11 @@ __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 (Note: discrete-time systems
-    are not yet supported by most functions.)
+    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.
 
     Parameters
     ----------
@@ -63,20 +64,13 @@ class NonlinearIOSystem(InputOutputSystem):
 
         where the arguments are the same as for `upfcn`.
 
-    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.
+    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.
 
-    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`.
+    params : dict, optional
+        Parameter values for the systems.  Passed to the evaluation functions
+        for the system as default values, overriding internal defaults.
 
     dt : timebase, optional
         The timebase for the system, used to specify whether the system is
@@ -88,13 +82,16 @@ 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 (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.
+        System name.
 
     See Also
     --------
@@ -713,6 +710,11 @@ 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)
@@ -1220,8 +1222,7 @@ class InterconnectedSystem(NonlinearIOSystem):
         return dropped_inputs, dropped_outputs
 
 
-def nlsys(
-        updfcn, outfcn=None, inputs=None, outputs=None, states=None, **kwargs):
+def nlsys(updfcn, outfcn=None, **kwargs):
     """Create a nonlinear input/output system.
 
     Creates an :class:`~control.InputOutputSystem` for a nonlinear system by
@@ -1230,7 +1231,7 @@ def nlsys(
 
     Parameters
     ----------
-    updfcn : callable
+    updfcn : callable (or StateSpace)
         Function returning the state update function
 
             `updfcn(t, x, u, params) -> array`
@@ -1240,6 +1241,10 @@ def nlsys(
         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
 
@@ -1285,6 +1290,12 @@ def nlsys(
     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
@@ -1308,9 +1319,34 @@ def nlsys(
     >>> timepts = np.linspace(0, 10)
     >>> response = ct.input_output_response(
     ...     kincar, timepts, [10, 0.05 * np.sin(timepts)])
+
     """
-    return NonlinearIOSystem(
-        updfcn, outfcn, inputs=inputs, outputs=outputs, states=states, **kwargs)
+    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)
 
 
 def input_output_response(
@@ -2240,7 +2276,8 @@ def interconnect(
 
     params : dict, optional
         Parameter values for the systems.  Passed to the evaluation functions
-        for the system as default values, overriding internal defaults.
+        for the system as default values, overriding internal defaults.  If
+        not specified, defaults to parameters from subsystems.
 
     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 0eb49c82..77cfd370 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 1dd0be32..c5c5a803 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -586,7 +586,8 @@ 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, **kwargs):
+        name=None, inputs=None, outputs=None, states=None, params=None,
+        **kwargs):
     r"""Create an I/O system using a (full) state feedback controller.
 
     This function creates an input/output system that implements a
@@ -751,6 +752,10 @@ 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
@@ -773,7 +778,8 @@ def create_statefbk_iosystem(
     if not isinstance(sys, NonlinearIOSystem):
         raise ControlArgument("Input system must be I/O system")
 
-    # Process (legacy) keywords
+    # Process keywords
+    params = sys.params if params is None else params
     controller_type = _process_legacy_keyword(
         kwargs, 'type', 'controller_type', controller_type)
     if kwargs:
@@ -970,10 +976,10 @@ def create_statefbk_iosystem(
 
             return u
 
-        params = {} if gainsched else {'K': K}
+        ctrl_params = {} if gainsched else {'K': K}
         ctrl = NonlinearIOSystem(
             _control_update, _control_output, name=name, inputs=inputs,
-            outputs=outputs, states=states, params=params)
+            outputs=outputs, states=states, params=ctrl_params)
 
     elif controller_type == 'iosystem' and feedfwd_pattern == 'trajgen':
         # Use the passed system to compute feedback compensation
@@ -1061,7 +1067,8 @@ 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
+        outlist=outlist, outputs=output_labels,
+        params= ctrl.params | params
     )
     return ctrl, closed
 
diff --git a/control/statesp.py b/control/statesp.py
index 7af9008f..3f53777e 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -1,52 +1,18 @@
-"""statesp.py
+# 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
 
-State space representation and functions.
+"""State space representation and functions.
 
-This file contains the StateSpace class, which is used to represent linear
-systems in state space.  This is the primary representation for the
+This module 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
@@ -101,30 +67,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.
+    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.
 
     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.
-    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).
+    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.
 
     Notes
     -----
@@ -164,13 +130,12 @@ 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 `'{:'`
@@ -194,12 +159,7 @@ class StateSpace(NonlinearIOSystem, LTI):
         True for unspecified sampling time).  To call the copy constructor,
         call StateSpace(sys), where sys is a StateSpace object.
 
-        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).
+        See :class:`StateSpace` and :func:`ss` for more information.
 
         """
         #
@@ -1541,7 +1501,7 @@ def ss(*args, **kwargs):
 
     Create a state space system.
 
-    The function accepts either 1, 2, 4 or 5 parameters:
+    The function accepts either 1, 4 or 5 positional parameters:
 
     ``ss(sys)``
         Convert a linear system into space system form. Always creates a
@@ -1565,11 +1525,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 *array like* data types or strings.
-        Everything that the constructor of :class:`numpy.matrix` accepts is
-        permissible here too.
+        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.
 
-    ``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
@@ -1584,23 +1544,33 @@ 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
     ------
@@ -1609,7 +1579,7 @@ def ss(*args, **kwargs):
 
     See Also
     --------
-    tf, ss2tf, tf2ss
+    tf, ss2tf, tf2ss, zpk
 
     Notes
     -----
@@ -2291,7 +2261,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[i][j][0] / sys.den[i][j][0]
+                    D[i, j] = sys.num_array[i, j][0] / sys.den_array[i, j][0]
                 newsys = StateSpace([], [], [], D, sys.dt)
             else:
                 if not issiso(sys):
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 5304ea63..9ebbaf4f 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -1,64 +1,29 @@
-"""xferfcn.py
+# 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
 
-Transfer function representation and functions.
+"""Transfer function representation and functions.
 
-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$
+This module contains the TransferFunction class and also functions that
+operate on transfer functions.  This is the primary representation for the
+python-control library.
 
 """
 
 from collections.abc import Iterable
 from copy import deepcopy
-from itertools import chain
+from itertools import chain, product
 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, ndarray, nonzero, \
-    ones, pi, poly, polyadd, polymul, polyval, real, roots, sqrt, squeeze, \
-    where, zeros
+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 scipy.signal import TransferFunction as signalTransferFunction
 from scipy.signal import cont2discrete, tf2zpk, zpk2tf
 
@@ -85,46 +50,68 @@ class TransferFunction(LTI):
     A class for representing transfer functions.
 
     The TransferFunction class is used to represent systems in transfer
-    function form.
+    function form.  Transfer functions are usually created with the
+    :func:`~control.tf` factory function.
 
     Parameters
     ----------
-    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
+    num : 2D list of coefficient arrays
+        Polynomial coefficients of the numerator.
+    den : 2D list of coefficient arrays
+        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).
-    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'].
+        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.
-    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).
+    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.
 
     Notes
     -----
-    The attribues 'num' and 'den' are 2-D lists of arrays containing MIMO
-    numerator and denominator coefficients.  For example,
+    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,
 
-    >>> num[2][5] = numpy.array([1., 4., 8.])                   # doctest: +SKIP
+    >>> sys.num_array[2, 5]         # doctest: +SKIP
 
-    means that the numerator of the transfer function from the 6th input to
-    the 3rd output is set to s^2 + 4s + 8.
+    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``.
 
     A discrete time transfer function is created by specifying a nonzero
     'timebase' dt when the system is constructed:
@@ -172,13 +159,15 @@ class TransferFunction(LTI):
         Construct a transfer function.
 
         The default constructor is TransferFunction(num, den), where num and
-        den are lists of lists of arrays containing polynomial coefficients.
+        den are 2D arrays 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
@@ -210,8 +199,8 @@ class TransferFunction(LTI):
             raise TypeError("Needs 1, 2 or 3 arguments; received %i."
                              % len(args))
 
-        num = _clean_part(num)
-        den = _clean_part(den)
+        num = _clean_part(num, "numerator")
+        den = _clean_part(den, "denominator")
 
         #
         # Process keyword arguments
@@ -230,13 +219,16 @@ class TransferFunction(LTI):
 
         # Determine if the transfer function is static (needed for dt)
         static = True
-        for col in num + den:
-            for poly in col:
-                if len(poly) > 1:
+        for arr in [num, den]:
+            for poly in np.nditer(arr, flags=['refs_ok']):
+                if poly.item().size > 1:
                     static = False
+                    break
+            if not static:
+                break
 
         defaults = args[0] if len(args) == 1 else \
-            {'inputs': len(num[0]), 'outputs': len(num)}
+            {'inputs': num.shape[1], 'outputs': num.shape[0]}
 
         name, inputs, outputs, states, dt = _process_iosys_keywords(
                 kwargs, defaults, static=static)
@@ -252,27 +244,17 @@ class TransferFunction(LTI):
         # Check to make sure everything is consistent
         #
         # Make sure numerator and denominator matrices have consistent sizes
-        if self.ninputs != len(den[0]):
+        if self.ninputs != den.shape[1]:
             raise ValueError(
                 "The numerator has %i input(s), but the denominator has "
-                "%i input(s)." % (self.ninputs, len(den[0])))
-        if self.noutputs != len(den):
+                "%i input(s)." % (self.ninputs, den.shape[1]))
+        if self.noutputs != den.shape[0]:
             raise ValueError(
                 "The numerator has %i output(s), but the denominator has "
-                "%i output(s)." % (self.noutputs, len(den)))
+                "%i output(s)." % (self.noutputs, den.shape[0]))
 
         # 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
@@ -280,8 +262,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 k:
+                for k in den[i, j]:
+                    if np.any(k):
                         zeroden = False
                         break
                 if zeroden:
@@ -291,16 +273,16 @@ class TransferFunction(LTI):
 
                 # If we have zero numerators, set the denominator to 1.
                 zeronum = True
-                for k in num[i][j]:
-                    if k:
+                for k in num[i, j]:
+                    if np.any(k):
                         zeronum = False
                         break
                 if zeronum:
                     den[i][j] = ones(1)
 
         # Store the numerator and denominator
-        self.num = num
-        self.den = den
+        self.num_array = num
+        self.den_array = den
 
         #
         # Final processing
@@ -325,29 +307,22 @@ class TransferFunction(LTI):
     #: :meta hide-value:
     noutputs = 1
 
-    #: 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]]
+    # 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()
 
-    #: 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]]
+    @property
+    def den_list(self):
+        """Denominator polynomial (as 2D nested lists of 1D arrays)."""
+        return self.den_array.tolist()
+
+    # Legacy versions (TODO: add DeprecationWarning in a later release?)
+    num, den = num_list, den_list
 
     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
@@ -427,8 +402,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[i][j], x_arr) /
-                                 polyval(self.den[i][j], x_arr))
+                    out[i][j] = (polyval(self.num_array[i, j], x_arr) /
+                                 polyval(self.den_array[i, j], x_arr))
         return out
 
     def _truncatecoeff(self):
@@ -441,14 +416,14 @@ class TransferFunction(LTI):
         """
 
         # Beware: this is a shallow copy.  This should be okay.
-        data = [self.num, self.den]
+        data = [self.num_array, self.den_array]
         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
 
@@ -458,7 +433,7 @@ class TransferFunction(LTI):
                     else:
                         # Truncate the trivial coefficients.
                         data[p][i][j] = data[p][i][j][nonzero:]
-        [self.num, self.den] = data
+        [self.num_array, self.den_array] = data
 
     def __str__(self, var=None):
         """String representation of the transfer function.
@@ -478,16 +453,18 @@ class TransferFunction(LTI):
 
                 # Convert the numerator and denominator polynomials to strings.
                 if self.display_format == 'poly':
-                    numstr = _tf_polynomial_to_string(self.num[no][ni], var=var)
-                    denstr = _tf_polynomial_to_string(self.den[no][ni], var=var)
+                    numstr = _tf_polynomial_to_string(
+                        self.num_array[no, ni], var=var)
+                    denstr = _tf_polynomial_to_string(
+                        self.den_array[no, ni], var=var)
                 elif self.display_format == 'zpk':
-                    num = self.num[no][ni]
+                    num = self.num_array[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[no][ni])
+                        z, p, k = tf2zpk([1.], self.den_array[no, ni])
                         k = 0
                     else:
-                        z, p, k = tf2zpk(self.num[no][ni], self.den[no][ni])
+                        z, p, k = tf2zpk(self.num[no][ni], self.den_array[no, ni])
                     numstr = _tf_factorized_polynomial_to_string(
                         z, gain=k, var=var)
                     denstr = _tf_factorized_polynomial_to_string(p, var=var)
@@ -512,20 +489,33 @@ 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[0][0].__repr__(), den=self.den[0][0].__repr__(),
+                num=self.num_array[0, 0].__repr__(),
+                den=self.den_array[0, 0].__repr__(),
                 dt=', {}'.format(self.dt) if isdtime(self, strict=True)
                 else '')
         else:
-            return "TransferFunction({num}, {den}{dt})".format(
-                num=self.num.__repr__(), den=self.den.__repr__(),
-                dt=', {}'.format(self.dt) if isdtime(self, strict=True)
-                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
 
     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()
 
@@ -541,10 +531,13 @@ 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[no][ni], var=var)
-                    denstr = _tf_polynomial_to_string(self.den[no][ni], var=var)
+                    numstr = _tf_polynomial_to_string(
+                        self.num_array[no, ni], var=var)
+                    denstr = _tf_polynomial_to_string(
+                        self.den_array[no, ni], var=var)
                 elif self.display_format == 'zpk':
-                    z, p, k = tf2zpk(self.num[no][ni], self.den[no][ni])
+                    z, p, k = tf2zpk(
+                        self.num_array[no, ni], self.den_array[no, ni])
                     numstr = _tf_factorized_polynomial_to_string(
                         z, gain=k, var=var)
                     denstr = _tf_factorized_polynomial_to_string(p, var=var)
@@ -573,10 +566,10 @@ class TransferFunction(LTI):
 
     def __neg__(self):
         """Negate a transfer function."""
-        num = deepcopy(self.num)
+        num = deepcopy(self.num_array)
         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):
@@ -606,14 +599,14 @@ class TransferFunction(LTI):
         dt = common_timebase(self.dt, other.dt)
 
         # Preallocate the numerator and denominator of the sum.
-        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)]
+        num = _create_poly_array((self.noutputs, self.ninputs))
+        den = _create_poly_array((self.noutputs, self.ninputs))
 
         for i in range(self.noutputs):
             for j in range(self.ninputs):
-                num[i][j], den[i][j] = _add_siso(
-                    self.num[i][j], self.den[i][j],
-                    other.num[i][j], other.den[i][j])
+                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])
 
         return TransferFunction(num, den, dt)
 
@@ -648,14 +641,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))
 
-        inputs = other.ninputs
-        outputs = self.noutputs
+        ninputs = other.ninputs
+        noutputs = self.noutputs
 
         dt = common_timebase(self.dt, other.dt)
 
         # Preallocate the numerator and denominator of the sum.
-        num = [[[0] for j in range(inputs)] for i in range(outputs)]
-        den = [[[1] for j in range(inputs)] for i in range(outputs)]
+        num = _create_poly_array((noutputs, ninputs), [0])
+        den = _create_poly_array((noutputs, ninputs), [1])
 
         # Temporary storage for the summands needed to find the (i, j)th
         # element of the product.
@@ -663,17 +656,16 @@ class TransferFunction(LTI):
         den_summand = [[] for k in range(self.ninputs)]
 
         # Multiply & add.
-        for row in range(outputs):
-            for col in range(inputs):
+        for row in range(noutputs):
+            for col in range(ninputs):
                 for k in range(self.ninputs):
                     num_summand[k] = polymul(
-                        self.num[row][k], other.num[k][col])
+                        self.num_array[row, k], other.num_array[k, col])
                     den_summand[k] = polymul(
-                        self.den[row][k], other.den[k][col])
-                    num[row][col], den[row][col] = _add_siso(
-                        num[row][col], den[row][col],
+                        self.den_array[row, k], other.den_array[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):
@@ -692,14 +684,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))
 
-        inputs = self.ninputs
-        outputs = other.noutputs
+        ninputs = self.ninputs
+        noutputs = other.noutputs
 
         dt = common_timebase(self.dt, other.dt)
 
         # Preallocate the numerator and denominator of the sum.
-        num = [[[0] for j in range(inputs)] for i in range(outputs)]
-        den = [[[1] for j in range(inputs)] for i in range(outputs)]
+        num = _create_poly_array((noutputs, ninputs), [0])
+        den = _create_poly_array((noutputs, ninputs), [1])
 
         # Temporary storage for the summands needed to find the
         # (i, j)th element
@@ -707,13 +699,15 @@ class TransferFunction(LTI):
         num_summand = [[] for k in range(other.ninputs)]
         den_summand = [[] for k in range(other.ninputs)]
 
-        for i in range(outputs):  # Iterate through rows of product.
-            for j in range(inputs):  # Iterate through columns of product.
+        for i in range(noutputs):  # Iterate through rows of product.
+            for j in range(ninputs):  # Iterate through columns of product.
                 for k in range(other.ninputs):  # Multiply & add.
-                    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_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[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)
@@ -736,8 +730,8 @@ class TransferFunction(LTI):
 
         dt = common_timebase(self.dt, other.dt)
 
-        num = polymul(self.num[0][0], other.den[0][0])
-        den = polymul(self.den[0][0], other.num[0][0])
+        num = polymul(self.num_array[0, 0], other.den_array[0, 0])
+        den = polymul(self.den_array[0, 0], other.num_array[0, 0])
 
         return TransferFunction(num, den, dt)
 
@@ -785,15 +779,14 @@ class TransferFunction(LTI):
             indices[1], self.input_labels, slice_to_list=True)
         
         # Construct the transfer function for the subsyste
-        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)
+        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
 
         # Create the system name
         sysname = config.defaults['iosys.indexed_system_name_prefix'] + \
@@ -832,7 +825,7 @@ class TransferFunction(LTI):
                 "for SISO systems.")
         else:
             # for now, just give zeros of a SISO tf
-            return roots(self.num[0][0]).astype(complex)
+            return roots(self.num_array[0, 0]).astype(complex)
 
     def feedback(self, other=1, sign=-1):
         """Feedback interconnection between two LTI objects."""
@@ -846,10 +839,10 @@ class TransferFunction(LTI):
                 "MIMO systems.")
         dt = common_timebase(self.dt, other.dt)
 
-        num1 = self.num[0][0]
-        den1 = self.den[0][0]
-        num2 = other.num[0][0]
-        den2 = other.den[0][0]
+        num1 = self.num_array[0, 0]
+        den1 = self.den_array[0, 0]
+        num2 = other.num_array[0, 0]
+        den2 = other.den_array[0, 0]
 
         num = polymul(num1, den2)
         den = polyadd(polymul(den2, den1), -sign * polymul(num2, num1))
@@ -862,7 +855,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
@@ -870,17 +863,17 @@ class TransferFunction(LTI):
         sqrt_eps = sqrt(float_info.epsilon)
 
         # pre-allocate arrays
-        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)]
+        num = _create_poly_array((self.noutputs, self.ninputs))
+        den = _create_poly_array((self.noutputs, self.ninputs))
 
         for i in range(self.noutputs):
             for j in range(self.ninputs):
 
                 # split up in zeros, poles and gain
                 newzeros = []
-                zeros = roots(self.num[i][j])
-                poles = roots(self.den[i][j])
-                gain = self.num[i][j][0] / self.den[i][j][0]
+                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]
 
                 # check all zeros
                 for z in zeros:
@@ -895,19 +888,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 list of a list of :class:`scipy.signal.lti` objects.
+        """Return a 2D array 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.
@@ -1116,7 +1109,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.
@@ -1263,7 +1256,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.
@@ -1276,7 +1269,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.
@@ -1320,7 +1313,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"
 
@@ -1367,7 +1360,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)
@@ -1410,9 +1403,11 @@ def _tf_factorized_polynomial_to_string(roots, gain=1, var='s'):
 
 
 def _tf_string_to_latex(thestr, var='s'):
-    """ make sure to superscript all digits in a polynomial string
-        and convert float coefficients in scientific notation
-        to prettier LaTeX representation """
+    """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)
@@ -1555,15 +1550,20 @@ 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 nested lists
-        of array_like objects. (A 3 dimensional data structure in total.)
-        (For details see note below.)
+        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.
 
     ``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').
@@ -1571,20 +1571,27 @@ def tf(*args, **kwargs):
     Parameters
     ----------
     sys : LTI (StateSpace or TransferFunction)
-        A linear system
+        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.
     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
+        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).
     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
     -------
-    out: :class:`TransferFunction`
-        The new linear system
+    sys : TransferFunction
+        The new linear system.
 
     Other Parameters
     ----------------
@@ -1592,6 +1599,8 @@ 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.
@@ -1599,9 +1608,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
     --------
@@ -1612,9 +1621,10 @@ 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.
-    ``den[i][j]`` works the same way.
+    for the transfer function from the (j+1)st input to the (i+1)st output,
+    and ``den[i][j]`` works the same way.
 
     The list ``[2, 3, 4]`` denotes the polynomial :math:`2s^2 + 3s + 4`.
 
@@ -1639,11 +1649,7 @@ def tf(*args, **kwargs):
     >>> sys_tf = ct.tf(sys_ss)
 
     """
-
-    if len(args) == 2 or len(args) == 3:
-        return TransferFunction(*args, **kwargs)
-
-    elif len(args) == 1 and isinstance(args[0], str):
+    if len(args) == 1 and isinstance(args[0], str):
         # Make sure there were no extraneous keywords
         if kwargs:
             raise TypeError("unrecognized keywords: ", str(kwargs))
@@ -1654,19 +1660,53 @@ 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
-        sys = args[0]
-        if isinstance(sys, StateSpace):
+        if isinstance(sys := args[0], 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))
-    else:
-        raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
+
+    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)
 
 
 def zpk(zeros, poles, gain, *args, **kwargs):
@@ -1688,7 +1728,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
@@ -1710,7 +1750,7 @@ def zpk(zeros, poles, gain, *args, **kwargs):
 
     Returns
     -------
-    out: :class:`TransferFunction`
+    out: `TransferFunction`
         Transfer function with given zeros, poles, and gain.
 
     Examples
@@ -1846,7 +1886,7 @@ def tfdata(sys):
     return tf.num, tf.den
 
 
-def _clean_part(data):
+def _clean_part(data, name="<unknown>"):
     """
     Return a valid, cleaned up numerator or denominator
     for the TransferFunction class.
@@ -1862,27 +1902,37 @@ def _clean_part(data):
     valid_types = (int, float, complex, np.number)
     valid_collection = (list, tuple, ndarray)
 
-    if (isinstance(data, valid_types) or
+    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
             (isinstance(data, ndarray) and data.ndim == 0)):
         # Data is a scalar (including 0d ndarray)
-        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])]
+        out = np.empty((1,1), dtype=np.ndarray)
+        out[0, 0] = array([data])
     elif (isinstance(data, valid_collection) and
             all([isinstance(d, valid_types) for d in data])):
-        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])
+        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])
     else:
         # If the user passed in anything else, then it's unclear what
         # the meaning is.
@@ -1891,20 +1941,36 @@ def _clean_part(data):
             "(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(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])
+    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
 
-    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 85e8d8bd..36dafd61 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.NonlinearIOSystem(
+vehicle = ct.nlsys(
     vehicle_update, vehicle_output, states=3, name='vehicle',
-    inputs=('v', 'phi'),
-    outputs=('x', 'y', 'theta'))
+    inputs=('v', 'phi'), outputs=('x', 'y', 'theta'),
+    params={'wheelbase': 3, 'maxsteer': 0.5})
 
 #
 # Gain scheduled controller
@@ -89,10 +89,12 @@ def control_output(t, x, u, params):
     return  np.array([v, phi])
 
 # Define the controller as an input/output system
-controller = ct.NonlinearIOSystem(
+controller = ct.nlsys(
     None, control_output, name='controller',        # static system
     inputs=('ex', 'ey', 'etheta', 'vd', 'phid'),    # system inputs
-    outputs=('v', 'phi')                            # system outputs
+    outputs=('v', 'phi'),                           # system outputs
+    params={'longpole': -2, 'latpole1': -1/2 + sqrt(-7)/2,
+            'latpole2': -1/2 - sqrt(-7)/2, 'wheelbase': 3}
 )
 
 #
@@ -113,7 +115,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.NonlinearIOSystem(
+trajgen = ct.nlsys(
     None, trajgen_output, name='trajgen',
     inputs=('vref', 'yref'),
     outputs=('xd', 'yd', 'thetad', 'vd', 'phid'))
@@ -156,10 +158,13 @@ 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']
+    outputs=['x', 'y', 'theta', 'v', 'phi'],
+
+    # Parameters
+    params=trajgen.params | vehicle.params | controller.params,
 )
 
 # Set up the simulation conditions
@@ -220,7 +225,8 @@ 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')
+    gainsched_indices=['vd', 'theta'], gainsched_method='linear',
+    params=vehicle.params | controller.params)
 
 # Connect everything together (note that controller inputs are different)
 steering = ct.interconnect(
@@ -245,10 +251,13 @@ 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']
+    outputs=['x', 'y', 'theta', 'v', 'phi'],
+
+    # Parameters
+    params=steering.params
 )
 
 # 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 8ea67e0f..e0d64f01 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[i][j],
-                tf_b.num[i][j],
+                tf_a.num_array[i, j],
+                tf_b.num_array[i, j],
                 rtol=rtol,
                 atol=atol,
             ):
                 return False
             if not np.allclose(
-                tf_a.den[i][j],
-                tf_b.den[i][j],
+                tf_a.den_array[i, j],
+                tf_b.den_array[i, j],
                 rtol=rtol,
                 atol=atol,
             ):
diff --git a/control/tests/docstrings_test.py b/control/tests/docstrings_test.py
index 991ead3e..16647895 100644
--- a/control/tests/docstrings_test.py
+++ b/control/tests/docstrings_test.py
@@ -21,6 +21,9 @@ 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
@@ -31,13 +34,10 @@ 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,14 +260,244 @@ 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=""):
+def _check_parameter_docs(
+        funcname, argname, docstring, prefix="", fail_if_missing=True):
     funcname = prefix + funcname
 
     # Find the "Parameters" section of docstring, where we start searching
-    if not (match := re.search(r"\nParameters\n----", docstring)):
-        pytest.fail(f"{funcname} docstring missing Parameters section")
+    if not (match := re.search(
+            r"\nParameters\n----", docstring)):
+        if fail_if_missing:
+            pytest.fail(f"{funcname} docstring missing Parameters section")
+        else:
+            return False
     else:
         start = match.start()
 
@@ -299,7 +529,10 @@ def _check_parameter_docs(funcname, argname, docstring, prefix=""):
             docstring)):
         if verbose:
             print(f"      {funcname}: {argname} not documented")
-        pytest.fail(f"{funcname} '{argname}' not documented")
+        if fail_if_missing:
+            pytest.fail(f"{funcname} '{argname}' not documented")
+        else:
+            return False
 
     # Make sure there isn't another instance
     second_match = re.search(
@@ -307,3 +540,5 @@ def _check_parameter_docs(funcname, argname, docstring, prefix=""):
             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 c2a29ee2..b08cd826 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 604488ca..d124859f 100644
--- a/control/tests/interconnect_test.py
+++ b/control/tests/interconnect_test.py
@@ -666,15 +666,29 @@ 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': 1})
+        states=1, inputs='u', outputs='y', params={'a': 2, 'c':2})
 
     # Simple system for serial interconnection
     sys2 = ct.nlsys(
         None, lambda t, x, u, params: u[0],
-        inputs='r', outputs='u')
+        inputs='r', outputs='u', params={'a': 4, 'b': 3})
 
-    # Create a series interconnection
+    # Make sure default parameters get set as expected
     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 0216235b..54d6d56c 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -17,7 +17,8 @@ import pytest
 import scipy
 
 import control as ct
-
+import control.flatsys as fs
+from control.tests.conftest import slycotonly
 
 class TestIOSys:
 
@@ -930,6 +931,8 @@ 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])
@@ -959,7 +962,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(
-            np.sort(w), np.sort(np.concatenate((w_default, [2j, -2j]))))
+             w, np.concatenate([w_update, w_update]))
 
         # Show that we can change the parameters at linearization
         lin_series_override = ct.linearize(
@@ -2284,3 +2287,55 @@ 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/nlsys_test.py b/control/tests/nlsys_test.py
index 7f649e0c..926ca436 100644
--- a/control/tests/nlsys_test.py
+++ b/control/tests/nlsys_test.py
@@ -7,11 +7,15 @@ The main test functions are contained in iosys_test.py.
 
 """
 
-import pytest
-import numpy as np
 import math
+import re
+
+import numpy as np
+import pytest
+
 import control as ct
 
+
 # Basic test of nlsys()
 def test_nlsys_basic():
     def kincar_update(t, x, u, params):
@@ -154,3 +158,44 @@ 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 3928fb72..22a946fe 100644
--- a/control/tests/statefbk_test.py
+++ b/control/tests/statefbk_test.py
@@ -930,9 +930,10 @@ def unicycle():
 
     return ct.NonlinearIOSystem(
         unicycle_update, None,
-        inputs = ['v', 'phi'],
-        outputs = ['x', 'y', 'theta'],
-        states = ['x_', 'y_', 'theta_'])
+        inputs=['v', 'phi'],
+        outputs=['x', 'y', 'theta'],
+        states=['x_', 'y_', 'theta_'],
+        params={'a': 1})        # only used for testing params
 
 from math import pi
 
@@ -1194,3 +1195,40 @@ 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 1798c524..2bb0badc 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,
-          "A is the wrong shape; expected \(3, 3\)"),
+          r"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,
-          "B is the wrong shape; expected \(3, 2\)"),
+          r"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,
-          "C is the wrong shape; expected \(2, 3\)"),
+          r"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,
-          "D is the wrong shape; expected \(2, 2\)"),
+          r"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,
-          "D is the wrong shape; expected \(2, 2\)"),
+          r"D is the wrong shape; expected \(2, 2\)"),
         ])
     def test_constructor_invalid(self, args, exc, errmsg):
         """Test invalid input to StateSpace() constructor"""
diff --git a/control/tests/xferfcn_input_test.py b/control/tests/xferfcn_input_test.py
index 46efbd25..c375d768 100644
--- a/control/tests/xferfcn_input_test.py
+++ b/control/tests/xferfcn_input_test.py
@@ -64,15 +64,18 @@ def test_clean_part(num, fun, dtype):
     num_ = _clean_part(numa)
     ref_ = np.array(num, dtype=float, ndmin=3)
 
-    assert isinstance(num_, list)
-    assert np.all([isinstance(part, list) for part in num_])
+    assert isinstance(num_, np.ndarray)
+    assert num_.ndim == 2
     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.]], "a"])
+@pytest.mark.parametrize("badinput", [
+    # [[0., 1.], [2., 3.]],             # OK: treated as static array
+    np.ones((2, 2, 2, 2)),
+    "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 d480cef6..3f87ef1d 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -30,13 +30,6 @@ 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])
@@ -190,8 +183,10 @@ class TestXferFcn:
 
         for i in range(sys3.noutputs):
             for j in range(sys3.ninputs):
-                np.testing.assert_allclose(sys2.num[i][j], sys3.num[i][j])
-                np.testing.assert_allclose(sys2.den[i][j], sys3.den[i][j])
+                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])
 
     # Tests for TransferFunction.__add__
 
@@ -236,8 +231,8 @@ class TestXferFcn:
 
         for i in range(sys3.noutputs):
             for j in range(sys3.ninputs):
-                np.testing.assert_allclose(sys3.num[i][j], num3[i][j])
-                np.testing.assert_allclose(sys3.den[i][j], den3[i][j])
+                np.testing.assert_allclose(sys3.num_array[i, j], num3[i][j])
+                np.testing.assert_allclose(sys3.den_array[i, j], den3[i][j])
 
     # Tests for TransferFunction.__sub__
 
@@ -284,8 +279,8 @@ class TestXferFcn:
 
         for i in range(sys3.noutputs):
             for j in range(sys3.ninputs):
-                np.testing.assert_allclose(sys3.num[i][j], num3[i][j])
-                np.testing.assert_allclose(sys3.den[i][j], den3[i][j])
+                np.testing.assert_allclose(sys3.num_array[i, j], num3[i][j])
+                np.testing.assert_allclose(sys3.den_array[i, j], den3[i][j])
 
     # Tests for TransferFunction.__mul__
 
@@ -340,8 +335,8 @@ class TestXferFcn:
 
         for i in range(sys3.noutputs):
             for j in range(sys3.ninputs):
-                np.testing.assert_allclose(sys3.num[i][j], num3[i][j])
-                np.testing.assert_allclose(sys3.den[i][j], den3[i][j])
+                np.testing.assert_allclose(sys3.num_array[i, j], num3[i][j])
+                np.testing.assert_allclose(sys3.den_array[i, j], den3[i][j])
 
     # Tests for TransferFunction.__div__
 
@@ -662,10 +657,10 @@ class TestXferFcn:
 
         for i in range(sys.noutputs):
             for j in range(sys.ninputs):
-                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])
+                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])
 
     def test_minreal(self):
         """Try the minreal function, and also test easy entry by creation
@@ -1121,8 +1116,10 @@ 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[p][m], H2.num[p][m])
-                np.testing.assert_array_almost_equal(H.den[p][m], H2.den[p][m])
+                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])
             assert H.dt == H2.dt
 
     def test_sample_named_signals(self):
@@ -1180,8 +1177,10 @@ 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[i][j])
-                np.testing.assert_allclose(sslti[i][j].den, mimotf.den[i][j])
+                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])
                 if mimotf.dt == 0:
                     assert sslti[i][j].dt is None
                 else:
@@ -1286,3 +1285,35 @@ 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)