diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 10001dd96..20d715483 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -10,7 +10,6 @@ import unittest import numpy as np -# import scipy as sp from control.timeresp import * from control.statesp import * from control.xferfcn import TransferFunction, _convert_to_transfer_function @@ -309,7 +308,7 @@ def check(u, x0, xtrue): def test_discrete_initial(self): h1 = TransferFunction([1.], [1., 0.], 1.) t, yout = impulse_response(h1, np.arange(4)) - np.testing.assert_array_equal(yout[0], [0., 1., 0., 0.]) + np.testing.assert_array_equal(yout, [0., 1., 0., 0.]) @unittest.skipIf(not slycot_check(), "slycot not installed") def test_step_robustness(self): @@ -345,27 +344,32 @@ def test_time_vector(self): # No timebase in system => output should match input # # Initial response - tout, yout = initial_response(self.siso_dtf1, Tin2, siso_x0) + tout, yout = initial_response(self.siso_dtf1, Tin2, siso_x0, + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin2) # Impulse response - tout, yout = impulse_response(self.siso_dtf1, Tin2) + tout, yout = impulse_response(self.siso_dtf1, Tin2, + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin2) # Step response - tout, yout = step_response(self.siso_dtf1, Tin2) + tout, yout = step_response(self.siso_dtf1, Tin2, + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin2) # Forced response with specified time vector - tout, yout, xout = forced_response(self.siso_dtf1, Tin2, np.sin(Tin2)) + tout, yout, xout = forced_response(self.siso_dtf1, Tin2, np.sin(Tin2), + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin2) # Forced response with no time vector, no sample time (should use 1) - tout, yout, xout = forced_response(self.siso_dtf1, None, np.sin(Tin1)) + tout, yout, xout = forced_response(self.siso_dtf1, None, np.sin(Tin1), + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin1) @@ -380,49 +384,58 @@ def test_time_vector(self): # Matching timebase in system => output should match input # # Initial response - tout, yout = initial_response(self.siso_dtf2, Tin2, siso_x0) + tout, yout = initial_response(self.siso_dtf2, Tin2, siso_x0, + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin2) # Impulse response - tout, yout = impulse_response(self.siso_dtf2, Tin2) + tout, yout = impulse_response(self.siso_dtf2, Tin2, + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin2) # Step response - tout, yout = step_response(self.siso_dtf2, Tin2) + tout, yout = step_response(self.siso_dtf2, Tin2, + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin2) # Forced response - tout, yout, xout = forced_response(self.siso_dtf2, Tin2, np.sin(Tin2)) + tout, yout, xout = forced_response(self.siso_dtf2, Tin2, np.sin(Tin2), + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin2) # Forced response with no time vector, use sample time - tout, yout, xout = forced_response(self.siso_dtf2, None, np.sin(Tin2)) + tout, yout, xout = forced_response(self.siso_dtf2, None, np.sin(Tin2), + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin2) # Compatible timebase in system => output should match input # # Initial response - tout, yout = initial_response(self.siso_dtf2, Tin1, siso_x0) + tout, yout = initial_response(self.siso_dtf2, Tin1, siso_x0, + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin1) # Impulse response - tout, yout = impulse_response(self.siso_dtf2, Tin1) + tout, yout = impulse_response(self.siso_dtf2, Tin1, + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin1) # Step response - tout, yout = step_response(self.siso_dtf2, Tin1) + tout, yout = step_response(self.siso_dtf2, Tin1, + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin1) # Forced response - tout, yout, xout = forced_response(self.siso_dtf2, Tin1, np.sin(Tin1)) + tout, yout, xout = forced_response(self.siso_dtf2, Tin1, np.sin(Tin1), + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) np.testing.assert_array_equal(tout, Tin1) @@ -431,7 +444,8 @@ def test_time_vector(self): # # Initial response tout, yout, xout = forced_response(self.siso_dtf2, Tin1, - np.sin(Tin1), interpolate=True) + np.sin(Tin1), interpolate=True, + squeeze=False) self.assertEqual(np.shape(tout), np.shape(yout[0,:])) self.assertTrue(np.allclose(tout[1:] - tout[:-1], self.siso_dtf2.dt)) @@ -442,9 +456,62 @@ def test_time_vector(self): # # Initial response with self.assertRaises(Exception) as context: - tout, yout = initial_response(self.siso_dtf2, Tin3, siso_x0) + tout, yout = initial_response(self.siso_dtf2, Tin3, siso_x0, + squeeze=False) self.assertTrue(isinstance(context.exception, ValueError)) + def test_time_series_data_convention(self): + """Make sure time series data matches documentation conventions""" + # SISO continuous time + t, y = step_response(self.siso_ss1) + self.assertTrue(isinstance(t, np.ndarray) + and not isinstance(t, np.matrix)) + self.assertTrue(len(t.shape) == 1) + self.assertTrue(len(y.shape) == 1) # SISO returns "scalar" output + self.assertTrue(len(t) == len(y)) # Allows direct plotting of output + + # SISO discrete time + t, y = step_response(self.siso_dss1) + self.assertTrue(isinstance(t, np.ndarray) + and not isinstance(t, np.matrix)) + self.assertTrue(len(t.shape) == 1) + self.assertTrue(len(y.shape) == 1) # SISO returns "scalar" output + self.assertTrue(len(t) == len(y)) # Allows direct plotting of output + + # MIMO continuous time + tin = np.linspace(0, 10, 100) + uin = [np.sin(tin), np.cos(tin)] + t, y, x = forced_response(self.mimo_ss1, tin, uin) + self.assertTrue(isinstance(t, np.ndarray) + and not isinstance(t, np.matrix)) + self.assertTrue(len(t.shape) == 1) + self.assertTrue(len(y[0].shape) == 1) + self.assertTrue(len(y[1].shape) == 1) + self.assertTrue(len(t) == len(y[0])) + self.assertTrue(len(t) == len(y[1])) + + # MIMO discrete time + tin = np.linspace(0, 10, 100) + uin = [np.sin(tin), np.cos(tin)] + t, y, x = forced_response(self.mimo_dss1, tin, uin) + self.assertTrue(isinstance(t, np.ndarray) + and not isinstance(t, np.matrix)) + self.assertTrue(len(t.shape) == 1) + self.assertTrue(len(y[0].shape) == 1) + self.assertTrue(len(y[1].shape) == 1) + self.assertTrue(len(t) == len(y[0])) + self.assertTrue(len(t) == len(y[1])) + + # Allow input time as 2D array (output should be 1D) + tin = np.array(np.linspace(0, 10, 100), ndmin=2) + t, y = step_response(self.siso_ss1, tin) + self.assertTrue(isinstance(t, np.ndarray) + and not isinstance(t, np.matrix)) + self.assertTrue(len(t.shape) == 1) + self.assertTrue(len(y.shape) == 1) # SISO returns "scalar" output + self.assertTrue(len(t) == len(y)) # Allows direct plotting of output + + def suite(): return unittest.TestLoader().loadTestsFromTestCase(TestTimeresp) diff --git a/control/timeresp.py b/control/timeresp.py index 14e7072d2..2fc794cb5 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -61,6 +61,7 @@ __all__ = ['forced_response', 'step_response', 'step_info', 'initial_response', 'impulse_response'] + # Helper function for checking array-like parameters def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False, transpose=False): @@ -170,7 +171,7 @@ def shape_matches(s_legal, s_actual): # Forced response of a linear system def forced_response(sys, T=None, U=0., X0=0., transpose=False, - interpolate=False): + interpolate=False, squeeze=True): """Simulate the output of a linear system. As a convenience for parameters `U`, `X0`: @@ -192,22 +193,28 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, Input array giving input at each time `T` (default = 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. + algorithm is faster than the general algorithm, which is used + otherwise. X0: array-like or number, optional Initial condition (default = 0). - transpose: bool + transpose: bool, optional (default=False) If True, transpose all input and output arrays (for backward compatibility with MATLAB and scipy.signal.lsim) - interpolate:bool + interpolate: bool, optional (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 (default = False). + squeeze: bool, optional (default=True) + If True, remove single-dimensional entries from the shape of + the output. For single output systems, this converts the + output response to a 1D array. + Returns ------- T: array @@ -226,6 +233,7 @@ 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`. + """ if not isinstance(sys, LTI): raise TypeError('Parameter ``sys``: must be a ``LTI`` object. ' @@ -245,7 +253,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, 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 - T = np.array(range(len(U))) * (1 if sys.dt == True else sys.dt) + T = np.array(range(len(U))) * (1 if sys.dt is True else sys.dt) else: # Make sure the input vector and time vector have same length # TODO: allow interpolation of the input vector @@ -264,7 +272,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, transpose=transpose) dt = T[1] - T[0] if not np.allclose(T[1:] - T[:-1], dt): - raise ValueError('Parameter ``T``: time values must be equally spaced.') + raise ValueError("Parameter ``T``: time values must be " + "equally spaced.") n_steps = len(T) # number of simulation steps # create X0 if not given, test if X0 has correct shape @@ -323,14 +332,11 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, xout[:, i] = (dot(Ad, xout[:, i-1]) + dot(Bd0, U[:, i-1]) + dot(Bd1, U[:, i])) yout = dot(C, xout) + dot(D, U) - tout = T - yout = squeeze(yout) - xout = squeeze(xout) else: # Discrete type system => use SciPy signal processing toolbox - if (sys.dt != True): + if sys.dt is not True: # Make sure that the time increment is a multiple of sampling time # First make sure that time increment is bigger than sampling time @@ -341,7 +347,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, # 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)): - raise ValueError("Time steps ``T`` must be multiples of " \ + raise ValueError("Time steps ``T`` must be multiples of " "sampling time") else: sys.dt = dt # For unspecified sampling time, use time incr @@ -350,28 +356,34 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, dsys = (A, B, C, D, sys.dt) # Use signal processing toolbox for the discrete time simulation - # Transpose the input to match toolbox convention + # Transpose the input to match toolbox convention tout, yout, xout = sp.signal.dlsim(dsys, np.transpose(U), T, X0) if not interpolate: # If dt is different from sys.dt, resample the output inc = int(round(dt / sys.dt)) tout = T # Return exact list of time steps - yout = yout[::inc,:] - xout = xout[::inc,:] + yout = yout[::inc, :] + xout = xout[::inc, :] # Transpose the output and state vectors to match local convention xout = sp.transpose(xout) yout = sp.transpose(yout) + # Get rid of unneeded dimensions + if squeeze: + yout = np.squeeze(yout) + xout = np.squeeze(xout) + # See if we need to transpose the data back into MATLAB form - if (transpose): + if transpose: tout = np.transpose(tout) yout = np.transpose(yout) xout = np.transpose(xout) return tout, yout, xout + def _get_ss_simo(sys, input=None, output=None): """Return a SISO or SIMO state-space version of sys @@ -390,8 +402,9 @@ def _get_ss_simo(sys, input=None, output=None): else: return _mimo2siso(sys_ss, input, output, warn_conversion=warn) + def step_response(sys, T=None, X0=0., input=None, output=None, - transpose=False, return_x=False): + transpose=False, return_x=False, squeeze=True): # pylint: disable=W0622 """Step response of a linear system @@ -430,6 +443,11 @@ def step_response(sys, T=None, X0=0., input=None, output=None, return_x: bool If True, return the state vector (default = False). + squeeze: bool, optional (default=True) + If True, remove single-dimensional entries from the shape of + the output. For single output systems, this converts the + output response to a 1D array. + Returns ------- T: array @@ -460,15 +478,17 @@ def step_response(sys, T=None, X0=0., input=None, output=None, U = np.ones_like(T) - T, yout, xout = forced_response(sys, T, U, X0, - transpose=transpose) + T, yout, xout = forced_response(sys, T, U, X0, transpose=transpose, + squeeze=squeeze) if return_x: return T, yout, xout return T, yout -def step_info(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1,0.9)): + +def step_info(sys, T=None, SettlingTimeThreshold=0.02, + RiseTimeLimits=(0.1, 0.9)): ''' Step response characteristics (Rise time, Settling Time, Peak and others). @@ -561,8 +581,9 @@ def step_info(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1,0.9)) return S + def initial_response(sys, T=None, X0=0., input=0, output=None, - transpose=False, return_x=False): + transpose=False, return_x=False, squeeze=True): # pylint: disable=W0622 """Initial condition response of a linear system @@ -601,6 +622,11 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, return_x: bool If True, return the state vector (default = False). + squeeze: bool, optional (default=True) + If True, remove single-dimensional entries from the shape of + the output. For single output systems, this converts the + output response to a 1D array. + Returns ------- T: array @@ -631,7 +657,8 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T = range(int(np.ceil(max(tvec)))) U = np.zeros_like(T) - T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose) + T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose, + squeeze=squeeze) if return_x: return T, yout, _xout @@ -640,7 +667,7 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, def impulse_response(sys, T=None, X0=0., input=0, output=None, - transpose=False, return_x=False): + transpose=False, return_x=False, squeeze=True): # pylint: disable=W0622 """Impulse response of a linear system @@ -679,6 +706,11 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, return_x: bool If True, return the state vector (default = False). + squeeze: bool, optional (default=True) + If True, remove single-dimensional entries from the shape of + the output. For single output systems, this converts the + output response to a 1D array. + Returns ------- T: array @@ -698,11 +730,13 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, """ sys = _get_ss_simo(sys, input, output) - # System has direct feedthrough, can't simulate impulse response numerically + # 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 infinite ' - 'impulse at ``t=0`` does not appear in the output. \n' - 'Results may be meaningless!') + warnings.warn("System has direct feedthrough: ``D != 0``. The " + "infinite impulse at ``t=0`` does not appear in the " + "output.\n" + "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. @@ -732,9 +766,8 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, new_X0 = X0 U[0] = 1. - T, yout, _xout = forced_response( - sys, T, U, new_X0, - transpose=transpose) + T, yout, _xout = forced_response(sys, T, U, new_X0, transpose=transpose, + squeeze=squeeze) if return_x: return T, yout, _xout diff --git a/doc/conventions.rst b/doc/conventions.rst index 7bdf3c628..1d5fe2029 100644 --- a/doc/conventions.rst +++ b/doc/conventions.rst @@ -162,12 +162,12 @@ The initial conditions are either 1D, or 2D with shape (j, 1):: As all simulation functions return *arrays*, plotting is convenient:: - t, y = step(sys) + t, y = step_response(sys) plot(t, y) The output of a MIMO system can be plotted like this:: - t, y, x = lsim(sys, u, t) + t, y, x = forced_response(sys, u, t) plot(t, y[0], label='y_0') plot(t, y[1], label='y_1')