fix discete time simulation time step issue #332 (#356)

murrayrm · web-flow · commit ccbd7acf4ed8 · 2020-01-02T09:09:14.000-08:00
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
@@ -471,6 +471,45 @@ def test_time_vector(self):
                                           squeeze=False)
         self.assertTrue(isinstance(context.exception, ValueError))
 
+    def test_discrete_time_steps(self):
+        """Make sure rounding errors in sample time are handled properly"""
+        # See https://github.com/python-control/python-control/issues/332)
+        #
+        # These tests play around with the input time vector to make sure that
+        # small rounding errors don't generate spurious errors.
+
+        # Discrete time system to use for simulation
+        # self.siso_dtf2 = TransferFunction([1], [1, 1, 0.25], 0.2)
+
+        # Set up a time range and simulate
+        T = np.arange(0, 100, 0.2)
+        tout1, yout1 = step_response(self.siso_dtf2, T)
+
+        # Simulate every other time step
+        T = np.arange(0, 100, 0.4)
+        tout2, yout2 = step_response(self.siso_dtf2, T)
+        np.testing.assert_array_almost_equal(tout1[::2], tout2)
+        np.testing.assert_array_almost_equal(yout1[::2], yout2)
+
+        # Add a small error into some of the time steps
+        T = np.arange(0, 100, 0.2)
+        T[1:-2:2] -= 1e-12      # tweak second value and a few others
+        tout3, yout3 = step_response(self.siso_dtf2, T)
+        np.testing.assert_array_almost_equal(tout1, tout3)
+        np.testing.assert_array_almost_equal(yout1, yout3)
+
+        # Add a small error into some of the time steps (w/ skipping)
+        T = np.arange(0, 100, 0.4)
+        T[1:-2:2] -= 1e-12      # tweak second value and a few others
+        tout4, yout4 = step_response(self.siso_dtf2, T)
+        np.testing.assert_array_almost_equal(tout2, tout4)
+        np.testing.assert_array_almost_equal(yout2, yout4)
+
+        # Make sure larger errors *do* generate an error
+        T = np.arange(0, 100, 0.2)
+        T[1:-2:2] -= 1e-3      # change second value and a few others
+        self.assertRaises(ValueError, step_response, self.siso_dtf2, T)
+
     def test_time_series_data_convention(self):
         """Make sure time series data matches documentation conventions"""
         # SISO continuous time
diff --git a/control/timeresp.py b/control/timeresp.py
@@ -365,7 +365,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
             # Make sure that the time increment is a multiple of sampling time
 
             # First make sure that time increment is bigger than sampling time
-            if dt < sys.dt:
+            # (with allowance for small precision errors)
+            if dt < sys.dt and not np.isclose(dt, sys.dt):
                 raise ValueError("Time steps ``T`` must match sampling time")
 
             # Now check to make sure it is a multiple (with check against
@@ -374,19 +375,21 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
                       np.isclose(dt % sys.dt, sys.dt)):
                 raise ValueError("Time steps ``T`` must be multiples of "
                                  "sampling time")
+            sys_dt = sys.dt
+
         else:
-            sys.dt = dt         # For unspecified sampling time, use time incr
+            sys_dt = dt         # For unspecified sampling time, use time incr
 
         # Discrete time simulation using signal processing toolbox
-        dsys = (A, B, C, D, sys.dt)
+        dsys = (A, B, C, D, sys_dt)
 
         # Use signal processing toolbox for the discrete time simulation
         # Transpose the input to match toolbox convention
         tout, yout, xout = sp.signal.dlsim(dsys, np.transpose(U), T, X0)
 
         if not interpolate:
             # If dt is different from sys.dt, resample the output
-            inc = int(round(dt / sys.dt))
+            inc = int(round(dt / sys_dt))
             tout = T            # Return exact list of time steps
             yout = yout[::inc, :]
             xout = xout[::inc, :]
