Improved default time vector and handling for time response functions step, impulse, and initial (#420)

sawyerbfuller · web-flow · commit 155d8ebca7f7 · 2020-07-11T08:59:53.000-07:00
* fixed default response time for time response of discrete-time functions step, impulse, and initial

* to pass tests, added convenient ability to specify simulation time and number of steps rather than complete time vector in timeresponse functions

* eliminated deprecation warnings by importing certain functions from numpy instead of scipy

* eliminated deprecation warnings by importing certain functions from numpy instead of scipy

* small fix to pass unit tests

* adjusted sisotool test so tests pass with new default step response time window

* added functionality to automatically choose dt in timeresp.py based on system poles. and unit tests.

* removed some leftover code and comments

* explanation in docstrings for how time vector T is auto-computed in time response functions
diff --git a/control/freqplot.py b/control/freqplot.py
@@ -822,10 +822,10 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None,
 
     # Set the range to be an order of magnitude beyond any features
     if number_of_samples:
-        omega = sp.logspace(
+        omega = np.logspace(
             lsp_min, lsp_max, num=number_of_samples, endpoint=True)
     else:
-        omega = sp.logspace(lsp_min, lsp_max, endpoint=True)
+        omega = np.logspace(lsp_min, lsp_max, endpoint=True)
     return omega
 
 
diff --git a/control/matlab/timeresp.py b/control/matlab/timeresp.py
@@ -21,8 +21,9 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
     sys: StateSpace, or TransferFunction
         LTI system to simulate
 
-    T: array-like object, optional
-        Time vector (argument is autocomputed if not given)
+    T: array-like or number, optional
+        Time vector, or simulation time duration if a number (time vector is 
+        autocomputed if not given)
 
     X0: array-like or number, optional
         Initial condition (default = 0)
@@ -59,7 +60,7 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
     from ..timeresp import step_response
 
     T, yout, xout = step_response(sys, T, X0, input, output,
-                                  transpose = True, return_x=True)
+                                  transpose=True, return_x=True)
 
     if return_x:
         return yout, T, xout
@@ -75,8 +76,9 @@ def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1,0.9)):
     sys: StateSpace, or TransferFunction
         LTI system to simulate
 
-    T: array-like object, optional
-        Time vector (argument is autocomputed if not given)
+    T: array-like or number, optional
+        Time vector, or simulation time duration if a number (time vector is 
+        autocomputed if not given)
 
     SettlingTimeThreshold: float value, optional
         Defines the error to compute settling time (default = 0.02)
@@ -127,9 +129,10 @@ def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False):
     sys: StateSpace, TransferFunction
         LTI system to simulate
 
-    T: array-like object, optional
-        Time vector (argument is autocomputed if not given)
-
+    T: array-like or number, optional
+        Time vector, or simulation time duration if a number (time vector is 
+        autocomputed if not given)
+    
     X0: array-like or number, optional
         Initial condition (default = 0)
 
@@ -182,9 +185,10 @@ def initial(sys, T=None, X0=0., input=None, output=None, return_x=False):
     sys: StateSpace, or TransferFunction
         LTI system to simulate
 
-    T: array-like object, optional
-        Time vector (argument is autocomputed if not given)
-
+    T: array-like or number, optional
+        Time vector, or simulation time duration if a number (time vector is 
+        autocomputed if not given)
+    
     X0: array-like object or number, optional
         Initial condition (default = 0)
 
@@ -245,9 +249,8 @@ def lsim(sys, U=0., T=None, X0=0.):
         If `U` is ``None`` or ``0``, a special algorithm is used. This special
         algorithm is faster than the general algorithm, which is used otherwise.
 
-    T: array-like
-        Time steps at which the input is defined, numbers must be (strictly
-        monotonic) increasing.
+    T: array-like, optional for discrete LTI `sys`
+        Time steps at which the input is defined; values must be evenly spaced.
 
     X0: array-like or number, optional
         Initial condition (default = 0).
diff --git a/control/rlocus.py b/control/rlocus.py
@@ -50,7 +50,7 @@
 import numpy as np
 import matplotlib
 import matplotlib.pyplot as plt
-from scipy import array, poly1d, row_stack, zeros_like, real, imag
+from numpy import array, poly1d, row_stack, zeros_like, real, imag
 import scipy.signal             # signal processing toolbox
 import pylab                    # plotting routines
 from .xferfcn import _convert_to_transfer_function
diff --git a/control/sisotool.py b/control/sisotool.py
@@ -136,7 +136,7 @@ def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect=None):
     # Generate the step response and plot it
     sys_closed = (K*sys).feedback(1)
     if tvect is None:
-        tvect, yout = step_response(sys_closed)
+        tvect, yout = step_response(sys_closed, T_num=100)
     else:
         tvect, yout = step_response(sys_closed,tvect)
     ax_step.plot(tvect, yout)
diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py
@@ -33,10 +33,10 @@ def test_sisotool(self):
 
         # Check the step response before moving the point
         step_response_original = np.array(
-            [0., 0.02233651, 0.13118374, 0.33078542, 0.5907113, 0.87041549,
-             1.13038536, 1.33851053, 1.47374666, 1.52757114])
-        assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],
-                                  step_response_original, 4)
+            [0., 0.0217, 0.1281, 0.3237, 0.5797, 0.8566, 1.116,
+             1.3261, 1.4659, 1.526])
+        assert_array_almost_equal(
+            ax_step.lines[0].get_data()[1][:10], step_response_original, 4)
 
         bode_plot_params = {
             'omega': None,
@@ -78,10 +78,10 @@ def test_sisotool(self):
 
         # Check if the step response has changed
         step_response_moved = np.array(
-            [0., 0.02458187, 0.16529784, 0.46602716, 0.91012035, 1.43364313,
-             1.93996334, 2.3190105, 2.47041552, 2.32724853])
-        assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],
-                                  step_response_moved, 4)
+            [0., 0.0239, 0.161 , 0.4547, 0.8903, 1.407,
+             1.9121, 2.2989, 2.4686, 2.353])
+        assert_array_almost_equal(
+            ax_step.lines[0].get_data()[1][:10], step_response_moved, 4)
 
 
 if __name__ == "__main__":
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
@@ -11,6 +11,7 @@
 import unittest
 import numpy as np
 from control.timeresp import *
+from control.timeresp import _ideal_tfinal_and_dt, _default_time_vector
 from control.statesp import *
 from control.xferfcn import TransferFunction, _convert_to_transfer_function
 from control.dtime import c2d
@@ -94,6 +95,7 @@ def test_step_response(self):
         np.testing.assert_array_equal(Tc.shape, Td.shape)
         np.testing.assert_array_equal(youtc.shape, youtd.shape)
 
+
     # Recreate issue #374 ("Bug in step_response()")
     def test_step_nostates(self):
         # Continuous time, constant system
@@ -346,10 +348,75 @@ def test_step_robustness(self):
         sys2 = TransferFunction(num, den2)
 
         # Compute step response from input 1 to output 1, 2
-        t1, y1 = step_response(sys1, input=0)
-        t2, y2 = step_response(sys2, input=0)
+        t1, y1 = step_response(sys1, input=0, T_num=100)
+        t2, y2 = step_response(sys2, input=0, T_num=100)
         np.testing.assert_array_almost_equal(y1, y2)
 
+    def test_auto_generated_time_vector(self):
+        # confirm a TF with a pole at p simulates for 7.0/p seconds
+        p = 0.5
+        np.testing.assert_array_almost_equal(
+            _ideal_tfinal_and_dt(TransferFunction(1, [1, .5]))[0], 
+            (7/p))
+        np.testing.assert_array_almost_equal(
+            _ideal_tfinal_and_dt(TransferFunction(1, [1, .5]).sample(.1))[0], 
+            (7/p))
+        # confirm a TF with poles at 0 and p simulates for 7.0/p seconds
+        np.testing.assert_array_almost_equal(
+            _ideal_tfinal_and_dt(TransferFunction(1, [1, .5, 0]))[0], 
+            (7/p))
+        # confirm a TF with a natural frequency of wn rad/s gets a 
+        # dt of 1/(7.0*wn)
+        wn = 10
+        np.testing.assert_array_almost_equal(
+            _ideal_tfinal_and_dt(TransferFunction(1, [1, 0, wn**2]))[1], 
+            1/(7.0*wn))
+        zeta = .1
+        np.testing.assert_array_almost_equal(
+            _ideal_tfinal_and_dt(TransferFunction(1, [1, 2*zeta*wn, wn**2]))[1], 
+            1/(7.0*wn))
+        # but a smapled one keeps its dt
+        np.testing.assert_array_almost_equal(
+            _ideal_tfinal_and_dt(TransferFunction(1, [1, 2*zeta*wn, wn**2]).sample(.1))[1], 
+            .1)
+        np.testing.assert_array_almost_equal(
+            np.diff(initial_response(TransferFunction(1, [1, 2*zeta*wn, wn**2]).sample(.1))[0][0:2]), 
+            .1)
+        np.testing.assert_array_almost_equal(
+            _ideal_tfinal_and_dt(TransferFunction(1, [1, 2*zeta*wn, wn**2]))[1], 
+            1/(7.0*wn))
+        # TF with fast oscillations simulates only 5000 time steps even with long tfinal
+        self.assertEqual(5000, 
+            len(_default_time_vector(TransferFunction(1, [1, 0, wn**2]),tfinal=100)))
+        # and simulates for 7.0/dt time steps 
+        self.assertEqual(
+            len(_default_time_vector(TransferFunction(1, [1, 0, wn**2]))), 
+            int(7.0/(1/(7.0*wn))))
+        
+        sys = TransferFunction(1, [1, .5, 0])
+        sysdt = TransferFunction(1, [1, .5, 0], .1)
+        # test impose number of time steps
+        self.assertEqual(10, len(step_response(sys, T_num=10)[0]))
+        self.assertEqual(10, len(step_response(sysdt, T_num=10)[0]))
+        # test impose final time
+        np.testing.assert_array_almost_equal(
+            100, 
+            step_response(sys, 100)[0][-1], 
+            decimal=.5)
+        np.testing.assert_array_almost_equal(
+            100, 
+            step_response(sysdt, 100)[0][-1], 
+            decimal=.5)
+        np.testing.assert_array_almost_equal(
+            100, 
+            impulse_response(sys, 100)[0][-1], 
+            decimal=.5)
+        np.testing.assert_array_almost_equal(
+            100, 
+            initial_response(sys, 100)[0][-1], 
+            decimal=.5)
+
+
     def test_time_vector(self):
         "Unit test: https://github.com/python-control/python-control/issues/239"
         # Discrete time simulations with specified time vectors
diff --git a/control/timeresp.py b/control/timeresp.py
