ensure dlsim output length

bnavigator · bnavigator · commit 9cac9fa0daa4 · 2021-04-01T13:11:06.000+02:00
diff --git a/control/timeresp.py b/control/timeresp.py
@@ -209,7 +209,9 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         LTI system to simulate
 
     T : array_like, optional for discrete LTI `sys`
-        Time steps at which the input is defined; values must be evenly spaced.
+        Time steps at which the input is defined; values must be evenly spaced
+        and start with 0.
+
         If None, `U` must be given and `len(U)` time steps of sys.dt are
         simulated. If sys.dt is None or True (undetermined time step), a time
         step of 1.0 is assumed.
@@ -355,13 +357,14 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
                              'Parameter ``T``: ', squeeze=True,
                              transpose=transpose)
 
-    # equally spaced also implies strictly monotonic increase
-    dt = T[1] - T[0]
-    if not np.allclose(np.diff(T), dt):
-        raise ValueError("Parameter ``T``: time values must be "
-                         "equally spaced.")
     n_steps = T.shape[0]            # number of simulation steps
 
+    # equally spaced also implies strictly monotonic increase,
+    dt = T[-1] / (n_steps - 1)
+    if not np.allclose(np.diff(T), dt):
+        raise ValueError("Parameter ``T`` must start with 0 and time values "
+                         "must be equally spaced.")
+
     # create X0 if not given, test if X0 has correct shape
     X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)],
                               'Parameter ``X0``: ', squeeze=True)
@@ -432,12 +435,21 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
             # Now check to make sure it is a multiple (with check against
             # sys.dt because floating point mod can have small errors
-            elif not (np.isclose(dt % sys.dt, 0) or
-                      np.isclose(dt % sys.dt, sys.dt)):
+            if not (np.isclose(dt % sys.dt, 0) or
+                    np.isclose(dt % sys.dt, sys.dt)):
                 raise ValueError("Time steps ``T`` must be multiples of "
                                  "sampling time")
             sys_dt = sys.dt
 
+            # sp.signal.dlsim returns not enough samples if
+            # T[-1] - T[0] < sys_dt * decimation * (n_steps - 1)
+            # due to rounding errors.
+            # https://github.com/scipyscipy/blob/v1.6.1/scipy/signal/ltisys.py#L3462
+            scipy_out_samples = int(np.floor(T[-1] / sys_dt)) + 1
+            if scipy_out_samples < n_steps:
+                # parantheses: order of evaluation is important
+                T[-1] = T[-1] * (n_steps / (T[-1] / sys_dt + 1))
+
         else:
             sys_dt = dt         # For unspecified sampling time, use time incr
 
@@ -459,7 +471,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         xout = np.transpose(xout)
         yout = np.transpose(yout)
 
-    return _process_time_response(sys, tout, yout, xout, transpose=transpose,
+    return _process_time_response(sys, tout[:n_steps], yout[:, :n_steps],
+                                  xout[:, :n_steps], transpose=transpose,
                                   return_x=return_x, squeeze=squeeze)
 
 
