fix computation of default response time for constant systems (#383)

murrayrm · web-flow · commit b9dad5e8bbfe · 2020-05-16T19:08:59.000-07:00
* fix computation of default response time for constant systems

* add single function to compute response times
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
@@ -94,9 +94,21 @@ 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
+        sys = TransferFunction([1], [1])
+        t, y = step_response(sys)
+        np.testing.assert_array_equal(y, np.ones(len(t)))
+
+        # Discrete time, constant system
+        sys = TransferFunction([1], [1], 1)
+        t, y = step_response(sys)
+        np.testing.assert_array_equal(y, np.ones(len(t)))
+
     def test_step_info(self):
         # From matlab docs:
-        sys = TransferFunction([1,5,5],[1,1.65,5,6.5,2])
+        sys = TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2])
         Strue = {
             'RiseTime': 3.8456,
             'SettlingTime': 27.9762,
diff --git a/control/timeresp.py b/control/timeresp.py
@@ -512,13 +512,7 @@ def step_response(sys, T=None, X0=0., input=None, output=None,
     """
     sys = _get_ss_simo(sys, input, output)
     if T is None:
-        if isctime(sys):
-            T = _default_response_times(sys.A, 100)
-        else:
-            # For discrete time, use integers
-            tvec = _default_response_times(sys.A, 100)
-            T = range(int(np.ceil(max(tvec))))
-
+        T = _get_response_times(sys, N=100)
     U = np.ones_like(T)
 
     T, yout, xout = forced_response(sys, T, U, X0, transpose=transpose,
@@ -573,12 +567,7 @@ def step_info(sys, T=None, SettlingTimeThreshold=0.02,
     '''
     sys = _get_ss_simo(sys)
     if T is None:
-        if isctime(sys):
-            T = _default_response_times(sys.A, 1000)
-        else:
-            # For discrete time, use integers
-            tvec = _default_response_times(sys.A, 1000)
-            T = range(int(np.ceil(max(tvec))))
+        T = _get_response_times(sys, N=1000)
 
     T, yout = step_response(sys, T)
 
@@ -697,12 +686,8 @@ def initial_response(sys, T=None, X0=0., input=0, output=None,
     # Create time and input vectors; checking is done in forced_response(...)
     # The initial vector X0 is created in forced_response(...) if necessary
     if T is None:
-        if isctime(sys):
-            T = _default_response_times(sys.A, 1000)
-        else:
-            # For discrete time, use integers
-            tvec = _default_response_times(sys.A, 1000)
-            T = range(int(np.ceil(max(tvec))))
+        # TODO: default step size inconsistent with step/impulse_response()
+        T = _get_response_times(sys, N=1000)
     U = np.zeros_like(T)
 
     T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,
@@ -801,13 +786,7 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None,
 
     # Compute T and U, no checks necessary, they will be checked in lsim
     if T is None:
-        if isctime(sys):
-            T = _default_response_times(sys.A, 100)
-        else:
-            # For discrete time, use integers
-            tvec = _default_response_times(sys.A, 100)
-            T = range(int(np.ceil(max(tvec))))
-
+        T = _get_response_times(sys, N=100)
     U = np.zeros_like(T)
 
     # Compute new X0 that contains the impulse
@@ -828,3 +807,22 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None,
         return T, yout, _xout
 
     return T, yout
+
+
+# Utility function to get response times
+def _get_response_times(sys, N=100):
+    if isctime(sys):
+        if sys.A.shape == (0, 0):
+            # No dynamics; use the unit time interval
+            T = np.linspace(0, 1, N, endpoint=False)
+        else:
+            T = _default_response_times(sys.A, N)
+    else:
+        # For discrete time, use integers
+        if sys.A.shape == (0, 0):
+            # No dynamics; use N time steps
+            T = range(N)
+        else:
+            tvec = _default_response_times(sys.A, N)
+            T = range(int(np.ceil(max(tvec))))
+    return T
