diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 37fcff763..8705f3805 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -193,6 +193,35 @@ def pole_cancellation(self): def no_pole_cancellation(self): return 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) + return TransferFunction(1, [1, 1, 0]) + + @pytest.fixture + def siso_tf_kpos(self): + # SISO under shoot response and positive final value G(s)=(-s+1)/(s²+s+1) + return TransferFunction([-1, 1], [1, 1, 1]) + + @pytest.fixture + def siso_tf_kneg(self): + # SISO under shoot response and negative final value k=-1 G(s)=-(-s+1)/(s²+s+1) + return TransferFunction([1, -1], [1, 1, 1]) + + @pytest.fixture + def tf1_matlab_help(self): + # example from matlab online help https://www.mathworks.com/help/control/ref/stepinfo.html + return TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2]) + + @pytest.fixture + def tf2_matlab_help(self): + A = [[0.68, - 0.34], [0.34, 0.68]] + B = [[0.18], [0.04]] + C = [-1.12, - 1.10] + D = [0.06] + sys = StateSpace(A, B, C, D, 0.2) + return sys @pytest.fixture def tsystem(self, @@ -202,7 +231,9 @@ def tsystem(self, siso_dtf0, siso_dtf1, siso_dtf2, siso_dss1, siso_dss2, mimo_dss1, mimo_dss2, mimo_dtf1, - pole_cancellation, no_pole_cancellation): + pole_cancellation, no_pole_cancellation, siso_tf_type1, + siso_tf_kpos, siso_tf_kneg, tf1_matlab_help, + tf2_matlab_help): systems = {"siso_ss1": siso_ss1, "siso_ss2": siso_ss2, "siso_tf1": siso_tf1, @@ -220,6 +251,11 @@ def tsystem(self, "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, + "tf1_matlab_help": tf1_matlab_help, + "tf2_matlab_help": tf2_matlab_help, } return systems[request.param] @@ -303,6 +339,73 @@ def test_step_info(self): [Strue[k] for k in Sktrue], rtol=rtol) + # tolerance for all parameters could be wrong for some systems + # discrete systems time parameters tolerance could be +/-dt + @pytest.mark.parametrize( + "tsystem, info_true, tolerance", + [("tf1_matlab_help", { + 'RiseTime': 3.8456, + 'SettlingTime': 27.9762, + 'SettlingMin': 2.0689, + 'SettlingMax': 2.6873, + 'Overshoot': 7.4915, + 'Undershoot': 0, + 'Peak': 2.6873, + 'PeakTime': 8.0530, + 'SteadyStateValue': 2.5}, 2e-2), + ("tf2_matlab_help", { + '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}, .2), + ("siso_tf_kpos", { + 'RiseTime': 1.242, + 'SettlingTime': 9.110, + 'SettlingMin': 0.950, + 'SettlingMax': 1.208, + 'Overshoot': 20.840, + 'Undershoot': 27.840, + 'Peak': 1.208, + 'PeakTime': 4.282, + 'SteadyStateValue': 1.0}, 2e-2), + ("siso_tf_kneg", { + 'RiseTime': 1.242, + 'SettlingTime': 9.110, + 'SettlingMin': -1.208, + 'SettlingMax': -0.950, + 'Overshoot': 20.840, + 'Undershoot': 27.840, + 'Peak': 1.208, + 'PeakTime': 4.282, + 'SteadyStateValue': -1.0}, 2e-2), + ("siso_tf_type1", {'RiseTime': np.NaN, + 'SettlingTime': np.NaN, + 'SettlingMin': np.NaN, + 'SettlingMax': np.NaN, + 'Overshoot': np.NaN, + 'Undershoot': np.NaN, + 'Peak': np.Inf, + 'PeakTime': np.Inf, + 'SteadyStateValue': np.NaN}, 2e-2)], + indirect=["tsystem"]) + def test_step_info(self, tsystem, info_true, tolerance): + """ """ + info = step_info(tsystem) + + info_true_sorted = sorted(info_true.keys()) + info_sorted = sorted(info.keys()) + + assert info_sorted == info_true_sorted + + np.testing.assert_allclose([info_true[k] for k in info_true_sorted], + [info[k] for k in info_sorted], + rtol=tolerance) + def test_step_pole_cancellation(self, pole_cancellation, no_pole_cancellation): # confirm that pole-zero cancellation doesn't perturb results diff --git a/control/timeresp.py b/control/timeresp.py index 3df225de9..9c7b7a990 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -746,7 +746,8 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, Parameters ---------- - sys : StateSpace or TransferFunction + sys : SISO dynamic system model. Dynamic systems that you can use include: + StateSpace or TransferFunction LTI system to simulate T : array_like or float, optional @@ -785,43 +786,82 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, -------- >>> info = step_info(sys, T) ''' - _, sys = _get_ss_simo(sys) + + if not sys.issiso(): + sys = _mimo2siso(sys,0,0) + warnings.warn(" Internal conversion from a MIMO system to a SISO system," + " the first input and the first output were used (u1 -> y1);" + " it may not be the result you are looking for") + if T is None or np.asarray(T).size == 1: T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True) T, yout = step_response(sys, T) # Steady state value - InfValue = sys.dcgain() - - # 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 + InfValue = sys.dcgain().real + + # TODO: this could be a function step_info4data(t,y,yfinal) + rise_time: float = np.NaN + settling_time: float = np.NaN + settling_min: float = np.NaN + settling_max: float = np.NaN + peak_value: float = np.Inf + peak_time: float = np.Inf + undershoot: float = np.NaN + overshoot: float = np.NaN + steady_state_value: float = np.NaN + + if not np.isnan(InfValue) and not np.isinf(InfValue): + # SteadyStateValue + steady_state_value = InfValue + # Peak + peak_index = np.abs(yout).argmax() + peak_value = np.abs(yout[peak_index]) + peak_time = T[peak_index] + + sup_margin = (1. + SettlingTimeThreshold) * InfValue + inf_margin = (1. - SettlingTimeThreshold) * InfValue + + # RiseTime + tr_lower_index = (np.where(np.sign(InfValue.real) * (yout- RiseTimeLimits[0] * InfValue) >= 0 )[0])[0] + tr_upper_index = (np.where(np.sign(InfValue.real) * yout >= np.sign(InfValue.real) * RiseTimeLimits[1] * InfValue)[0])[0] + + # SettlingTime + settling_time = T[np.where(np.abs(yout-InfValue) >= np.abs(SettlingTimeThreshold*InfValue))[0][-1]+1] + # Overshoot and Undershoot + y_os = (np.sign(InfValue.real)*yout).max() + dy_os = np.abs(y_os) - np.abs(InfValue) + if dy_os > 0: + overshoot = np.abs(100. * dy_os / InfValue) + else: + overshoot = 0 + + y_us = (np.sign(InfValue.real)*yout).min() + dy_us = np.abs(y_us) + if dy_us > 0: + undershoot = np.abs(100. * dy_us / InfValue) + else: + undershoot = 0 + + # RiseTime + rise_time = T[tr_upper_index] - T[tr_lower_index] + + settling_max = (yout[tr_upper_index:]).max() + settling_min = (yout[tr_upper_index:]).min() - PeakIndex = np.abs(yout).argmax() return { - 'RiseTime': RiseTime, - 'SettlingTime': SettlingTime, - 'SettlingMin': yout[tr_upper_index:].min(), - 'SettlingMax': yout.max(), - 'Overshoot': 100. * (yout.max() - InfValue) / InfValue, - 'Undershoot': yout.min(), # not very confident about this - 'Peak': yout[PeakIndex], - 'PeakTime': T[PeakIndex], - 'SteadyStateValue': InfValue + 'RiseTime': rise_time, + 'SettlingTime': settling_time, + 'SettlingMin': settling_min, + 'SettlingMax': settling_max, + 'Overshoot': overshoot, + 'Undershoot': undershoot, + 'Peak': peak_value, + 'PeakTime': peak_time, + 'SteadyStateValue': steady_state_value } - def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None, transpose=False, return_x=False, squeeze=None): # pylint: disable=W0622