From b39184be78d3ae32481250b2b4953e66bf6e95b8 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Thu, 18 Jun 2020 15:57:34 -0700
Subject: [PATCH 01/10] fixed default response time for time response of
 discrete-time functions step, impulse, and initial

---
 control/timeresp.py | 48 +++++++++++++++++++++++----------------------
 1 file changed, 25 insertions(+), 23 deletions(-)

diff --git a/control/timeresp.py b/control/timeresp.py
index 4c0fbd940..41f2f443b 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -69,7 +69,7 @@
 from scipy.signal.ltisys import _default_response_times
 import warnings
 from .lti import LTI     # base class of StateSpace, TransferFunction
-from .statesp import _convertToStateSpace, _mimo2simo, _mimo2siso
+from .statesp import _convertToStateSpace, _mimo2simo, _mimo2siso, ssdata
 from .lti import isdtime, isctime
 
 __all__ = ['forced_response', 'step_response', 'step_info', 'initial_response',
@@ -512,7 +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:
-        T = _get_response_times(sys, N=100)
+        T = _get_response_times(sys)
     U = np.ones_like(T)
 
     T, yout, xout = forced_response(sys, T, U, X0, transpose=transpose,
@@ -567,7 +567,7 @@ def step_info(sys, T=None, SettlingTimeThreshold=0.02,
     '''
     sys = _get_ss_simo(sys)
     if T is None:
-        T = _get_response_times(sys, N=1000)
+        T = _get_response_times(sys)
 
     T, yout = step_response(sys, T)
 
@@ -686,8 +686,7 @@ 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:
-        # TODO: default step size inconsistent with step/impulse_response()
-        T = _get_response_times(sys, N=1000)
+        T = _get_response_times(sys)
     U = np.zeros_like(T)
 
     T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,
@@ -786,7 +785,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:
-        T = _get_response_times(sys, N=100)
+        T = _get_response_times(sys)
     U = np.zeros_like(T)
 
     # Compute new X0 that contains the impulse
@@ -808,21 +807,24 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None,
 
     return T, yout
 
-
+# utility function to find time period for time responses if not provided
+def _get_default_tfinal(sys):
+    A = ssdata(sys)[0]
+    if A.shape == (0,0): 
+        return 1.0 # no dynamics, use unit time interval
+    if isdtime(sys, strict=True):
+        A = 1.0/sys.dt * (A - np.eye(A.shape[0])) # zoh approximation
+    poles = sp.linalg.eigvals(A)
+    tfinal = 7.0 / min(abs(poles.real))
+    return tfinal if np.isfinite(tfinal) else 1.0
+    
 # 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
+def _get_response_times(sys, N=1000, tfinal=None):
+    """Returns a time vector suitable for observing the response of the 
+    slowest poles of a system. if system is discrete-time, N is ignored"""
+    if tfinal == None: 
+        tfinal = _get_default_tfinal(sys)
+    if isdtime(sys, strict=True):
+        return np.arange(0, tfinal, sys.dt)
+    else: 
+        return np.linspace(0, tfinal, N, endpoint=False)
\ No newline at end of file

From 32b610bb88503c6714ea69e46ae6d44ecdf1846a Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Fri, 19 Jun 2020 11:01:24 -0700
Subject: [PATCH 02/10] to pass tests, added convenient ability to specify
 simulation time and number of steps rather than complete time vector in
 timeresponse functions

---
 control/matlab/timeresp.py |  31 +++++-----
 control/sisotool.py        |   2 +-
 control/timeresp.py        | 117 ++++++++++++++++++++-----------------
 3 files changed, 82 insertions(+), 68 deletions(-)

diff --git a/control/matlab/timeresp.py b/control/matlab/timeresp.py
index 647210a9c..b9d4004ca 100644
--- 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/sisotool.py b/control/sisotool.py
index e700875ca..c2db4b5ab 100644
--- 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/timeresp.py b/control/timeresp.py
index 41f2f443b..97afc87eb 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -440,7 +440,7 @@ def _get_ss_simo(sys, input=None, output=None):
         return _mimo2siso(sys_ss, input, output, warn_conversion=warn)
 
 
-def step_response(sys, T=None, X0=0., input=None, output=None,
+def step_response(sys, T=None, X0=0., input=None, output=None, T_num=1000,
                   transpose=False, return_x=False, squeeze=True):
     # pylint: disable=W0622
     """Step response of a linear system
@@ -458,8 +458,9 @@ def step_response(sys, T=None, X0=0., input=None, output=None,
     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)
@@ -472,6 +473,10 @@ def step_response(sys, T=None, X0=0., input=None, output=None,
     output: int
         Index of the output that will be used in this simulation. Set to None
         to not trim outputs
+    
+    T_num: number, optional (default=1000)
+        Number of time steps to use in simulation if T is not provided as an 
+        array; ignored if sys is discrete-time.  
 
     transpose: bool
         If True, transpose all input and output arrays (for backward
@@ -511,8 +516,8 @@ def step_response(sys, T=None, X0=0., input=None, output=None,
 
     """
     sys = _get_ss_simo(sys, input, output)
-    if T is None:
-        T = _get_response_times(sys)
+    if T is None or np.asarray(T).size == 1:
+        T = _get_response_times(sys, N=T_num, tfinal=T)
     U = np.ones_like(T)
 
     T, yout, xout = forced_response(sys, T, U, X0, transpose=transpose,
@@ -524,7 +529,7 @@ def step_response(sys, T=None, X0=0., input=None, output=None,
     return T, yout
 
 
-def step_info(sys, T=None, SettlingTimeThreshold=0.02,
+def step_info(sys, T=None, T_num=1000, SettlingTimeThreshold=0.02,
               RiseTimeLimits=(0.1, 0.9)):
     '''
     Step response characteristics (Rise time, Settling Time, Peak and others).
@@ -534,8 +539,13 @@ def step_info(sys, T=None, SettlingTimeThreshold=0.02,
     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)
+
+    T_num: number, optional (default=1000)
+        Number of time steps to use in simulation if T is not provided as an 
+        array; ignored if sys is discrete-time.  
 
     SettlingTimeThreshold: float value, optional
         Defines the error to compute settling time (default = 0.02)
@@ -566,9 +576,9 @@ def step_info(sys, T=None, SettlingTimeThreshold=0.02,
     >>> info = step_info(sys, T)
     '''
     sys = _get_ss_simo(sys)
-    if T is None:
-        T = _get_response_times(sys)
-
+    if T is None or np.asarray(T).size == 1:
+        T = _get_response_times(sys, N=T_num, tfinal=T)
+    
     T, yout = step_response(sys, T)
 
     # Steady state value
@@ -588,33 +598,21 @@ def step_info(sys, T=None, SettlingTimeThreshold=0.02,
             SettlingTime = T[i + 1]
             break
 
-    # Peak
     PeakIndex = np.abs(yout).argmax()
-    PeakValue = yout[PeakIndex]
-    PeakTime = T[PeakIndex]
-    SettlingMax = (yout).max()
-    SettlingMin = (yout[tr_upper_index:]).min()
-    # I'm really not very confident about UnderShoot:
-    UnderShoot = yout.min()
-    OverShoot = 100. * (yout.max() - InfValue) / (InfValue - yout[0])
-
-    # Return as a dictionary
-    S = {
+    return {
         'RiseTime': RiseTime,
         'SettlingTime': SettlingTime,
-        'SettlingMin': SettlingMin,
-        'SettlingMax': SettlingMax,
-        'Overshoot': OverShoot,
-        'Undershoot': UnderShoot,
-        'Peak': PeakValue,
-        'PeakTime': PeakTime,
+        'SettlingMin': yout[tr_upper_index:].min(),
+        'SettlingMax': yout.max(),
+        'Overshoot': 100. * (yout.max() - InfValue) / (InfValue - yout[0]),
+        'Undershoot': yout.min(), # not very confident about this
+        'Peak': yout[PeakIndex],
+        'PeakTime':  T[PeakIndex],
         'SteadyStateValue': InfValue
-    }
+        }
 
-    return S
 
-
-def initial_response(sys, T=None, X0=0., input=0, output=None,
+def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
                      transpose=False, return_x=False, squeeze=True):
     # pylint: disable=W0622
     """Initial condition response of a linear system
@@ -631,10 +629,11 @@ def initial_response(sys, T=None, X0=0., input=0, output=None,
     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
+    X0: array-like or number, optional
         Initial condition (default = 0)
 
         Numbers are converted to constant arrays with the correct shape.
@@ -646,6 +645,10 @@ def initial_response(sys, T=None, X0=0., input=0, output=None,
     output: int
         Index of the output that will be used in this simulation. Set to None
         to not trim outputs
+    
+    T_num: number, optional (default=1000)
+        Number of time steps to use in simulation if T is not provided as an 
+        array; ignored if sys is discrete-time.  
 
     transpose: bool
         If True, transpose all input and output arrays (for backward
@@ -685,9 +688,9 @@ 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:
-        T = _get_response_times(sys)
-    U = np.zeros_like(T)
+    if T is None or np.asarray(T).size == 1:
+        T = _get_response_times(sys, N=T_num, tfinal=T)
+    U = 0 #np.zeros_like(T)
 
     T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,
                                      squeeze=squeeze)
@@ -698,7 +701,7 @@ def initial_response(sys, T=None, X0=0., input=0, output=None,
     return T, yout
 
 
-def impulse_response(sys, T=None, X0=0., input=0, output=None,
+def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
                      transpose=False, return_x=False, squeeze=True):
     # pylint: disable=W0622
     """Impulse response of a linear system
@@ -716,10 +719,11 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None,
     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 object or number, optional
+    X0: array-like or number, optional
         Initial condition (default = 0)
 
         Numbers are converted to constant arrays with the correct shape.
@@ -731,6 +735,10 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None,
         Index of the output that will be used in this simulation. Set to None
         to not trim outputs
 
+    T_num: number, optional (default=1000)
+        Number of time steps to use in simulation if T is not provided as an 
+        array; ignored if sys is discrete-time.  
+
     transpose: bool
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and scipy.signal.lsim)
@@ -783,9 +791,9 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None,
     X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)],
                               'Parameter ``X0``: \n', squeeze=True)
 
-    # Compute T and U, no checks necessary, they will be checked in lsim
-    if T is None:
-        T = _get_response_times(sys)
+    # Compute T and U, no checks necessary, will be checked in forced_response
+    if T is None or np.asarray(T).size == 1:
+        T = _get_response_times(sys, N=T_num, tfinal=T)
     U = np.zeros_like(T)
 
     # Compute new X0 that contains the impulse
@@ -807,22 +815,25 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None,
 
     return T, yout
 
-# utility function to find time period for time responses if not provided
+# utility function to find time period for time responses if not provided. 
 def _get_default_tfinal(sys):
     A = ssdata(sys)[0]
     if A.shape == (0,0): 
-        return 1.0 # no dynamics, use unit time interval
-    if isdtime(sys, strict=True):
-        A = 1.0/sys.dt * (A - np.eye(A.shape[0])) # zoh approximation
-    poles = sp.linalg.eigvals(A)
-    tfinal = 7.0 / min(abs(poles.real))
-    return tfinal if np.isfinite(tfinal) else 1.0
+        slowest_time_constant = 1.0 # no dynamics
+    else:
+        if isdtime(sys, strict=True):
+            A = 1.0/sys.dt * (A - np.eye(A.shape[0])) # zoh approximation
+        poles = sp.linalg.eigvals(A)
+        slowest_time_constant = 1.0 / min(abs(poles.real))
+        if np.isinf(slowest_time_constant): 
+            slowest_time_constant = 1.0
+    return 7.0 * slowest_time_constant
     
 # Utility function to get response times
 def _get_response_times(sys, N=1000, tfinal=None):
     """Returns a time vector suitable for observing the response of the 
     slowest poles of a system. if system is discrete-time, N is ignored"""
-    if tfinal == None: 
+    if tfinal is None: 
         tfinal = _get_default_tfinal(sys)
     if isdtime(sys, strict=True):
         return np.arange(0, tfinal, sys.dt)

From 9a8da65352bd7b8e96d9d4f21c23259272f6a799 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Fri, 19 Jun 2020 11:51:05 -0700
Subject: [PATCH 03/10] eliminated deprecation warnings by importing certain
 functions from numpy instead of scipy

---
 control/freqplot.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/control/freqplot.py b/control/freqplot.py
index a1772fea7..7b296c111 100644
--- 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
 
 

From db9f1977247e8347a3e61eff06df49ac6ec63755 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Fri, 19 Jun 2020 11:51:41 -0700
Subject: [PATCH 04/10] eliminated deprecation warnings by importing certain
 functions from numpy instead of scipy

---
 control/rlocus.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/control/rlocus.py b/control/rlocus.py
index 955c5c56d..56e0c55d1 100644
--- 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

From 6d629106be8a3109801b4ef3cd6378a911856ab7 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Fri, 19 Jun 2020 11:52:31 -0700
Subject: [PATCH 05/10] small fix to pass unit tests

---
 control/timeresp.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/control/timeresp.py b/control/timeresp.py
index 97afc87eb..7f606f4e2 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -690,7 +690,7 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
     # The initial vector X0 is created in forced_response(...) if necessary
     if T is None or np.asarray(T).size == 1:
         T = _get_response_times(sys, N=T_num, tfinal=T)
-    U = 0 #np.zeros_like(T)
+    U = np.zeros_like(T)
 
     T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,
                                      squeeze=squeeze)
@@ -777,7 +777,7 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
     """
     sys = _get_ss_simo(sys, input, output)
 
-    # System has direct feedthrough, can't simulate impulse response
+    # if system has direct feedthrough, can't simulate impulse response
     # numerically
     if np.any(sys.D != 0) and isctime(sys):
         warnings.warn("System has direct feedthrough: ``D != 0``. The "
@@ -786,7 +786,7 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
                       "Results may be meaningless!")
 
     # create X0 if not given, test if X0 has correct shape.
-    # Must be done here because it is used for computations here.
+    # Must be done here because it is used for computations below.
     n_states = sys.A.shape[0]
     X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)],
                               'Parameter ``X0``: \n', squeeze=True)

From 9f021192d670a789041ab4fa40941d84be787da5 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Fri, 19 Jun 2020 15:20:55 -0700
Subject: [PATCH 06/10] adjusted sisotool test so tests pass with new default
 step response time window

---
 control/tests/sisotool_test.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py
index f2cdf9106..34ea32414 100644
--- a/control/tests/sisotool_test.py
+++ b/control/tests/sisotool_test.py
@@ -27,7 +27,7 @@ def test_sisotool(self):
         assert_array_almost_equal(ax_rlocus.lines[2].get_data(),initial_point_2)
 
         # 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])
+        step_response_original = np.array([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)
 
         bode_plot_params = {
@@ -59,7 +59,7 @@ def test_sisotool(self):
         assert_array_almost_equal(ax_mag.lines[0].get_data()[1][10:20],bode_mag_moved)
 
         # 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] ])
+        step_response_moved = np.array([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)
 
 

From 1bc60baf5ed9deb8fd8d8bd794eb9fff3be4c91c Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Sun, 21 Jun 2020 02:50:48 -0700
Subject: [PATCH 07/10] added functionality to automatically choose dt in
 timeresp.py based on system poles. and unit tests.

---
 control/tests/timeresp_test.py |  75 ++++++++++++++-
 control/timeresp.py            | 168 ++++++++++++++++++++++++++-------
 2 files changed, 209 insertions(+), 34 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index b208e70d2..90eae5b50 100644
--- 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,79 @@ 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)
+        
+        
+        # slow system max 
+        # to add: impose N, impose tfinal
+
+
     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
index 7f606f4e2..03cb3df01 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -60,13 +60,17 @@
 
 Initial Author: Eike Welk
 Date: 12 May 2011
+
+Modified: Sawyer B. Fuller (minster@uw.edu) to add discrete-time 
+capability and better automatic time vector creation
+Date: June 2020
+
 $Id$
 """
 
 # Libraries that we make use of
 import scipy as sp              # SciPy library (used all over)
 import numpy as np              # NumPy library
-from scipy.signal.ltisys import _default_response_times
 import warnings
 from .lti import LTI     # base class of StateSpace, TransferFunction
 from .statesp import _convertToStateSpace, _mimo2simo, _mimo2siso, ssdata
@@ -440,7 +444,7 @@ def _get_ss_simo(sys, input=None, output=None):
         return _mimo2siso(sys_ss, input, output, warn_conversion=warn)
 
 
-def step_response(sys, T=None, X0=0., input=None, output=None, T_num=1000,
+def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
                   transpose=False, return_x=False, squeeze=True):
     # pylint: disable=W0622
     """Step response of a linear system
@@ -474,9 +478,9 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=1000,
         Index of the output that will be used in this simulation. Set to None
         to not trim outputs
     
-    T_num: number, optional (default=1000)
+    T_num: number, optional
         Number of time steps to use in simulation if T is not provided as an 
-        array; ignored if sys is discrete-time.  
+        array (autocomputed if not given); ignored if sys is discrete-time.  
 
     transpose: bool
         If True, transpose all input and output arrays (for backward
@@ -517,7 +521,7 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=1000,
     """
     sys = _get_ss_simo(sys, input, output)
     if T is None or np.asarray(T).size == 1:
-        T = _get_response_times(sys, N=T_num, tfinal=T)
+        T = _default_time_vector(sys, N=T_num, tfinal=T)
     U = np.ones_like(T)
 
     T, yout, xout = forced_response(sys, T, U, X0, transpose=transpose,
@@ -529,7 +533,7 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=1000,
     return T, yout
 
 
-def step_info(sys, T=None, T_num=1000, SettlingTimeThreshold=0.02,
+def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
               RiseTimeLimits=(0.1, 0.9)):
     '''
     Step response characteristics (Rise time, Settling Time, Peak and others).
@@ -543,9 +547,9 @@ def step_info(sys, T=None, T_num=1000, SettlingTimeThreshold=0.02,
         Time vector, or simulation time duration if a number (time vector is 
         autocomputed if not given)
 
-    T_num: number, optional (default=1000)
+    T_num: number, optional
         Number of time steps to use in simulation if T is not provided as an 
-        array; ignored if sys is discrete-time.  
+        array (autocomputed if not given); ignored if sys is discrete-time.  
 
     SettlingTimeThreshold: float value, optional
         Defines the error to compute settling time (default = 0.02)
@@ -577,7 +581,7 @@ def step_info(sys, T=None, T_num=1000, SettlingTimeThreshold=0.02,
     '''
     sys = _get_ss_simo(sys)
     if T is None or np.asarray(T).size == 1:
-        T = _get_response_times(sys, N=T_num, tfinal=T)
+        T = _default_time_vector(sys, N=T_num, tfinal=T)
     
     T, yout = step_response(sys, T)
 
@@ -612,7 +616,7 @@ def step_info(sys, T=None, T_num=1000, SettlingTimeThreshold=0.02,
         }
 
 
-def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
+def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
                      transpose=False, return_x=False, squeeze=True):
     # pylint: disable=W0622
     """Initial condition response of a linear system
@@ -646,9 +650,9 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
         Index of the output that will be used in this simulation. Set to None
         to not trim outputs
     
-    T_num: number, optional (default=1000)
+    T_num: number, optional
         Number of time steps to use in simulation if T is not provided as an 
-        array; ignored if sys is discrete-time.  
+        array (autocomputed if not given); ignored if sys is discrete-time.  
 
     transpose: bool
         If True, transpose all input and output arrays (for backward
@@ -689,7 +693,7 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
     # 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 or np.asarray(T).size == 1:
-        T = _get_response_times(sys, N=T_num, tfinal=T)
+        T = _default_time_vector(sys, N=T_num, tfinal=T)
     U = np.zeros_like(T)
 
     T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,
@@ -701,7 +705,7 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
     return T, yout
 
 
-def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
+def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
                      transpose=False, return_x=False, squeeze=True):
     # pylint: disable=W0622
     """Impulse response of a linear system
@@ -735,9 +739,9 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
         Index of the output that will be used in this simulation. Set to None
         to not trim outputs
 
-    T_num: number, optional (default=1000)
+    T_num: number, optional
         Number of time steps to use in simulation if T is not provided as an 
-        array; ignored if sys is discrete-time.  
+        array (autocomputed if not given); ignored if sys is discrete-time.  
 
     transpose: bool
         If True, transpose all input and output arrays (for backward
@@ -793,7 +797,7 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
 
     # Compute T and U, no checks necessary, will be checked in forced_response
     if T is None or np.asarray(T).size == 1:
-        T = _get_response_times(sys, N=T_num, tfinal=T)
+        T = _default_time_vector(sys, N=T_num, tfinal=T)
     U = np.zeros_like(T)
 
     # Compute new X0 that contains the impulse
@@ -815,27 +819,127 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=1000,
 
     return T, yout
 
-# utility function to find time period for time responses if not provided. 
-def _get_default_tfinal(sys):
+# utility function to find time period and time increment using pole locations
+def _ideal_tfinal_and_dt(sys):
+    constant = 7.0
+    tolerance = 1e-10
     A = ssdata(sys)[0]
     if A.shape == (0,0): 
-        slowest_time_constant = 1.0 # no dynamics
+        # no dynamics
+        tfinal = constant * 1.0 
+        dt = sys.dt if isdtime(sys, strict=True) else 1.0
     else:
-        if isdtime(sys, strict=True):
-            A = 1.0/sys.dt * (A - np.eye(A.shape[0])) # zoh approximation
         poles = sp.linalg.eigvals(A)
-        slowest_time_constant = 1.0 / min(abs(poles.real))
-        if np.isinf(slowest_time_constant): 
-            slowest_time_constant = 1.0
-    return 7.0 * slowest_time_constant
-    
-# Utility function to get response times
-def _get_response_times(sys, N=1000, tfinal=None):
+        if isdtime(sys, strict=True):
+            poles = np.log(poles)/sys.dt # z-poles to s-plane using s=(lnz)/dt
+        
+        # calculate ideal dt
+        if isdtime(sys, strict=True):
+            dt = sys.dt
+        else:
+            fastest_natural_frequency = max(abs(poles))
+            dt = 1/constant / fastest_natural_frequency
+        
+        # calculate ideal tfinal
+        poles = poles[abs(poles.real) > tolerance] # ignore poles near im axis
+        if poles.size == 0:        
+            slowest_decay_rate = 1.0
+        else:
+            slowest_decay_rate = min(abs(poles.real))
+        tfinal = constant / slowest_decay_rate  
+
+    return tfinal, dt
+
+# test below: ct with pole at the origin is 7 seconds, ct with pole at .5 is 14 s long, 
+def _default_time_vector(sys, N=None, tfinal=None):
     """Returns a time vector suitable for observing the response of the 
-    slowest poles of a system. if system is discrete-time, N is ignored"""
+    both the slowest poles and fastest resonant modes. if system is 
+    discrete-time, N is ignored """
+    
+    N_max = 5000
+    N_min_ct = 100
+    N_min_dt = 7 # more common to see just a few samples in discrete-time
+    
+    ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys)
+    
+    if isdtime(sys, strict=True):
+        if tfinal is None: 
+            # for discrete time, change from ideal_tfinal if N too large/small
+            N = int(np.clip(ideal_tfinal/sys.dt, N_min_dt, N_max))# [N_min, N_max]
+            tfinal = sys.dt * N
+        else: 
+            N = int(tfinal/sys.dt)
+    else: 
+        if tfinal is None: 
+            # for continuous time, simulate to ideal_tfinal but limit N
+            tfinal = ideal_tfinal
+        if N is None:     
+            N = int(np.clip(tfinal/ideal_dt, N_min_ct, N_max)) # N<-[N_min, N_max]
+            
+    return np.linspace(0, tfinal, N, endpoint=False)
+
+def junk(sys, N, tfinal):
+    N_max = 5000
+    N_min = 10
+    if tfinal is None: 
+        ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys)
+        N = int(np.clip(ideal_tfinal/ideal_dt, N_min, N_max))# N->[N_min, N_max]
+        if isdtime(sys, strict=True):
+            # for discrete time, change tfinal if N gets too large/small
+            tfinal = sys.dt * N
+        else: 
+            # for continuous time, simulate to tfinal but change dt if N large
+            tfinal = ideal_tfinal
+    else: 
+        if isdtime(sys, strict=True):
+            N = int(tfinal/sys.dt)
+        else:
+            _, ideal_dt = _ideal_tfinal_and_dt(sys)
+            N = int(np.clip(tfinal/sys.dt, N_min, N_max)) # N->[N_min, N_max]
+            
+         
+    if isdtime(sys, strict=True):
+        if tfinal is None: 
+            tfinal, dt = _ideal_tfinal_and_dt(sys)
+            # sanity bounds on tfinal
+            N = int(np.clip(tfinal/sys.dt, N_min, N_max)) # [N_min, N_max]
+            tfinal = sys.dt * N
+        else:
+            N = int(tfinal/sys.dt)
+    else: 
+        if tfinal is None: 
+            tfinal, dt = _ideal_tfinal_and_dt(sys)
+        N = int(np.clip(tfinal/dt, N_min, N_max)) # [N_min, N_max]
+        dt = tfinal/N
+    return np.linspace(0, tfinal, N, endpoint=False)
+    
+    if 1: 
+        # run to N steps, or to ideal tfinal if N not given 
+        if N is None:
+            # run to tfinal
+            if tfinal is None: 
+                tfinal, dt = _ideal_tfinal_and_dt(sys)
+            N = int(np.clip(tfinal/sys.dt, N_min, N_max)) # [N_min, N_max]
+            tfinal = sys.dt * N
+        else: 
+            tfinal = sys.dt * N
+    else: 
+        # continuous-time 
+        ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys)
+ 
+
+    # tfinal 
     if tfinal is None: 
-        tfinal = _get_default_tfinal(sys)
+        ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys)
+        
+    if N is None: 
+        ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys)
+        if tfinal is None: 
+            tfinal = ideal_tfinal 
+        N = int(np.clip(tfinal/dt, N_min, N_max)) # limit to [N_min, N_max]
     if isdtime(sys, strict=True):
-        return np.arange(0, tfinal, sys.dt)
+        tfinal = sys.dt * N
     else: 
+        dt = tfinal/N
+        return np.arange(0, tfinal, dt)
         return np.linspace(0, tfinal, N, endpoint=False)
\ No newline at end of file

From 512f822cef56f57f2e268a117d24511dfcdab948 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Sun, 21 Jun 2020 14:41:00 -0700
Subject: [PATCH 08/10] removed some leftover code and comments

---
 control/tests/timeresp_test.py |  4 --
 control/timeresp.py            | 68 +---------------------------------
 2 files changed, 1 insertion(+), 71 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 90eae5b50..5549b2a88 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -415,10 +415,6 @@ def test_auto_generated_time_vector(self):
             100, 
             initial_response(sys, 100)[0][-1], 
             decimal=.5)
-        
-        
-        # slow system max 
-        # to add: impose N, impose tfinal
 
 
     def test_time_vector(self):
diff --git a/control/timeresp.py b/control/timeresp.py
index 03cb3df01..eb1b9372a 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -876,70 +876,4 @@ def _default_time_vector(sys, N=None, tfinal=None):
         if N is None:     
             N = int(np.clip(tfinal/ideal_dt, N_min_ct, N_max)) # N<-[N_min, N_max]
             
-    return np.linspace(0, tfinal, N, endpoint=False)
-
-def junk(sys, N, tfinal):
-    N_max = 5000
-    N_min = 10
-    if tfinal is None: 
-        ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys)
-        N = int(np.clip(ideal_tfinal/ideal_dt, N_min, N_max))# N->[N_min, N_max]
-        if isdtime(sys, strict=True):
-            # for discrete time, change tfinal if N gets too large/small
-            tfinal = sys.dt * N
-        else: 
-            # for continuous time, simulate to tfinal but change dt if N large
-            tfinal = ideal_tfinal
-    else: 
-        if isdtime(sys, strict=True):
-            N = int(tfinal/sys.dt)
-        else:
-            _, ideal_dt = _ideal_tfinal_and_dt(sys)
-            N = int(np.clip(tfinal/sys.dt, N_min, N_max)) # N->[N_min, N_max]
-            
-         
-    if isdtime(sys, strict=True):
-        if tfinal is None: 
-            tfinal, dt = _ideal_tfinal_and_dt(sys)
-            # sanity bounds on tfinal
-            N = int(np.clip(tfinal/sys.dt, N_min, N_max)) # [N_min, N_max]
-            tfinal = sys.dt * N
-        else:
-            N = int(tfinal/sys.dt)
-    else: 
-        if tfinal is None: 
-            tfinal, dt = _ideal_tfinal_and_dt(sys)
-        N = int(np.clip(tfinal/dt, N_min, N_max)) # [N_min, N_max]
-        dt = tfinal/N
-    return np.linspace(0, tfinal, N, endpoint=False)
-    
-    if 1: 
-        # run to N steps, or to ideal tfinal if N not given 
-        if N is None:
-            # run to tfinal
-            if tfinal is None: 
-                tfinal, dt = _ideal_tfinal_and_dt(sys)
-            N = int(np.clip(tfinal/sys.dt, N_min, N_max)) # [N_min, N_max]
-            tfinal = sys.dt * N
-        else: 
-            tfinal = sys.dt * N
-    else: 
-        # continuous-time 
-        ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys)
- 
-
-    # tfinal 
-    if tfinal is None: 
-        ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys)
-        
-    if N is None: 
-        ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys)
-        if tfinal is None: 
-            tfinal = ideal_tfinal 
-        N = int(np.clip(tfinal/dt, N_min, N_max)) # limit to [N_min, N_max]
-    if isdtime(sys, strict=True):
-        tfinal = sys.dt * N
-    else: 
-        dt = tfinal/N
-        return np.arange(0, tfinal, dt)
-        return np.linspace(0, tfinal, N, endpoint=False)
\ No newline at end of file
+    return np.linspace(0, tfinal, N, endpoint=False)
\ No newline at end of file

From 79951a10bfe7df0c9f60589b2b33bf1c911c5fdc Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Sat, 27 Jun 2020 20:49:25 -0700
Subject: [PATCH 09/10] explanation in docstrings for how time vector T is
 auto-computed in time response functions

---
 control/timeresp.py | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/control/timeresp.py b/control/timeresp.py
index eb1b9372a..8670c180d 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -463,8 +463,14 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
         LTI system to simulate
 
     T: array-like or number, optional
-        Time vector, or simulation time duration if a number (time vector is 
-        autocomputed if not given)
+        Time vector, or simulation time duration if a number. If T is not 
+        provided, an attempt is made to create it automatically from the 
+        dynamics of sys. If sys is continuous-time, the time increment dt 
+        is chosen small enough to show the fastest mode, and the simulation 
+        time period tfinal long enough to show the slowest mode, excluding 
+        poles at the origin. If this results in too many time steps (>5000),
+        dt is reduced. If sys is discrete-time, only tfinal is computed, and 
+        tfinal is reduced if it requires too many simulation steps. 
 
     X0: array-like or number, optional
         Initial condition (default = 0)
@@ -545,7 +551,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
 
     T: array-like or number, optional
         Time vector, or simulation time duration if a number (time vector is 
-        autocomputed if not given)
+        autocomputed if not given, see :func:`step_response` for more detail)
 
     T_num: number, optional
         Number of time steps to use in simulation if T is not provided as an 
@@ -635,7 +641,7 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
 
     T: array-like or number, optional
         Time vector, or simulation time duration if a number (time vector is 
-        autocomputed if not given)
+        autocomputed if not given; see  :func:`step_response` for more detail)
 
     X0: array-like or number, optional
         Initial condition (default = 0)
@@ -725,7 +731,7 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
 
     T: array-like or number, optional
         Time vector, or simulation time duration if a number (time vector is 
-        autocomputed if not given)
+        autocomputed if not given; see :func:`step_response` for more detail)
 
     X0: array-like or number, optional
         Initial condition (default = 0)

From e155a15e5330eabe4ca2fdac354e692ffb6d6979 Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Sat, 11 Jul 2020 08:47:52 -0700
Subject: [PATCH 10/10] fix unit test precision in sisotool_test.py

---
 control/tests/sisotool_test.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py
index 234d6d28d..f93de54f8 100644
--- a/control/tests/sisotool_test.py
+++ b/control/tests/sisotool_test.py
@@ -33,8 +33,10 @@ def test_sisotool(self):
 
         # Check the step response before moving the point
         step_response_original = np.array(
-            [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)
+            [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,
@@ -76,8 +78,10 @@ def test_sisotool(self):
 
         # Check if the step response has changed
         step_response_moved = np.array(
-            [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)
+            [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__":