From 1f515268bbe43971759139412666c9f8f87adeca Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Wed, 24 Mar 2021 21:46:10 +0100
Subject: [PATCH 01/12] update parameter doc and error message for
 forced_response

---
 control/tests/timeresp_test.py | 10 +--------
 control/timeresp.py            | 37 ++++++++++++++++++++--------------
 2 files changed, 23 insertions(+), 24 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 1c97b9385..88a18483d 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -1,12 +1,4 @@
-"""timeresp_test.py - test time response functions
-
-RMM, 17 Jun 2011 (based on TestMatlab from v0.4c)
-
-This test suite just goes through and calls all of the MATLAB
-functions using different systems and arguments to make sure that
-nothing crashes.  It doesn't test actual functionality; the module
-specific unit tests will do that.
-"""
+"""timeresp_test.py - test time response functions"""
 
 from copy import copy
 from distutils.version import StrictVersion
diff --git a/control/timeresp.py b/control/timeresp.py
index fd2e19c91..479bc3372 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -212,30 +212,32 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         Time steps at which the input is defined; values must be evenly spaced.
 
     U : array_like or float, optional
-        Input array giving input at each time `T` (default = 0).
+        Input array giving input at each time `T`
 
         If `U` is ``None`` or ``0``, a special algorithm is used. This special
         algorithm is faster than the general algorithm, which is used
         otherwise.
 
     X0 : array_like or float, optional
-        Initial condition (default = 0).
+        Initial condition.
 
     transpose : bool, optional
         If True, transpose all input and output arrays (for backward
-        compatibility with MATLAB and :func:`scipy.signal.lsim`).  Default
-        value is False.
+        compatibility with MATLAB and :func:`scipy.signal.lsim`).
 
-    interpolate : bool, optional (default=False)
+    interpolate : bool, optional
         If True and system is a discrete time system, the input will
         be interpolated between the given time steps and the output
         will be given at system sampling rate.  Otherwise, only return
         the output at the times given in `T`.  No effect on continuous
-        time simulations (default = False).
+        time simulations.
 
     return_x : bool, optional
-        If True (default), return the the state vector.  Set to False to
-        return only the time and output vectors.
+        - If False, return only the time and output vectors.
+        - If True, also return the the state vector.
+        - If None, determine the returned variables by
+          config.defaults['forced_response.return_x'], which was True
+          before version 0.9 and is False since then.
 
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then
@@ -243,7 +245,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         squeeze=True, remove single-dimensional entries from the shape of
         the output even if the system is not SISO. If squeeze=False, keep
         the output as a 2D array (indexed by the output number and time)
-        even if the system is SISO. The default value can be set using
+        even if the system is SISO. The default value can be overruled by
         config.defaults['control.squeeze_time_response'].
 
     Returns
@@ -252,13 +254,15 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         Time values of the output.
 
     yout : array
-        Response of the system.  If the system is SISO and squeeze is not
+        Response of the system.  If the system is SISO and `squeeze` is not
         True, the array is 1D (indexed by time).  If the system is not SISO or
-        squeeze is False, the array is 2D (indexed by the output number and
+        `squeeze` is False, the array is 2D (indexed by the output number and
         time).
 
     xout : array
-        Time evolution of the state vector. Not affected by squeeze.
+        Time evolution of the state vector. Not affected by squeeze. Only
+        returned if `return_x` is True, or `return_x` is None and
+        config.defaults['forced_response.return_x'] is True.
 
     See Also
     --------
@@ -297,7 +301,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     sys = _convert_to_statespace(sys)
     A, B, C, D = np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), \
         np.asarray(sys.D)
-#    d_type = A.dtype
+    # d_type = A.dtype
     n_states = A.shape[0]
     n_inputs = B.shape[1]
     n_outputs = C.shape[0]
@@ -332,8 +336,11 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     # T must be array-like and values must be increasing.
     # The length of T determines the length of the input vector.
     if T is None:
-        raise ValueError('Parameter ``T``: must be array-like, and contain '
-                         '(strictly monotonic) increasing numbers.')
+        if not isdtime(sys, strict=True):
+            errmsg_ctime = 'is mandatory for continuous time systems, '
+        raise ValueError('Parameter ``T`` ' + errmsg_ctime + 'must be '
+                         'array-like, and contain (strictly monotonic) '
+                         'increasing numbers.')
     T = _check_convert_array(T, [('any',), (1, 'any')],
                              'Parameter ``T``: ', squeeze=True,
                              transpose=transpose)

From 2870432e48fbba476bc6c04dfd77a791fb80afae Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Thu, 25 Mar 2021 00:03:59 +0100
Subject: [PATCH 02/12] allow T=None with sys.dt=None in forced_response

---
 control/timeresp.py | 27 +++++++++++++++------------
 1 file changed, 15 insertions(+), 12 deletions(-)

diff --git a/control/timeresp.py b/control/timeresp.py
index 479bc3372..4233eee88 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -210,6 +210,9 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
     T : array_like, optional for discrete LTI `sys`
         Time steps at which the input is defined; values must be evenly spaced.
+        If None, `U` must be given and and `len(U)` time steps of sys.dt are
+        simulated. If sys.dt is None or True (undetermined time step), a dt
+        of 1.0 is assumed.
 
     U : array_like or float, optional
         Input array giving input at each time `T`
@@ -245,7 +248,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         squeeze=True, remove single-dimensional entries from the shape of
         the output even if the system is not SISO. If squeeze=False, keep
         the output as a 2D array (indexed by the output number and time)
-        even if the system is SISO. The default value can be overruled by
+        even if the system is SISO. The default value can be overridden by
         config.defaults['control.squeeze_time_response'].
 
     Returns
@@ -310,10 +313,11 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     if U is not None:
         U = np.asarray(U)
     if T is not None:
+        # T must be array-like
         T = np.asarray(T)
 
     # Set and/or check time vector in discrete time case
-    if isdtime(sys, strict=True):
+    if isdtime(sys):
         if T is None:
             if U is None:
                 raise ValueError('Parameters ``T`` and ``U`` can\'t both be'
@@ -323,7 +327,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
                 n_steps = U.shape[0]
             else:
                 n_steps = U.shape[1]
-            T = np.array(range(n_steps)) * (1 if sys.dt is True else sys.dt)
+            dt = 1. if sys.dt in [True, None] else sys.dt
+            T = np.array(range(n_steps)) * dt
         else:
             # Make sure the input vector and time vector have same length
             # TODO: allow interpolation of the input vector
@@ -331,21 +336,19 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
                     (U.ndim > 1 and U.shape[1] != T.shape[0]):
                 ValueError('Pamameter ``T`` must have same elements as'
                            ' the number of columns in input array ``U``')
+    else:
+        if T is None:
+            raise ValueError('Parameter ``T`` is mandatory for continuous '
+                             'time systems.')
 
     # Test if T has shape (n,) or (1, n);
-    # T must be array-like and values must be increasing.
-    # The length of T determines the length of the input vector.
-    if T is None:
-        if not isdtime(sys, strict=True):
-            errmsg_ctime = 'is mandatory for continuous time systems, '
-        raise ValueError('Parameter ``T`` ' + errmsg_ctime + 'must be '
-                         'array-like, and contain (strictly monotonic) '
-                         'increasing numbers.')
     T = _check_convert_array(T, [('any',), (1, 'any')],
                              'Parameter ``T``: ', squeeze=True,
                              transpose=transpose)
+
+    # equally spaced also implies strictly monotonic increase
     dt = T[1] - T[0]
-    if not np.allclose(T[1:] - T[:-1], dt):
+    if not np.allclose(np.diff(T), np.full_like(T[:-1], dt)):
         raise ValueError("Parameter ``T``: time values must be "
                          "equally spaced.")
     n_steps = T.shape[0]            # number of simulation steps

From 517b1eb8d3733f91684286d2155b5630e2525730 Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Thu, 25 Mar 2021 01:32:48 +0100
Subject: [PATCH 03/12] timeresp: add test system with dt=None

---
 control/tests/timeresp_test.py | 24 +++++++++++++++++++-----
 1 file changed, 19 insertions(+), 5 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 88a18483d..a42f71383 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -60,6 +60,14 @@ def siso_ss2(self, siso_ss1):
 
         return T
 
+    @pytest.fixture
+    def siso_ss2_dtnone(self, siso_ss2):
+        """System with unspecified timebase"""
+        ss2 = siso_ss2.sys
+        T = TSys(StateSpace(ss2.A, ss2.B, ss2.C, 0, None))
+        T.t = np.arange(0, 10, 1.)
+        return T
+
     @pytest.fixture
     def siso_tf1(self):
         # Create some transfer functions
@@ -349,7 +357,7 @@ def mimo_tf_step_info(self,
     @pytest.fixture
     def tsystem(self,
                 request,
-                siso_ss1, siso_ss2, siso_tf1, siso_tf2,
+                siso_ss1, siso_ss2, siso_ss2_dtnone, siso_tf1, siso_tf2,
                 mimo_ss1, mimo_ss2, mimo_tf2,
                 siso_dtf0, siso_dtf1, siso_dtf2,
                 siso_dss1, siso_dss2,
@@ -361,6 +369,7 @@ def tsystem(self,
                 siso_tf_asymptotic_from_neg1):
         systems = {"siso_ss1": siso_ss1,
                    "siso_ss2": siso_ss2,
+                   "siso_ss2_dtnone": siso_ss2_dtnone,
                    "siso_tf1": siso_tf1,
                    "siso_tf2": siso_tf2,
                    "mimo_ss1": mimo_ss1,
@@ -840,10 +849,11 @@ def test_default_timevector_functions_d(self, fun, dt):
     @pytest.mark.parametrize("tsystem",
                              ["siso_ss2",   # continuous
                               "siso_tf1",
-                              "siso_dss1",  # no timebase
+                              "siso_dss1",  # unspecified sampling time
                               "siso_dtf1",
                               "siso_dss2",  # matching timebase
                               "siso_dtf2",
+                              "siso_ss2_dtnone",  # undetermined timebase
                               "mimo_ss2",   # MIMO
                               pytest.param("mimo_tf2", marks=slycotonly),
                               "mimo_dss1",
@@ -868,9 +878,9 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout):
             kw['T'] = t
             if fun == forced_response:
                 kw['U'] = np.vstack([np.sin(t) for i in range(sys.ninputs)])
-        elif fun == forced_response and isctime(sys):
+        elif fun == forced_response and isctime(sys, strict=True):
             pytest.skip("No continuous forced_response without time vector.")
-        if hasattr(tsystem.sys, "nstates"):
+        if hasattr(sys, "nstates"):
             kw['X0'] = np.arange(sys.nstates) + 1
         if sys.ninputs > 1 and fun in [step_response, impulse_response]:
             kw['input'] = 1
@@ -884,6 +894,9 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout):
         if hasattr(tsystem, 't'):
             # tout should always match t, which has shape (n, )
             np.testing.assert_allclose(tout, tsystem.t)
+        elif fun == forced_response and sys.dt in [None, True]:
+            np.testing.assert_allclose(
+                np.diff(tout), np.full_like(tout[:-1], 1.))
 
         if squeeze is False or not sys.issiso():
             assert yout.shape[0] == sys.noutputs
@@ -891,7 +904,8 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout):
         else:
             assert yout.shape == tout.shape
 
-        if sys.dt > 0 and sys.dt is not True and not np.isclose(sys.dt, 0.5):
+        if sys.isdtime(strict=True) and sys.dt is not True and not \
+                np.isclose(sys.dt, 0.5):
             kw['T'] = np.arange(0, 5, 0.5)  # incompatible timebase
             with pytest.raises(ValueError):
                 fun(sys, **kw)

From 52a3fd4547fd0b39fb80a057b3d038db765e3110 Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Thu, 25 Mar 2021 13:16:00 +0100
Subject: [PATCH 04/12] more docstring fixups

---
 control/tests/timeresp_test.py |  3 +-
 control/timeresp.py            | 51 +++++++++++++++++-----------------
 doc/conventions.rst            | 28 +++++++++++--------
 3 files changed, 43 insertions(+), 39 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index a42f71383..29e2fd49c 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -895,8 +895,7 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout):
             # tout should always match t, which has shape (n, )
             np.testing.assert_allclose(tout, tsystem.t)
         elif fun == forced_response and sys.dt in [None, True]:
-            np.testing.assert_allclose(
-                np.diff(tout), np.full_like(tout[:-1], 1.))
+            np.testing.assert_allclose(np.diff(tout), 1.)
 
         if squeeze is False or not sys.issiso():
             assert yout.shape[0] == sys.noutputs
diff --git a/control/timeresp.py b/control/timeresp.py
index 4233eee88..ed05b7dae 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -210,32 +210,32 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
     T : array_like, optional for discrete LTI `sys`
         Time steps at which the input is defined; values must be evenly spaced.
-        If None, `U` must be given and and `len(U)` time steps of sys.dt are
-        simulated. If sys.dt is None or True (undetermined time step), a dt
-        of 1.0 is assumed.
+        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.
 
     U : array_like or float, optional
-        Input array giving input at each time `T`
+        Input array giving input at each time `T`.
+        If `U` is None or 0, `T` must be given, even for discrete
+        time systems. In this case, for continuous time systems, a direct
+        calculation of the matrix exponential is used, which is faster than the
+        general interpolating algorithm used otherwise.
 
-        If `U` is ``None`` or ``0``, a special algorithm is used. This special
-        algorithm is faster than the general algorithm, which is used
-        otherwise.
-
-    X0 : array_like or float, optional
+    X0 : array_like or float, default=0.
         Initial condition.
 
-    transpose : bool, optional
+    transpose : bool, default=False
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and :func:`scipy.signal.lsim`).
 
-    interpolate : bool, optional
+    interpolate : bool, default=False
         If True and system is a discrete time system, the input will
         be interpolated between the given time steps and the output
         will be given at system sampling rate.  Otherwise, only return
         the output at the times given in `T`.  No effect on continuous
         time simulations.
 
-    return_x : bool, optional
+    return_x : bool, default=None
         - If False, return only the time and output vectors.
         - If True, also return the the state vector.
         - If None, determine the returned variables by
@@ -245,10 +245,10 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then
         the output response is returned as a 1D array (indexed by time).  If
-        squeeze=True, remove single-dimensional entries from the shape of
-        the output even if the system is not SISO. If squeeze=False, keep
+        `squeeze` is True, remove single-dimensional entries from the shape of
+        the output even if the system is not SISO. If `squeeze` is False, keep
         the output as a 2D array (indexed by the output number and time)
-        even if the system is SISO. The default value can be overridden by
+        even if the system is SISO. The default behavior can be overridden by
         config.defaults['control.squeeze_time_response'].
 
     Returns
@@ -263,7 +263,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         time).
 
     xout : array
-        Time evolution of the state vector. Not affected by squeeze. Only
+        Time evolution of the state vector. Not affected by `squeeze`. Only
         returned if `return_x` is True, or `return_x` is None and
         config.defaults['forced_response.return_x'] is True.
 
@@ -284,7 +284,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     --------
     >>> T, yout, xout = forced_response(sys, T, u, X0)
 
-    See :ref:`time-series-convention`.
+    See :ref:`time-series-convention` and
+    :ref:`package-configuration-parameters`.
 
     """
     if not isinstance(sys, (StateSpace, TransferFunction)):
@@ -301,6 +302,13 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
             "return_x specified for a transfer function system. Internal "
             "conversion to state space used; results may meaningless.")
 
+    # If we are passed a transfer function and X0 is non-zero, warn the user
+    if isinstance(sys, TransferFunction) and np.any(X0 != 0):
+        warnings.warn(
+            "Non-zero initial condition given for transfer function system. "
+            "Internal conversion to state space used; may not be consistent "
+            "with given X0.")
+
     sys = _convert_to_statespace(sys)
     A, B, C, D = np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), \
         np.asarray(sys.D)
@@ -348,7 +356,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
     # equally spaced also implies strictly monotonic increase
     dt = T[1] - T[0]
-    if not np.allclose(np.diff(T), np.full_like(T[:-1], dt)):
+    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
@@ -357,13 +365,6 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)],
                               'Parameter ``X0``: ', squeeze=True)
 
-    # If we are passed a transfer function and X0 is non-zero, warn the user
-    if isinstance(sys, TransferFunction) and np.any(X0 != 0):
-        warnings.warn(
-            "Non-zero initial condition given for transfer function system. "
-            "Internal conversion to state space used; may not be consistent "
-            "with given X0.")
-
     xout = np.zeros((n_states, n_steps))
     xout[:, 0] = X0
     yout = np.zeros((n_outputs, n_steps))
diff --git a/doc/conventions.rst b/doc/conventions.rst
index 4a3d78926..8a4d11992 100644
--- a/doc/conventions.rst
+++ b/doc/conventions.rst
@@ -83,17 +83,17 @@ The timebase argument can be given when a system is constructed:
 * dt = 0: continuous time system (default)
 * dt > 0: discrete time system with sampling period 'dt'
 * dt = True: discrete time with unspecified sampling period
-* dt = None: no timebase specified 
+* dt = None: no timebase specified
 
 Only the :class:`StateSpace`, :class:`TransferFunction`, and
 :class:`InputOutputSystem` classes allow explicit representation of
 discrete time systems.
 
-Systems must have compatible timebases in order to be combined. A discrete time 
-system with unspecified sampling time (`dt = True`) can be combined with a system 
+Systems must have compatible timebases in order to be combined. A discrete time
+system with unspecified sampling time (`dt = True`) can be combined with a system
 having a specified sampling time; the result will be a discrete time system with the sample time of the latter
 system.  Similarly, a system with timebase `None` can be combined with a system having a specified
-timebase; the result will have the timebase of the latter system. For continuous 
+timebase; the result will have the timebase of the latter system. For continuous
 time systems, the :func:`sample_system` function or the :meth:`StateSpace.sample` and :meth:`TransferFunction.sample` methods
 can be used to create a discrete time system from a continuous time system.
 See :ref:`utility-and-conversions`. The default value of 'dt' can be changed by
@@ -179,6 +179,10 @@ can be computed like this::
 
     ft = D * U
 
+
+.. currentmodule:: control
+.. _package-configuration-parameters:
+
 Package configuration parameters
 ================================
 
@@ -207,25 +211,25 @@ on standard configurations.
 Selected variables that can be configured, along with their default values:
 
   * bode.dB (False): Bode plot magnitude plotted in dB (otherwise powers of 10)
-    
+
   * bode.deg (True): Bode plot phase plotted in degrees (otherwise radians)
-    
+
   * bode.Hz (False): Bode plot frequency plotted in Hertz (otherwise rad/sec)
-    
+
   * bode.grid (True): Include grids for magnitude and phase plots
-    
+
   * freqplot.number_of_samples (None): Number of frequency points in Bode plots
-    
+
   * freqplot.feature_periphery_decade (1.0): How many decades to include in the
     frequency range on both sides of features (poles, zeros).
-    
+
   * statesp.use_numpy_matrix (True): set the return type for state space matrices to
     `numpy.matrix` (verus numpy.ndarray)
 
-  * statesp.default_dt and xferfcn.default_dt (None): set the default value of dt when 
+  * statesp.default_dt and xferfcn.default_dt (None): set the default value of dt when
     constructing new LTI systems
 
-  * statesp.remove_useless_states (True): remove states that have no effect on the 
+  * statesp.remove_useless_states (True): remove states that have no effect on the
     input-output dynamics of the system
 
 Additional parameter variables are documented in individual functions

From f0d764f2dd34158bb44a23127c56a64fbb21e208 Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Thu, 25 Mar 2021 16:10:40 +0100
Subject: [PATCH 05/12] test forced_response omitting parameters as documented

---
 control/tests/timeresp_test.py | 60 ++++++++++++++++++++++++++++++++++
 control/timeresp.py            | 15 +++++----
 2 files changed, 68 insertions(+), 7 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 29e2fd49c..e6232dc13 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -66,6 +66,8 @@ def siso_ss2_dtnone(self, siso_ss2):
         ss2 = siso_ss2.sys
         T = TSys(StateSpace(ss2.A, ss2.B, ss2.C, 0, None))
         T.t = np.arange(0, 10, 1.)
+        T.ystep = np.array([    0.,    86.,   -72.,   230.,  -360.,   806.,
+                            -1512.,  3110., -6120., 12326.])
         return T
 
     @pytest.fixture
@@ -128,12 +130,18 @@ def siso_dtf0(self):
     def siso_dtf1(self):
         T =  TSys(TransferFunction([1], [1, 1, 0.25], True))
         T.t = np.arange(0, 5, 1)
+        T.ystep = np.array([0.  , 0.  , 1.  , 0.  , 0.75])
         return T
 
     @pytest.fixture
     def siso_dtf2(self):
         T = TSys(TransferFunction([1], [1, 1, 0.25], 0.2))
         T.t = np.arange(0, 5, 0.2)
+        T.ystep =np.array([0.    , 0.    , 1.    , 0.    , 0.75  , 0.25  ,
+                           0.5625, 0.375 , 0.4844, 0.4219, 0.457 , 0.4375,
+                           0.4482, 0.4424, 0.4456, 0.4438, 0.4448, 0.4443,
+                           0.4445, 0.4444, 0.4445, 0.4444, 0.4445, 0.4444,
+                           0.4444])
         return T
 
     @pytest.fixture
@@ -719,6 +727,58 @@ def test_forced_response_legacy(self):
         t, y = ct.forced_response(sys, T, U)
         t, y, x = ct.forced_response(sys, T, U, return_x=True)
 
+    @pytest.mark.parametrize(
+        "tsystem, fr_kwargs, refattr",
+        [pytest.param("siso_ss1",
+                      {'X0': [0.5, 1], 'T':  np.linspace(0, 1, 10)},
+                      'yinitial',
+                      id="ctime no T"),
+         pytest.param("siso_dtf1",
+                      {'U': np.ones(5,)}, 'ystep',
+                      id="dt=True, no U"),
+         pytest.param("siso_dtf2",
+                      {'U': np.ones(25,)}, 'ystep',
+                      id="dt=0.2, no U"),
+         pytest.param("siso_ss2_dtnone",
+                      {'U': np.ones(10,)}, 'ystep',
+                      id="dt=None, no U")],
+        indirect=["tsystem"])
+    def test_forced_response_T_U(self, tsystem, fr_kwargs, refattr):
+        """Test documented forced_response behavior for parameters T and U."""
+        t, y = forced_response(tsystem.sys, **fr_kwargs)
+        np.testing.assert_allclose(t, tsystem.t)
+        np.testing.assert_allclose(y, getattr(tsystem, refattr), rtol=1e-3)
+
+    def test_forced_response_invalid(self, siso_ss1, siso_dss2):
+        """Test invalid parameters."""
+        with pytest.raises(TypeError,
+                           match="StateSpace.*or.*TransferFunction"):
+            forced_response("not a system")
+
+        # ctime
+        with pytest.raises(ValueError, match="T.*is mandatory for continuous"):
+            forced_response(siso_ss1.sys)
+        with pytest.raises(ValueError, match="time values must be equally "
+                                             "spaced"):
+            forced_response(siso_ss1.sys, [0, 0.1, 0.12, 0.4])
+
+        # dtime with sys.dt > 0
+        with pytest.raises(ValueError, match="can't both be zero"):
+            forced_response(siso_dss2.sys)
+        with pytest.raises(ValueError, match="must have same elements"):
+            forced_response(siso_dss2.sys,
+                            T=siso_dss2.t, U=np.random.randn(1, 12))
+        with pytest.raises(ValueError, match="must have same elements"):
+            forced_response(siso_dss2.sys,
+                            T=siso_dss2.t, U=np.random.randn(12))
+        with pytest.raises(ValueError, match="must match sampling time"):
+            forced_response(siso_dss2.sys, T=siso_dss2.t*0.9)
+        with pytest.raises(ValueError, match="must be multiples of "
+                                             "sampling time"):
+            forced_response(siso_dss2.sys, T=siso_dss2.t*1.1)
+        # but this is ok
+        forced_response(siso_dss2.sys, T=siso_dss2.t*2)
+
 
     @pytest.mark.parametrize("u, x0, xtrue",
                              [(np.zeros((10,)),
diff --git a/control/timeresp.py b/control/timeresp.py
index ed05b7dae..89202ca79 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -327,8 +327,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     # Set and/or check time vector in discrete time case
     if isdtime(sys):
         if T is None:
-            if U is None:
-                raise ValueError('Parameters ``T`` and ``U`` can\'t both be'
+            if U is None or (U.ndim == 0 and U == 0.):
+                raise ValueError('Parameters ``T`` and ``U`` can\'t both be '
                                  'zero for discrete-time simulation')
             # Set T to equally spaced samples with same length as U
             if U.ndim == 1:
@@ -339,11 +339,12 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
             T = np.array(range(n_steps)) * dt
         else:
             # Make sure the input vector and time vector have same length
-            # TODO: allow interpolation of the input vector
             if (U.ndim == 1 and U.shape[0] != T.shape[0]) or \
                     (U.ndim > 1 and U.shape[1] != T.shape[0]):
-                ValueError('Pamameter ``T`` must have same elements as'
-                           ' the number of columns in input array ``U``')
+                raise ValueError('Pamameter ``T`` must have same elements as'
+                                 ' the number of columns in input array ``U``')
+            if U.ndim == 0:
+                U = np.full((n_inputs, T.shape[0]), U)
     else:
         if T is None:
             raise ValueError('Parameter ``T`` is mandatory for continuous '
@@ -370,7 +371,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     yout = np.zeros((n_outputs, n_steps))
 
     # Separate out the discrete and continuous time cases
-    if isctime(sys):
+    if isctime(sys, strict=True):
         # Solve the differential equation, copied from scipy.signal.ltisys.
         dot = np.dot  # Faster and shorter code
 
@@ -421,7 +422,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
     else:
         # Discrete type system => use SciPy signal processing toolbox
-        if sys.dt is not True:
+        if sys.dt is not True and sys.dt is not None:
             # Make sure that the time increment is a multiple of sampling time
 
             # First make sure that time increment is bigger than sampling time

From 9cac9fa0daa4479d2ecd3ee27b92203590d3ed27 Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Fri, 26 Mar 2021 14:08:54 +0100
Subject: [PATCH 06/12] ensure dlsim output length

---
 control/timeresp.py | 31 ++++++++++++++++++++++---------
 1 file changed, 22 insertions(+), 9 deletions(-)

diff --git a/control/timeresp.py b/control/timeresp.py
index 89202ca79..8c0d59718 100644
--- 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)
 
 

From 7be508b05c68b22b229ebe19fe6e480ccc27198c Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Fri, 26 Mar 2021 23:48:32 +0100
Subject: [PATCH 07/12] rework fixture setup for timeresp_test

---
 control/tests/timeresp_test.py | 562 ++++++++++++++-------------------
 1 file changed, 240 insertions(+), 322 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index e6232dc13..b4a3598e5 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -30,60 +30,45 @@ def __repr__(self):
 class TestTimeresp:
 
     @pytest.fixture
-    def siso_ss1(self):
+    def tsystem(self, request):
+        """Define some test systems"""
 
+        """continuous"""
         A = np.array([[1., -2.], [3., -4.]])
         B = np.array([[5.], [7.]])
         C = np.array([[6., 8.]])
         D = np.array([[9.]])
-        T = TSys(StateSpace(A, B, C, D, 0))
-
-        T.t = np.linspace(0, 1, 10)
-        T.ystep = np.array([9., 17.6457, 24.7072, 30.4855, 35.2234,
-                            39.1165, 42.3227, 44.9694, 47.1599, 48.9776])
-
-        T.yinitial = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118,
-                               1.9092, 1.1508, 0.5833, 0.1645, -0.1391])
-
-        return T
-
-    @pytest.fixture
-    def siso_ss2(self, siso_ss1):
-        """System siso_ss2 with D=0"""
+        siso_ss1 = TSys(StateSpace(A, B, C, D, 0))
+        siso_ss1.t = np.linspace(0, 1, 10)
+        siso_ss1.ystep = np.array([9., 17.6457, 24.7072, 30.4855, 35.2234,
+                                   39.1165, 42.3227, 44.9694, 47.1599,
+                                   48.9776])
+        siso_ss1.yinitial = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118,
+                                      1.9092, 1.1508, 0.5833, 0.1645, -0.1391])
         ss1 = siso_ss1.sys
-        T = TSys(StateSpace(ss1.A, ss1.B, ss1.C, 0, 0))
-        T.t = siso_ss1.t
-        T.ystep = siso_ss1.ystep - 9
-        T.initial = siso_ss1.yinitial - 9
-        T.yimpulse = np.array([86., 70.1808, 57.3753, 46.9975, 38.5766,
-                               31.7344, 26.1668, 21.6292, 17.9245, 14.8945])
 
-        return T
+        """D=0, continuous"""
+        siso_ss2 = TSys(StateSpace(ss1.A, ss1.B, ss1.C, 0, 0))
+        siso_ss2.t = siso_ss1.t
+        siso_ss2.ystep = siso_ss1.ystep - 9
+        siso_ss2.initial = siso_ss1.yinitial - 9
+        siso_ss2.yimpulse = np.array([86., 70.1808, 57.3753, 46.9975, 38.5766,
+                                      31.7344, 26.1668, 21.6292, 17.9245,
+                                      14.8945])
 
-    @pytest.fixture
-    def siso_ss2_dtnone(self, siso_ss2):
         """System with unspecified timebase"""
-        ss2 = siso_ss2.sys
-        T = TSys(StateSpace(ss2.A, ss2.B, ss2.C, 0, None))
-        T.t = np.arange(0, 10, 1.)
-        T.ystep = np.array([    0.,    86.,   -72.,   230.,  -360.,   806.,
-                            -1512.,  3110., -6120., 12326.])
-        return T
+        siso_ss2_dtnone = TSys(StateSpace(ss1.A, ss1.B, ss1.C, 0, None))
+        siso_ss2_dtnone.t = np.arange(0, 10, 1.)
+        siso_ss2_dtnone.ystep = np.array([0., 86., -72., 230., -360.,  806.,
+                                          -1512.,  3110., -6120., 12326.])
 
-    @pytest.fixture
-    def siso_tf1(self):
-        # Create some transfer functions
-        return TSys(TransferFunction([1], [1, 2, 1], 0))
+        siso_tf1 = TSys(TransferFunction([1], [1, 2, 1], 0))
 
-    @pytest.fixture
-    def siso_tf2(self, siso_ss1):
-        T = copy(siso_ss1)
-        T.sys = ss2tf(siso_ss1.sys)
-        return T
+        siso_tf2 = copy(siso_ss1)
+        siso_tf2.sys = ss2tf(siso_ss1.sys)
 
-    @pytest.fixture
-    def mimo_ss1(self, siso_ss1):
-        # Create MIMO system, contains ``siso_ss1`` twice
+        """MIMO system, contains ``siso_ss1`` twice"""
+        mimo_ss1 = copy(siso_ss1)
         A = np.zeros((4, 4))
         A[:2, :2] = siso_ss1.sys.A
         A[2:, 2:] = siso_ss1.sys.A
@@ -96,13 +81,10 @@ def mimo_ss1(self, siso_ss1):
         D = np.zeros((2, 2))
         D[:1, :1] = siso_ss1.sys.D
         D[1:, 1:] = siso_ss1.sys.D
-        T = copy(siso_ss1)
-        T.sys = StateSpace(A, B, C, D)
-        return T
+        mimo_ss1.sys = StateSpace(A, B, C, D)
 
-    @pytest.fixture
-    def mimo_ss2(self, siso_ss2):
-        # Create MIMO system, contains ``siso_ss2`` twice
+        """MIMO system, contains ``siso_ss2`` twice"""
+        mimo_ss2 = copy(siso_ss2)
         A = np.zeros((4, 4))
         A[:2, :2] = siso_ss2.sys.A
         A[2:, 2:] = siso_ss2.sys.A
@@ -113,99 +95,74 @@ def mimo_ss2(self, siso_ss2):
         C[:1, :2] = siso_ss2.sys.C
         C[1:, 2:] = siso_ss2.sys.C
         D = np.zeros((2, 2))
-        T = copy(siso_ss2)
-        T.sys = StateSpace(A, B, C, D, 0)
-        return T
-
-    # Create discrete time systems
-
-    @pytest.fixture
-    def siso_dtf0(self):
-        T = TSys(TransferFunction([1.], [1., 0.], 1.))
-        T.t = np.arange(4)
-        T.yimpulse = [0., 1., 0., 0.]
-        return T
-
-    @pytest.fixture
-    def siso_dtf1(self):
-        T =  TSys(TransferFunction([1], [1, 1, 0.25], True))
-        T.t = np.arange(0, 5, 1)
-        T.ystep = np.array([0.  , 0.  , 1.  , 0.  , 0.75])
-        return T
-
-    @pytest.fixture
-    def siso_dtf2(self):
-        T = TSys(TransferFunction([1], [1, 1, 0.25], 0.2))
-        T.t = np.arange(0, 5, 0.2)
-        T.ystep =np.array([0.    , 0.    , 1.    , 0.    , 0.75  , 0.25  ,
-                           0.5625, 0.375 , 0.4844, 0.4219, 0.457 , 0.4375,
-                           0.4482, 0.4424, 0.4456, 0.4438, 0.4448, 0.4443,
-                           0.4445, 0.4444, 0.4445, 0.4444, 0.4445, 0.4444,
-                           0.4444])
-        return T
-
-    @pytest.fixture
-    def siso_dss1(self, siso_dtf1):
-        T = copy(siso_dtf1)
-        T.sys = tf2ss(siso_dtf1.sys)
-        return T
-
-    @pytest.fixture
-    def siso_dss2(self, siso_dtf2):
-        T = copy(siso_dtf2)
-        T.sys = tf2ss(siso_dtf2.sys)
-        return T
-
-    @pytest.fixture
-    def mimo_dss1(self, mimo_ss1):
-        ss1 = mimo_ss1.sys
-        T = TSys(
-            StateSpace(ss1.A, ss1.B, ss1.C, ss1.D, True))
-        T.t = np.arange(0, 5, 0.2)
-        return T
-
-    @pytest.fixture
-    def mimo_dss2(self, mimo_ss1):
-        T = copy(mimo_ss1)
-        T.sys = c2d(mimo_ss1.sys, T.t[1]-T.t[0])
-        return T
-
-    @pytest.fixture
-    def mimo_tf2(self, siso_ss2, mimo_ss2):
-        T = copy(mimo_ss2)
-        # construct from siso to avoid slycot during fixture setup
+        mimo_ss2.sys = StateSpace(A, B, C, D, 0)
+
+        """discrete"""
+        siso_dtf0 = TSys(TransferFunction([1.], [1., 0.], 1.))
+        siso_dtf0.t = np.arange(4)
+        siso_dtf0.yimpulse = [0., 1., 0., 0.]
+
+        siso_dtf1 =  TSys(TransferFunction([1], [1, 1, 0.25], True))
+        siso_dtf1.t = np.arange(0, 5, 1)
+        siso_dtf1.ystep = np.array([0.  , 0.  , 1.  , 0.  , 0.75])
+
+        siso_dtf2 = TSys(TransferFunction([1], [1, 1, 0.25], 0.2))
+        siso_dtf2.t = np.arange(0, 5, 0.2)
+        siso_dtf2.ystep = np.array(
+            [0.    , 0.    , 1.    , 0.    , 0.75  , 0.25  ,
+             0.5625, 0.375 , 0.4844, 0.4219, 0.457 , 0.4375,
+             0.4482, 0.4424, 0.4456, 0.4438, 0.4448, 0.4443,
+             0.4445, 0.4444, 0.4445, 0.4444, 0.4445, 0.4444,
+             0.4444])
+
+        """Time step which leads to rounding errors for time vector length"""
+        num = [-0.10966442, 0.12431949]
+        den = [1., -1.86789511, 0.88255018]
+        dt = 0.12493963338370018
+        siso_dtf3 = TSys(TransferFunction(num, den, dt))
+        siso_dtf3.t = np.linspace(0, 9*dt, 10)
+        siso_dtf3.ystep = np.array(
+            [ 0.    , -0.1097, -0.1902, -0.2438, -0.2729,
+             -0.2799, -0.2674, -0.2377, -0.1934, -0.1368])
+
+        siso_dss1 = copy(siso_dtf1)
+        siso_dss1.sys = tf2ss(siso_dtf1.sys)
+
+        siso_dss2 = copy(siso_dtf2)
+        siso_dss2.sys = tf2ss(siso_dtf2.sys)
+
+        mimo_dss1 = TSys(StateSpace(ss1.A, ss1.B, ss1.C, ss1.D, True))
+        mimo_dss1.t = np.arange(0, 5, 0.2)
+
+        mimo_dss2 = copy(mimo_ss1)
+        mimo_dss2.sys = c2d(mimo_ss1.sys, mimo_ss1.t[1]-mimo_ss1.t[0])
+
+        mimo_tf2 = copy(mimo_ss2)
         tf_ = ss2tf(siso_ss2.sys)
-        T.sys = TransferFunction([[tf_.num[0][0], [0]], [[0], tf_.num[0][0]]],
-                                 [[tf_.den[0][0], [1]], [[1], tf_.den[0][0]]],
-                                 0)
-        return T
+        mimo_tf2.sys = TransferFunction(
+            [[tf_.num[0][0], [0]], [[0], tf_.num[0][0]]],
+            [[tf_.den[0][0], [1]], [[1], tf_.den[0][0]]],
+            0)
 
-    @pytest.fixture
-    def mimo_dtf1(self, siso_dtf1):
-        T = copy(siso_dtf1)
-        # construct from siso to avoid slycot during fixture setup
+        mimo_dtf1 = copy(siso_dtf1)
         tf_ = siso_dtf1.sys
-        T.sys = TransferFunction([[tf_.num[0][0], [0]], [[0], tf_.num[0][0]]],
-                                 [[tf_.den[0][0], [1]], [[1], tf_.den[0][0]]],
-                                 True)
-        return T
+        mimo_dtf1.sys = TransferFunction(
+            [[tf_.num[0][0], [0]], [[0], tf_.num[0][0]]],
+            [[tf_.den[0][0], [1]], [[1], tf_.den[0][0]]],
+            True)
 
-    @pytest.fixture
-    def pole_cancellation(self):
         # for pole cancellation tests
-        return TransferFunction([1.067e+05, 5.791e+04],
-                                [10.67, 1.067e+05, 5.791e+04])
+        pole_cancellation = TSys(TransferFunction(
+            [1.067e+05, 5.791e+04],
+            [10.67, 1.067e+05, 5.791e+04]))
 
-    @pytest.fixture
-    def no_pole_cancellation(self):
-        return TransferFunction([1.881e+06],
-                                [188.1, 1.881e+06])
+        no_pole_cancellation = TSys(TransferFunction(
+            [1.881e+06],
+            [188.1, 1.881e+06]))
 
-    @pytest.fixture
-    def siso_tf_type1(self):
         # System Type 1 - Step response not stationary:  G(s)=1/s(s+1)
-        T = TSys(TransferFunction(1, [1, 1, 0]))
-        T.step_info = {
+        siso_tf_type1 = TSys(TransferFunction(1, [1, 1, 0]))
+        siso_tf_type1.step_info = {
              'RiseTime': np.NaN,
              'SettlingTime': np.NaN,
              'SettlingMin': np.NaN,
@@ -215,14 +172,11 @@ def siso_tf_type1(self):
              'Peak': np.Inf,
              'PeakTime': np.Inf,
              'SteadyStateValue': np.NaN}
-        return T
 
-    @pytest.fixture
-    def siso_tf_kpos(self):
         # SISO under shoot response and positive final value
         # G(s)=(-s+1)/(s²+s+1)
-        T = TSys(TransferFunction([-1, 1], [1, 1, 1]))
-        T.step_info = {
+        siso_tf_kpos = TSys(TransferFunction([-1, 1], [1, 1, 1]))
+        siso_tf_kpos.step_info = {
              'RiseTime': 1.242,
              'SettlingTime': 9.110,
              'SettlingMin': 0.90,
@@ -232,14 +186,11 @@ def siso_tf_kpos(self):
              'Peak': 1.208,
              'PeakTime': 4.282,
              'SteadyStateValue': 1.0}
-        return T
 
-    @pytest.fixture
-    def siso_tf_kneg(self):
         # SISO under shoot response and negative final value
         # k=-1 G(s)=-(-s+1)/(s²+s+1)
-        T = TSys(TransferFunction([1, -1], [1, 1, 1]))
-        T.step_info = {
+        siso_tf_kneg = TSys(TransferFunction([1, -1], [1, 1, 1]))
+        siso_tf_kneg.step_info = {
              'RiseTime': 1.242,
              'SettlingTime': 9.110,
              'SettlingMin': -1.208,
@@ -249,13 +200,9 @@ def siso_tf_kneg(self):
              'Peak': 1.208,
              'PeakTime': 4.282,
              'SteadyStateValue': -1.0}
-        return T
 
-    @pytest.fixture
-    def siso_tf_asymptotic_from_neg1(self):
-        # Peak_value = Undershoot = y_final(y(t=inf))
-        T = TSys(TransferFunction([-1, 1], [1, 1]))
-        T.step_info = {
+        siso_tf_asymptotic_from_neg1 = TSys(TransferFunction([-1, 1], [1, 1]))
+        siso_tf_asymptotic_from_neg1.step_info = {
             'RiseTime': 2.197,
             'SettlingTime': 4.605,
             'SettlingMin': 0.9,
@@ -265,15 +212,14 @@ def siso_tf_asymptotic_from_neg1(self):
             'Peak': 1.0,
             'PeakTime': 0.0,
             'SteadyStateValue': 1.0}
-        T.kwargs = {'step_info': {'T': np.arange(0, 5, 1e-3)}}
-        return T
+        siso_tf_asymptotic_from_neg1.kwargs = {
+            'step_info': {'T': np.arange(0, 5, 1e-3)}}
 
-    @pytest.fixture
-    def siso_tf_step_matlab(self):
         # example from matlab online help
         # https://www.mathworks.com/help/control/ref/stepinfo.html
-        T = TSys(TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2]))
-        T.step_info = {
+        siso_tf_step_matlab = TSys(TransferFunction([1, 5, 5],
+                                                    [1, 1.65, 5, 6.5, 2]))
+        siso_tf_step_matlab.step_info = {
             'RiseTime': 3.8456,
             'SettlingTime': 27.9762,
             'SettlingMin': 2.0689,
@@ -283,10 +229,7 @@ def siso_tf_step_matlab(self):
             'Peak': 2.6873,
             'PeakTime': 8.0530,
             'SteadyStateValue': 2.5}
-        return T
 
-    @pytest.fixture
-    def mimo_ss_step_matlab(self):
         A = [[0.68, -0.34],
              [0.34, 0.68]]
         B = [[0.18, -0.05],
@@ -295,114 +238,70 @@ def mimo_ss_step_matlab(self):
              [-1.12, -1.10]]
         D = [[0, 0],
              [0.06, -0.37]]
-        T = TSys(StateSpace(A, B, C, D, 0.2))
-        T.kwargs['step_info'] = {'T': 4.6}
-        T.step_info = [[{'RiseTime': 0.6000,
-                         'SettlingTime': 3.0000,
-                         'SettlingMin': -0.5999,
-                         'SettlingMax': -0.4689,
-                         'Overshoot': 15.5072,
-                         'Undershoot': 0.,
-                         'Peak': 0.5999,
-                         'PeakTime': 1.4000,
-                         'SteadyStateValue': -0.5193},
-                        {'RiseTime': 0.,
-                         'SettlingTime': 3.6000,
-                         'SettlingMin': -0.2797,
-                         'SettlingMax': -0.1043,
-                         'Overshoot': 118.9918,
-                         'Undershoot': 0,
-                         'Peak': 0.2797,
-                         'PeakTime': .6000,
-                         'SteadyStateValue': -0.1277}],
-                       [{'RiseTime': 0.4000,
-                         'SettlingTime': 2.8000,
-                         'SettlingMin': -0.6724,
-                         'SettlingMax': -0.5188,
-                         'Overshoot': 24.6476,
-                         'Undershoot': 11.1224,
-                         'Peak': 0.6724,
-                         'PeakTime': 1,
-                         'SteadyStateValue': -0.5394},
-                         {'RiseTime': 0.0000, # (*)
-                         'SettlingTime': 3.4000,
-                         'SettlingMin': -0.4350,  # (*)
-                         'SettlingMax': -0.1485,
-                         'Overshoot': 132.0170,
-                         'Undershoot': 0.,
-                         'Peak': 0.4350,
-                         'PeakTime': .2,
-                         'SteadyStateValue': -0.1875}]]
-                         # (*): MATLAB gives 0.4 for RiseTime and -0.1034 for
-                         # SettlingMin, but it is unclear what 10% and 90% of
-                         # the steady state response mean, when the step for
-                         # this channel does not start a 0.
-        return T
-
-    @pytest.fixture
-    def siso_ss_step_matlab(self, mimo_ss_step_matlab):
-        T = copy(mimo_ss_step_matlab)
-        T.sys = T.sys[1, 0]
-        T.step_info = T.step_info[1][0]
-        return T
+        mimo_ss_step_matlab = TSys(StateSpace(A, B, C, D, 0.2))
+        mimo_ss_step_matlab.kwargs['step_info'] = {'T': 4.6}
+        mimo_ss_step_matlab.step_info = [[
+             {'RiseTime': 0.6000,
+              'SettlingTime': 3.0000,
+              'SettlingMin': -0.5999,
+              'SettlingMax': -0.4689,
+              'Overshoot': 15.5072,
+              'Undershoot': 0.,
+              'Peak': 0.5999,
+              'PeakTime': 1.4000,
+              'SteadyStateValue': -0.5193},
+             {'RiseTime': 0.,
+              'SettlingTime': 3.6000,
+              'SettlingMin': -0.2797,
+              'SettlingMax': -0.1043,
+              'Overshoot': 118.9918,
+              'Undershoot': 0,
+              'Peak': 0.2797,
+              'PeakTime': .6000,
+              'SteadyStateValue': -0.1277}],
+            [{'RiseTime': 0.4000,
+              'SettlingTime': 2.8000,
+              'SettlingMin': -0.6724,
+              'SettlingMax': -0.5188,
+              'Overshoot': 24.6476,
+              'Undershoot': 11.1224,
+              'Peak': 0.6724,
+              'PeakTime': 1,
+              'SteadyStateValue': -0.5394},
+              {'RiseTime': 0.0000, # (*)
+              'SettlingTime': 3.4000,
+              'SettlingMin': -0.4350,  # (*)
+              'SettlingMax': -0.1485,
+              'Overshoot': 132.0170,
+              'Undershoot': 0.,
+              'Peak': 0.4350,
+              'PeakTime': .2,
+              'SteadyStateValue': -0.1875}]]
+              # (*): MATLAB gives 0.4 for RiseTime and -0.1034 for
+              # SettlingMin, but it is unclear what 10% and 90% of
+              # the steady state response mean, when the step for
+              # this channel does not start a 0.
+
+        siso_ss_step_matlab = copy(mimo_ss_step_matlab)
+        siso_ss_step_matlab.sys = siso_ss_step_matlab.sys[1, 0]
+        siso_ss_step_matlab.step_info = siso_ss_step_matlab.step_info[1][0]
 
-    @pytest.fixture
-    def mimo_tf_step_info(self,
-                          siso_tf_kpos, siso_tf_kneg,
-                          siso_tf_step_matlab):
         Ta = [[siso_tf_kpos, siso_tf_kneg, siso_tf_step_matlab],
-               [siso_tf_step_matlab, siso_tf_kpos, siso_tf_kneg]]
-        T = TSys(TransferFunction(
+              [siso_tf_step_matlab, siso_tf_kpos, siso_tf_kneg]]
+        mimo_tf_step_info = TSys(TransferFunction(
             [[Ti.sys.num[0][0] for Ti in Tr] for Tr in Ta],
             [[Ti.sys.den[0][0] for Ti in Tr] for Tr in Ta]))
-        T.step_info = [[Ti.step_info for Ti in Tr] for Tr in Ta]
+        mimo_tf_step_info.step_info = [[Ti.step_info for Ti in Tr]
+                                       for Tr in Ta]
         # enforce enough sample points for all channels (they have different
         # characteristics)
-        T.kwargs['step_info'] = {'T_num': 2000}
-        return T
-
+        mimo_tf_step_info.kwargs['step_info'] = {'T_num': 2000}
 
-    @pytest.fixture
-    def tsystem(self,
-                request,
-                siso_ss1, siso_ss2, siso_ss2_dtnone, siso_tf1, siso_tf2,
-                mimo_ss1, mimo_ss2, mimo_tf2,
-                siso_dtf0, siso_dtf1, siso_dtf2,
-                siso_dss1, siso_dss2,
-                mimo_dss1, mimo_dss2, mimo_dtf1,
-                pole_cancellation, no_pole_cancellation, siso_tf_type1,
-                siso_tf_kpos, siso_tf_kneg,
-                siso_tf_step_matlab, siso_ss_step_matlab,
-                mimo_ss_step_matlab, mimo_tf_step_info,
-                siso_tf_asymptotic_from_neg1):
-        systems = {"siso_ss1": siso_ss1,
-                   "siso_ss2": siso_ss2,
-                   "siso_ss2_dtnone": siso_ss2_dtnone,
-                   "siso_tf1": siso_tf1,
-                   "siso_tf2": siso_tf2,
-                   "mimo_ss1": mimo_ss1,
-                   "mimo_ss2": mimo_ss2,
-                   "mimo_tf2": mimo_tf2,
-                   "siso_dtf0": siso_dtf0,
-                   "siso_dtf1": siso_dtf1,
-                   "siso_dtf2": siso_dtf2,
-                   "siso_dss1": siso_dss1,
-                   "siso_dss2": siso_dss2,
-                   "mimo_dss1": mimo_dss1,
-                   "mimo_dss2": mimo_dss2,
-                   "mimo_dtf1": mimo_dtf1,
-                   "pole_cancellation": pole_cancellation,
-                   "no_pole_cancellation": no_pole_cancellation,
-                   "siso_tf_type1": siso_tf_type1,
-                   "siso_tf_kpos": siso_tf_kpos,
-                   "siso_tf_kneg": siso_tf_kneg,
-                   "siso_tf_step_matlab": siso_tf_step_matlab,
-                   "siso_ss_step_matlab": siso_ss_step_matlab,
-                   "mimo_ss_step_matlab": mimo_ss_step_matlab,
-                   "mimo_tf_step": mimo_tf_step_info,
-                   "siso_tf_asymptotic_from_neg1": siso_tf_asymptotic_from_neg1,
-                   }
-        return systems[request.param]
+        systems = locals()
+        if isinstance(request.param, str):
+            return systems[request.param]
+        else:
+            return [systems[sys] for sys in request.param]
 
     @pytest.mark.parametrize(
         "kwargs",
@@ -411,11 +310,12 @@ def tsystem(self,
          {'X0': np.array([0, 0])},
          {'X0': 0, 'return_x': True},
          ])
-    def test_step_response_siso(self, siso_ss1, kwargs):
+    @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)
+    def test_step_response_siso(self, tsystem, kwargs):
         """Test SISO system step response"""
-        sys = siso_ss1.sys
-        t = siso_ss1.t
-        yref = siso_ss1.ystep
+        sys = tsystem.sys
+        t = tsystem.t
+        yref = tsystem.ystep
         # SISO call
         out = step_response(sys, T=t, **kwargs)
         tout, yout = out[:2]
@@ -423,19 +323,21 @@ def test_step_response_siso(self, siso_ss1, kwargs):
         np.testing.assert_array_almost_equal(tout, t)
         np.testing.assert_array_almost_equal(yout, yref, decimal=4)
 
-    def test_step_response_mimo(self, mimo_ss1):
-        """Test MIMO system, which contains ``siso_ss1`` twice"""
-        sys = mimo_ss1.sys
-        t = mimo_ss1.t
-        yref = mimo_ss1.ystep
+    @pytest.mark.parametrize("tsystem", ["mimo_ss1"], indirect=True)
+    def test_step_response_mimo(self, tsystem):
+        """Test MIMO system, which contains ``siso_ss1`` twice."""
+        sys = tsystem.sys
+        t = tsystem.t
+        yref = tsystem.ystep
         _t, y_00 = step_response(sys, T=t, input=0, output=0)
         _t, y_11 = step_response(sys, T=t, input=1, output=1)
         np.testing.assert_array_almost_equal(y_00, yref, decimal=4)
         np.testing.assert_array_almost_equal(y_11, yref, decimal=4)
 
-    def test_step_response_return(self, mimo_ss1):
-        """Verify continuous and discrete time use same return conventions"""
-        sysc = mimo_ss1.sys
+    @pytest.mark.parametrize("tsystem", ["mimo_ss1"], indirect=True)
+    def test_step_response_return(self, tsystem):
+        """Verify continuous and discrete time use same return conventions."""
+        sysc = tsystem.sys
         sysd = c2d(sysc, 1)            # discrete time system
         Tvec = np.linspace(0, 10, 11)  # make sure to use integer times 0..10
         Tc, youtc = step_response(sysc, Tvec, input=0)
@@ -443,10 +345,9 @@ def test_step_response_return(self, mimo_ss1):
         np.testing.assert_array_equal(Tc.shape, Td.shape)
         np.testing.assert_array_equal(youtc.shape, youtd.shape)
 
-
     @pytest.mark.parametrize("dt", [0, 1], ids=["continuous", "discrete"])
     def test_step_nostates(self, dt):
-        """Constant system, continuous and discrete time
+        """Constant system, continuous and discrete time.
 
         gh-374 "Bug in step_response()"
         """
@@ -524,7 +425,7 @@ def test_step_info(self, tsystem, systype, time_2d, yfinal):
     @pytest.mark.parametrize(
         "tsystem",
         ['mimo_ss_step_matlab',
-         pytest.param('mimo_tf_step', marks=slycotonly)],
+         pytest.param('mimo_tf_step_info', marks=slycotonly)],
         indirect=["tsystem"])
     def test_step_info_mimo(self, tsystem, systype, yfinal):
         """Test step info for MIMO systems."""
@@ -560,13 +461,15 @@ def test_step_info_invalid(self):
         with pytest.raises(ValueError, match="matching time vector"):
             step_info(np.ones((2, 2, 15)))                   # no time vector
 
-    def test_step_pole_cancellation(self, pole_cancellation,
-                                    no_pole_cancellation):
+    @pytest.mark.parametrize("tsystem",
+                             [("no_pole_cancellation", "pole_cancellation")],
+                             indirect=True)
+    def test_step_pole_cancellation(self, tsystem):
         # confirm that pole-zero cancellation doesn't perturb results
         # https://github.com/python-control/python-control/issues/440
-        step_info_no_cancellation = step_info(no_pole_cancellation)
-        step_info_cancellation = step_info(pole_cancellation)
-        self.assert_step_info_match(no_pole_cancellation,
+        step_info_no_cancellation = step_info(tsystem[0].sys)
+        step_info_cancellation = step_info(tsystem[1].sys)
+        self.assert_step_info_match(tsystem[0].sys,
                                     step_info_no_cancellation,
                                     step_info_cancellation)
 
@@ -579,7 +482,7 @@ def test_step_pole_cancellation(self, pole_cancellation,
          ("siso_dtf0", {})],
         indirect=["tsystem"])
     def test_impulse_response_siso(self, tsystem, kwargs):
-        """Test impulse response of SISO systems"""
+        """Test impulse response of SISO systems."""
         sys = tsystem.sys
         t = tsystem.t
         yref = tsystem.yimpulse
@@ -590,12 +493,13 @@ def test_impulse_response_siso(self, tsystem, kwargs):
         np.testing.assert_array_almost_equal(tout, t)
         np.testing.assert_array_almost_equal(yout, yref, decimal=4)
 
-    def test_impulse_response_mimo(self, mimo_ss2):
-        """"Test impulse response of MIMO systems"""
-        sys = mimo_ss2.sys
-        t = mimo_ss2.t
+    @pytest.mark.parametrize("tsystem", ["mimo_ss2"], indirect=True)
+    def test_impulse_response_mimo(self, tsystem):
+        """"Test impulse response of MIMO systems."""
+        sys = tsystem.sys
+        t = tsystem.t
 
-        yref = mimo_ss2.yimpulse
+        yref = tsystem.yimpulse
         _t, y_00 = impulse_response(sys, T=t, input=0, output=0)
         np.testing.assert_array_almost_equal(y_00, yref, decimal=4)
         _t, y_11 = impulse_response(sys, T=t, input=1, output=1)
@@ -608,19 +512,21 @@ def test_impulse_response_mimo(self, mimo_ss2):
 
     @pytest.mark.skipif(StrictVersion(sp.__version__) < "1.3",
                         reason="requires SciPy 1.3 or greater")
-    def test_discrete_time_impulse(self, siso_tf1):
+    @pytest.mark.parametrize("tsystem", ["siso_tf1"], indirect=True)
+    def test_discrete_time_impulse(self, tsystem):
         # discrete time impulse sampled version should match cont time
         dt = 0.1
         t = np.arange(0, 3, dt)
-        sys = siso_tf1.sys
+        sys = tsystem.sys
         sysdt = sys.sample(dt, 'impulse')
         np.testing.assert_array_almost_equal(impulse_response(sys, t)[1],
                                              impulse_response(sysdt, t)[1])
 
-    def test_impulse_response_warnD(self, siso_ss1):
+    @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)
+    def test_impulse_response_warnD(self, tsystem):
         """Test warning about direct feedthrough"""
         with pytest.warns(UserWarning, match="System has direct feedthrough"):
-            _ = impulse_response(siso_ss1.sys, siso_ss1.t)
+            _ = impulse_response(tsystem.sys, tsystem.t)
 
     @pytest.mark.parametrize(
         "kwargs",
@@ -630,12 +536,13 @@ def test_impulse_response_warnD(self, siso_ss1):
          {'X0': np.array([[0.5], [1]])},
          {'X0': np.array([0.5, 1]), 'return_x': True},
          ])
-    def test_initial_response(self, siso_ss1, kwargs):
+    @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)
+    def test_initial_response(self, tsystem, kwargs):
         """Test initial response of SISO system"""
-        sys = siso_ss1.sys
-        t = siso_ss1.t
+        sys = tsystem.sys
+        t = tsystem.t
         x0 = kwargs.get('X0', 0)
-        yref = siso_ss1.yinitial if np.any(x0) else np.zeros_like(t)
+        yref = tsystem.yinitial if np.any(x0) else np.zeros_like(t)
 
         out = initial_response(sys, T=t, **kwargs)
         tout, yout = out[:2]
@@ -643,12 +550,13 @@ def test_initial_response(self, siso_ss1, kwargs):
         np.testing.assert_array_almost_equal(tout, t)
         np.testing.assert_array_almost_equal(yout, yref, decimal=4)
 
-    def test_initial_response_mimo(self, mimo_ss1):
+    @pytest.mark.parametrize("tsystem", ["mimo_ss1"], indirect=True)
+    def test_initial_response_mimo(self, tsystem):
         """Test initial response of MIMO system"""
-        sys = mimo_ss1.sys
-        t = mimo_ss1.t
+        sys = tsystem.sys
+        t = tsystem.t
         x0 = np.array([[.5], [1.], [.5], [1.]])
-        yref = mimo_ss1.yinitial
+        yref = tsystem.yinitial
         yref_notrim = np.broadcast_to(yref, (2, len(t)))
 
         _t, y_00 = initial_response(sys, T=t, X0=x0, input=0, output=0)
@@ -676,12 +584,13 @@ def test_forced_response_step(self, tsystem):
                              [np.zeros((10,), dtype=float),
                               0]  # special algorithm
                              )
-    def test_forced_response_initial(self, siso_ss1, u):
+    @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)
+    def test_forced_response_initial(self, tsystem, u):
         """Test forced response of SISO system as intitial response"""
-        sys = siso_ss1.sys
-        t = siso_ss1.t
+        sys = tsystem.sys
+        t = tsystem.t
         x0 = np.array([[.5], [1.]])
-        yref = siso_ss1.yinitial
+        yref = tsystem.yinitial
 
         tout, yout = forced_response(sys, t, u, X0=x0)
         np.testing.assert_array_almost_equal(tout, t)
@@ -741,7 +650,11 @@ def test_forced_response_legacy(self):
                       id="dt=0.2, no U"),
          pytest.param("siso_ss2_dtnone",
                       {'U': np.ones(10,)}, 'ystep',
-                      id="dt=None, no U")],
+                      id="dt=None, no U"),
+         pytest.param("siso_dtf3",
+                      {'U': np.ones(10,)}, 'ystep',
+                      id="dt with rounding error"),
+         ],
         indirect=["tsystem"])
     def test_forced_response_T_U(self, tsystem, fr_kwargs, refattr):
         """Test documented forced_response behavior for parameters T and U."""
@@ -749,7 +662,8 @@ def test_forced_response_T_U(self, tsystem, fr_kwargs, refattr):
         np.testing.assert_allclose(t, tsystem.t)
         np.testing.assert_allclose(y, getattr(tsystem, refattr), rtol=1e-3)
 
-    def test_forced_response_invalid(self, siso_ss1, siso_dss2):
+    @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)
+    def test_forced_response_invalid_c(self, tsystem):
         """Test invalid parameters."""
         with pytest.raises(TypeError,
                            match="StateSpace.*or.*TransferFunction"):
@@ -757,28 +671,29 @@ def test_forced_response_invalid(self, siso_ss1, siso_dss2):
 
         # ctime
         with pytest.raises(ValueError, match="T.*is mandatory for continuous"):
-            forced_response(siso_ss1.sys)
+            forced_response(tsystem.sys)
         with pytest.raises(ValueError, match="time values must be equally "
                                              "spaced"):
-            forced_response(siso_ss1.sys, [0, 0.1, 0.12, 0.4])
+            forced_response(tsystem.sys, [0, 0.1, 0.12, 0.4])
 
-        # dtime with sys.dt > 0
+    @pytest.mark.parametrize("tsystem", ["siso_dss2"], indirect=True)
+    def test_forced_response_invalid_d(self, tsystem):
+        """Test invalid parameters dtime with sys.dt > 0."""
         with pytest.raises(ValueError, match="can't both be zero"):
-            forced_response(siso_dss2.sys)
+            forced_response(tsystem.sys)
         with pytest.raises(ValueError, match="must have same elements"):
-            forced_response(siso_dss2.sys,
-                            T=siso_dss2.t, U=np.random.randn(1, 12))
+            forced_response(tsystem.sys,
+                            T=tsystem.t, U=np.random.randn(1, 12))
         with pytest.raises(ValueError, match="must have same elements"):
-            forced_response(siso_dss2.sys,
-                            T=siso_dss2.t, U=np.random.randn(12))
+            forced_response(tsystem.sys,
+                            T=tsystem.t, U=np.random.randn(12))
         with pytest.raises(ValueError, match="must match sampling time"):
-            forced_response(siso_dss2.sys, T=siso_dss2.t*0.9)
+            forced_response(tsystem.sys, T=tsystem.t*0.9)
         with pytest.raises(ValueError, match="must be multiples of "
                                              "sampling time"):
-            forced_response(siso_dss2.sys, T=siso_dss2.t*1.1)
+            forced_response(tsystem.sys, T=tsystem.t*1.1)
         # but this is ok
-        forced_response(siso_dss2.sys, T=siso_dss2.t*2)
-
+        forced_response(tsystem.sys, T=tsystem.t*2)
 
     @pytest.mark.parametrize("u, x0, xtrue",
                              [(np.zeros((10,)),
@@ -970,14 +885,15 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout):
                 fun(sys, **kw)
 
     @pytest.mark.parametrize("squeeze", [None, True, False])
-    def test_time_vector_interpolation(self, siso_dtf2, squeeze):
-        """Test time vector handling in case of interpolation
+    @pytest.mark.parametrize("tsystem", ["siso_dtf2"], indirect=True)
+    def test_time_vector_interpolation(self, tsystem, squeeze):
+        """Test time vector handling in case of interpolation.
 
         Interpolation of the input (to match scipy.signal.dlsim)
 
         gh-239, gh-295
         """
-        sys = siso_dtf2.sys
+        sys = tsystem.sys
         t = np.arange(0, 10, 1.)
         u = np.sin(t)
         x0 = 0
@@ -993,7 +909,8 @@ def test_time_vector_interpolation(self, siso_dtf2, squeeze):
             assert yout.shape == tout.shape
         assert np.allclose(tout[1:] - tout[:-1], sys.dt)
 
-    def test_discrete_time_steps(self, siso_dtf2):
+    @pytest.mark.parametrize("tsystem", ["siso_dtf2"], indirect=True)
+    def test_discrete_time_steps(self, tsystem):
         """Make sure rounding errors in sample time are handled properly
 
         These tests play around with the input time vector to make sure that
@@ -1001,7 +918,7 @@ def test_discrete_time_steps(self, siso_dtf2):
 
         gh-332
         """
-        sys = siso_dtf2.sys
+        sys = tsystem.sys
 
         # Set up a time range and simulate
         T = np.arange(0, 100, 0.2)
@@ -1033,10 +950,11 @@ def test_discrete_time_steps(self, siso_dtf2):
         with pytest.raises(ValueError):
             step_response(sys, T)
 
-    def test_time_series_data_convention_2D(self, siso_ss1):
+    @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)
+    def test_time_series_data_convention_2D(self, tsystem):
         """Allow input time as 2D array (output should be 1D)"""
         tin = np.array(np.linspace(0, 10, 100), ndmin=2)
-        t, y = step_response(siso_ss1.sys, tin)
+        t, y = step_response(tsystem.sys, tin)
         assert isinstance(t, np.ndarray) and not isinstance(t, np.matrix)
         assert t.ndim == 1
         assert y.ndim == 1  # SISO returns "scalar" output

From a72441b9f629b0f0de03fe5a740128a4df15cc57 Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Fri, 26 Mar 2021 23:48:52 +0100
Subject: [PATCH 08/12] specify tolerance and add comment for markov test

---
 control/tests/modelsimp_test.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
index 70607419e..fd474f9d0 100644
--- a/control/tests/modelsimp_test.py
+++ b/control/tests/modelsimp_test.py
@@ -65,8 +65,6 @@ def testMarkovSignature(self, matarrayout, matarrayin):
             markov(Y, U, m)
 
     # Make sure markov() returns the right answer
-    # forced response can return wrong shape until gh-488 is merged
-    @pytest.mark.xfail
     @pytest.mark.parametrize("k, m, n",
                              [(2, 2, 2),
                               (2, 5, 5),
@@ -81,7 +79,7 @@ def testMarkovResults(self, k, m, n):
         # m = number of Markov parameters
         # n = size of the data vector
         #
-        # Values should match exactly for n = m, otherewise you get a
+        # Values *should* match exactly for n = m, otherewise you get a
         # close match but errors due to the assumption that C A^k B =
         # 0 for k > m-2 (see modelsimp.py).
         #
@@ -108,7 +106,10 @@ def testMarkovResults(self, k, m, n):
         Mcomp = markov(Y, U, m)
 
         # Compare to results from markov()
-        np.testing.assert_array_almost_equal(Mtrue, Mcomp)
+        # experimentally determined probability to get non matching results
+        # with rtot=1e-6 and atol=1e-8 due to numerical errors
+        # for k=5, m=n=10: 0.015 %
+        np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
 
     def testModredMatchDC(self, matarrayin):
         #balanced realization computed in matlab for the transfer function:
@@ -219,4 +220,3 @@ def testBalredMatchDC(self, matarrayin):
         np.testing.assert_array_almost_equal(rsys.B, Brtrue, decimal=4)
         np.testing.assert_array_almost_equal(rsys.C, Crtrue, decimal=4)
         np.testing.assert_array_almost_equal(rsys.D, Drtrue, decimal=4)
-

From 177290c9d157aa7fd5b9fb2ab35efafc1aba6424 Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Sat, 27 Mar 2021 00:40:46 +0100
Subject: [PATCH 09/12] reactivate the the fast forced_response(U=0) algorithm
 and test it

---
 control/tests/timeresp_test.py | 20 +++++++++++++-------
 control/timeresp.py            |  3 ++-
 2 files changed, 15 insertions(+), 8 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index b4a3598e5..73048a2e7 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -43,6 +43,7 @@ def tsystem(self, request):
         siso_ss1.ystep = np.array([9., 17.6457, 24.7072, 30.4855, 35.2234,
                                    39.1165, 42.3227, 44.9694, 47.1599,
                                    48.9776])
+        # X0 = [0.5, 1]
         siso_ss1.yinitial = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118,
                                       1.9092, 1.1508, 0.5833, 0.1645, -0.1391])
         ss1 = siso_ss1.sys
@@ -127,6 +128,7 @@ def tsystem(self, request):
 
         siso_dss1 = copy(siso_dtf1)
         siso_dss1.sys = tf2ss(siso_dtf1.sys)
+        siso_dss1.yinitial = np.array([-1., -0.5, 0.75, -0.625, 0.4375])
 
         siso_dss2 = copy(siso_dtf2)
         siso_dss2.sys = tf2ss(siso_dtf2.sys)
@@ -641,19 +643,23 @@ def test_forced_response_legacy(self):
         [pytest.param("siso_ss1",
                       {'X0': [0.5, 1], 'T':  np.linspace(0, 1, 10)},
                       'yinitial',
-                      id="ctime no T"),
+                      id="ctime no U"),
+         pytest.param("siso_dss1",
+                      {'T': np.arange(0, 5, 1,),
+                       'X0': [0.5, 1]}, 'yinitial',
+                      id="dt=True, no U"),
          pytest.param("siso_dtf1",
                       {'U': np.ones(5,)}, 'ystep',
-                      id="dt=True, no U"),
+                      id="dt=True, no T"),
          pytest.param("siso_dtf2",
                       {'U': np.ones(25,)}, 'ystep',
-                      id="dt=0.2, no U"),
+                      id="dt=0.2, no T"),
          pytest.param("siso_ss2_dtnone",
                       {'U': np.ones(10,)}, 'ystep',
-                      id="dt=None, no U"),
+                      id="dt=None, no T"),
          pytest.param("siso_dtf3",
                       {'U': np.ones(10,)}, 'ystep',
-                      id="dt with rounding error"),
+                      id="dt with rounding error, no T"),
          ],
         indirect=["tsystem"])
     def test_forced_response_T_U(self, tsystem, fr_kwargs, refattr):
@@ -668,13 +674,13 @@ def test_forced_response_invalid_c(self, tsystem):
         with pytest.raises(TypeError,
                            match="StateSpace.*or.*TransferFunction"):
             forced_response("not a system")
-
-        # ctime
         with pytest.raises(ValueError, match="T.*is mandatory for continuous"):
             forced_response(tsystem.sys)
         with pytest.raises(ValueError, match="time values must be equally "
                                              "spaced"):
             forced_response(tsystem.sys, [0, 0.1, 0.12, 0.4])
+        with pytest.raises(ValueError, match="must start with 0"):
+            forced_response(tsystem.sys, [1, 1.1, 1.2, 1.3])
 
     @pytest.mark.parametrize("tsystem", ["siso_dss2"], indirect=True)
     def test_forced_response_invalid_d(self, tsystem):
diff --git a/control/timeresp.py b/control/timeresp.py
index 8c0d59718..a30df18b5 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -379,7 +379,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         dot = np.dot  # Faster and shorter code
 
         # Faster algorithm if U is zero
-        if U is None or (isinstance(U, (int, float)) and U == 0):
+        # (if not None, it was converted to array above)
+        if U is None or np.all(U == 0):
             # Solve using matrix exponential
             expAdt = sp.linalg.expm(A * dt)
             for i in range(1, n_steps):

From ce51d34e1be57596c26a5661e0adc2bf728b25ac Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Sat, 27 Mar 2021 01:03:11 +0100
Subject: [PATCH 10/12] cover TF initial condition warning in forced_response

---
 control/tests/timeresp_test.py | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 73048a2e7..9a69756dd 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -43,7 +43,7 @@ def tsystem(self, request):
         siso_ss1.ystep = np.array([9., 17.6457, 24.7072, 30.4855, 35.2234,
                                    39.1165, 42.3227, 44.9694, 47.1599,
                                    48.9776])
-        # X0 = [0.5, 1]
+        siso_ss1.X0 = np.array([[.5], [1.]])
         siso_ss1.yinitial = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118,
                                       1.9092, 1.1508, 0.5833, 0.1645, -0.1391])
         ss1 = siso_ss1.sys
@@ -586,17 +586,23 @@ def test_forced_response_step(self, tsystem):
                              [np.zeros((10,), dtype=float),
                               0]  # special algorithm
                              )
-    @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)
+    @pytest.mark.parametrize("tsystem", ["siso_ss1", "siso_tf2"],
+                             indirect=True)
     def test_forced_response_initial(self, tsystem, u):
-        """Test forced response of SISO system as intitial response"""
+        """Test forced response of SISO system as intitial response."""
         sys = tsystem.sys
         t = tsystem.t
-        x0 = np.array([[.5], [1.]])
+        x0 = tsystem.X0
         yref = tsystem.yinitial
 
-        tout, yout = forced_response(sys, t, u, X0=x0)
-        np.testing.assert_array_almost_equal(tout, t)
-        np.testing.assert_array_almost_equal(yout, yref, decimal=4)
+        if isinstance(sys, StateSpace):
+            tout, yout = forced_response(sys, t, u, X0=x0)
+            np.testing.assert_array_almost_equal(tout, t)
+            np.testing.assert_array_almost_equal(yout, yref, decimal=4)
+        else:
+            with pytest.warns(UserWarning, match="Non-zero initial condition "
+                              "given for transfer function"):
+                tout, yout = forced_response(sys, t, u, X0=x0)
 
     @pytest.mark.parametrize("tsystem, useT",
                              [("mimo_ss1", True),

From 44041692eac186a7ac9df31ae98e9f19e7a69067 Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Sat, 27 Mar 2021 01:20:40 +0100
Subject: [PATCH 11/12] reallow non zero timevector start

---
 control/tests/timeresp_test.py |  2 --
 control/timeresp.py            | 23 +++++++++++++----------
 2 files changed, 13 insertions(+), 12 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 9a69756dd..e6261533d 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -685,8 +685,6 @@ def test_forced_response_invalid_c(self, tsystem):
         with pytest.raises(ValueError, match="time values must be equally "
                                              "spaced"):
             forced_response(tsystem.sys, [0, 0.1, 0.12, 0.4])
-        with pytest.raises(ValueError, match="must start with 0"):
-            forced_response(tsystem.sys, [1, 1.1, 1.2, 1.3])
 
     @pytest.mark.parametrize("tsystem", ["siso_dss2"], indirect=True)
     def test_forced_response_invalid_d(self, tsystem):
diff --git a/control/timeresp.py b/control/timeresp.py
index a30df18b5..989a832cb 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -209,8 +209,7 @@ 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
-        and start with 0.
+        Time steps at which the input is defined; values must be evenly spaced.
 
         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
@@ -360,10 +359,10 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     n_steps = T.shape[0]            # number of simulation steps
 
     # equally spaced also implies strictly monotonic increase,
-    dt = T[-1] / (n_steps - 1)
+    dt = (T[-1] - T[0]) / (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.")
+        raise ValueError("Parameter ``T``: 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)],
@@ -426,6 +425,10 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
     else:
         # Discrete type system => use SciPy signal processing toolbox
+
+        # sp.signal.dlsim assumes T[0] == 0
+        spT = T - T[0]
+
         if sys.dt is not True and sys.dt is not None:
             # Make sure that the time increment is a multiple of sampling time
 
@@ -446,10 +449,10 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
             # 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
+            scipy_out_samples = int(np.floor(spT[-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))
+                spT[-1] = spT[-1] * (n_steps / (spT[-1] / sys_dt + 1))
 
         else:
             sys_dt = dt         # For unspecified sampling time, use time incr
@@ -459,7 +462,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
         # 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)
+        tout, yout, xout = sp.signal.dlsim(dsys, np.transpose(U), spT, X0)
+        tout = tout + T[0]
 
         if not interpolate:
             # If dt is different from sys.dt, resample the output
@@ -472,8 +476,7 @@ 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[:n_steps], yout[:, :n_steps],
-                                  xout[:, :n_steps], transpose=transpose,
+    return _process_time_response(sys, tout, yout, xout, transpose=transpose,
                                   return_x=return_x, squeeze=squeeze)
 
 

From a8b72f5f6ac5317b956f305514dd7fce1628b12f Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Sat, 27 Mar 2021 01:42:46 +0100
Subject: [PATCH 12/12] avoid different realizations with and without slycot

---
 control/tests/timeresp_test.py | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index e6261533d..a91507a83 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -126,9 +126,18 @@ def tsystem(self, request):
             [ 0.    , -0.1097, -0.1902, -0.2438, -0.2729,
              -0.2799, -0.2674, -0.2377, -0.1934, -0.1368])
 
+        """dtf1 converted statically, because Slycot and Scipy produce
+        different realizations, wich means different initial condtions,"""
         siso_dss1 = copy(siso_dtf1)
-        siso_dss1.sys = tf2ss(siso_dtf1.sys)
-        siso_dss1.yinitial = np.array([-1., -0.5, 0.75, -0.625, 0.4375])
+        siso_dss1.sys = StateSpace([[-1., -0.25],
+                                    [ 1.,  0.]],
+                                   [[1.],
+                                    [0.]],
+                                   [[0., 1.]],
+                                   [[0.]],
+                                   True)
+        siso_dss1.X0 = [0.5, 1.]
+        siso_dss1.yinitial = np.array([1., 0.5, -0.75, 0.625, -0.4375])
 
         siso_dss2 = copy(siso_dtf2)
         siso_dss2.sys = tf2ss(siso_dtf2.sys)
@@ -647,12 +656,10 @@ def test_forced_response_legacy(self):
     @pytest.mark.parametrize(
         "tsystem, fr_kwargs, refattr",
         [pytest.param("siso_ss1",
-                      {'X0': [0.5, 1], 'T':  np.linspace(0, 1, 10)},
-                      'yinitial',
+                      {'T':  np.linspace(0, 1, 10)}, 'yinitial',
                       id="ctime no U"),
          pytest.param("siso_dss1",
-                      {'T': np.arange(0, 5, 1,),
-                       'X0': [0.5, 1]}, 'yinitial',
+                      {'T': np.arange(0, 5, 1,)}, 'yinitial',
                       id="dt=True, no U"),
          pytest.param("siso_dtf1",
                       {'U': np.ones(5,)}, 'ystep',
@@ -670,6 +677,8 @@ def test_forced_response_legacy(self):
         indirect=["tsystem"])
     def test_forced_response_T_U(self, tsystem, fr_kwargs, refattr):
         """Test documented forced_response behavior for parameters T and U."""
+        if refattr == 'yinitial':
+            fr_kwargs['X0'] = tsystem.X0
         t, y = forced_response(tsystem.sys, **fr_kwargs)
         np.testing.assert_allclose(t, tsystem.t)
         np.testing.assert_allclose(y, getattr(tsystem, refattr), rtol=1e-3)