python-control · murrayrm · Jan 20, 2025 · Dec 7, 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)

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, _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/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/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
@@ -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',