diff --git a/control/matlab/timeresp.py b/control/matlab/timeresp.py index ffdf9dd6f..647210a9c 100644 --- a/control/matlab/timeresp.py +++ b/control/matlab/timeresp.py @@ -4,7 +4,7 @@ Note that the return arguments are different than in the standard control package. """ -__all__ = ['step', 'impulse', 'initial', 'lsim'] +__all__ = ['step', 'stepinfo', 'impulse', 'initial', 'lsim'] def step(sys, T=None, X0=0., input=0, output=None, return_x=False): ''' @@ -66,6 +66,52 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False): return yout, T +def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1,0.9)): + ''' + Step response characteristics (Rise time, Settling Time, Peak and others). + + Parameters + ---------- + sys: StateSpace, or TransferFunction + LTI system to simulate + + T: array-like object, optional + Time vector (argument is autocomputed if not given) + + SettlingTimeThreshold: float value, optional + Defines the error to compute settling time (default = 0.02) + + RiseTimeLimits: tuple (lower_threshold, upper_theshold) + Defines the lower and upper threshold for RiseTime computation + + Returns + ------- + S: a dictionary containing: + RiseTime: Time from 10% to 90% of the steady-state value. + SettlingTime: Time to enter inside a default error of 2% + SettlingMin: Minimum value after RiseTime + SettlingMax: Maximum value after RiseTime + Overshoot: Percentage of the Peak relative to steady value + Undershoot: Percentage of undershoot + Peak: Absolute peak value + PeakTime: time of the Peak + SteadyStateValue: Steady-state value + + + See Also + -------- + step, lsim, initial, impulse + + Examples + -------- + >>> S = stepinfo(sys, T) + ''' + from ..timeresp import step_info + + S = step_info(sys, T, SettlingTimeThreshold, RiseTimeLimits) + + return S + def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False): ''' Impulse response of a linear system diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index f0372f734..3a6c6c665 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -87,6 +87,65 @@ def test_step_response(self): np.testing.assert_array_equal(Tc.shape, Td.shape) np.testing.assert_array_equal(youtc.shape, youtd.shape) + def test_step_info(self): + # From matlab docs: + sys = TransferFunction([1,5,5],[1,1.65,5,6.5,2]) + Strue = { + 'RiseTime': 3.8456, + 'SettlingTime': 27.9762, + 'SettlingMin': 2.0689, + 'SettlingMax': 2.6873, + 'Overshoot': 7.4915, + 'Undershoot': 0, + 'Peak': 2.6873, + 'PeakTime': 8.0530 + } + + S = step_info(sys) + + # Very arbitrary tolerance because I don't know if the + # response from the MATLAB is really that accurate. + # maybe it is a good idea to change the Strue to match + # but I didn't do it because I don't know if it is + # accurate either... + rtol = 2e-2 + np.testing.assert_allclose( + S.get('RiseTime'), + Strue.get('RiseTime'), + rtol=rtol) + np.testing.assert_allclose( + S.get('SettlingTime'), + Strue.get('SettlingTime'), + rtol=rtol) + np.testing.assert_allclose( + S.get('SettlingMin'), + Strue.get('SettlingMin'), + rtol=rtol) + np.testing.assert_allclose( + S.get('SettlingMax'), + Strue.get('SettlingMax'), + rtol=rtol) + np.testing.assert_allclose( + S.get('Overshoot'), + Strue.get('Overshoot'), + rtol=rtol) + np.testing.assert_allclose( + S.get('Undershoot'), + Strue.get('Undershoot'), + rtol=rtol) + np.testing.assert_allclose( + S.get('Peak'), + Strue.get('Peak'), + rtol=rtol) + np.testing.assert_allclose( + S.get('PeakTime'), + Strue.get('PeakTime'), + rtol=rtol) + np.testing.assert_allclose( + S.get('SteadyStateValue'), + 2.50, + rtol=rtol) + def test_impulse_response(self): # Test SISO system sys = self.siso_ss1 diff --git a/control/timeresp.py b/control/timeresp.py index a0ecf9d93..059e42061 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -58,7 +58,7 @@ from .statesp import _convertToStateSpace, _mimo2simo, _mimo2siso from .lti import isdtime, isctime -__all__ = ['forced_response', 'step_response', 'initial_response', +__all__ = ['forced_response', 'step_response', 'step_info', 'initial_response', 'impulse_response'] # Helper function for checking array-like parameters @@ -432,6 +432,98 @@ def step_response(sys, T=None, X0=0., input=None, output=None, return T, yout +def step_info(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1,0.9)): + ''' + Step response characteristics (Rise time, Settling Time, Peak and others). + + Parameters + ---------- + sys: StateSpace, or TransferFunction + LTI system to simulate + + T: array-like object, optional + Time vector (argument is autocomputed if not given) + + SettlingTimeThreshold: float value, optional + Defines the error to compute settling time (default = 0.02) + + RiseTimeLimits: tuple (lower_threshold, upper_theshold) + Defines the lower and upper threshold for RiseTime computation + + Returns + ------- + S: a dictionary containing: + RiseTime: Time from 10% to 90% of the steady-state value. + SettlingTime: Time to enter inside a default error of 2% + SettlingMin: Minimum value after RiseTime + SettlingMax: Maximum value after RiseTime + Overshoot: Percentage of the Peak relative to steady value + Undershoot: Percentage of undershoot + Peak: Absolute peak value + PeakTime: time of the Peak + SteadyStateValue: Steady-state value + + + See Also + -------- + step, lsim, initial, impulse + + Examples + -------- + >>> info = step_info(sys, T) + ''' + sys = _get_ss_simo(sys) + if T is None: + if isctime(sys): + T = _default_response_times(sys.A, 1000) + else: + # For discrete time, use integers + tvec = _default_response_times(sys.A, 1000) + T = range(int(np.ceil(max(tvec)))) + + T, yout = step_response(sys, T) + + # Steady state value + InfValue = yout[-1] + + # RiseTime + tr_lower_index = (np.where(yout >= RiseTimeLimits[0] * InfValue)[0])[0] + tr_upper_index = (np.where(yout >= RiseTimeLimits[1] * InfValue)[0])[0] + RiseTime = T[tr_upper_index] - T[tr_lower_index] + + # SettlingTime + sup_margin = (1. + SettlingTimeThreshold) * InfValue + inf_margin = (1. - SettlingTimeThreshold) * InfValue + # find Steady State looking for the first point out of specified limits + for i in reversed(range(T.size)): + if((yout[i] <= inf_margin) | (yout[i] >= sup_margin)): + SettlingTime = T[i + 1] + break + + # Peak + PeakIndex = np.abs(yout).argmax() + PeakValue = yout[PeakIndex] + PeakTime = T[PeakIndex] + SettlingMax = (yout).max() + SettlingMin = (yout[tr_upper_index:]).min() + # I'm really not very confident about UnderShoot: + UnderShoot = yout.min() + OverShoot = 100. * (yout.max() - InfValue) / (InfValue - yout[0]) + + # Return as a dictionary + S = { + 'RiseTime': RiseTime, + 'SettlingTime': SettlingTime, + 'SettlingMin': SettlingMin, + 'SettlingMax': SettlingMax, + 'Overshoot': OverShoot, + 'Undershoot': UnderShoot, + 'Peak': PeakValue, + 'PeakTime': PeakTime, + 'SteadyStateValue': InfValue + } + + return S def initial_response(sys, T=None, X0=0., input=0, output=None, transpose=False, return_x=False):