Fixing bug in timeresp.forced_response that caused scipy interp1d ValueError when applied to discrete sys with T=None.

Andrew D. McGuire · Andrew D. McGuire · commit efc42c8ab132 · 2019-08-08T14:13:42.000-04:00
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
@@ -259,7 +259,7 @@ def test_forced_response(self):
                               42.3227, 44.9694, 47.1599, 48.9776]])
         _t, yout, _xout = forced_response(self.mimo_ss1, t, u, x0)
         np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
-        
+
         # Test discrete MIMO system to use correct convention for input
         sysc = self.mimo_ss1
         dt=t[1]-t[0]
@@ -270,6 +270,17 @@ def test_forced_response(self):
         np.testing.assert_array_equal(youtc.shape, youtd.shape)
         np.testing.assert_array_almost_equal(youtc, youtd, decimal=4)
 
+        # Test discrete MIMO system without default T argument
+        u = np.array([[0., 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                      [1., 1, 1, 1, 1, 1, 1, 1, 1, 1]])
+        x0 = np.matrix(".5; 1; 0; 0")
+        youttrue = np.array([[11., 8.1494, 5.9361, 4.2258, 2.9118, 1.9092,
+                              1.1508, 0.5833, 0.1645, -0.1391],
+                             [9., 17.6457, 24.7072, 30.4855, 35.2234, 39.1165,
+                              42.3227, 44.9694, 47.1599, 48.9776]])
+        _t, yout, _xout = forced_response(sysd, U=u, X0=x0)
+        np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
+
     def test_lsim_double_integrator(self):
         # Note: scipy.signal.lsim fails if A is not invertible
         A = np.mat("0. 1.;0. 0.")
diff --git a/control/timeresp.py b/control/timeresp.py
@@ -200,7 +200,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     sys: LTI (StateSpace, or TransferFunction)
         LTI system to simulate
 
-    T: array-like
+    T: array-like, optional for discrete LTI `sys`
         Time steps at which the input is defined; values must be evenly spaced.
 
     U: array-like or number, optional
@@ -260,20 +260,29 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     n_inputs = B.shape[1]
     n_outputs = C.shape[0]
 
+    # Convert inputs to numpy arrays for easier shape checking
+    if U is not None:
+        U = np.asarray(U)
+    if T is not None:
+        T = np.asarray(T)
+
     # Set and/or check time vector in discrete time case
     if isdtime(sys, strict=True):
         if T is None:
             if U is None:
                 raise ValueError('Parameters ``T`` and ``U`` can\'t both be'
                                  'zero for discrete-time simulation')
             # Set T to equally spaced samples with same length as U
-            T = np.array(range(len(U))) * (1 if sys.dt is True else sys.dt)
+            if len(U.shape) == 1:
+                n_steps = U.shape[0]
+            else:
+                n_steps = U.shape[1]
+            T = np.array(range(n_steps)) * (1 if sys.dt is True else sys.dt)
         else:
             # Make sure the input vector and time vector have same length
-            # TODO: allow interpolation of the input vector
-            if len(U) != len(T):
-                ValueError('Pamameter ``T`` must have same length as'
-                           'input vector ``U``')
+            if (len(U.shape) == 1 and U.shape[0] != T.shape[0]) or (len(U.shape) > 1 and U.shape[1] != T.shape[0]):
+                ValueError('Pamameter ``T`` must have same elements as'
+                           ' the number of columns in input array ``U``')
 
     # Test if T has shape (n,) or (1, n);
     # T must be array-like and values must be increasing.
@@ -288,7 +297,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     if not np.allclose(T[1:] - T[:-1], dt):
         raise ValueError("Parameter ``T``: time values must be "
                          "equally spaced.")
-    n_steps = len(T)            # number of simulation steps
+    n_steps = T.shape[0]            # number of simulation steps
 
     # create X0 if not given, test if X0 has correct shape
     X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)],
