python-control · murrayrm · Apr 2, 2022 · Mar 12, 2022 · Mar 23, 2022 · Mar 16, 2022
diff --git a/control/__init__.py b/control/__init__.py
@@ -61,6 +61,7 @@
 from .rlocus import *
 from .statefbk import *
 from .statesp import *
+from .stochsys import *
 from .timeresp import *
 from .xferfcn import *
 from .ctrlutil import *

diff --git a/control/iosys.py b/control/iosys.py
@@ -1585,11 +1585,17 @@ def input_output_response(
     T : array-like
         Time steps at which the input is defined; values must be evenly spaced.
 
-    U : array-like or number, optional
-        Input array giving input at each time `T` (default = 0).
-
-    X0 : array-like or number, optional
-        Initial condition (default = 0).
+    U : array-like, list, or number, optional
+        Input array giving input at each time `T` (default = 0).  If a list
+        is specified, each element in the list will be treated as a portion
+        of the input and broadcast (if necessary) to match the time vector.
+
+    X0 : array-like, list, or number, optional
+        Initial condition (default = 0).  If a list is given, each element
+        in the list will be flattened and stacked into the initial
+        condition.  If a smaller number of elements are given that the
+        number of states in the system, the initial condition will be padded
+        with zeros.
 
     return_x : bool, optional
         If True, return the state vector when assigning to a tuple (default =
@@ -1641,6 +1647,16 @@ def input_output_response(
     ValueError
         If time step does not match sampling time (for discrete time systems).
 
+    Notes
+    -----
+    1. If a smaller number of initial conditions are given than the number of
+       states in the system, the initial conditions will be padded with
+       zeros.  This is often useful for interconnected control systems where
+       the process dynamics are the first system and all other components
+       start with zero initial condition since this can be specified as
+       [xsys_0, 0].  A warning is issued if the initial conditions are padded
+       and and the final listed initial state is not zero.
+
     """
     #
     # Process keyword arguments
@@ -1656,14 +1672,14 @@ def input_output_response(
             raise ValueError("ivp_method specified more than once")
         solve_ivp_kwargs['method'] = kwargs.pop('solve_ivp_method')
 
-    # Make sure there were no extraneous keywords
-    if kwargs:
-        raise TypeError("unrecognized keywords: ", str(kwargs))
-
     # Set the default method to 'RK45'
     if solve_ivp_kwargs.get('method', None) is None:
         solve_ivp_kwargs['method'] = 'RK45'
 
+    # Make sure there were no extraneous keywords
+    if kwargs:
+        raise TypeError("unrecognized keyword(s): ", str(kwargs))
+
     # Sanity checking on the input
     if not isinstance(sys, InputOutputSystem):
         raise TypeError("System of type ", type(sys), " not valid")
@@ -1683,19 +1699,75 @@ def input_output_response(
         # Use the input time points as the output time points
         t_eval = T
 
-    # Check and convert the input, if needed
-    # TODO: improve MIMO ninputs check (choose from U)
+    # If we were passed a list of input, concatenate them (w/ broadcast)
+    if isinstance(U, (tuple, list)) and len(U) != ntimepts:
+        U_elements = []
+        for i, u in enumerate(U):
+            u = np.array(u)     # convert everyting to an array
+            # Process this input
+            if u.ndim == 0 or (u.ndim == 1 and u.shape[0] != T.shape[0]):
+                # Broadcast array to the length of the time input
+                u = np.outer(u, np.ones_like(T))
+
+            elif (u.ndim == 1 and u.shape[0] == T.shape[0]) or \
+                 (u.ndim == 2 and u.shape[1] == T.shape[0]):
+                # No processing necessary; just stack
+                pass
+
+            else:
+                raise ValueError(f"Input element {i} has inconsistent shape")
+
+            # Append this input to our list
+            U_elements.append(u)
+
+        # Save the newly created input vector
+        U = np.vstack(U_elements)
+
+    # Make sure the input has the right shape
     if sys.ninputs is None or sys.ninputs == 1:
         legal_shapes = [(ntimepts,), (1, ntimepts)]
     else:
         legal_shapes = [(sys.ninputs, ntimepts)]
+
     U = _check_convert_array(U, legal_shapes,
                              'Parameter ``U``: ', squeeze=False)
+
+    # Always store the input as a 2D array
     U = U.reshape(-1, ntimepts)
     ninputs = U.shape[0]
 
-    # create X0 if not given, test if X0 has correct shape
+    # If we were passed a list of initial states, concatenate them
+    if isinstance(X0, (tuple, list)):
+        X0_list = []
+        for i, x0 in enumerate(X0):
+            x0 = np.array(x0).reshape(-1)       # convert everyting to 1D array
+            X0_list += x0.tolist()              # add elements to initial state
+
+        # Save the newly created input vector
+        X0 = np.array(X0_list)
+
+    # If the initial state is too short, make it longer (NB: sys.nstates
+    # could be None if nstates comes from size of initial condition)
+    if sys.nstates and isinstance(X0, np.ndarray) and X0.size < sys.nstates:
+        if X0[-1] != 0:
+            warn("initial state too short; padding with zeros")
+        X0 = np.hstack([X0, np.zeros(sys.nstates - X0.size)])
+
+    # Check to make sure this is not a static function
     nstates = _find_size(sys.nstates, X0)
+    if nstates == 0:
+        # No states => map input to output
+        u = U[0] if len(U.shape) == 1 else U[:, 0]
+        y = np.zeros((np.shape(sys._out(T[0], X0, u))[0], len(T)))
+        for i in range(len(T)):
+            u = U[i] if len(U.shape) == 1 else U[:, i]
+            y[:, i] = sys._out(T[i], [], u)
+        return TimeResponseData(
+            T, y, None, U, issiso=sys.issiso(),
+            output_labels=sys.output_index, input_labels=sys.input_index,
+            transpose=transpose, return_x=return_x, squeeze=squeeze)
+
+    # create X0 if not given, test if X0 has correct shape
     X0 = _check_convert_array(X0, [(nstates,), (nstates, 1)],
                               'Parameter ``X0``: ', squeeze=True)
 

diff --git a/control/optimal.py b/control/optimal.py
@@ -776,7 +776,7 @@ def create_mpc_iosystem(self):
         """Create an I/O system implementing an MPC controller"""
         # Check to make sure we are in discrete time
         if self.system.dt == 0:
-            raise ControlNotImplemented(
+            raise ct.ControlNotImplemented(
                 "MPC for continuous time systems not implemented")
 
         def _update(t, x, u, params={}):

diff --git a/control/sisotool.py b/control/sisotool.py
@@ -215,14 +215,14 @@ def rootlocus_pid_designer(plant, gain='P', sign=+1, input_signal='r',
     derivative terms are given instead by Kp, Ki*dt/2*(z+1)/(z-1), and
     Kd/dt*(z-1)/z, respectively.
 
-        ------> C_ff ------    d
-        |                 |    |
-    r   |     e           V    V  u         y
-    ------->O---> C_f --->O--->O---> plant --->
-            ^-            ^-                |
-            |             |                 |
-            |             ----- C_b <-------|
-            ---------------------------------
+          ------> C_ff ------    d
+          |                 |    |
+      r   |     e           V    V  u         y
+      ------->O---> C_f --->O--->O---> plant --->
+              ^-            ^-                |
+              |             |                 |
+              |             ----- C_b <-------|
+              ---------------------------------
 
     It is also possible to move the derivative term into the feedback path
     `C_b` using `derivative_in_feedback_path=True`. This may be desired to
@@ -234,8 +234,8 @@ def rootlocus_pid_designer(plant, gain='P', sign=+1, input_signal='r',
 
     Remark: It may be helpful to zoom in using the magnifying glass on the
     plot. Just ake sure to deactivate magnification mode when you are done by
-    clicking the magnifying glass. Otherwise you will not be able to be able to choose
-    a gain on the root locus plot.
+    clicking the magnifying glass. Otherwise you will not be able to be able
+    to choose a gain on the root locus plot.
 
     Parameters
     ----------
@@ -269,6 +269,7 @@ def rootlocus_pid_designer(plant, gain='P', sign=+1, input_signal='r',
     ----------
     closedloop : class:`StateSpace` system
         The closed-loop system using initial gains.
+
     """
 
     plant = _convert_to_statespace(plant)