From f7d3db5468329acde058bd9e8b03ee92f809f3e1 Mon Sep 17 00:00:00 2001 From: jpp Date: Thu, 18 Mar 2021 15:01:04 -0300 Subject: [PATCH 1/2] new feature and bug fixes: Now can compute the step response characteristics for MIMO systems. Fix undershoot caculus for sistems with negative response. --- control/tests/timeresp_test.py | 122 ++++++++++++++++++++++- control/timeresp.py | 172 +++++++++++++++++++-------------- 2 files changed, 222 insertions(+), 72 deletions(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 8705f3805..3adeaa120 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -222,6 +222,15 @@ def tf2_matlab_help(self): D = [0.06] sys = StateSpace(A, B, C, D, 0.2) return sys + + @pytest.fixture + def mimo_matlab_help(self): + A = [[0.68, -0.34], [0.34, 0.68]]; + B = [[0.18, -0.05], [0.04, 0.11]]; + C = [[0, -1.53], [-1.12, -1.10]]; + D = [[0, 0], [0.06, -0.37]]; + sys = StateSpace(A,B,C,D,0.2); + return sys @pytest.fixture def tsystem(self, @@ -233,7 +242,7 @@ def tsystem(self, mimo_dss1, mimo_dss2, mimo_dtf1, pole_cancellation, no_pole_cancellation, siso_tf_type1, siso_tf_kpos, siso_tf_kneg, tf1_matlab_help, - tf2_matlab_help): + tf2_matlab_help, mimo_matlab_help): systems = {"siso_ss1": siso_ss1, "siso_ss2": siso_ss2, "siso_tf1": siso_tf1, @@ -256,6 +265,7 @@ def tsystem(self, "siso_tf_kneg": siso_tf_kneg, "tf1_matlab_help": tf1_matlab_help, "tf2_matlab_help": tf2_matlab_help, + "mimo_matlab_help":mimo_matlab_help, } return systems[request.param] @@ -341,6 +351,116 @@ def test_step_info(self): # tolerance for all parameters could be wrong for some systems # discrete systems time parameters tolerance could be +/-dt + 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, + 'SteadyStateValue': 2.50 + } + + S = step_info(sys) + + Sk = sorted(S.keys()) + Sktrue = sorted(Strue.keys()) + assert Sk == Sktrue + # 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[k] for k in Sk], + [Strue[k] for k in Sktrue], + rtol=rtol) + + def test_step_info_mimo(self): + # From matlab docs: + A = [[0.68, -0.34], [0.34, 0.68]]; + B = [[0.18, -0.05], [0.04, 0.11]]; + C = [[0, -1.53], [-1.12, -1.10]]; + D = [[0, 0], [0.06, -0.37]]; + sys = StateSpace(A,B,C,D,0.2); + + Strue = [[{'RiseTime': 0.6000, + 'SettlingTime': 3.0, + 'SettlingMin': -0.5999, + 'SettlingMax': -0.4689, + 'Overshoot': 15.5072, + 'Undershoot': 0, + 'Peak': 0.5999, + 'PeakTime': 1.4000, + 'SteadyStateValue': -0.51935}, + {'RiseTime': 0.0, + 'SettlingTime': 3.6000, + 'SettlingMin': -0.2797, + 'SettlingMax': -0.1043, + 'Overshoot': 118.9918, + 'Undershoot': 0, + 'Peak': 0.2797, + 'PeakTime': 0.60000, + 'SteadyStateValue': -0.1277}], + [{'RiseTime': 0.40000, + 'SettlingTime': 2.8, + 'SettlingMin': -0.6724, + 'SettlingMax': -0.5188, + 'Overshoot': 24.6476, + 'Undershoot': 11.1224, + 'Peak': 0.6724, + 'PeakTime': 1.0, + 'SteadyStateValue': -0.5394}, + {'RiseTime': 0.0, + 'SettlingTime': 3.4, + 'SettlingMin': -0.435, #this value is not a result of step_info + 'SettlingMax': -0.1485, + 'Overshoot': 132.0170, + 'Undershoot': 0, + 'Peak': 0.4350, + 'PeakTime': 0.2, + 'SteadyStateValue': -0.18748}]] + + S = step_info(sys) + + S1k = sorted(S[0][0].keys()) + S1ktrue = sorted(Strue[0][0].keys()) + assert S1k == S1ktrue + + S2k = sorted(S[0][1].keys()) + S2ktrue = sorted(Strue[0][1].keys()) + assert S2k == S2ktrue + + S3k = sorted(S[1][0].keys()) + S3ktrue = sorted(Strue[1][0].keys()) + assert S3k == S3ktrue + + S4k = sorted(S[1][1].keys()) + S4ktrue = sorted(Strue[1][1].keys()) + assert S4k == S4ktrue + + rtol = 0.2 # tolerance could be better was chosen to meet the discret times parameters + np.testing.assert_allclose([S[0][0][k] for k in S1k], + [Strue[0][0][k] for k in S1ktrue], + rtol=rtol) + np.testing.assert_allclose([S[0][1][k] for k in S2k], + [Strue[0][1][k] for k in S2ktrue], + rtol=rtol) + np.testing.assert_allclose([S[1][0][k] for k in S3k], + [Strue[1][0][k] for k in S3ktrue], + rtol=rtol) + np.testing.assert_allclose([S[1][1][k] for k in S4k], + [Strue[1][1][k] for k in S4ktrue], + 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", { diff --git a/control/timeresp.py b/control/timeresp.py index 9c7b7a990..9c6968836 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -787,80 +787,110 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, >>> info = step_info(sys, T) ''' - 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: + if T is None: + T, y = step_response(sys) + elif np.asarray(T).size == 1: T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True) - - T, yout = step_response(sys, T) + T, y = step_response(sys, T) + else: + T, y = step_response(sys, T) # Steady state value - 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() - - return { - '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 - } + yfinal = sys.dcgain().real + + # TODO: this could be a function step_info4data(t,y,yfinal) + # y: Step-response data, specified as: + # For SISO response data, a vector of length Ns, where Ns is the number of samples in the response data. + # For MIMO response data, an Ny-by-Nu-by-Ns array, where Ny is the number of system outputs, and Nu is the number of system inputs. + # T: Time vector corresponding to the response data in y, specified as a vector of length Ns. + + # TODO: check dimentions for time vs step output and final output + if len(y.shape) == 3: # MIMO + Ny,Nu,Ns = y.shape + elif len(y.shape) == 1: # SISO + Ns = y.shape + Ny = 1 + Nu = 1 + else: + raise TypeError("Poner un mensaje") + + S = [[[] for x in range(Ny)] for x in range(Nu)] + + # Struct for each couple of inputs + for i in range(Ny): # Ny + for j in range(Nu): # Nu + 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 len(y.shape) > 1: + yout = y[i,j,:] + InfValue = yfinal[i,j] + else: + yout = y + InfValue = yfinal + + InfValue_sign = np.sign(InfValue) + + 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(InfValue_sign * (yout- RiseTimeLimits[0] * InfValue) >= 0 )[0])[0] + tr_upper_index = (np.where(InfValue_sign * yout >= InfValue_sign * 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 = (InfValue_sign*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 = (InfValue_sign*yout).min() + y_us_index = (InfValue_sign*yout).argmin() + if (InfValue_sign * yout[y_us_index]) < 0: # must have oposite sign + undershoot = np.abs(100. * np.abs(y_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() + + S[i][j] = { + '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 + } + if Ny > 1 or Nu > 1: + return S + else: + return S[0][0] def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None, transpose=False, return_x=False, squeeze=None): From 8d5810f05627fd426534de08217b1c5f71c0379a Mon Sep 17 00:00:00 2001 From: jpp Date: Thu, 18 Mar 2021 15:06:45 -0300 Subject: [PATCH 2/2] modify the help according to the new characteristics and add a working example too, finaly I corrected the see also to referring to functions of this library --- control/timeresp.py | 36 ++++++++++++++++++++++++++---------- 1 file changed, 26 insertions(+), 10 deletions(-) diff --git a/control/timeresp.py b/control/timeresp.py index 9c6968836..3cbbccf56 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -746,7 +746,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, Parameters ---------- - sys : SISO dynamic system model. Dynamic systems that you can use include: + sys : SISO or MIMO dynamic system model. Dynamic systems that you can use include: StateSpace or TransferFunction LTI system to simulate @@ -766,25 +766,41 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, Returns ------- - S: a dictionary containing: - RiseTime: Time from 10% to 90% of the steady-state value. + S: a dictionary for SISO systems or a list of Ny-by-Nu of dictionaries + for MIMO systems, each dictionary containing: + + RiseTime: Rise time (by default 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 + SettlingMin: Minimum value of `y(t)` after RiseTime + SettlingMax: Maximum value of `y(t)` after RiseTime Overshoot: Percentage of the Peak relative to steady value - Undershoot: Percentage of undershoot - Peak: Absolute peak value - PeakTime: time of the Peak + Undershoot: Percentage of undershoot + Peak: Peak absolute value of `y(t)` + PeakTime: Time of the peak value SteadyStateValue: Steady-state value See Also -------- - step, lsim, initial, impulse + step_response, initial_response, impulse_response, forced_response + Examples -------- - >>> info = step_info(sys, T) + >>> from control import tf, step_info + >>> sys = tf([3.32, 0, -162.8],[1, 24.56, 186.5, 457.8, 116.2]) + >>> step_info(sys) + >>> + >>> {'RiseTime': 7.70245002940221, + >>> 'SettlingTime': 14.13151114265325, + >>> 'SettlingMin': -1.3994264225116284, + >>> 'SettlingMax': -1.2612640212267976, + >>> 'Overshoot': 0, + >>> 'Undershoot': 0.6864433321178778, + >>> 'Peak': 1.3994264225116284, + >>> 'PeakTime': 24.13227287437709, + >>> 'SteadyStateValue': -1.4010327022375215} + ''' if T is None: