murrayrm · Jan 20, 2025 · Dec 7, 2024 · Dec 9, 2024 · Dec 9, 2024 · Dec 9, 2024
diff --git 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)
@@ -541,26 +542,20 @@ def combine_tf(tf_array):
     --------
     Combine two transfer functions
 
-    >>> s = control.TransferFunction.s
-    >>> control.combine_tf([
+    >>> s = ct.TransferFunction.s
+    >>> ct.combine_tf([
     ...     [1 / (s + 1)],
     ...     [s / (s + 2)],
     ... ])
-    TransferFunction([[array([1])], [array([1, 0])]],
-                     [[array([1, 1])], [array([1, 2])]])
+    TransferFunction([[array([1])], [array([1, 0])]], [[array([1, 1])], [array([1, 2])]])
 
     Combine NumPy arrays with transfer functions
 
-    >>> control.combine_tf([
+    >>> ct.combine_tf([
     ...     [np.eye(2), np.zeros((2, 1))],
-    ...     [np.zeros((1, 2)), control.TransferFunction([1], [1, 0])],
+    ...     [np.zeros((1, 2)), ct.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])]])
+    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 = []
@@ -635,7 +630,7 @@ def split_tf(transfer_function):
     --------
     Split a MIMO transfer function
 
-    >>> G = control.TransferFunction(
+    >>> G = ct.TransferFunction(
     ...     [
     ...         [[87.8], [-86.4]],
     ...         [[108.2], [-109.6]],
@@ -645,7 +640,7 @@ def split_tf(transfer_function):
     ...         [[1, 1], [1, 1]],
     ...     ],
     ... )
-    >>> control.split_tf(G)
+    >>> ct.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])),

diff --git 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, _process_iosys_keywords, \
     _process_subsys_index, common_timebase
@@ -405,6 +406,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(
@@ -442,6 +449,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(
@@ -469,6 +482,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(
@@ -497,11 +516,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
@@ -516,11 +533,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
 
@@ -753,6 +768,35 @@ 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/statesp.py b/control/statesp.py
@@ -63,6 +63,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, \
@@ -611,6 +612,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
@@ -621,6 +625,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)):
@@ -675,6 +685,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
@@ -686,6 +700,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(
@@ -725,23 +745,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.
 
@@ -1146,7 +1210,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/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[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