Repository: python-control/python-control

Issue:
Proposal for multiplying a NumPy array by a scalar `TransferFunction`
Currently, multiplying a SISO `TransferFunction` object by a NumPy array results in a dimension error:

```python
import control
import numpy as np

tf = control.TransferFunction([1], [1, 0])
arr = np.array([
    [1, 0],
    [2, 1],
])

print(tf * arr)
print(arr * tf)
```

```
ValueError: C = A * B: A has 1 column(s) (input(s)), but B has 2 row(s) (output(s)).
```

I think it would be convenient if output was a MIMO `TransferFunction` object of the array's dimension. So, for example,

```python
G = control.TransferFunction(
    [
        [[1], [0]],
        [[2], [1]],
    ],
    [
        [[1, 0], [1, 0]],
        [[1, 0], [1, 0]],
    ],
)
````

I would be happy to implement this and do a PR, but I want to make sure it's something the maintainers actually want. I would think all SISO systems (state-space, etc) should behave in the same way.

This appears to have been discussed in #459 , but I'm not sure what the conclusion was.

Hints from issue comments:
@murrayrm: I think it makes sense that if a "scalar" I/O system (= SISO system) is used in an operation involving either an ndarray or a MIMO transfer function, it should broadcast to the right shape.

Specific examples/comments:
* <2x2 MIMO system> * <1x1 SISO system> should broadcast the dimensions so that each element in the 2x2 system is multiplied by the SISO system.  Same thing would happen with + or - (each element has the SISO system added or subtracted) and everything should also work if the arguments are switched.
* <mxn ndarray> * <1x1 SISO system> should generate an <mxn MIMO system> with each entry being the SISO system.
* <mxn MIMO system> * <nxr MIMO system> becomes a <mxr> MIMO system (as currently happens).  I bring this one up because this is *not* what happens with ndarrays (you need to use `@` to get matrix multiplication) => maybe we want to implement `__matmul__`?
* If we do implement `__matmul__`, then I think <2x2 MIMO system> @ <1x1 SISO system> should probably generate an error, as it does with ndarrays of those sizes.

For now, perhaps we just implement the SISO → MIMO broadcast functionality and leave out `_matmul__`?  Note that all of this should work for all I/O systems, not just transfer functions.  That might be tricky for `NonlinearIOSystem`'s, so perhaps we just implement for `LTI` systems (`FrequencyResponseData`, `StateSpace`, `TransferFuncion`).@sdahdah: > I think it makes sense that if a "scalar" I/O system (= SISO system) is used in an operation involving either an ndarray or a MIMO transfer function, it should broadcast to the right shape.
> 
> Specific examples/comments:
> 
>     * <2x2 MIMO system> * <1x1 SISO system> should broadcast the dimensions so that each element in the 2x2 system is multiplied by the SISO system.  Same thing would happen with + or - (each element has the SISO system added or subtracted) and everything should also work if the arguments are switched.

^ I very much agree with this

>     * * <1x1 SISO system> should generate an  with each entry being the SISO system.
> 
>     * *  becomes a  MIMO system (as currently happens).  I bring this one up because this is _not_ what happens with ndarrays (you need to use `@` to get matrix multiplication) => maybe we want to implement `__matmul__`?

^ I'm not sure I follow these two points
 
>     * If we do implement `__matmul__`, then I think <2x2 MIMO system> @ <1x1 SISO system> should probably generate an error, as it does with ndarrays of those sizes.

^ I do think that this should raise an error if `__matmul__` is implemented. My gut feeling is that implementing `__matmul__` would be confusing though, since `series()` is essentially "multiplying transfer matrices." Is the idea that `*` would allow broadcasting and `@` would not, but otherwise they would function identically?
 
> For now, perhaps we just implement the SISO → MIMO broadcast functionality and leave out `_matmul__`? Note that all of this should work for all I/O systems, not just transfer functions. That might be tricky for `NonlinearIOSystem`'s, so perhaps we just implement for `LTI` systems (`FrequencyResponseData`, `StateSpace`, `TransferFuncion`).

^ I will take a look at this next week and see if I can find a nice way to do it for all I/O systems. If not, I will start small with just LTI.

Should I also look at #1075 at the same time?@murrayrm: > Should I also look at https://github.com/python-control/python-control/issues/1075 at the same time?

I'll probably take a look at this over the week and if the fix is simple, will submit a PR that you can build on.  But if I haven't gotten to it by the time you start, it would be great to take care of #1075 at the same time (it will be in the same block of code).@murrayrm: #1075 is fixed.  Make sure to build on PR #1078 when you make the changes for this issue.@sdahdah: Thanks. I started working on it earlier this week, and I will probably be ready to do a PR next week.

Solution patch:
diff --git a/control/bdalg.py b/control/bdalg.py
index 2be239177..f066b72b5 100644
--- a/control/bdalg.py
+++ b/control/bdalg.py
@@ -356,14 +356,14 @@ def feedback(sys1, sys2=1, sign=-1, **kwargs):
 def append(*sys, **kwargs):
     """append(sys1, sys2, [..., sysn])
 
-    Group LTI state space models by appending their inputs and outputs.
+    Group LTI models by appending their inputs and outputs.
 
     Forms an augmented system model, and appends the inputs and
     outputs together.
 
     Parameters
     ----------
-    sys1, sys2, ..., sysn: scalar, array, or :class:`StateSpace`
+    sys1, sys2, ..., sysn : scalar, array, or :class:`LTI`
         I/O systems to combine.
 
     Other Parameters
@@ -382,9 +382,10 @@ def append(*sys, **kwargs):
 
     Returns
     -------
-    out: :class:`StateSpace`
+    out : :class:`LTI`
         Combined system, with input/output vectors consisting of all
-        input/output vectors appended.
+        input/output vectors appended. Specific type returned is the type of
+        the first argument.
 
     See Also
     --------
@@ -405,7 +406,7 @@ def append(*sys, **kwargs):
     (3, 8, 7)
 
     """
-    s1 = ss._convert_to_statespace(sys[0])
+    s1 = sys[0]
     for s in sys[1:]:
         s1 = s1.append(s)
     s1.update_names(**kwargs)
diff --git a/control/frdata.py b/control/frdata.py
index cb6925661..1200bfffa 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -20,6 +20,7 @@
 from scipy.interpolate import splev, splprep
 
 from . import config
+from . import bdalg
 from .exception import pandas_check
 from .iosys import InputOutputSystem, NamedSignal, _extended_system_name, \
     _process_iosys_keywords, _process_subsys_index, common_timebase
@@ -455,6 +456,12 @@ def __add__(self, other):
         else:
             other = _convert_to_frd(other, omega=self.omega)
 
+        # Promote SISO object to compatible dimension
+        if self.issiso() and not other.issiso():
+            self = np.ones((other.noutputs, other.ninputs)) * self
+        elif not self.issiso() and other.issiso():
+            other = np.ones((self.noutputs, self.ninputs)) * other
+
         # Check that the input-output sizes are consistent.
         if self.ninputs != other.ninputs:
             raise ValueError(
@@ -492,6 +499,12 @@ def __mul__(self, other):
         else:
             other = _convert_to_frd(other, omega=self.omega)
 
+        # Promote SISO object to compatible dimension
+        if self.issiso() and not other.issiso():
+            self = bdalg.append(*([self] * other.noutputs))
+        elif not self.issiso() and other.issiso():
+            other = bdalg.append(*([other] * self.ninputs))
+
         # Check that the input-output sizes are consistent.
         if self.ninputs != other.noutputs:
             raise ValueError(
@@ -519,6 +532,12 @@ def __rmul__(self, other):
         else:
             other = _convert_to_frd(other, omega=self.omega)
 
+        # Promote SISO object to compatible dimension
+        if self.issiso() and not other.issiso():
+            self = bdalg.append(*([self] * other.ninputs))
+        elif not self.issiso() and other.issiso():
+            other = bdalg.append(*([other] * self.noutputs))
+
         # Check that the input-output sizes are consistent.
         if self.noutputs != other.ninputs:
             raise ValueError(
@@ -547,11 +566,9 @@ def __truediv__(self, other):
         else:
             other = _convert_to_frd(other, omega=self.omega)
 
-        if (self.ninputs > 1 or self.noutputs > 1 or
-            other.ninputs > 1 or other.noutputs > 1):
-            raise NotImplementedError(
-                "FRD.__truediv__ is currently only implemented for SISO "
-                "systems.")
+        if (other.ninputs > 1 or other.noutputs > 1):
+            # FRD.__truediv__ is currently only implemented for SISO systems
+            return NotImplemented
 
         return FRD(self.fresp/other.fresp, self.omega,
                    smooth=(self._ifunc is not None) and
@@ -566,11 +583,9 @@ def __rtruediv__(self, other):
         else:
             other = _convert_to_frd(other, omega=self.omega)
 
-        if (self.ninputs > 1 or self.noutputs > 1 or
-            other.ninputs > 1 or other.noutputs > 1):
-            raise NotImplementedError(
-                "FRD.__rtruediv__ is currently only implemented for "
-                "SISO systems.")
+        if (self.ninputs > 1 or self.noutputs > 1):
+            # FRD.__rtruediv__ is currently only implemented for SISO systems
+            return NotImplemented
 
         return other / self
 
@@ -803,6 +818,26 @@ def feedback(self, other=1, sign=-1):
 
         return FRD(fresp, other.omega, smooth=(self._ifunc is not None))
 
+    def append(self, other):
+        """Append a second model to the present model.
+
+        The second model is converted to FRD if necessary, inputs and
+        outputs are appended and their order is preserved"""
+        other = _convert_to_frd(other, omega=self.omega, inputs=other.ninputs,
+                                outputs=other.noutputs)
+
+        # TODO: handle omega re-mapping
+
+        new_fresp = np.zeros(
+            (self.noutputs + other.noutputs, self.ninputs + other.ninputs,
+             self.omega.shape[-1]), dtype=complex)
+        new_fresp[:self.noutputs, :self.ninputs, :] = np.reshape(
+            self.fresp, (self.noutputs, self.ninputs, -1))
+        new_fresp[self.noutputs:, self.ninputs:, :] = np.reshape(
+            other.fresp, (other.noutputs, other.ninputs, -1))
+
+        return FRD(new_fresp, self.omega, smooth=(self._ifunc is not None))
+
     # Plotting interface
     def plot(self, plot_type=None, *args, **kwargs):
         """Plot the frequency response using a Bode plot.
diff --git a/control/statefbk.py b/control/statefbk.py
index c5c5a8030..57df28d45 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -47,12 +47,12 @@
 from . import statesp
 from .config import _process_legacy_keyword
 from .exception import ControlArgument, ControlDimension, \
-    ControlNotImplemented, ControlSlycot
+    ControlSlycot
 from .iosys import _process_indices, _process_labels, isctime, isdtime
 from .lti import LTI
 from .mateqn import _check_shape, care, dare
 from .nlsys import NonlinearIOSystem, interconnect
-from .statesp import StateSpace, _convert_to_statespace, _ssmatrix, ss
+from .statesp import StateSpace, _ssmatrix, ss
 
 # Make sure we have access to the right slycot routines
 try:
@@ -1098,19 +1098,26 @@ def ctrb(A, B, t=None):
     """
 
     # Convert input parameters to matrices (if they aren't already)
-    amat = _ssmatrix(A)
-    bmat = _ssmatrix(B)
-    n = np.shape(amat)[0]
-    m = np.shape(bmat)[1]
+    A = _ssmatrix(A)
+    if np.asarray(B).ndim == 1 and len(B) == A.shape[0]:
+        B = _ssmatrix(B, axis=0)
+    else:
+        B = _ssmatrix(B)
+
+    n = A.shape[0]
+    m = B.shape[1]
+
+    _check_shape('A', A, n, n, square=True)
+    _check_shape('B', B, n, m)
 
     if t is None or t > n:
         t = n
 
     # Construct the controllability matrix
     ctrb = np.zeros((n, t * m))
-    ctrb[:, :m] = bmat
+    ctrb[:, :m] = B
     for k in range(1, t):
-        ctrb[:, k * m:(k + 1) * m] = np.dot(amat, ctrb[:, (k - 1) * m:k * m])
+        ctrb[:, k * m:(k + 1) * m] = np.dot(A, ctrb[:, (k - 1) * m:k * m])
 
     return _ssmatrix(ctrb)
 
@@ -1140,20 +1147,24 @@ def obsv(A, C, t=None):
     """
 
     # Convert input parameters to matrices (if they aren't already)
-    amat = _ssmatrix(A)
-    cmat = _ssmatrix(C)
-    n = np.shape(amat)[0]
-    p = np.shape(cmat)[0]
+    A = _ssmatrix(A)
+    C = _ssmatrix(C)
+
+    n = np.shape(A)[0]
+    p = np.shape(C)[0]
+
+    _check_shape('A', A, n, n, square=True)
+    _check_shape('C', C, p, n)
 
     if t is None or t > n:
         t = n
 
     # Construct the observability matrix
     obsv = np.zeros((t * p, n))
-    obsv[:p, :] = cmat
+    obsv[:p, :] = C
 
     for k in range(1, t):
-        obsv[k * p:(k + 1) * p, :] = np.dot(obsv[(k - 1) * p:k * p, :], amat)
+        obsv[k * p:(k + 1) * p, :] = np.dot(obsv[(k - 1) * p:k * p, :], A)
 
     return _ssmatrix(obsv)
 
diff --git a/control/statesp.py b/control/statesp.py
index 98adc942f..44fe8b605 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -30,6 +30,7 @@
 from scipy.signal import cont2discrete
 
 from . import config
+from . import bdalg
 from .exception import ControlMIMONotImplemented, ControlSlycot, slycot_check
 from .frdata import FrequencyResponseData
 from .iosys import InputOutputSystem, NamedSignal, _process_dt_keyword, \
@@ -572,6 +573,9 @@ def __add__(self, other):
 
         elif isinstance(other, np.ndarray):
             other = np.atleast_2d(other)
+            # Special case for SISO
+            if self.issiso():
+                self = np.ones_like(other) * self
             if self.ninputs != other.shape[0]:
                 raise ValueError("array has incompatible shape")
             A, B, C = self.A, self.B, self.C
@@ -582,6 +586,12 @@ def __add__(self, other):
             return NotImplemented       # let other.__rmul__ handle it
 
         else:
+            # Promote SISO object to compatible dimension
+            if self.issiso() and not other.issiso():
+                self = np.ones((other.noutputs, other.ninputs)) * self
+            elif not self.issiso() and other.issiso():
+                other = np.ones((self.noutputs, self.ninputs)) * other
+
             # Check to make sure the dimensions are OK
             if ((self.ninputs != other.ninputs) or
                     (self.noutputs != other.noutputs)):
@@ -636,6 +646,10 @@ def __mul__(self, other):
 
         elif isinstance(other, np.ndarray):
             other = np.atleast_2d(other)
+            # Special case for SISO
+            if self.issiso():
+                self = bdalg.append(*([self] * other.shape[0]))
+            # Dimension check after broadcasting
             if self.ninputs != other.shape[0]:
                 raise ValueError("array has incompatible shape")
             A, C = self.A, self.C
@@ -647,6 +661,12 @@ def __mul__(self, other):
             return NotImplemented       # let other.__rmul__ handle it
 
         else:
+            # Promote SISO object to compatible dimension
+            if self.issiso() and not other.issiso():
+                self = bdalg.append(*([self] * other.noutputs))
+            elif not self.issiso() and other.issiso():
+                other = bdalg.append(*([other] * self.ninputs))
+
             # Check to make sure the dimensions are OK
             if self.ninputs != other.noutputs:
                 raise ValueError(
@@ -686,23 +706,67 @@ def __rmul__(self, other):
             return StateSpace(self.A, B, self.C, D, self.dt)
 
         elif isinstance(other, np.ndarray):
-            C = np.atleast_2d(other) @ self.C
-            D = np.atleast_2d(other) @ self.D
+            other = np.atleast_2d(other)
+            # Special case for SISO transfer function
+            if self.issiso():
+                self = bdalg.append(*([self] * other.shape[1]))
+            # Dimension check after broadcasting
+            if self.noutputs != other.shape[1]:
+                raise ValueError("array has incompatible shape")
+            C = other @ self.C
+            D = other @ self.D
             return StateSpace(self.A, self.B, C, D, self.dt)
 
         if not isinstance(other, StateSpace):
             return NotImplemented
 
+        # Promote SISO object to compatible dimension
+        if self.issiso() and not other.issiso():
+            self = bdalg.append(*([self] * other.ninputs))
+        elif not self.issiso() and other.issiso():
+            other = bdalg.append(*([other] * self.noutputs))
+
         return other * self
 
     # TODO: general __truediv__ requires descriptor system support
     def __truediv__(self, other):
         """Division of state space systems by TFs, FRDs, scalars, and arrays"""
-        if not isinstance(other, (LTI, InputOutputSystem)):
-            return self * (1/other)
-        else:
+        # Let ``other.__rtruediv__`` handle it
+        try:
+            return self * (1 / other)
+        except ValueError:
             return NotImplemented
 
+    def __rtruediv__(self, other):
+        """Division by state space system"""
+        return other * self**-1
+
+    def __pow__(self, other):
+        """Power of a state space system"""
+        if not type(other) == int:
+            raise ValueError("Exponent must be an integer")
+        if self.ninputs != self.noutputs:
+            # System must have same number of inputs and outputs
+            return NotImplemented
+        if other < -1:
+            return (self**-1)**(-other)
+        elif other == -1:
+            try:
+                Di = scipy.linalg.inv(self.D)
+            except scipy.linalg.LinAlgError:
+                # D matrix must be nonsingular
+                return NotImplemented
+            Ai = self.A - self.B @ Di @ self.C
+            Bi = self.B @ Di
+            Ci = -Di @ self.C
+            return StateSpace(Ai, Bi, Ci, Di, self.dt)
+        elif other == 0:
+            return StateSpace([], [], [], np.eye(self.ninputs), self.dt)
+        elif other == 1:
+            return self
+        elif other > 1:
+            return self * (self**(other - 1))
+
     def __call__(self, x, squeeze=None, warn_infinite=True):
         """Evaluate system's frequency response at complex frequencies.
 
@@ -1107,7 +1171,7 @@ def minreal(self, tol=0.0):
                 A, B, C, nr = tb01pd(self.nstates, self.ninputs, self.noutputs,
                                      self.A, B, C, tol=tol)
                 return StateSpace(A[:nr, :nr], B[:nr, :self.ninputs],
-                                  C[:self.noutputs, :nr], self.D)
+                                  C[:self.noutputs, :nr], self.D, self.dt)
             except ImportError:
                 raise TypeError("minreal requires slycot tb01pd")
         else:
diff --git a/control/xferfcn.py b/control/xferfcn.py
index dea61f9f3..c9043bd0e 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -29,6 +29,7 @@
 from scipy.signal import cont2discrete, tf2zpk, zpk2tf
 
 from . import config
+from . import bdalg
 from .exception import ControlMIMONotImplemented
 from .frdata import FrequencyResponseData
 from .iosys import InputOutputSystem, NamedSignal, _process_iosys_keywords, \
@@ -582,6 +583,12 @@ def __add__(self, other):
         if not isinstance(other, TransferFunction):
             return NotImplemented
 
+        # Promote SISO object to compatible dimension
+        if self.issiso() and not other.issiso():
+            self = np.ones((other.noutputs, other.ninputs)) * self
+        elif not self.issiso() and other.issiso():
+            other = np.ones((self.noutputs, self.ninputs)) * other
+
         # Check that the input-output sizes are consistent.
         if self.ninputs != other.ninputs:
             raise ValueError(
@@ -631,6 +638,12 @@ def __mul__(self, other):
         if not isinstance(other, TransferFunction):
             return NotImplemented
 
+        # Promote SISO object to compatible dimension
+        if self.issiso() and not other.issiso():
+            self = bdalg.append(*([self] * other.noutputs))
+        elif not self.issiso() and other.issiso():
+            other = bdalg.append(*([other] * self.ninputs))
+
         # Check that the input-output sizes are consistent.
         if self.ninputs != other.noutputs:
             raise ValueError(
@@ -674,6 +687,12 @@ def __rmul__(self, other):
         else:
             other = _convert_to_transfer_function(other)
 
+        # Promote SISO object to compatible dimension
+        if self.issiso() and not other.issiso():
+            self = bdalg.append(*([self] * other.ninputs))
+        elif not self.issiso() and other.issiso():
+            other = bdalg.append(*([other] * self.noutputs))
+
         # Check that the input-output sizes are consistent.
         if other.ninputs != self.noutputs:
             raise ValueError(
@@ -718,12 +737,16 @@ def __truediv__(self, other):
         else:
             other = _convert_to_transfer_function(other)
 
+        # Special case for SISO ``other``
+        if not self.issiso() and other.issiso():
+            other = bdalg.append(*([other**-1] * self.noutputs))
+            return self * other
+
         if (self.ninputs > 1 or self.noutputs > 1 or
                 other.ninputs > 1 or other.noutputs > 1):
-            raise NotImplementedError(
-                "TransferFunction.__truediv__ is currently \
-                implemented only for SISO systems.")
-
+            # TransferFunction.__truediv__ is currently implemented only for
+            # SISO systems.
+            return NotImplemented
         dt = common_timebase(self.dt, other.dt)
 
         num = polymul(self.num_array[0, 0], other.den_array[0, 0])
@@ -741,11 +764,16 @@ def __rtruediv__(self, other):
         else:
             other = _convert_to_transfer_function(other)
 
+        # Special case for SISO ``self``
+        if self.issiso() and not other.issiso():
+            self = bdalg.append(*([self**-1] * other.ninputs))
+            return other * self
+
         if (self.ninputs > 1 or self.noutputs > 1 or
                 other.ninputs > 1 or other.noutputs > 1):
-            raise NotImplementedError(
-                "TransferFunction.__rtruediv__ is currently implemented only "
-                "for SISO systems.")
+            # TransferFunction.__rtruediv__ is currently implemented only for
+            # SISO systems
+            return NotImplemented
 
         return other / self
 
@@ -850,6 +878,20 @@ def feedback(self, other=1, sign=-1):
         # But this does not work correctly because the state size will be too
         # large.
 
+    def append(self, other):
+        """Append a second model to the present model.
+
+        The second model is converted to a transfer function if necessary,
+        inputs and outputs are appended and their order is preserved"""
+        other = _convert_to_transfer_function(other)
+
+        new_tf = bdalg.combine_tf([
+            [self, np.zeros((self.noutputs, other.ninputs))],
+            [np.zeros((other.noutputs, self.ninputs)), other],
+        ])
+
+        return new_tf
+
     def minreal(self, tol=None):
         """Remove cancelling pole/zero pairs from a transfer function."""
         # based on octave minreal
@@ -1751,7 +1793,7 @@ def zpk(zeros, poles, gain, *args, **kwargs):
 
     Returns
     -------
-    out: `TransferFunction`
+    out : `TransferFunction`
         Transfer function with given zeros, poles, and gain.
 
     Examples
@@ -1803,7 +1845,7 @@ def ss2tf(*args, **kwargs):
 
     Returns
     -------
-    out: TransferFunction
+    out : TransferFunction
         New linear system in transfer function form
 
     Other Parameters
@@ -1968,6 +2010,53 @@ def _float2str(value):
     return f"{value:{_num_format}}"
 
 
+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
+
+
 def _create_poly_array(shape, default=None):
     out = np.empty(shape, dtype=np.ndarray)
     if default is not None:
diff --git a/doc/control.rst b/doc/control.rst
index dd418f2af..766e593d8 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -37,6 +37,8 @@ System interconnections
     parallel
     series
     connection_table
+    combine_tf
+    split_tf
 
 
 Frequency domain plotting


Test patch:
diff --git a/control/tests/bdalg_test.py b/control/tests/bdalg_test.py
index e0d64f018..f69574d9a 100644
--- a/control/tests/bdalg_test.py
+++ b/control/tests/bdalg_test.py
@@ -8,7 +8,7 @@
 import pytest
 
 import control as ctrl
-from control.xferfcn import TransferFunction
+from control.xferfcn import TransferFunction, _tf_close_coeff
 from control.statesp import StateSpace
 from control.bdalg import feedback, append, connect
 from control.lti import zeros, poles
@@ -870,50 +870,3 @@ 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_array[i, j],
-                tf_b.num_array[i, j],
-                rtol=rtol,
-                atol=atol,
-            ):
-                return False
-            if not np.allclose(
-                tf_a.den_array[i, j],
-                tf_b.den_array[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 0c0a7904f..b1fce53e0 100644
--- a/control/tests/docstrings_test.py
+++ b/control/tests/docstrings_test.py
@@ -28,7 +28,7 @@
 
 # Checksums to use for checking whether a docstring has changed
 function_docstring_hash = {
-    control.append:                     '48548c4c4e0083312b3ea9e56174b0b5',
+    control.append:                     '1bddbac0fe932755c85e9fb0bfb97d88',
     control.describing_function_plot:   '95f894706b1d3eeb3b854934596af09f',
     control.dlqe:                       '9db995ed95c2214ce97074b0616a3191',
     control.dlqr:                       '896cfa651dbbd80e417635904d13c9d6',
@@ -37,7 +37,7 @@
     control.margin:                     'f02b3034f5f1d44ce26f916cc3e51600',
     control.parallel:                   'bfc470aef75dbb923f9c6fb8bf3c9b43',
     control.series:                     '9aede1459667738f05cf4fc46603a4f6',
-    control.ss2tf:                      '48ff25d22d28e7b396e686dd5eb58831',
+    control.ss2tf:                      'e779b8d70205bc1218cc2a4556a66e4b',
     control.tf2ss:                      '086a3692659b7321c2af126f79f4bc11',
     control.markov:                     'a4199c54cb50f07c0163d3790739eafe',
     control.gangof4:                    '0e52eb6cf7ce024f9a41f3ae3ebf04f7',
diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index c63d9e217..5c0acd4c8 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -187,6 +187,56 @@ def testFeedback2(self):
                         [[1.0, 0], [0, 1]], [[0.0], [0.0]])
         # h2.feedback([[0.3, 0.2], [0.1, 0.1]])
 
+    def testAppendSiso(self):
+        # Create frequency responses
+        d1 = np.array([1 + 2j, 1 - 2j, 1 + 4j, 1 - 4j, 1 + 6j, 1 - 6j])
+        d2 = d1 + 2
+        d3 = d1 - 1j
+        w = np.arange(d1.shape[-1])
+        frd1 = FrequencyResponseData(d1, w)
+        frd2 = FrequencyResponseData(d2, w)
+        frd3 = FrequencyResponseData(d3, w)
+        # Create appended frequency responses
+        d_app_1 = np.zeros((2, 2, d1.shape[-1]), dtype=complex)
+        d_app_1[0, 0, :] = d1
+        d_app_1[1, 1, :] = d2
+        d_app_2 = np.zeros((3, 3, d1.shape[-1]), dtype=complex)
+        d_app_2[0, 0, :] = d1
+        d_app_2[1, 1, :] = d2
+        d_app_2[2, 2, :] = d3
+        # Test appending two FRDs
+        frd_app_1 = frd1.append(frd2)
+        np.testing.assert_allclose(d_app_1, frd_app_1.fresp)
+        # Test appending three FRDs
+        frd_app_2 = frd1.append(frd2).append(frd3)
+        np.testing.assert_allclose(d_app_2, frd_app_2.fresp)
+
+    def testAppendMimo(self):
+        # Create frequency responses
+        rng = np.random.default_rng(1234)
+        n = 100
+        w = np.arange(n)
+        d1 = rng.uniform(size=(2, 2, n)) + 1j * rng.uniform(size=(2, 2, n))
+        d2 = rng.uniform(size=(3, 1, n)) + 1j * rng.uniform(size=(3, 1, n))
+        d3 = rng.uniform(size=(1, 2, n)) + 1j * rng.uniform(size=(1, 2, n))
+        frd1 = FrequencyResponseData(d1, w)
+        frd2 = FrequencyResponseData(d2, w)
+        frd3 = FrequencyResponseData(d3, w)
+        # Create appended frequency responses
+        d_app_1 = np.zeros((5, 3, d1.shape[-1]), dtype=complex)
+        d_app_1[:2, :2, :] = d1
+        d_app_1[2:, 2:, :] = d2
+        d_app_2 = np.zeros((6, 5, d1.shape[-1]), dtype=complex)
+        d_app_2[:2, :2, :] = d1
+        d_app_2[2:5, 2:3, :] = d2
+        d_app_2[5:, 3:, :] = d3
+        # Test appending two FRDs
+        frd_app_1 = frd1.append(frd2)
+        np.testing.assert_allclose(d_app_1, frd_app_1.fresp)
+        # Test appending three FRDs
+        frd_app_2 = frd1.append(frd2).append(frd3)
+        np.testing.assert_allclose(d_app_2, frd_app_2.fresp)
+
     def testAuto(self):
         omega = np.logspace(-1, 2, 10)
         f1 = _convert_to_frd(1, omega)
@@ -424,15 +474,238 @@ def test_operator_conversion(self):
         np.testing.assert_array_almost_equal(sys_add.omega, chk_add.omega)
         np.testing.assert_array_almost_equal(sys_add.fresp, chk_add.fresp)
 
+        # Test broadcasting with SISO system
+        sys_tf_mimo = TransferFunction([1], [1, 0]) * np.eye(2)
+        frd_tf_mimo = frd(sys_tf_mimo, np.logspace(-1, 1, 10))
+        result = FrequencyResponseData.__rmul__(frd_tf, frd_tf_mimo)
+        expected = frd(sys_tf_mimo * sys_tf, np.logspace(-1, 1, 10))
+        np.testing.assert_array_almost_equal(expected.omega, result.omega)
+        np.testing.assert_array_almost_equal(expected.fresp, result.fresp)
+
         # Input/output mismatch size mismatch in rmul
         sys1 = frd(ct.rss(2, 2, 2), np.logspace(-1, 1, 10))
+        sys2 = frd(ct.rss(3, 3, 3), np.logspace(-1, 1, 10))
         with pytest.raises(ValueError):
-            FrequencyResponseData.__rmul__(frd_2, sys1)
+            FrequencyResponseData.__rmul__(sys2, sys1)
 
         # Make sure conversion of something random generates exception
         with pytest.raises(TypeError):
             FrequencyResponseData.__add__(frd_tf, 'string')
 
+    def test_add_sub_mimo_siso(self):
+        omega = np.logspace(-1, 1, 10)
+        sys_mimo = frd(ct.rss(2, 2, 2), omega)
+        sys_siso = frd(ct.rss(2, 1, 1), omega)
+
+        for op, expected_fresp in [
+            (FrequencyResponseData.__add__, sys_mimo.fresp + sys_siso.fresp),
+            (FrequencyResponseData.__radd__, sys_mimo.fresp + sys_siso.fresp),
+            (FrequencyResponseData.__sub__, sys_mimo.fresp - sys_siso.fresp),
+            (FrequencyResponseData.__rsub__, -sys_mimo.fresp + sys_siso.fresp),
+        ]:
+            result = op(sys_mimo, sys_siso)
+            np.testing.assert_array_almost_equal(omega, result.omega)
+            np.testing.assert_array_almost_equal(expected_fresp, result.fresp)
+
+    @pytest.mark.parametrize(
+        "left, right, expected",
+        [
+            (
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                TransferFunction([2], [1, 0]),
+                np.eye(3),
+                TransferFunction(
+                    [
+                        [[2], [0], [0]],
+                        [[0], [2], [0]],
+                        [[0], [0], [2]],
+                    ],
+                    [
+                        [[1, 0], [1], [1]],
+                        [[1], [1, 0], [1]],
+                        [[1], [1], [1, 0]],
+                    ],
+                ),
+            ),
+        ]
+    )
+    def test_mul_mimo_siso(self, left, right, expected):
+        result = frd(left, np.logspace(-1, 1, 10)).__mul__(right)
+        expected_frd = frd(expected, np.logspace(-1, 1, 10))
+        np.testing.assert_array_almost_equal(expected_frd.omega, result.omega)
+        np.testing.assert_array_almost_equal(expected_frd.fresp, result.fresp)
+
+    @slycotonly
+    def test_truediv_mimo_siso(self):
+        omega = np.logspace(-1, 1, 10)
+        tf_mimo = TransferFunction([1], [1, 0]) * np.eye(2)
+        frd_mimo = frd(tf_mimo, omega)
+        ss_mimo = ct.tf2ss(tf_mimo)
+        tf_siso = TransferFunction([1], [1, 1])
+        frd_siso = frd(tf_siso, omega)
+        expected = frd(tf_mimo.__truediv__(tf_siso), omega)
+        ss_siso = ct.tf2ss(tf_siso)
+
+        # Test division of MIMO FRD by SISO FRD
+        result = frd_mimo.__truediv__(frd_siso)
+        np.testing.assert_array_almost_equal(expected.omega, result.omega)
+        np.testing.assert_array_almost_equal(expected.fresp, result.fresp)
+
+        # Test division of MIMO FRD by SISO TF
+        result = frd_mimo.__truediv__(tf_siso)
+        np.testing.assert_array_almost_equal(expected.omega, result.omega)
+        np.testing.assert_array_almost_equal(expected.fresp, result.fresp)
+
+        # Test division of MIMO FRD by SISO TF
+        result = frd_mimo.__truediv__(ss_siso)
+        np.testing.assert_array_almost_equal(expected.omega, result.omega)
+        np.testing.assert_array_almost_equal(expected.fresp, result.fresp)
+
+    @slycotonly
+    def test_rtruediv_mimo_siso(self):
+        omega = np.logspace(-1, 1, 10)
+        tf_mimo = TransferFunction([1], [1, 0]) * np.eye(2)
+        frd_mimo = frd(tf_mimo, omega)
+        ss_mimo = ct.tf2ss(tf_mimo)
+        tf_siso = TransferFunction([1], [1, 1])
+        frd_siso = frd(tf_siso, omega)
+        ss_siso = ct.tf2ss(tf_siso)
+        expected = frd(tf_siso.__rtruediv__(tf_mimo), omega)
+
+        # Test division of MIMO FRD by SISO FRD
+        result = frd_siso.__rtruediv__(frd_mimo)
+        np.testing.assert_array_almost_equal(expected.omega, result.omega)
+        np.testing.assert_array_almost_equal(expected.fresp, result.fresp)
+
+        # Test division of MIMO TF by SISO FRD
+        result = frd_siso.__rtruediv__(tf_mimo)
+        np.testing.assert_array_almost_equal(expected.omega, result.omega)
+        np.testing.assert_array_almost_equal(expected.fresp, result.fresp)
+
+        # Test division of MIMO SS by SISO FRD
+        result = frd_siso.__rtruediv__(ss_mimo)
+        np.testing.assert_array_almost_equal(expected.omega, result.omega)
+        np.testing.assert_array_almost_equal(expected.fresp, result.fresp)
+
+
+    @pytest.mark.parametrize(
+        "left, right, expected",
+        [
+            (
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                np.eye(3),
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[2], [0], [0]],
+                        [[0], [2], [0]],
+                        [[0], [0], [2]],
+                    ],
+                    [
+                        [[1, 0], [1], [1]],
+                        [[1], [1, 0], [1]],
+                        [[1], [1], [1, 0]],
+                    ],
+                ),
+            ),
+        ]
+    )
+    def test_rmul_mimo_siso(self, left, right, expected):
+        result = frd(right, np.logspace(-1, 1, 10)).__rmul__(left)
+        expected_frd = frd(expected, np.logspace(-1, 1, 10))
+        np.testing.assert_array_almost_equal(expected_frd.omega, result.omega)
+        np.testing.assert_array_almost_equal(expected_frd.fresp, result.fresp)
+
     def test_eval(self):
         sys_tf = ct.tf([1], [1, 2, 1])
         frd_tf = frd(sys_tf, np.logspace(-1, 1, 3))
diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py
index 22a946fe3..110c8dbc7 100644
--- a/control/tests/statefbk_test.py
+++ b/control/tests/statefbk_test.py
@@ -7,10 +7,10 @@
 import pytest
 import itertools
 import warnings
-from math import pi, atan
+from math import pi
 
 import control as ct
-from control import lqe, dlqe, poles, rss, ss, tf
+from control import poles, rss, ss, tf
 from control.exception import ControlDimension, ControlSlycot, \
     ControlArgument, slycot_check
 from control.mateqn import care, dare
@@ -57,6 +57,23 @@ def testCtrbT(self):
         Wc = ctrb(A, B, t=t)
         np.testing.assert_array_almost_equal(Wc, Wctrue)
 
+    def testCtrbNdim1(self):
+        # gh-1097: treat 1-dim B as nx1
+        A = np.array([[1., 2.], [3., 4.]])
+        B = np.array([5., 7.])
+        Wctrue = np.array([[5., 19.], [7., 43.]])
+        Wc = ctrb(A, B)
+        np.testing.assert_array_almost_equal(Wc, Wctrue)
+
+    def testCtrbRejectMismatch(self):
+        # gh-1097: check A, B for compatible shapes
+        with pytest.raises(ControlDimension, match='A must be a square matrix'):
+            ctrb([[1,2]],[1])
+        with pytest.raises(ControlDimension, match='Incompatible dimensions of B matrix'):
+            ctrb([[1,2],[2,3]], 1)
+        with pytest.raises(ControlDimension, match='Incompatible dimensions of B matrix'):
+            ctrb([[1,2],[2,3]], [[1,2]])
+
     def testObsvSISO(self):
         A = np.array([[1., 2.], [3., 4.]])
         C = np.array([[5., 7.]])
@@ -79,6 +96,23 @@ def testObsvT(self):
         Wo = obsv(A, C, t=t)
         np.testing.assert_array_almost_equal(Wo, Wotrue)
 
+    def testObsvNdim1(self):
+        # gh-1097: treat 1-dim C as 1xn
+        A = np.array([[1., 2.], [3., 4.]])
+        C = np.array([5., 7.])
+        Wotrue = np.array([[5., 7.], [26., 38.]])
+        Wo = obsv(A, C)
+        np.testing.assert_array_almost_equal(Wo, Wotrue)
+
+    def testObsvRejectMismatch(self):
+        # gh-1097: check A, C for compatible shapes
+        with pytest.raises(ControlDimension, match='A must be a square matrix'):
+            obsv([[1,2]],[1])
+        with pytest.raises(ControlDimension, match='Incompatible dimensions of C matrix'):
+            obsv([[1,2],[2,3]], 1)
+        with pytest.raises(ControlDimension, match='Incompatible dimensions of C matrix'):
+            obsv([[1,2],[2,3]], [[1],[2]])
+
     def testCtrbObsvDuality(self):
         A = np.array([[1.2, -2.3], [3.4, -4.5]])
         B = np.array([[5.8, 6.9], [8., 9.1]])
@@ -935,7 +969,6 @@ def unicycle_update(t, x, u, params):
         states=['x_', 'y_', 'theta_'],
         params={'a': 1})        # only used for testing params
 
-from math import pi
 
 @pytest.mark.parametrize("method", ['nearest', 'linear', 'cubic'])
 def test_gainsched_unicycle(unicycle, method):
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index a80168649..366c5ef3b 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -20,7 +20,7 @@
 from control.lti import evalfr
 from control.statesp import StateSpace, _convert_to_statespace, tf2ss, \
     _statesp_defaults, _rss_generate, linfnorm, ss, rss, drss
-from control.xferfcn import TransferFunction, ss2tf
+from control.xferfcn import TransferFunction, ss2tf, _tf_close_coeff
 
 from .conftest import editsdefaults, slycotonly
 
@@ -320,6 +320,343 @@ def test_multiply_ss(self, sys222, sys322):
         np.testing.assert_array_almost_equal(sys.C, C)
         np.testing.assert_array_almost_equal(sys.D, D)
 
+    @slycotonly
+    def test_add_sub_mimo_siso(self):
+        # Test SS with SS
+        ss_siso = StateSpace(
+            np.array([
+                [1, 2],
+                [3, 4],
+            ]),
+            np.array([
+                [1],
+                [4],
+            ]),
+            np.array([
+                [1, 1],
+            ]),
+            np.array([
+                [0],
+            ]),
+        )
+        ss_siso_1 = StateSpace(
+            np.array([
+                [1, 1],
+                [3, 1],
+            ]),
+            np.array([
+                [3],
+                [-4],
+            ]),
+            np.array([
+                [-1, 1],
+            ]),
+            np.array([
+                [0.1],
+            ]),
+        )
+        ss_siso_2 = StateSpace(
+            np.array([
+                [1, 0],
+                [0, 1],
+            ]),
+            np.array([
+                [0],
+                [2],
+            ]),
+            np.array([
+                [0, 1],
+            ]),
+            np.array([
+                [0],
+            ]),
+        )
+        ss_mimo = ss_siso_1.append(ss_siso_2)
+        expected_add = ct.combine_tf([
+            [ss2tf(ss_siso_1 + ss_siso), ss2tf(ss_siso)],
+            [ss2tf(ss_siso), ss2tf(ss_siso_2 + ss_siso)],
+        ])
+        expected_sub = ct.combine_tf([
+            [ss2tf(ss_siso_1 - ss_siso), -ss2tf(ss_siso)],
+            [-ss2tf(ss_siso), ss2tf(ss_siso_2 - ss_siso)],
+        ])
+        for op, expected in [
+            (StateSpace.__add__, expected_add),
+            (StateSpace.__radd__, expected_add),
+            (StateSpace.__sub__, expected_sub),
+            (StateSpace.__rsub__, -expected_sub),
+        ]:
+            result = op(ss_mimo, ss_siso)
+            assert _tf_close_coeff(
+                expected.minreal(),
+                ss2tf(result).minreal(),
+            )
+        # Test SS with array
+        expected_add = ct.combine_tf([
+            [ss2tf(1 + ss_siso), ss2tf(ss_siso)],
+            [ss2tf(ss_siso), ss2tf(1 + ss_siso)],
+        ])
+        expected_sub = ct.combine_tf([
+            [ss2tf(-1 + ss_siso), ss2tf(ss_siso)],
+            [ss2tf(ss_siso), ss2tf(-1 + ss_siso)],
+        ])
+        for op, expected in [
+            (StateSpace.__add__, expected_add),
+            (StateSpace.__radd__, expected_add),
+            (StateSpace.__sub__, expected_sub),
+            (StateSpace.__rsub__, -expected_sub),
+        ]:
+            result = op(ss_siso, np.eye(2))
+            assert _tf_close_coeff(
+                expected.minreal(),
+                ss2tf(result).minreal(),
+            )
+
+    @slycotonly
+    @pytest.mark.parametrize(
+        "left, right, expected",
+        [
+            (
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                TransferFunction([2], [1, 0]),
+                np.eye(3),
+                TransferFunction(
+                    [
+                        [[2], [0], [0]],
+                        [[0], [2], [0]],
+                        [[0], [0], [2]],
+                    ],
+                    [
+                        [[1, 0], [1], [1]],
+                        [[1], [1, 0], [1]],
+                        [[1], [1], [1, 0]],
+                    ],
+                ),
+            ),
+        ]
+    )
+    def test_mul_mimo_siso(self, left, right, expected):
+        result = tf2ss(left).__mul__(right)
+        assert _tf_close_coeff(
+            expected.minreal(),
+            ss2tf(result).minreal(),
+        )
+
+    @slycotonly
+    @pytest.mark.parametrize(
+        "left, right, expected",
+        [
+            (
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                np.eye(3),
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[2], [0], [0]],
+                        [[0], [2], [0]],
+                        [[0], [0], [2]],
+                    ],
+                    [
+                        [[1, 0], [1], [1]],
+                        [[1], [1, 0], [1]],
+                        [[1], [1], [1, 0]],
+                    ],
+                ),
+            ),
+        ]
+    )
+    def test_rmul_mimo_siso(self, left, right, expected):
+        result = tf2ss(right).__rmul__(left)
+        assert _tf_close_coeff(
+            expected.minreal(),
+            ss2tf(result).minreal(),
+        )
+
+    @slycotonly
+    def test_pow(self, sys222, sys322):
+        """Test state space powers."""
+        for sys in [sys222, sys322]:
+            # Power of 0
+            result = sys**0
+            expected = StateSpace([], [], [], np.eye(2), dt=0)
+            np.testing.assert_allclose(expected.A, result.A)
+            np.testing.assert_allclose(expected.B, result.B)
+            np.testing.assert_allclose(expected.C, result.C)
+            np.testing.assert_allclose(expected.D, result.D)
+            # Power of 1
+            result = sys**1
+            expected = sys
+            np.testing.assert_allclose(expected.A, result.A)
+            np.testing.assert_allclose(expected.B, result.B)
+            np.testing.assert_allclose(expected.C, result.C)
+            np.testing.assert_allclose(expected.D, result.D)
+            # Power of -1 (inverse of biproper system)
+            # Testing transfer function representations to avoid the
+            # non-uniqueness of the state-space representation. Once MIMO
+            # canonical forms are supported, can check canonical state-space
+            # matrices instead.
+            result = (sys * sys**-1).minreal()
+            expected = StateSpace([], [], [], np.eye(2), dt=0)
+            assert _tf_close_coeff(
+                ss2tf(expected).minreal(),
+                ss2tf(result).minreal(),
+            )
+            result = (sys**-1 * sys).minreal()
+            expected = StateSpace([], [], [], np.eye(2), dt=0)
+            assert _tf_close_coeff(
+                ss2tf(expected).minreal(),
+                ss2tf(result).minreal(),
+            )
+            # Power of 3
+            result = sys**3
+            expected = sys * sys * sys
+            np.testing.assert_allclose(expected.A, result.A)
+            np.testing.assert_allclose(expected.B, result.B)
+            np.testing.assert_allclose(expected.C, result.C)
+            np.testing.assert_allclose(expected.D, result.D)
+            # Power of -3
+            result = sys**-3
+            expected = sys**-1 * sys**-1 * sys**-1
+            np.testing.assert_allclose(expected.A, result.A)
+            np.testing.assert_allclose(expected.B, result.B)
+            np.testing.assert_allclose(expected.C, result.C)
+            np.testing.assert_allclose(expected.D, result.D)
+
+    @slycotonly
+    def test_truediv(self, sys222, sys322):
+        """Test state space truediv"""
+        for sys in [sys222, sys322]:
+            # Divide by self
+            result = (sys.__truediv__(sys)).minreal()
+            expected = StateSpace([], [], [], np.eye(2), dt=0)
+            assert _tf_close_coeff(
+                ss2tf(expected).minreal(),
+                ss2tf(result).minreal(),
+            )
+            # Divide by TF
+            result = sys.__truediv__(TransferFunction.s)
+            expected = ss2tf(sys) / TransferFunction.s
+            assert _tf_close_coeff(
+                expected.minreal(),
+                ss2tf(result).minreal(),
+            )
+
+    @slycotonly
+    def test_rtruediv(self, sys222, sys322):
+        """Test state space rtruediv"""
+        for sys in [sys222, sys322]:
+            result = (sys.__rtruediv__(sys)).minreal()
+            expected = StateSpace([], [], [], np.eye(2), dt=0)
+            assert _tf_close_coeff(
+                ss2tf(expected).minreal(),
+                ss2tf(result).minreal(),
+            )
+            # Divide TF by SS
+            result = sys.__rtruediv__(TransferFunction.s)
+            expected = TransferFunction.s / sys
+            assert _tf_close_coeff(
+                expected.minreal(),
+                result.minreal(),
+            )
+        # Divide array by SS
+        sys = tf2ss(TransferFunction([1, 2], [2, 1]))
+        result = sys.__rtruediv__(np.eye(2))
+        expected = TransferFunction([2, 1], [1, 2]) * np.eye(2)
+        assert _tf_close_coeff(
+            expected.minreal(),
+            ss2tf(result).minreal(),
+        )
+
     @pytest.mark.parametrize("k", [2, -3.141, np.float32(2.718), np.array([[4.321], [5.678]])])
     def test_truediv_ss_scalar(self, sys322, k):
         """Divide SS by scalar."""
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index c7f92379a..cab556cc2 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -15,7 +15,7 @@
                      ss, ss2tf, tf, tf2ss, zpk)
 from control.statesp import _convert_to_statespace
 from control.tests.conftest import slycotonly
-from control.xferfcn import _convert_to_transfer_function
+from control.xferfcn import _convert_to_transfer_function, _tf_close_coeff
 
 
 class TestXferFcn:
@@ -106,9 +106,28 @@ def test_constructor_double_dt(self):
 
     def test_add_inconsistent_dimension(self):
         """Add two transfer function matrices of different sizes."""
-        sys1 = TransferFunction([[[1., 2.]]], [[[4., 5.]]])
-        sys2 = TransferFunction([[[4., 3.]], [[1., 2.]]],
-                                [[[1., 6.]], [[2., 4.]]])
+        sys1 = TransferFunction(
+            [
+                [[1., 2.]],
+                [[2., -2.]],
+                [[2., 1.]],
+            ],
+            [
+                [[4., 5.]],
+                [[5., 2.]],
+                [[3., 2.]],
+            ],
+        )
+        sys2 = TransferFunction(
+            [
+                [[4., 3.]],
+                [[1., 2.]],
+            ],
+            [
+                [[1., 6.]],
+                [[2., 4.]],
+            ]
+        )
         with pytest.raises(ValueError):
             sys1.__add__(sys2)
         with pytest.raises(ValueError):
@@ -385,6 +404,239 @@ def test_pow(self):
         with pytest.raises(ValueError):
             TransferFunction.__pow__(sys1, 0.5)
 
+    def test_add_sub_mimo_siso(self):
+        for op in [
+            TransferFunction.__add__,
+            TransferFunction.__radd__,
+            TransferFunction.__sub__,
+            TransferFunction.__rsub__,
+        ]:
+            tf_mimo = TransferFunction(
+                [
+                    [[1], [1]],
+                    [[1], [1]],
+                ],
+                [
+                    [[10, 1], [20, 1]],
+                    [[20, 1], [30, 1]],
+                ],
+            )
+            tf_siso = TransferFunction([1], [5, 0])
+            tf_arr = ct.split_tf(tf_mimo)
+            expected = ct.combine_tf([
+                [op(tf_arr[0, 0], tf_siso), op(tf_arr[0, 1], tf_siso)],
+                [op(tf_arr[1, 0], tf_siso), op(tf_arr[1, 1], tf_siso)],
+            ])
+            result = op(tf_mimo, tf_siso)
+            assert _tf_close_coeff(expected.minreal(), result.minreal())
+
+    @pytest.mark.parametrize(
+        "left, right, expected",
+        [
+            (
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                TransferFunction([2], [1, 0]),
+                np.eye(3),
+                TransferFunction(
+                    [
+                        [[2], [0], [0]],
+                        [[0], [2], [0]],
+                        [[0], [0], [2]],
+                    ],
+                    [
+                        [[1, 0], [1], [1]],
+                        [[1], [1, 0], [1]],
+                        [[1], [1], [1, 0]],
+                    ],
+                ),
+            ),
+        ]
+    )
+    def test_mul_mimo_siso(self, left, right, expected):
+        """Test multiplication of a MIMO and a SISO system."""
+        result = left.__mul__(right)
+        assert _tf_close_coeff(expected.minreal(), result.minreal())
+
+    @pytest.mark.parametrize(
+        "left, right, expected",
+        [
+            (
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                TransferFunction(
+                    [
+                        [[2], [1]],
+                        [[-1], [4]],
+                    ],
+                    [
+                        [[10, 1], [20, 1]],
+                        [[20, 1], [30, 1]],
+                    ],
+                ),
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[4], [2]],
+                        [[-2], [8]],
+                    ],
+                    [
+                        [[10, 1, 0], [20, 1, 0]],
+                        [[20, 1, 0], [30, 1, 0]],
+                    ],
+                ),
+            ),
+            (
+                np.eye(3),
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[2], [0], [0]],
+                        [[0], [2], [0]],
+                        [[0], [0], [2]],
+                    ],
+                    [
+                        [[1, 0], [1], [1]],
+                        [[1], [1, 0], [1]],
+                        [[1], [1], [1, 0]],
+                    ],
+                ),
+            ),
+        ]
+    )
+    def test_rmul_mimo_siso(self, left, right, expected):
+        """Test right multiplication of a MIMO and a SISO system."""
+        result = right.__rmul__(left)
+        assert _tf_close_coeff(expected.minreal(), result.minreal())
+
+    @pytest.mark.parametrize(
+        "left, right, expected",
+        [
+            (
+                TransferFunction(
+                    [
+                        [[1], [0], [0]],
+                        [[0], [2], [0]],
+                        [[0], [0], [3]],
+                    ],
+                    [
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                        [[1], [1], [1]],
+                    ],
+                ),
+                TransferFunction([-2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[1, 0], [0], [0]],
+                        [[0], [2, 0], [0]],
+                        [[0], [0], [3, 0]],
+                    ],
+                    [
+                        [[-2], [1], [1]],
+                        [[1], [-2], [1]],
+                        [[1], [1], [-2]],
+                    ],
+                ),
+            ),
+        ]
+    )
+    def test_truediv_mimo_siso(self, left, right, expected):
+        """Test true division of a MIMO and a SISO system."""
+        result = left.__truediv__(right)
+        assert _tf_close_coeff(expected.minreal(), result.minreal())
+
+    @pytest.mark.parametrize(
+        "left, right, expected",
+        [
+            (
+                np.eye(3),
+                TransferFunction([2], [1, 0]),
+                TransferFunction(
+                    [
+                        [[1, 0], [0], [0]],
+                        [[0], [1, 0], [0]],
+                        [[0], [0], [1, 0]],
+                    ],
+                    [
+                        [[2], [1], [1]],
+                        [[1], [2], [1]],
+                        [[1], [1], [2]],
+                    ],
+                ),
+            ),
+        ]
+    )
+    def test_rtruediv_mimo_siso(self, left, right, expected):
+        """Test right true division of a MIMO and a SISO system."""
+        result = right.__rtruediv__(left)
+        assert _tf_close_coeff(expected.minreal(), result.minreal())
+
     @pytest.mark.parametrize("named", [False, True])
     def test_slice(self, named):
         sys = TransferFunction(
@@ -638,6 +890,52 @@ def test_feedback_siso(self):
         np.testing.assert_allclose(sys4.num, [[[-1., 7., -16., 16., 0.]]])
         np.testing.assert_allclose(sys4.den, [[[1., 0., 2., -8., 8., 0.]]])
 
+    def test_append(self):
+        """Test ``TransferFunction.append()``."""
+        tf1 = TransferFunction(
+            [
+                [[1], [1]]
+            ],
+            [
+                [[10, 1], [20, 1]]
+            ],
+        )
+        tf2 = TransferFunction(
+            [
+                [[2], [2]]
+            ],
+            [
+                [[10, 1], [1, 1]]
+            ],
+        )
+        tf3 = TransferFunction([100], [100, 1])
+        tf_exp_1 = TransferFunction(
+            [
+                [[1], [1], [0], [0]],
+                [[0], [0], [2], [2]],
+            ],
+            [
+                [[10, 1], [20, 1], [1], [1]],
+                [[1], [1], [10, 1], [1, 1]],
+            ],
+        )
+        tf_exp_2 = TransferFunction(
+            [
+                [[1], [1], [0], [0], [0]],
+                [[0], [0], [2], [2], [0]],
+                [[0], [0], [0], [0], [100]],
+            ],
+            [
+                [[10, 1], [20, 1], [1], [1], [1]],
+                [[1], [1], [10, 1], [1, 1], [1]],
+                [[1], [1], [1], [1], [100, 1]],
+            ],
+        )
+        tf_appended_1 = tf1.append(tf2)
+        assert _tf_close_coeff(tf_exp_1, tf_appended_1)
+        tf_appended_2 = tf1.append(tf2).append(tf3)
+        assert _tf_close_coeff(tf_exp_2, tf_appended_2)
+
     @slycotonly
     def test_convert_to_transfer_function(self):
         """Test for correct state space to transfer function conversion."""