Repository: python-control/python-control

Issue:
Multiplying a MIMO `TransferFunction` by `1.0` multiplies by a matrix of ones, not identity
Multiplying a MIMO `TransferFunction` object by the Python float `1.0` appears to multiply it by a matrix of ones, rather than the identity matrix. The same goes for any scalar. Is this the intended behaviour? I checked in the documentation and didn't see anything.

Example:
```python
import control

G = control.TransferFunction(
    [
        [[87.8], [-86.4]],
        [[108.2], [-109.6]],
    ],
    [
        [[1, 1], [1, 1]],
        [[1, 1], [1, 1]],
    ],
)

I = control.TransferFunction(
    [
        [[1], [1]],
        [[1], [1]],
    ],
    [
        [[1], [1]],
        [[1], [1]],
    ],
)

print(G)
print(1 * G)
print(I * G)
```

Discussion:

Solution patch:
diff --git a/control/bdalg.py b/control/bdalg.py
index 024d95fba..d907cd3c5 100644
--- a/control/bdalg.py
+++ b/control/bdalg.py
@@ -10,6 +10,8 @@
 negate
 feedback
 connect
+combine_tf
+split_tf
 
 """
 
@@ -63,7 +65,8 @@
 from . import xferfcn as tf
 from .iosys import InputOutputSystem
 
-__all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect']
+__all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect',
+           'combine_tf', 'split_tf']
 
 
 def series(sys1, *sysn, **kwargs):
@@ -507,3 +510,212 @@ def connect(sys, Q, inputv, outputv):
         Ytrim[i,y-1] = 1.
 
     return Ytrim * sys * Utrim
+
+def combine_tf(tf_array):
+    """Combine array-like of transfer functions into MIMO transfer function.
+
+    Parameters
+    ----------
+    tf_array : list of list of TransferFunction or array_like
+        Transfer matrix represented as a two-dimensional array or list-of-lists
+        containing TransferFunction objects. The TransferFunction objects can
+        have multiple outputs and inputs, as long as the dimensions are
+        compatible.
+
+    Returns
+    -------
+    TransferFunction
+        Transfer matrix represented as a single MIMO TransferFunction object.
+
+    Raises
+    ------
+    ValueError
+        If timesteps of transfer functions do not match.
+    ValueError
+        If ``tf_array`` has incorrect dimensions.
+    ValueError
+        If the transfer functions in a row have mismatched output or input
+        dimensions.
+
+    Examples
+    --------
+    Combine two transfer functions
+
+    >>> s = control.TransferFunction.s
+    >>> control.combine_tf([
+    ...     [1 / (s + 1)],
+    ...     [s / (s + 2)],
+    ... ])
+    TransferFunction([[array([1])], [array([1, 0])]],
+                     [[array([1, 1])], [array([1, 2])]])
+
+    Combine NumPy arrays with transfer functions
+
+    >>> control.combine_tf([
+    ...     [np.eye(2), np.zeros((2, 1))],
+    ...     [np.zeros((1, 2)), control.TransferFunction([1], [1, 0])],
+    ... ])
+    TransferFunction([[array([1.]), array([0.]), array([0.])],
+                      [array([0.]), array([1.]), array([0.])],
+                      [array([0.]), array([0.]), array([1])]],
+                     [[array([1.]), array([1.]), array([1.])],
+                      [array([1.]), array([1.]), array([1.])],
+                      [array([1.]), array([1.]), array([1, 0])]])
+    """
+    # Find common timebase or raise error
+    dt_list = []
+    try:
+        for row in tf_array:
+            for tfn in row:
+                dt_list.append(getattr(tfn, "dt", None))
+    except OSError:
+        raise ValueError("`tf_array` has too few dimensions.")
+    dt_set = set(dt_list)
+    dt_set.discard(None)
+    if len(dt_set) > 1:
+        raise ValueError("Timesteps of transfer functions are "
+                         f"mismatched: {dt_set}")
+    elif len(dt_set) == 0:
+        dt = None
+    else:
+        dt = dt_set.pop()
+    # Convert all entries to transfer function objects
+    ensured_tf_array = []
+    for row in tf_array:
+        ensured_row = []
+        for tfn in row:
+            ensured_row.append(_ensure_tf(tfn, dt))
+        ensured_tf_array.append(ensured_row)
+    # Iterate over
+    num = []
+    den = []
+    for row_index, row in enumerate(ensured_tf_array):
+        for j_out in range(row[0].noutputs):
+            num_row = []
+            den_row = []
+            for col in row:
+                if col.noutputs != row[0].noutputs:
+                    raise ValueError(
+                        "Mismatched number of transfer function outputs in "
+                        f"row {row_index}."
+                    )
+                for j_in in range(col.ninputs):
+                    num_row.append(col.num[j_out][j_in])
+                    den_row.append(col.den[j_out][j_in])
+            num.append(num_row)
+            den.append(den_row)
+    for row_index, row in enumerate(num):
+        if len(row) != len(num[0]):
+            raise ValueError(
+                "Mismatched number transfer function inputs in row "
+                f"{row_index} of numerator."
+            )
+    for row_index, row in enumerate(den):
+        if len(row) != len(den[0]):
+            raise ValueError(
+                "Mismatched number transfer function inputs in row "
+                f"{row_index} of denominator."
+            )
+    return tf.TransferFunction(num, den, dt=dt)
+
+def split_tf(transfer_function):
+    """Split MIMO transfer function into NumPy array of SISO tranfer functions.
+
+    Parameters
+    ----------
+    transfer_function : TransferFunction
+        MIMO transfer function to split.
+
+    Returns
+    -------
+    np.ndarray
+        NumPy array of SISO transfer functions.
+
+    Examples
+    --------
+    Split a MIMO transfer function
+
+    >>> G = control.TransferFunction(
+    ...     [
+    ...         [[87.8], [-86.4]],
+    ...         [[108.2], [-109.6]],
+    ...     ],
+    ...     [
+    ...         [[1, 1], [1, 1]],
+    ...         [[1, 1], [1, 1]],
+    ...     ],
+    ... )
+    >>> control.split_tf(G)
+    array([[TransferFunction(array([87.8]), array([1, 1])),
+            TransferFunction(array([-86.4]), array([1, 1]))],
+           [TransferFunction(array([108.2]), array([1, 1])),
+            TransferFunction(array([-109.6]), array([1, 1]))]], dtype=object)
+    """
+    tf_split_lst = []
+    for i_out in range(transfer_function.noutputs):
+        row = []
+        for i_in in range(transfer_function.ninputs):
+            row.append(
+                tf.TransferFunction(
+                    transfer_function.num[i_out][i_in],
+                    transfer_function.den[i_out][i_in],
+                    dt=transfer_function.dt,
+                )
+            )
+        tf_split_lst.append(row)
+    return np.array(tf_split_lst, dtype=object)
+
+def _ensure_tf(arraylike_or_tf, dt=None):
+    """Convert an array-like to a transfer function.
+
+    Parameters
+    ----------
+    arraylike_or_tf : TransferFunction or array_like
+        Array-like or transfer function.
+    dt : None, True or float, optional
+        System timebase. 0 (default) indicates continuous
+        time, True indicates discrete time with unspecified sampling
+        time, positive number is discrete time with specified
+        sampling time, None indicates unspecified timebase (either
+        continuous or discrete time). If None, timestep is not validated.
+
+    Returns
+    -------
+    TransferFunction
+        Transfer function.
+
+    Raises
+    ------
+    ValueError
+        If input cannot be converted to a transfer function.
+    ValueError
+        If the timesteps do not match.
+    """
+    # If the input is already a transfer function, return it right away
+    if isinstance(arraylike_or_tf, tf.TransferFunction):
+        # If timesteps don't match, raise an exception
+        if (dt is not None) and (arraylike_or_tf.dt != dt):
+            raise ValueError(
+                f"`arraylike_or_tf.dt={arraylike_or_tf.dt}` does not match "
+                f"argument `dt={dt}`."
+            )
+        return arraylike_or_tf
+    if np.ndim(arraylike_or_tf) > 2:
+        raise ValueError(
+            "Array-like must have less than two dimensions to be converted "
+            "into a transfer function."
+        )
+    # If it's not, then convert it to a transfer function
+    arraylike_3d = np.atleast_3d(arraylike_or_tf)
+    try:
+        tfn = tf.TransferFunction(
+            arraylike_3d,
+            np.ones_like(arraylike_3d),
+            dt,
+        )
+    except TypeError:
+        raise ValueError(
+            "`arraylike_or_tf` must only contain array-likes or transfer "
+            "functions."
+        )
+    return tfn
diff --git a/control/canonical.py b/control/canonical.py
index a62044322..7be7f88ad 100644
--- a/control/canonical.py
+++ b/control/canonical.py
@@ -205,7 +205,7 @@ def similarity_transform(xsys, T, timescale=1, inverse=False):
     timescale : float, optional
         If present, also rescale the time unit to tau = timescale * t
     inverse : bool, optional
-        If True (default), transform so z = T x.  If False, transform
+        If False (default), transform so z = T x.  If True, transform
         so x = T z.
 
     Returns
diff --git a/control/frdata.py b/control/frdata.py
index 1bdf28528..ac032d3f7 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -261,6 +261,11 @@ def __init__(self, *args, **kwargs):
 
         # create interpolation functions
         if smooth:
+            # Set the order of the fit
+            if self.omega.size < 2:
+                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]),
                                dtype=tuple)
             for i in range(self.fresp.shape[0]):
@@ -268,7 +273,8 @@ def __init__(self, *args, **kwargs):
                     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)
+                        w=1.0/(absolute(self.fresp[i, j, :]) + 0.001),
+                        s=0.0, k=degree)
         else:
             self.ifunc = None
 
@@ -392,7 +398,12 @@ def __add__(self, other):
 
         # Convert the second argument to a frequency response function.
         # or re-base the frd to the current omega (if needed)
-        other = _convert_to_frd(other, omega=self.omega)
+        if isinstance(other, (int, float, complex, np.number)):
+            other = _convert_to_frd(
+                other, omega=self.omega,
+                inputs=self.ninputs, outputs=self.noutputs)
+        else:
+            other = _convert_to_frd(other, omega=self.omega)
 
         # Check that the input-output sizes are consistent.
         if self.ninputs != other.ninputs:
diff --git a/control/modelsimp.py b/control/modelsimp.py
index 966da1ce7..fe519b82d 100644
--- a/control/modelsimp.py
+++ b/control/modelsimp.py
@@ -1,4 +1,3 @@
-#! TODO: add module docstring
 # modelsimp.py - tools for model simplification
 #
 # Author: Steve Brunton, Kevin Chen, Lauren Padilla
@@ -6,47 +5,16 @@
 #
 # This file contains routines for obtaining reduced order models
 #
-# 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.
-#
-# $Id$
+
+import warnings
 
 # External packages and modules
 import numpy as np
-import warnings
-from .exception import ControlSlycot, ControlArgument, ControlDimension
-from .iosys import isdtime, isctime
-from .statesp import StateSpace
+
+from .exception import ControlArgument, ControlDimension, ControlSlycot
+from .iosys import isctime, isdtime
 from .statefbk import gram
+from .statesp import StateSpace
 from .timeresp import TimeResponseData
 
 __all__ = ['hankel_singular_values', 'balanced_reduction', 'model_reduction',
@@ -108,89 +76,146 @@ def hankel_singular_values(sys):
     return hsv[::-1]
 
 
-def model_reduction(sys, ELIM, method='matchdc'):
-    """Model reduction by state elimination.
-    
-    Model reduction of `sys` by eliminating the states in `ELIM` using a given
-    method.
+def model_reduction(
+        sys, elim_states=None, method='matchdc', elim_inputs=None,
+        elim_outputs=None, keep_states=None, keep_inputs=None,
+        keep_outputs=None, warn_unstable=True):
+    """Model reduction by input, output, or state elimination.
+
+    This function produces a reduced-order model of a system by eliminating
+    specified inputs, outputs, and/or states from the original system.  The
+    specific states, inputs, or outputs that are eliminated can be
+    specified by either listing the states, inputs, or outputs to be
+    eliminated or those to be kept.
+
+    Two methods of state reduction are possible: 'truncate' removes the
+    states marked for elimination, while 'matchdc' replaces the eliminated
+    states with their equilibrium values (thereby keeping the input/output
+    gain unchanged at zero frequency ["DC"]).
 
     Parameters
     ----------
     sys : StateSpace
         Original system to reduce.
-    ELIM : array
-        Vector of states to eliminate.
+    elim_inputs, elim_outputs, elim_states : array of int or str, optional
+        Vector of inputs, outputs, or states to eliminate.  Can be specified
+        either as an offset into the appropriate vector or as a signal name.
+    keep_inputs, keep_outputs, keep_states : array, optional
+        Vector of inputs, outputs, or states to keep.  Can be specified
+        either as an offset into the appropriate vector or as a signal name.
     method : string
-        Method of removing states in `ELIM`: either 'truncate' or
-        'matchdc'.
+        Method of removing states: either 'truncate' or 'matchdc' (default).
+    warn_unstable : bool, option
+        If `False`, don't warn if system is unstable.
 
     Returns
     -------
     rsys : StateSpace
-        A reduced order model.
+        Reduced order model.
 
     Raises
     ------
     ValueError
-        Raised under the following conditions:
+        If `method` is not either 'matchdc' or 'truncate'.
+    NotImplementedError
+        If the 'matchdc' method is used for a discrete time system.
 
-            * if `method` is not either ``'matchdc'`` or ``'truncate'``
-
-            * if eigenvalues of `sys.A` are not all in left half plane
-              (`sys` must be stable)
+    Warns
+    -----
+    UserWarning
+        If eigenvalues of `sys.A` are not all stable.
 
     Examples
     --------
     >>> G = ct.rss(4)
-    >>> Gr = ct.modred(G, [0, 2], method='matchdc')
+    >>> Gr = ct.model_reduction(G, [0, 2], method='matchdc')
     >>> Gr.nstates
     2
 
-    """
+    See Also
+    --------
+    balanced_reduction : Eliminate states using Hankel singular values.
+    minimal_realization : Eliminate unreachable or unobservable states.
 
-    # Check for ss system object, need a utility for this?
+    Notes
+    -----
+    The model_reduction function issues a warning if the system has
+    unstable eigenvalues, since in those situations the stability of the
+    reduced order model may be different than the stability of the full
+    model.  No other checking is done, so users must to be careful not to
+    render a system unobservable or unreachable.
 
-    # TODO: Check for continous or discrete, only continuous supported for now
-    #   if isCont():
-    #       dico = 'C'
-    #   elif isDisc():
-    #       dico = 'D'
-    #   else:
-    if (isctime(sys)):
-        dico = 'C'
-    else:
-        raise NotImplementedError("Function not implemented in discrete time")
+    States, inputs, and outputs can be specified using integer offsets or
+    using signal names.  Slices can also be specified, but must use the
+    Python ``slice()`` function.
+
+    """
+    if not isinstance(sys, StateSpace):
+        raise TypeError("system must be a StateSpace system")
 
     # Check system is stable
-    if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
-        raise ValueError("Oops, the system is unstable!")
-
-    ELIM = np.sort(ELIM)
-    # Create list of elements not to eliminate (NELIM)
-    NELIM = [i for i in range(len(sys.A)) if i not in ELIM]
-    # A1 is a matrix of all columns of sys.A not to eliminate
-    A1 = sys.A[:, NELIM[0]].reshape(-1, 1)
-    for i in NELIM[1:]:
-        A1 = np.hstack((A1, sys.A[:, i].reshape(-1, 1)))
-    A11 = A1[NELIM, :]
-    A21 = A1[ELIM, :]
-    # A2 is a matrix of all columns of sys.A to eliminate
-    A2 = sys.A[:, ELIM[0]].reshape(-1, 1)
-    for i in ELIM[1:]:
-        A2 = np.hstack((A2, sys.A[:, i].reshape(-1, 1)))
-    A12 = A2[NELIM, :]
-    A22 = A2[ELIM, :]
-
-    C1 = sys.C[:, NELIM]
-    C2 = sys.C[:, ELIM]
-    B1 = sys.B[NELIM, :]
-    B2 = sys.B[ELIM, :]
-
-    if method == 'matchdc':
-        # if matchdc, residualize
+    if warn_unstable:
+        if isctime(sys) and np.any(np.linalg.eigvals(sys.A).real >= 0.0) or \
+           isdtime(sys) and np.any(np.abs(np.linalg.eigvals(sys.A)) >= 1):
+            warnings.warn("System is unstable; reduction may be meaningless")
+
+    # Utility function to process keep/elim keywords
+    def _process_elim_or_keep(elim, keep, labels):
+        def _expand_key(key):
+            if key is None:
+                return []
+            elif isinstance(key, str):
+                return labels.index(key)
+            elif isinstance(key, list):
+                return [_expand_key(k) for k in key]
+            elif isinstance(key, slice):
+                return range(len(labels))[key]
+            else:
+                return key
 
+        elim = np.atleast_1d(_expand_key(elim))
+        keep = np.atleast_1d(_expand_key(keep))
+            
+        if len(elim) > 0 and len(keep) > 0:
+            raise ValueError(
+                "can't provide both 'keep' and 'elim' for same variables")
+        elif len(keep) > 0:
+            keep = np.sort(keep).tolist()
+            elim = [i for i in range(len(labels)) if i not in keep]
+        else:
+            elim = [] if elim is None else np.sort(elim).tolist()
+            keep = [i for i in range(len(labels)) if i not in elim]
+        return elim, keep
+
+    # Determine which states to keep
+    elim_states, keep_states = _process_elim_or_keep(
+        elim_states, keep_states, sys.state_labels)
+    elim_inputs, keep_inputs = _process_elim_or_keep(
+        elim_inputs, keep_inputs, sys.input_labels)
+    elim_outputs, keep_outputs = _process_elim_or_keep(
+        elim_outputs, keep_outputs, sys.output_labels)
+
+    # Create submatrix of states we are keeping
+    A11 = sys.A[:, keep_states][keep_states, :]     # states we are keeping
+    A12 = sys.A[:, elim_states][keep_states, :]     # needed for 'matchdc'
+    A21 = sys.A[:, keep_states][elim_states, :]
+    A22 = sys.A[:, elim_states][elim_states, :]
+
+    B1 = sys.B[keep_states, :]
+    B2 = sys.B[elim_states, :]
+
+    C1 = sys.C[:, keep_states]
+    C2 = sys.C[:, elim_states]
+
+    # Figure out the new state space system
+    if method == 'matchdc' and A22.size > 0:
+        if sys.isdtime(strict=True):
+            raise NotImplementedError(
+                "'matchdc' not (yet) supported for discrete time systems")
+        
+        # if matchdc, residualize
         # Check if the matrix A22 is invertible
-        if np.linalg.matrix_rank(A22) != len(ELIM):
+        if np.linalg.matrix_rank(A22) != len(elim_states):
             raise ValueError("Matrix A22 is singular to working precision.")
 
         # Now precompute A22\A21 and A22\B2 (A22I = inv(A22))
@@ -206,22 +231,29 @@ def model_reduction(sys, ELIM, method='matchdc'):
         Br = B1 - A12 @ A22I_B2
         Cr = C1 - C2 @ A22I_A21
         Dr = sys.D - C2 @ A22I_B2
-    elif method == 'truncate':
-        # if truncate, simply discard state x2
+
+    elif method == 'truncate' or A22.size == 0:
+        # Get rid of unwanted states
         Ar = A11
         Br = B1
         Cr = C1
         Dr = sys.D
+
     else:
         raise ValueError("Oops, method is not supported!")
 
+    # Get rid of additional inputs and outputs
+    Br = Br[:, keep_inputs]
+    Cr = Cr[keep_outputs, :]
+    Dr = Dr[keep_outputs, :][:, keep_inputs]
+
     rsys = StateSpace(Ar, Br, Cr, Dr)
     return rsys
 
 
 def balanced_reduction(sys, orders, method='truncate', alpha=None):
     """Balanced reduced order model of sys of a given order.
-    
+
     States are eliminated based on Hankel singular value.
     If sys has unstable modes, they are removed, the
     balanced realization is done on the stable part, then
@@ -276,7 +308,7 @@ def balanced_reduction(sys, orders, method='truncate', alpha=None):
         raise ValueError("supported methods are 'truncate' or 'matchdc'")
     elif method == 'truncate':
         try:
-            from slycot import ab09md, ab09ad
+            from slycot import ab09ad, ab09md
         except ImportError:
             raise ControlSlycot(
                 "can't find slycot subroutine ab09md or ab09ad")
@@ -346,7 +378,7 @@ def balanced_reduction(sys, orders, method='truncate', alpha=None):
 
 def minimal_realization(sys, tol=None, verbose=True):
     """ Eliminate uncontrollable or unobservable states.
-    
+
     Eliminates uncontrollable or unobservable states in state-space
     models or cancelling pole-zero pairs in transfer functions. The
     output sysr has minimal order and the same response
@@ -378,14 +410,14 @@ def _block_hankel(Y, m, n):
     """Create a block Hankel matrix from impulse response."""
     q, p, _ = Y.shape
     YY = Y.transpose(0, 2, 1) # transpose for reshape
-    
+
     H = np.zeros((q*m, p*n))
-    
+
     for r in range(m):
         # shift and add row to Hankel matrix
         new_row = YY[:, r:r+n, :]
         H[q*r:q*(r+1), :] = new_row.reshape((q, p*n))
-            
+
     return H
 
 
@@ -431,7 +463,7 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
         unit-area impulse response of python-control. Default is True.
     transpose : bool, optional
         Assume that input data is transposed relative to the standard
-        :ref:`time-series-convention`. For TimeResponseData this parameter 
+        :ref:`time-series-convention`. For TimeResponseData this parameter
         is ignored. Default is False.
 
     Returns
@@ -466,7 +498,7 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
         YY = np.array(arg, ndmin=3)
         if transpose:
             YY = np.transpose(YY)
-    
+
     q, p, l = YY.shape
 
     if m is None:
@@ -476,14 +508,14 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
 
     if m*q < r or n*p < r:
         raise ValueError("Hankel parameters are to small")
-    
+
     if (l-1) < m+n:
         raise ValueError("not enough data for requested number of parameters")
-    
+
     H = _block_hankel(YY[:, :, 1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
     Hf = H[:, :-p] # first p*n columns of H
     Hl = H[:, p:] # last p*n columns of H
-    
+
     U,S,Vh = np.linalg.svd(Hf, True)
     Ur =U[:, 0:r]
     Vhr =Vh[0:r, :]
@@ -500,7 +532,7 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
 
 def markov(*args, m=None, transpose=False, dt=None, truncate=False):
     """markov(Y, U, [, m])
-    
+
     Calculate the first `m` Markov parameters [D CB CAB ...] from data.
 
     This function computes the Markov parameters for a discrete time
@@ -579,7 +611,7 @@ def markov(*args, m=None, transpose=False, dt=None, truncate=False):
     # Get the system description
     if len(args) < 1:
         raise ControlArgument("not enough input arguments")
-    
+
     if isinstance(args[0], TimeResponseData):
         data = args[0]
         Umat = np.array(data.inputs, ndmin=2)
@@ -639,7 +671,7 @@ def markov(*args, m=None, transpose=False, dt=None, truncate=False):
     # This algorithm sets up the following problem and solves it for
     # the Markov parameters
     #
-    # (l,q)   = (l,p*m) @ (p*m,q) 
+    # (l,q)   = (l,p*m) @ (p*m,q)
     # YY      = UU @ H.T
     #
     # [ y(0)   ]   [ u(0)    0       0                 ] [ D           ]
@@ -675,7 +707,7 @@ def markov(*args, m=None, transpose=False, dt=None, truncate=False):
     # Truncate first t=0 or t=m time steps, transpose the problem for lsq
     YY = Ymat[:, t:].T
     UU = UUT[:, t:].T
-    
+
     # Solve for the Markov parameters from  YY = UU @ H.T
     HT, _, _, _ = np.linalg.lstsq(UU, YY, rcond=None)
     H = HT.T/dt # scaling
diff --git a/control/statefbk.py b/control/statefbk.py
index 16eeb36ee..1dd0be325 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -39,24 +39,25 @@
 #
 # $Id$
 
-# External packages and modules
+import warnings
+
 import numpy as np
 import scipy as sp
-import warnings
 
 from . import statesp
-from .mateqn import care, dare, _check_shape
-from .statesp import StateSpace, _ssmatrix, _convert_to_statespace, ss
+from .config import _process_legacy_keyword
+from .exception import ControlArgument, ControlDimension, \
+    ControlNotImplemented, ControlSlycot
+from .iosys import _process_indices, _process_labels, isctime, isdtime
 from .lti import LTI
-from .iosys import isdtime, isctime, _process_indices, _process_labels
+from .mateqn import _check_shape, care, dare
 from .nlsys import NonlinearIOSystem, interconnect
-from .exception import ControlSlycot, ControlArgument, ControlDimension, \
-    ControlNotImplemented
-from .config import _process_legacy_keyword
+from .statesp import StateSpace, _convert_to_statespace, _ssmatrix, ss
 
 # Make sure we have access to the right slycot routines
 try:
     from slycot import sb03md57
+
     # wrap without the deprecation warning
     def sb03md(n, C, A, U, dico, job='X',fact='N',trana='N',ldwork=None):
         ret = sb03md57(A, U, C, dico, job, fact, trana, ldwork)
@@ -581,8 +582,9 @@ def dlqr(*args, **kwargs):
 
 # Function to create an I/O sytems representing a state feedback controller
 def create_statefbk_iosystem(
-        sys, gain, integral_action=None, estimator=None, controller_type=None,
-        xd_labels=None, ud_labels=None, gainsched_indices=None,
+        sys, gain, feedfwd_gain=None, integral_action=None, estimator=None,
+        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):
     r"""Create an I/O system using a (full) state feedback controller.
@@ -592,7 +594,7 @@ def create_statefbk_iosystem(
 
     .. math:: u = u_d - K_p (x - x_d) - K_i \int(C x - C x_d)
 
-    It can be called in the form::
+    by calling
 
         ctrl, clsys = ct.create_statefbk_iosystem(sys, K)
 
@@ -608,6 +610,18 @@ def create_statefbk_iosystem(
 
     where :math:`\mu` represents the scheduling variable.
 
+    Alternatively, a controller of the form
+
+    .. math:: u = k_f r - K_p x - K_i \int(C x - r)
+
+    can be created by calling
+
+        ctrl, clsys = ct.create_statefbk_iosystem(
+            sys, K, kf, feedfwd_pattern='refgain')
+
+    In either form, an estimator can also be used to compute the estimated
+    state from the input and output measurements.
+
     Parameters
     ----------
     sys : NonlinearIOSystem
@@ -631,14 +645,15 @@ def create_statefbk_iosystem(
         If an I/O system is given, the error e = x - xd is passed to the
         system and the output is used as the feedback compensation term.
 
-    xd_labels, ud_labels : str or list of str, optional
-        Set the name of the signals to use for the desired state and
-        inputs.  If a single string is specified, it should be a format
-        string using the variable `i` as an index.  Otherwise, a list of
-        strings matching the size of `x_d` and `u_d`, respectively, should
-        be used.  Default is "xd[{i}]" for xd_labels and "ud[{i}]" for
-        ud_labels.  These settings can also be overridden using the
-        `inputs` keyword.
+    feedfwd_gain : array_like, optional
+        Specify the feedforward gain, `k_f`.  Used only for the reference
+        gain design pattern.  If not given and if `sys` is a `StateSpace`
+        (linear) system, will be computed as -1/(C (A-BK)^{-1}) B.
+
+    feedfwd_pattern : str, optional
+        If set to 'refgain', the reference gain design pattern is used to
+        create the controller instead of the trajectory generation
+        ('trajgen') pattern.
 
     integral_action : ndarray, optional
         If this keyword is specified, the controller can include integral
@@ -680,17 +695,19 @@ def create_statefbk_iosystem(
     Returns
     -------
     ctrl : NonlinearIOSystem
-        Input/output system representing the controller.  This system
-        takes as inputs the desired state `x_d`, the desired input
-        `u_d`, and either the system state `x` or the estimated state
-        `xhat`.  It outputs the controller action `u` according to the
-        formula `u = u_d - K(x - x_d)`.  If the keyword
-        `integral_action` is specified, then an additional set of
-        integrators is included in the control system (with the gain
-        matrix `K` having the integral gains appended after the state
-        gains).  If a gain scheduled controller is specified, the gain
-        (proportional and integral) are evaluated using the scheduling
-        variables specified by `gainsched_indices`.
+        Input/output system representing the controller.  For the 'trajgen'
+        design pattern (default), this system takes as inputs the desired
+        state `x_d`, the desired input `u_d`, and either the system state
+        `x` or the estimated state `xhat`.  It outputs the controller
+        action `u` according to the formula `u = u_d - K(x - x_d)`.  For
+        the 'refgain' design pattern, the system takes as inputs the
+        reference input `r` and the system or estimated state. If the
+        keyword `integral_action` is specified, then an additional set of
+        integrators is included in the control system (with the gain matrix
+        `K` having the integral gains appended after the state gains).  If
+        a gain scheduled controller is specified, the gain (proportional
+        and integral) are evaluated using the scheduling variables
+        specified by `gainsched_indices`.
 
     clsys : NonlinearIOSystem
         Input/output system representing the closed loop system.  This
@@ -716,6 +733,15 @@ def create_statefbk_iosystem(
         specified as either integer offsets or as estimator/system output
         signal names.  If not specified, defaults to the system states.
 
+    xd_labels, ud_labels, ref_labels : str or list of str, optional
+        Set the name of the signals to use for the desired state and inputs
+        or the reference inputs (for the 'refgain' design pattern).  If a
+        single string is specified, it should be a format string using the
+        variable `i` as an index.  Otherwise, a list of strings matching
+        the size of `x_d` and `u_d`, respectively, should be used.  Default
+        is "xd[{i}]" for xd_labels and "ud[{i}]" for ud_labels.  These
+        settings can also be overridden using the `inputs` keyword.
+
     inputs, outputs, states : str, or list of str, optional
         List of strings that name the individual signals of the transformed
         system.  If not given, the inputs, outputs, and states are the same
@@ -753,6 +779,11 @@ def create_statefbk_iosystem(
     if kwargs:
         raise TypeError("unrecognized keywords: ", str(kwargs))
 
+    # Check for consistency of positional parameters
+    if feedfwd_gain is not None and feedfwd_pattern != 'refgain':
+        raise ControlArgument(
+            "feedfwd_gain specified but feedfwd_pattern != 'refgain'")
+
     # Figure out what inputs to the system to use
     control_indices = _process_indices(
         control_indices, 'control', sys.input_labels, sys.ninputs)
@@ -812,6 +843,10 @@ def create_statefbk_iosystem(
         # Check for gain scheduled controller
         if len(gain) != 2:
             raise ControlArgument("gain must be a 2-tuple for gain scheduling")
+        elif feedfwd_pattern != 'trajgen':
+            raise NotImplementedError(
+                "Gain scheduling is not implemented for pattern "
+                f"'{feedfwd_pattern}'")
         gains, points = gain[0:2]
 
         # Stack gains and points if past as a list
@@ -819,7 +854,7 @@ def create_statefbk_iosystem(
         points = np.stack(points)
         gainsched = True
 
-    elif isinstance(gain, NonlinearIOSystem):
+    elif isinstance(gain, NonlinearIOSystem) and feedfwd_pattern != 'refgain':
         if controller_type not in ['iosystem', None]:
             raise ControlArgument(
                 f"incompatible controller type '{controller_type}'")
@@ -841,20 +876,29 @@ def create_statefbk_iosystem(
         raise ControlArgument(f"unknown controller_type '{controller_type}'")
 
     # Figure out the labels to use
-    xd_labels = _process_labels(
-        xd_labels, 'xd', ['xd[{i}]'.format(i=i) for i in range(sys_nstates)])
-    ud_labels = _process_labels(
-        ud_labels, 'ud', ['ud[{i}]'.format(i=i) for i in range(sys_ninputs)])
-
-    # Create the signal and system names
-    if inputs is None:
-        inputs = xd_labels + ud_labels + estimator.output_labels
+    if feedfwd_pattern == 'trajgen':
+        xd_labels = _process_labels(xd_labels, 'xd', [
+            'xd[{i}]'.format(i=i) for i in range(sys_nstates)])
+        ud_labels = _process_labels(ud_labels, 'ud', [
+            'ud[{i}]'.format(i=i) for i in range(sys_ninputs)])
+
+        # Create the signal and system names
+        if inputs is None:
+            inputs = xd_labels + ud_labels + estimator.output_labels
+    elif feedfwd_pattern == 'refgain':
+        ref_labels = _process_labels(ref_labels, 'r', [
+            f'r[{i}]' for i in range(sys_ninputs)])
+        if inputs is None:
+            inputs = ref_labels + estimator.output_labels
+    else:
+        raise NotImplementedError(f"unknown pattern '{feedfwd_pattern}'")
+
     if outputs is None:
         outputs = [sys.input_labels[i] for i in control_indices]
     if states is None:
         states = nintegrators
 
-    # Process gainscheduling variables, if present
+    # Process gain scheduling variables, if present
     if gainsched:
         # Create a copy of the scheduling variable indices (default = xd)
         gainsched_indices = _process_indices(
@@ -897,7 +941,7 @@ def _compute_gain(mu):
             return K
 
     # Define the controller system
-    if controller_type == 'nonlinear':
+    if controller_type == 'nonlinear' and feedfwd_pattern == 'trajgen':
         # Create an I/O system for the state feedback gains
         def _control_update(t, states, inputs, params):
             # Split input into desired state, nominal input, and current state
@@ -931,7 +975,7 @@ def _control_output(t, states, inputs, params):
             _control_update, _control_output, name=name, inputs=inputs,
             outputs=outputs, states=states, params=params)
 
-    elif controller_type == 'iosystem':
+    elif controller_type == 'iosystem' and feedfwd_pattern == 'trajgen':
         # Use the passed system to compute feedback compensation
         def _control_update(t, states, inputs, params):
             # Split input into desired state, nominal input, and current state
@@ -955,7 +999,7 @@ def _control_output(t, states, inputs, params):
             _control_update, _control_output, name=name, inputs=inputs,
             outputs=outputs, states=fbkctrl.state_labels, dt=fbkctrl.dt)
 
-    elif controller_type == 'linear' or controller_type is None:
+    elif controller_type in 'linear' and feedfwd_pattern == 'trajgen':
         # Create the matrices implementing the controller
         if isctime(sys):
             # Continuous time: integrator
@@ -973,6 +1017,37 @@ def _control_output(t, states, inputs, params):
             A_lqr, B_lqr, C_lqr, D_lqr, dt=sys.dt, name=name,
             inputs=inputs, outputs=outputs, states=states)
 
+    elif feedfwd_pattern == 'refgain':
+        if controller_type not in ['linear', 'iosystem']:
+            raise ControlArgument(
+                "refgain design pattern only supports linear controllers")
+
+        if feedfwd_gain is None:
+            raise ControlArgument(
+                "'feedfwd_gain' required for reference gain pattern")
+
+        # Check to make sure the reference gain is valid
+        Kf = np.atleast_2d(feedfwd_gain)
+        if Kf.ndim != 2 or Kf.shape[0] != sys.ninputs or \
+           Kf.shape[1] != sys.ninputs:
+            raise ControlArgument("feedfwd_gain is not the right shape")
+
+        # Create the matrices implementing the controller
+        # [r, x]->[u]: u = k_f r - K_p x - K_i \int(C x - r)
+        if isctime(sys):
+            # Continuous time: integrator
+            A_lqr = np.zeros((C.shape[0], C.shape[0]))
+        else:
+            # Discrete time: summer
+            A_lqr = np.eye(C.shape[0])
+        B_lqr = np.hstack([-np.eye(C.shape[0], sys_ninputs), C])
+        C_lqr = -K[:, sys_nstates:]                     # integral gain (opt)
+        D_lqr = np.hstack([Kf, -K])
+
+        ctrl = ss(
+            A_lqr, B_lqr, C_lqr, D_lqr, dt=sys.dt, name=name,
+            inputs=inputs, outputs=outputs, states=states)
+
     else:
         raise ControlArgument(f"unknown controller_type '{controller_type}'")
 
@@ -1020,7 +1095,7 @@ def ctrb(A, B, t=None):
     bmat = _ssmatrix(B)
     n = np.shape(amat)[0]
     m = np.shape(bmat)[1]
-    
+
     if t is None or t > n:
         t = n
 
@@ -1042,7 +1117,7 @@ def obsv(A, C, t=None):
         Dynamics and output matrix of the system
     t : None or integer
         maximum time horizon of the controllability matrix, max = A.shape[0]
-        
+
     Returns
     -------
     O : 2D array (or matrix)
@@ -1062,14 +1137,14 @@ def obsv(A, C, t=None):
     cmat = _ssmatrix(C)
     n = np.shape(amat)[0]
     p = np.shape(cmat)[0]
-    
+
     if t is None or t > n:
         t = n
 
     # Construct the observability matrix
     obsv = np.zeros((t * p, n))
     obsv[:p, :] = cmat
-        
+
     for k in range(1, t):
         obsv[k * p:(k + 1) * p, :] = np.dot(obsv[(k - 1) * p:k * p, :], amat)
 
diff --git a/control/statesp.py b/control/statesp.py
index bfe5f996b..7af9008f4 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -252,6 +252,12 @@ def __init__(self, *args, **kwargs):
             D = np.zeros((C.shape[0], B.shape[1]))
         D = _ssmatrix(D)
 
+        # If only direct term is present, adjust sizes of C and D if needed
+        if D.size > 0 and B.size == 0:
+            B = np.zeros((0, D.shape[1]))
+        if D.size > 0 and C.size == 0:
+            C = np.zeros((D.shape[0], 0))
+
         # Matrices defining the linear system
         self.A = A
         self.B = B
@@ -268,7 +274,7 @@ def __init__(self, *args, **kwargs):
 
         # Process iosys keywords
         defaults = args[0] if len(args) == 1 else \
-            {'inputs': D.shape[1], 'outputs': D.shape[0],
+            {'inputs': B.shape[1], 'outputs': C.shape[0],
              'states': A.shape[0]}
         name, inputs, outputs, states, dt = _process_iosys_keywords(
             kwargs, defaults, static=(A.size == 0))
@@ -295,16 +301,15 @@ def __init__(self, *args, **kwargs):
         # Check to make sure everything is consistent
         #
         # Check that the matrix sizes are consistent
-        if A.shape[0] != A.shape[1] or self.nstates != A.shape[0]:
-            raise ValueError("A must be square.")
-        if self.nstates != B.shape[0]:
-            raise ValueError("A and B must have the same number of rows.")
-        if self.nstates != C.shape[1]:
-            raise ValueError("A and C must have the same number of columns.")
-        if self.ninputs != B.shape[1] or self.ninputs != D.shape[1]:
-            raise ValueError("B and D must have the same number of columns.")
-        if self.noutputs != C.shape[0] or self.noutputs != D.shape[0]:
-            raise ValueError("C and D must have the same number of rows.")
+        def _check_shape(matrix, expected, name):
+            if matrix.shape != expected:
+                raise ValueError(
+                    f"{name} is the wrong shape; "
+                    f"expected {expected} instead of {matrix.shape}")
+        _check_shape(A, (self.nstates, self.nstates), "A")
+        _check_shape(B, (self.nstates, self.ninputs), "B")
+        _check_shape(C, (self.noutputs, self.nstates), "C")
+        _check_shape(D, (self.noutputs, self.ninputs), "D")
 
         #
         # Final processing
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 56ec7395f..5304ea636 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -634,11 +634,11 @@ def __mul__(self, other):
         from .statesp import StateSpace
 
         # Convert the second argument to a transfer function.
-        if isinstance(other, StateSpace):
+        if isinstance(other, (StateSpace, np.ndarray)):
             other = _convert_to_transfer_function(other)
-        elif isinstance(other, (int, float, complex, np.number, np.ndarray)):
-            other = _convert_to_transfer_function(other, inputs=self.ninputs,
-                                                  outputs=self.noutputs)
+        elif isinstance(other, (int, float, complex, np.number)):
+            # Multiply by a scaled identity matrix (transfer function)
+            other = _convert_to_transfer_function(np.eye(self.ninputs) * other)
         if not isinstance(other, TransferFunction):
             return NotImplemented
 
@@ -681,8 +681,8 @@ def __rmul__(self, other):
 
         # Convert the second argument to a transfer function.
         if isinstance(other, (int, float, complex, np.number)):
-            other = _convert_to_transfer_function(other, inputs=self.ninputs,
-                                                  outputs=self.ninputs)
+            # Multiply by a scaled identity matrix (transfer function)
+            other = _convert_to_transfer_function(np.eye(self.noutputs) * other)
         else:
             other = _convert_to_transfer_function(other)
 
@@ -723,9 +723,8 @@ def __truediv__(self, other):
         """Divide two LTI objects."""
 
         if isinstance(other, (int, float, complex, np.number)):
-            other = _convert_to_transfer_function(
-                other, inputs=self.ninputs,
-                outputs=self.ninputs)
+            # Multiply by a scaled identity matrix (transfer function)
+            other = _convert_to_transfer_function(np.eye(self.ninputs) * other)
         else:
             other = _convert_to_transfer_function(other)
 
diff --git a/doc/iosys.rst b/doc/iosys.rst
index f2dfbff4d..0554683f3 100644
--- a/doc/iosys.rst
+++ b/doc/iosys.rst
@@ -57,7 +57,7 @@ Example
 
 To illustrate the use of the input/output systems module, we create a
 model for a predator/prey system, following the notation and parameter
-values in FBS2e.
+values in `Feedback Systems <http://fbsbook.org>`_.
 
 We begin by defining the dynamics of the system
 
@@ -129,7 +129,8 @@ system and computing the linearization about that point.
 
 We next compute a controller that stabilizes the equilibrium point using
 eigenvalue placement and computing the feedforward gain using the number of
-lynxes as the desired output (following FBS2e, Example 7.5):
+lynxes as the desired output (following `Feedback Systems
+<http://fbsbook.org>`_, Example 7.5):
 
 .. code-block:: python
 
@@ -470,6 +471,29 @@ closed loop systems `clsys`, both as I/O systems.  The input to the
 controller is the vector of desired states :math:`x_\text{d}`, desired
 inputs :math:`u_\text{d}`, and system states :math:`x`.
 
+The above design pattern is referred to as the "trajectory generation"
+('trajgen') pattern, since it assumes that the input to the controller is a
+feasible trajectory :math:`(x_\text{d}, u_\text{d})`.  Alternatively, a
+controller using the "reference gain" pattern can be created, which
+implements a state feedback controller of the form
+
+.. math::
+
+  u = k_\text{f}\, r - K x,
+
+where :math:`r` is the reference input and :math:`k_\text{f}` is the
+feedforward gain (normally chosen so that the steady state output
+:math:`y_\text{ss}` will be equal to :math:`r`).
+
+A reference gain controller can be created with the command::
+
+  ctrl, clsys = ct.create_statefbk_iosystem(sys, K, kf, feedfwd_pattern='refgain')
+
+This reference gain design pattern is described in more detail in
+Section 7.2 of `Feedback Systems <http://fbsbook.org>`_ (Stabilization
+by State Feedback) and the trajectory generation design pattern is
+described in Section 8.5 (State Space Controller Design).
+
 If the full system state is not available, the output of a state
 estimator can be used to construct the controller using the command::
 
@@ -507,10 +531,11 @@ must match the number of additional columns in the `K` matrix.  If an
 estimator is specified, :math:`\hat x` will be used in place of
 :math:`x`.
 
-Finally, gain scheduling on the desired state, desired input, or
-system state can be implemented by setting the gain to a 2-tuple
-consisting of a list of gains and a list of points at which the gains
-were computed, as well as a description of the scheduling variables::
+Finally, for the trajectory generation design pattern, gain scheduling on
+the desired state, desired input, or system state can be implemented by
+setting the gain to a 2-tuple consisting of a list of gains and a list of
+points at which the gains were computed, as well as a description of the
+scheduling variables::
 
   ctrl, clsys = ct.create_statefbk_iosystem(
       sys, ([g1, ..., gN], [p1, ..., pN]), gainsched_indices=[s1, ..., sq])


Test patch:
diff --git a/control/tests/bdalg_test.py b/control/tests/bdalg_test.py
index 5629f27f9..8ea67e0f7 100644
--- a/control/tests/bdalg_test.py
+++ b/control/tests/bdalg_test.py
@@ -12,6 +12,7 @@
 from control.statesp import StateSpace
 from control.bdalg import feedback, append, connect
 from control.lti import zeros, poles
+from control.bdalg import _ensure_tf
 
 
 class TestFeedback:
@@ -362,3 +363,557 @@ def test_bdalg_udpate_names_errors():
 
     with pytest.raises(TypeError, match="unrecognized keywords"):
         sys = ctrl.series(sys1, sys2, dt=1)
+
+
+class TestEnsureTf:
+    """Test ``_ensure_tf``."""
+
+    @pytest.mark.parametrize(
+        "arraylike_or_tf, dt, tf",
+        [
+            (
+                ctrl.TransferFunction([1], [1, 2, 3]),
+                None,
+                ctrl.TransferFunction([1], [1, 2, 3]),
+            ),
+            (
+                ctrl.TransferFunction([1], [1, 2, 3]),
+                0,
+                ctrl.TransferFunction([1], [1, 2, 3]),
+            ),
+            (
+                2,
+                None,
+                ctrl.TransferFunction([2], [1]),
+            ),
+            (
+                np.array([2]),
+                None,
+                ctrl.TransferFunction([2], [1]),
+            ),
+            (
+                np.array([[2]]),
+                None,
+                ctrl.TransferFunction([2], [1]),
+            ),
+            (
+                np.array(
+                    [
+                        [2, 0, 3],
+                        [1, 2, 3],
+                    ]
+                ),
+                None,
+                ctrl.TransferFunction(
+                    [
+                        [[2], [0], [3]],
+                        [[1], [2], [3]],
+                    ],
+                    [
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                    ],
+                ),
+            ),
+            (
+                np.array([2, 0, 3]),
+                None,
+                ctrl.TransferFunction(
+                    [
+                        [[2], [0], [3]],
+                    ],
+                    [
+                        [[1], [1], [1]],
+                    ],
+                ),
+            ),
+        ],
+    )
+    def test_ensure(self, arraylike_or_tf, dt, tf):
+        """Test nominal cases"""
+        ensured_tf = _ensure_tf(arraylike_or_tf, dt)
+        assert _tf_close_coeff(tf, ensured_tf)
+
+    @pytest.mark.parametrize(
+        "arraylike_or_tf, dt, exception",
+        [
+            (
+                ctrl.TransferFunction([1], [1, 2, 3]),
+                0.1,
+                ValueError,
+            ),
+            (
+                ctrl.TransferFunction([1], [1, 2, 3], 0.1),
+                0,
+                ValueError,
+            ),
+            (
+                np.ones((1, 1, 1)),
+                None,
+                ValueError,
+            ),
+            (
+                np.ones((1, 1, 1, 1)),
+                None,
+                ValueError,
+            ),
+        ],
+    )
+    def test_error_ensure(self, arraylike_or_tf, dt, exception):
+        """Test error cases"""
+        with pytest.raises(exception):
+            _ensure_tf(arraylike_or_tf, dt)
+
+
+class TestTfCombineSplit:
+    """Test ``combine_tf`` and ``split_tf``."""
+
+    @pytest.mark.parametrize(
+        "tf_array, tf",
+        [
+            # Continuous-time
+            (
+                [
+                    [ctrl.TransferFunction([1], [1, 1])],
+                    [ctrl.TransferFunction([2], [1, 0])],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[1]],
+                        [[2]],
+                    ],
+                    [
+                        [[1, 1]],
+                        [[1, 0]],
+                    ],
+                ),
+            ),
+            # Discrete-time
+            (
+                [
+                    [ctrl.TransferFunction([1], [1, 1], dt=1)],
+                    [ctrl.TransferFunction([2], [1, 0], dt=1)],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[1]],
+                        [[2]],
+                    ],
+                    [
+                        [[1, 1]],
+                        [[1, 0]],
+                    ],
+                    dt=1,
+                ),
+            ),
+            # Scalar
+            (
+                [
+                    [2],
+                    [ctrl.TransferFunction([2], [1, 0])],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[2]],
+                        [[2]],
+                    ],
+                    [
+                        [[1]],
+                        [[1, 0]],
+                    ],
+                ),
+            ),
+            # Matrix
+            (
+                [
+                    [np.eye(3)],
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [0], [3]],
+                                [[1], [2], [3]],
+                            ],
+                            [
+                                [[1], [1], [1]],
+                                [[1], [1], [1]],
+                            ],
+                        )
+                    ],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[1], [0], [0]],
+                        [[0], [1], [0]],
+                        [[0], [0], [1]],
+                        [[2], [0], [3]],
+                        [[1], [2], [3]],
+                    ],
+                    [
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                    ],
+                ),
+            ),
+            # Inhomogeneous
+            (
+                [
+                    [np.eye(3)],
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [0]],
+                                [[1], [2]],
+                            ],
+                            [
+                                [[1], [1]],
+                                [[1], [1]],
+                            ],
+                        ),
+                        ctrl.TransferFunction(
+                            [
+                                [[3]],
+                                [[3]],
+                            ],
+                            [
+                                [[1]],
+                                [[1]],
+                            ],
+                        ),
+                    ],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[1], [0], [0]],
+                        [[0], [1], [0]],
+                        [[0], [0], [1]],
+                        [[2], [0], [3]],
+                        [[1], [2], [3]],
+                    ],
+                    [
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                    ],
+                ),
+            ),
+            # Discrete-time
+            (
+                [
+                    [2],
+                    [ctrl.TransferFunction([2], [1, 0], dt=0.1)],
+                ],
+                ctrl.TransferFunction(
+                    [
+                        [[2]],
+                        [[2]],
+                    ],
+                    [
+                        [[1]],
+                        [[1, 0]],
+                    ],
+                    dt=0.1,
+                ),
+            ),
+        ],
+    )
+    def test_combine_tf(self, tf_array, tf):
+        """Test combining transfer functions."""
+        tf_combined = ctrl.combine_tf(tf_array)
+        assert _tf_close_coeff(tf_combined, tf)
+
+    @pytest.mark.parametrize(
+        "tf_array, tf",
+        [
+            (
+                np.array(
+                    [
+                        [ctrl.TransferFunction([1], [1, 1])],
+                    ],
+                    dtype=object,
+                ),
+                ctrl.TransferFunction(
+                    [
+                        [[1]],
+                    ],
+                    [
+                        [[1, 1]],
+                    ],
+                ),
+            ),
+            (
+                np.array(
+                    [
+                        [ctrl.TransferFunction([1], [1, 1])],
+                        [ctrl.TransferFunction([2], [1, 0])],
+                    ],
+                    dtype=object,
+                ),
+                ctrl.TransferFunction(
+                    [
+                        [[1]],
+                        [[2]],
+                    ],
+                    [
+                        [[1, 1]],
+                        [[1, 0]],
+                    ],
+                ),
+            ),
+            (
+                np.array(
+                    [
+                        [ctrl.TransferFunction([1], [1, 1], dt=1)],
+                        [ctrl.TransferFunction([2], [1, 0], dt=1)],
+                    ],
+                    dtype=object,
+                ),
+                ctrl.TransferFunction(
+                    [
+                        [[1]],
+                        [[2]],
+                    ],
+                    [
+                        [[1, 1]],
+                        [[1, 0]],
+                    ],
+                    dt=1,
+                ),
+            ),
+            (
+                np.array(
+                    [
+                        [ctrl.TransferFunction([2], [1], dt=0.1)],
+                        [ctrl.TransferFunction([2], [1, 0], dt=0.1)],
+                    ],
+                    dtype=object,
+                ),
+                ctrl.TransferFunction(
+                    [
+                        [[2]],
+                        [[2]],
+                    ],
+                    [
+                        [[1]],
+                        [[1, 0]],
+                    ],
+                    dt=0.1,
+                ),
+            ),
+        ],
+    )
+    def test_split_tf(self, tf_array, tf):
+        """Test splitting transfer functions."""
+        tf_split = ctrl.split_tf(tf)
+        # Test entry-by-entry
+        for i in range(tf_split.shape[0]):
+            for j in range(tf_split.shape[1]):
+                assert _tf_close_coeff(
+                    tf_split[i, j],
+                    tf_array[i, j],
+                )
+        # Test combined
+        assert _tf_close_coeff(
+            ctrl.combine_tf(tf_split),
+            ctrl.combine_tf(tf_array),
+        )
+
+    @pytest.mark.parametrize(
+        "tf_array, exception",
+        [
+            # Wrong timesteps
+            (
+                [
+                    [ctrl.TransferFunction([1], [1, 1], 0.1)],
+                    [ctrl.TransferFunction([2], [1, 0], 0.2)],
+                ],
+                ValueError,
+            ),
+            (
+                [
+                    [ctrl.TransferFunction([1], [1, 1], 0.1)],
+                    [ctrl.TransferFunction([2], [1, 0], 0)],
+                ],
+                ValueError,
+            ),
+            # Too few dimensions
+            (
+                [
+                    ctrl.TransferFunction([1], [1, 1]),
+                    ctrl.TransferFunction([2], [1, 0]),
+                ],
+                ValueError,
+            ),
+            # Too many dimensions
+            (
+                [
+                    [[ctrl.TransferFunction([1], [1, 1], 0.1)]],
+                    [[ctrl.TransferFunction([2], [1, 0], 0)]],
+                ],
+                ValueError,
+            ),
+            # Incompatible dimensions
+            (
+                [
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [
+                                    [1],
+                                ]
+                            ],
+                            [
+                                [
+                                    [1, 1],
+                                ]
+                            ],
+                        ),
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1]],
+                                [[1], [3]],
+                            ],
+                            [
+                                [[1, 0], [1, 0]],
+                                [[1, 0], [1, 0]],
+                            ],
+                        ),
+                    ],
+                ],
+                ValueError,
+            ),
+            (
+                [
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1]],
+                                [[1], [3]],
+                            ],
+                            [
+                                [[1, 0], [1, 0]],
+                                [[1, 0], [1, 0]],
+                            ],
+                        ),
+                        ctrl.TransferFunction(
+                            [
+                                [
+                                    [1],
+                                ]
+                            ],
+                            [
+                                [
+                                    [1, 1],
+                                ]
+                            ],
+                        ),
+                    ],
+                ],
+                ValueError,
+            ),
+            (
+                [
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1]],
+                                [[1], [3]],
+                            ],
+                            [
+                                [[1, 0], [1, 0]],
+                                [[1, 0], [1, 0]],
+                            ],
+                        ),
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1]],
+                                [[1], [3]],
+                            ],
+                            [
+                                [[1, 0], [1, 0]],
+                                [[1, 0], [1, 0]],
+                            ],
+                        ),
+                    ],
+                    [
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1], [1]],
+                                [[1], [3], [2]],
+                            ],
+                            [
+                                [[1, 0], [1, 0], [1, 0]],
+                                [[1, 0], [1, 0], [1, 0]],
+                            ],
+                        ),
+                        ctrl.TransferFunction(
+                            [
+                                [[2], [1]],
+                                [[1], [3]],
+                            ],
+                            [
+                                [[1, 0], [1, 0]],
+                                [[1, 0], [1, 0]],
+                            ],
+                        ),
+                    ],
+                ],
+                ValueError,
+            ),
+        ],
+    )
+    def test_error_combine_tf(self, tf_array, exception):
+        """Test error cases."""
+        with pytest.raises(exception):
+            ctrl.combine_tf(tf_array)
+
+
+def _tf_close_coeff(tf_a, tf_b, rtol=1e-5, atol=1e-8):
+    """Check if two transfer functions have close coefficients.
+
+    Parameters
+    ----------
+    tf_a : TransferFunction
+        First transfer function.
+    tf_b : TransferFunction
+        Second transfer function.
+    rtol : float
+        Relative tolerance for ``np.allclose``.
+    atol : float
+        Absolute tolerance for ``np.allclose``.
+
+    Returns
+    -------
+    bool
+        True if transfer function cofficients are all close.
+    """
+    # Check number of outputs and inputs
+    if tf_a.noutputs != tf_b.noutputs:
+        return False
+    if tf_a.ninputs != tf_b.ninputs:
+        return False
+    # Check timestep
+    if tf_a.dt != tf_b.dt:
+        return False
+    # Check coefficient arrays
+    for i in range(tf_a.noutputs):
+        for j in range(tf_a.ninputs):
+            if not np.allclose(
+                tf_a.num[i][j],
+                tf_b.num[i][j],
+                rtol=rtol,
+                atol=atol,
+            ):
+                return False
+            if not np.allclose(
+                tf_a.den[i][j],
+                tf_b.den[i][j],
+                rtol=rtol,
+                atol=atol,
+            ):
+                return False
+    return True
diff --git a/control/tests/docstrings_test.py b/control/tests/docstrings_test.py
index 27ced105f..991ead3e5 100644
--- a/control/tests/docstrings_test.py
+++ b/control/tests/docstrings_test.py
@@ -39,7 +39,7 @@
     control.ss2tf:                      '48ff25d22d28e7b396e686dd5eb58831',
     control.tf:                         '53a13f4a7f75a31c81800e10c88730ef',
     control.tf2ss:                      '086a3692659b7321c2af126f79f4bc11',
-    control.markov:                     'eda7c4635bbb863ae6659e574285d356',
+    control.markov:                     'a4199c54cb50f07c0163d3790739eafe',
     control.gangof4:                    '0e52eb6cf7ce024f9a41f3ae3ebf04f7',
 }
 
diff --git a/control/tests/lti_test.py b/control/tests/lti_test.py
index 3f001c17b..5359ceea3 100644
--- a/control/tests/lti_test.py
+++ b/control/tests/lti_test.py
@@ -350,3 +350,21 @@ def test_subsys_indexing(fcn, outdx, inpdx, key):
             np.testing.assert_almost_equal(
                 subsys_fcn.frequency_response(omega).response,
                 subsys_chk.frequency_response(omega).response)
+
+
+@slycotonly
+@pytest.mark.parametrize("op", [
+    '__mul__', '__rmul__', '__add__', '__radd__', '__sub__', '__rsub__'])
+@pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.frd])
+def test_scalar_algebra(op, fcn):
+    sys_ss = ct.rss(4, 2, 2)
+    match fcn:
+        case ct.ss:
+            sys = sys_ss
+        case ct.tf:
+            sys = ct.tf(sys_ss)
+        case ct.frd:
+            sys = ct.frd(sys_ss, [0.1, 1, 10])
+
+    scaled = getattr(sys, op)(2)
+    np.testing.assert_almost_equal(getattr(sys(1j), op)(2), scaled(1j))
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
index 616ef5f09..dead4eb75 100644
--- a/control/tests/modelsimp_test.py
+++ b/control/tests/modelsimp_test.py
@@ -3,11 +3,15 @@
 RMM, 30 Mar 2011 (based on TestModelSimp from v0.4a)
 """
 
+import math
+import warnings
+
 import numpy as np
 import pytest
 
+import control as ct
 from control import StateSpace, TimeResponseData, c2d, forced_response, \
-    impulse_response, step_response, rss, tf
+    impulse_response, rss, step_response, tf
 from control.exception import ControlArgument, ControlDimension
 from control.modelsimp import balred, eigensys_realization, hsvd, markov, \
     modred
@@ -42,7 +46,7 @@ def testMarkovSignature(self):
                                     inputs=U,
                                     input_labels='u',
                                     )
-        
+
         # setup
         m = 3
         Htrue = np.array([1., 0., 0.])
@@ -101,7 +105,6 @@ def testMarkovSignature(self):
         HT = markov(response, m)
         np.testing.assert_array_almost_equal(HT, np.transpose(Htrue))
         response.transpose=False
-        
 
         # Test example from docstring
         # TODO: There is a problem here, last markov parameter does not fit
@@ -212,7 +215,7 @@ def testMarkovResults(self, k, m, n):
         Mtrue = np.hstack([Hd.D] + [
             Hd.C @ np.linalg.matrix_power(Hd.A, i) @ Hd.B
             for i in range(m-1)])
-        
+
         Mtrue = np.squeeze(Mtrue)
 
         # Generate input/output data
@@ -238,8 +241,8 @@ def testMarkovResults(self, k, m, n):
         Mcomp_scaled = markov(response, m, dt=Ts)
 
         np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
-        np.testing.assert_allclose(Mtrue_scaled, Mcomp_scaled, rtol=1e-6, atol=1e-8)
-        
+        np.testing.assert_allclose(
+            Mtrue_scaled, Mcomp_scaled, rtol=1e-6, atol=1e-8)
 
     def testERASignature(self):
 
@@ -253,7 +256,7 @@ def testERASignature(self):
         B = np.array([[1.],[1.,]])
         C = np.array([[1., 0.,]])
         D = np.array([[0.,]])
-        
+
         T = np.arange(0,10,1)
         sysd_true = StateSpace(A,B,C,D,True)
         ir_true = impulse_response(sysd_true,T=T)
@@ -262,7 +265,7 @@ def testERASignature(self):
         sysd_est, _  = eigensys_realization(ir_true,r=2)
         ir_est = impulse_response(sysd_est, T=T)
         _, H_est = ir_est
-    
+
         np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
 
         # test ndarray
@@ -270,7 +273,7 @@ def testERASignature(self):
         sysd_est, _  = eigensys_realization(YY_true,r=2)
         ir_est = impulse_response(sysd_est, T=T)
         _, H_est = ir_est
-    
+
         np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
 
         # test mimo
@@ -304,10 +307,10 @@ def testERASignature(self):
         step_true = step_response(sysd_true)
         step_est = step_response(sysd_est)
 
-        np.testing.assert_allclose(step_true.outputs, 
+        np.testing.assert_allclose(step_true.outputs,
                                    step_est.outputs,
                                    rtol=1e-6, atol=1e-8)
-        
+
         # test ndarray
         _, YY_true = ir_true
         sysd_est, _  = eigensys_realization(YY_true,r=4,dt=dt)
@@ -315,7 +318,7 @@ def testERASignature(self):
         step_true = step_response(sysd_true, T=T)
         step_est = step_response(sysd_est, T=T)
 
-        np.testing.assert_allclose(step_true.outputs, 
+        np.testing.assert_allclose(step_true.outputs,
                                    step_est.outputs,
                                    rtol=1e-6, atol=1e-8)
 
@@ -343,7 +346,7 @@ def testModredMatchDC(self):
         np.testing.assert_array_almost_equal(rsys.D, Drtrue, decimal=2)
 
     def testModredUnstable(self):
-        """Check if an error is thrown when an unstable system is given"""
+        """Check if warning is issued when an unstable system is given"""
         A = np.array(
             [[4.5418, 3.3999, 5.0342, 4.3808],
              [0.3890, 0.3599, 0.4195, 0.1760],
@@ -353,7 +356,16 @@ def testModredUnstable(self):
         C = np.array([[1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 4.0]])
         D = np.array([[0.0, 0.0], [0.0, 0.0]])
         sys = StateSpace(A, B, C, D)
-        np.testing.assert_raises(ValueError, modred, sys, [2, 3])
+
+        # Make sure we get a warning message
+        with pytest.warns(UserWarning, match="System is unstable"):
+            newsys1 = modred(sys, [2, 3])
+
+        # Make sure we can turn the warning off
+        with warnings.catch_warnings():
+            warnings.simplefilter('error')
+            newsys2 = ct.model_reduction(sys, [2, 3], warn_unstable=False)
+            np.testing.assert_equal(newsys1.A, newsys2.A)
 
     def testModredTruncate(self):
         #balanced realization computed in matlab for the transfer function:
@@ -391,7 +403,7 @@ def testBalredTruncate(self):
         B = np.array([[2.], [0.], [0.], [0.]])
         C = np.array([[0.5, 0.6875, 0.7031, 0.5]])
         D = np.array([[0.]])
-        
+
         sys = StateSpace(A, B, C, D)
         orders = 2
         rsys = balred(sys, orders, method='truncate')
@@ -412,7 +424,7 @@ def testBalredTruncate(self):
                 # Apply a similarity transformation
                 Ar, Br, Cr = T @ Ar @ T, T @ Br, Cr @ T
                 break
-            
+
         # Make sure we got the correct answer
         np.testing.assert_array_almost_equal(Ar, Artrue, decimal=2)
         np.testing.assert_array_almost_equal(Br, Brtrue, decimal=4)
@@ -432,12 +444,12 @@ def testBalredMatchDC(self):
         B = np.array([[2.], [0.], [0.], [0.]])
         C = np.array([[0.5, 0.6875, 0.7031, 0.5]])
         D = np.array([[0.]])
-        
+
         sys = StateSpace(A, B, C, D)
         orders = 2
         rsys = balred(sys,orders,method='matchdc')
         Ar, Br, Cr, Dr = rsys.A, rsys.B, rsys.C, rsys.D
-        
+
         # Result from MATLAB
         Artrue = np.array(
             [[-4.43094773, -4.55232904],
@@ -445,7 +457,7 @@ def testBalredMatchDC(self):
         Brtrue = np.array([[1.36235673], [1.03114388]])
         Crtrue = np.array([[1.36235673, 1.03114388]])
         Drtrue = np.array([[-0.08383902]])
-        
+
         # Look for possible changes in state in slycot
         T1 = np.array([[1, 0], [0, -1]])
         T2 = np.array([[-1, 0], [0, 1]])
@@ -455,9 +467,46 @@ def testBalredMatchDC(self):
                 # Apply a similarity transformation
                 Ar, Br, Cr = T @ Ar @ T, T @ Br, Cr @ T
                 break
-            
+
         # Make sure we got the correct answer
         np.testing.assert_array_almost_equal(Ar, Artrue, decimal=2)
         np.testing.assert_array_almost_equal(Br, Brtrue, decimal=4)
         np.testing.assert_array_almost_equal(Cr, Crtrue, decimal=4)
         np.testing.assert_array_almost_equal(Dr, Drtrue, decimal=4)
+
+
+@pytest.mark.parametrize("kwargs, nstates, noutputs, ninputs", [
+    ({'elim_states': [1, 3]}, 3, 3, 3),
+    ({'elim_inputs': [1, 2], 'keep_states': [1, 3]}, 2, 3, 1),
+    ({'elim_outputs': [1, 2], 'keep_inputs': [0, 1],}, 5, 1, 2),
+    ({'keep_states': [2, 0], 'keep_outputs': [0, 1]}, 2, 2, 3),
+    ({'keep_states': slice(0, 4, 2), 'keep_outputs': slice(None, 2)}, 2, 2, 3),
+    ({'keep_states': ['x[0]', 'x[3]'], 'keep_inputs': 'u[0]'}, 2, 3, 1),
+    ({'elim_inputs': [0, 1, 2]}, 5, 3, 0),              # no inputs
+    ({'elim_outputs': [0, 1, 2]}, 5, 0, 3),             # no outputs
+    ({'elim_states': [0, 1, 2, 3, 4]}, 0, 3, 3),        # no states
+    ({'elim_states': [0, 1], 'keep_states': [1, 2]}, None, None, None),
+])
+@pytest.mark.parametrize("method", ['truncate', 'matchdc'])
+def test_model_reduction(method, kwargs, nstates, noutputs, ninputs):
+    sys = ct.rss(5, 3, 3)
+
+    if nstates is None:
+        # Arguments should generate an error
+        with pytest.raises(ValueError, match="can't provide both"):
+            red = ct.model_reduction(sys, **kwargs, method=method)
+        return
+    else:
+        red = ct.model_reduction(sys, **kwargs, method=method)
+
+    assert red.nstates == nstates
+    assert red.ninputs == ninputs
+    assert red.noutputs == noutputs
+
+    if method == 'matchdc':
+        # Define a new system with truncated inputs and outputs
+        # (assumes we always keep the initial inputs and outputs)
+        chk = ct.ss(
+            sys.A, sys.B[:, :ninputs], sys.C[:noutputs, :],
+            sys.D[:noutputs, :][:, :ninputs])
+        np.testing.assert_allclose(red(0), chk(0))
diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py
index 2d96ad225..3928fb725 100644
--- a/control/tests/statefbk_test.py
+++ b/control/tests/statefbk_test.py
@@ -56,7 +56,7 @@ def testCtrbT(self):
         Wctrue = np.array([[5., 6.], [7., 8.]])
         Wc = ctrb(A, B, t=t)
         np.testing.assert_array_almost_equal(Wc, Wctrue)
-        
+
     def testObsvSISO(self):
         A = np.array([[1., 2.], [3., 4.]])
         C = np.array([[5., 7.]])
@@ -70,7 +70,7 @@ def testObsvMIMO(self):
         Wotrue = np.array([[5., 6.], [7., 8.], [23., 34.], [31., 46.]])
         Wo = obsv(A, C)
         np.testing.assert_array_almost_equal(Wo, Wotrue)
-        
+
     def testObsvT(self):
         A = np.array([[1., 2.], [3., 4.]])
         C = np.array([[5., 6.], [7., 8.]])
@@ -128,15 +128,14 @@ def testGramRc(self):
         C = np.array([[4., 5.], [6., 7.]])
         D = np.array([[13., 14.], [15., 16.]])
         sys = ss(A, B, C, D)
-        Rctrue = np.array([[4.30116263, 5.6961343], 
-                           [0., 0.23249528]])
+        Rctrue = np.array([[4.30116263, 5.6961343], [0., 0.23249528]])
         Rc = gram(sys, 'cf')
         np.testing.assert_array_almost_equal(Rc, Rctrue)
         sysd = ct.c2d(sys, 0.2)
         Rctrue = np.array([[1.91488054, 2.53468814],
                            [0.        , 0.10290372]])
         Rc = gram(sysd, 'cf')
-        np.testing.assert_array_almost_equal(Rc, Rctrue)        
+        np.testing.assert_array_almost_equal(Rc, Rctrue)
 
     @slycotonly
     def testGramWo(self):
@@ -149,7 +148,7 @@ def testGramWo(self):
         Wo = gram(sys, 'o')
         np.testing.assert_array_almost_equal(Wo, Wotrue)
         sysd = ct.c2d(sys, 0.2)
-        Wotrue = np.array([[ 1305.369179, -440.046414], 
+        Wotrue = np.array([[ 1305.369179, -440.046414],
                            [ -440.046414,  333.034844]])
         Wo = gram(sysd, 'o')
         np.testing.assert_array_almost_equal(Wo, Wotrue)
@@ -184,7 +183,7 @@ def testGramRo(self):
         Rotrue = np.array([[ 36.12989315, -12.17956588],
                            [  0.        ,  13.59018097]])
         Ro = gram(sysd, 'of')
-        np.testing.assert_array_almost_equal(Ro, Rotrue) 
+        np.testing.assert_array_almost_equal(Ro, Rotrue)
 
     def testGramsys(self):
         sys = tf([1.], [1., 1., 1.])
@@ -1143,3 +1142,55 @@ def test_gainsched_errors(unicycle):
         ctrl, clsys = ct.create_statefbk_iosystem(
             unicycle, (gains, points),
             gainsched_indices=[3, 2], gainsched_method='unknown')
+
+
+@pytest.mark.parametrize("ninputs, Kf", [
+    (1, 1),
+    (1, None),
+    (2, np.diag([1, 1])),
+    (2, None),
+])
+def test_refgain_pattern(ninputs, Kf):
+    sys = ct.rss(2, 2, ninputs, strictly_proper=True)
+    sys.C = np.eye(2)
+
+    K, _, _ = ct.lqr(sys.A, sys.B, np.eye(sys.nstates), np.eye(sys.ninputs))
+    if Kf is None:
+        # Make sure we get an error if we don't specify Kf
+        with pytest.raises(ControlArgument, match="'feedfwd_gain' required"):
+            ctrl, clsys = ct.create_statefbk_iosystem(
+                sys, K, Kf, feedfwd_pattern='refgain')
+
+        # Now compute the gain to give unity zero frequency gain
+        C = np.eye(ninputs, sys.nstates)
+        Kf = -np.linalg.inv(
+            C @ np.linalg.inv(sys.A - sys.B @ K) @ sys.B)
+        ctrl, clsys = ct.create_statefbk_iosystem(
+            sys, K, Kf, feedfwd_pattern='refgain')
+
+        np.testing.assert_almost_equal(
+            C @ clsys(0)[0:sys.nstates], np.eye(ninputs))
+
+    else:
+        ctrl, clsys = ct.create_statefbk_iosystem(
+            sys, K, Kf, feedfwd_pattern='refgain')
+
+    manual = ct.feedback(sys, K) * Kf
+    np.testing.assert_almost_equal(clsys.A, manual.A)
+    np.testing.assert_almost_equal(clsys.B, manual.B)
+    np.testing.assert_almost_equal(clsys.C[:sys.nstates, :], manual.C)
+    np.testing.assert_almost_equal(clsys.D[:sys.nstates, :], manual.D)
+
+
+def test_create_statefbk_errors():
+    sys = ct.rss(2, 2, 1, strictly_proper=True)
+    sys.C = np.eye(2)
+    K = -np.ones((1, 4))
+    Kf = 1
+
+    K, _, _ = ct.lqr(sys.A, sys.B, np.eye(sys.nstates), np.eye(sys.ninputs))
+    with pytest.raises(NotImplementedError, match="unknown pattern"):
+        ct.create_statefbk_iosystem(sys, K, feedfwd_pattern='mypattern')
+
+    with pytest.raises(ControlArgument, match="feedfwd_pattern != 'refgain'"):
+        ct.create_statefbk_iosystem(sys, K, Kf, feedfwd_pattern='trajgen')
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index cb200c4ab..1798c524c 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -121,29 +121,30 @@ def test_constructor(self, sys322ABCD, dt, argfun):
         np.testing.assert_almost_equal(sys.D, sys322ABCD[3])
         assert sys.dt == dtref
 
-    @pytest.mark.parametrize("args, exc, errmsg",
-                             [((True, ), TypeError,
-                               "(can only take in|sys must be) a StateSpace"),
-                              ((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 must be square"),
-                              ((np.ones((3, 3)), np.ones((2, 2)),
-                                np.ones((2, 3)), np.ones((2, 2))),
-                               ValueError, "A and B"),
-                              ((np.ones((3, 3)), np.ones((3, 2)),
-                                np.ones((2, 2)), np.ones((2, 2))),
-                               ValueError, "A and C"),
-                              ((np.ones((3, 3)), np.ones((3, 2)),
-                                np.ones((2, 3)), np.ones((2, 3))),
-                               ValueError, "B and D"),
-                              ((np.ones((3, 3)), np.ones((3, 2)),
-                                np.ones((2, 3)), np.ones((3, 2))),
-                               ValueError, "C and D"),
-                              ])
+    @pytest.mark.parametrize(
+        "args, exc, errmsg",
+        [((True, ), TypeError, "(can only take in|sys must be) a StateSpace"),
+         ((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\)"),
+         ((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\)"),
+         ((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\)"),
+         ((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\)"),
+         ((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\)"),
+        ])
     def test_constructor_invalid(self, args, exc, errmsg):
         """Test invalid input to StateSpace() constructor"""
-        with pytest.raises(exc, match=errmsg):
+
+        with pytest.raises(exc, match=errmsg) as w:
             StateSpace(*args)
         with pytest.raises(exc, match=errmsg):
             ss(*args)