From 12635500e42d89cb4abfbb0fed88d2d6b5594e72 Mon Sep 17 00:00:00 2001 From: jpp Date: Fri, 12 Mar 2021 17:41:03 -0300 Subject: [PATCH 1/4] =?UTF-8?q?Added=205=20test=20for=20step=5Finfo:=201)?= =?UTF-8?q?=20System=20Type=201=20-=20Step=20response=20not=20stationary:?= =?UTF-8?q?=C2=A0=20G(s)=3D1/s(s+1)=202)=20SISO=20system=20with=20under=20?= =?UTF-8?q?shoot=20response=20and=20positive=20final=20value=20G(s)=3D(-s+?= =?UTF-8?q?1)/(s=C2=B2+s+1)=203)=20Same=20system=20that=202)=20with=20k=3D?= =?UTF-8?q?-1=204)=20example=20from=20matlab=20online=20help=20https://www?= =?UTF-8?q?.mathworks.com/help/control/ref/stepinfo.html=20G(s)=3D(s=C2=B2?= =?UTF-8?q?+5s+5)/(s^4+1.65s^3+6.5s+2)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit with stepinfo output:         RiseTime: 3.8456     SettlingTime: 27.9762      SettlingMin: 2.0689      SettlingMax: 2.6873        Overshoot: 7.4915       Undershoot: 0             Peak: 2.6873         PeakTime: 8.0530 5) example from matlab online help https://www.mathworks.com/help/control/ref/stepinfo.html 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); examine the response characteristics for the response from the first input to the second output of sys. with stepinfo output:         RiseTime: 0.4000     SettlingTime: 2.8000      SettlingMin: -0.6724      SettlingMax: -0.5188        Overshoot: 24.6476       Undershoot: 11.1224             Peak: 0.6724         PeakTime: 1 --- control/tests/timeresp_test.py | 105 ++++++++++++++++++++++++++++++++- 1 file changed, 104 insertions(+), 1 deletion(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 37fcff763..54a349523 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 ej: 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 From f9641709f8f4b647715ad1e71a2533a8d3757ec0 Mon Sep 17 00:00:00 2001 From: jpp Date: Fri, 12 Mar 2021 17:45:47 -0300 Subject: [PATCH 2/4] Solve issue #337, #565 and #564 --- control/timeresp.py | 79 ++++++++++++++++++++++++++++++++------------- 1 file changed, 56 insertions(+), 23 deletions(-) diff --git a/control/timeresp.py b/control/timeresp.py index 3df225de9..7df7f34d6 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -794,34 +794,67 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, # 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 + # 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 + + if not np.isnan(InfValue) and not np.isinf(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 + + if InfValue < 0.0: + # RiseTime + tr_lower_index = (np.where(yout <= RiseTimeLimits[0] * InfValue)[0])[0] + tr_upper_index = (np.where(yout <= RiseTimeLimits[1] * InfValue)[0])[0] + # SettlingTime + for i in reversed(range(T.size - 1)): + if (-yout[i] <= np.abs(inf_margin)) | (-yout[i] >= np.abs(sup_margin)): + settling_time = T[i + 1] + break + # Overshoot and Undershoot + overshoot = np.abs(100. * ((-yout).max() - np.abs(InfValue)) / np.abs(InfValue)) + undershoot = np.abs(100. * (-yout).min() / np.abs(InfValue)) + else: + tr_lower_index = (np.where(yout >= RiseTimeLimits[0] * InfValue)[0])[0] + tr_upper_index = (np.where(yout >= RiseTimeLimits[1] * InfValue)[0])[0] + # SettlingTime + for i in reversed(range(T.size - 1)): + if (yout[i] <= inf_margin) | (yout[i] >= sup_margin): + settling_time = T[i + 1] + break + # Overshoot and Undershoot + overshoot = np.abs(100. * (yout.max() - InfValue) / InfValue) + undershoot = np.abs(100. * yout.min() / InfValue) + + # 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], + 'RiseTime': rise_time, + 'SettlingTime': settling_time, + 'SettlingMin': settling_min, + 'SettlingMax': settling_max, + 'Overshoot': overshoot, + 'Undershoot': undershoot, + 'Peak': peak_value, + 'PeakTime': peak_time, 'SteadyStateValue': InfValue } - 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 From 018d12818b7ab19924eca82a71124d9dbb9835b1 Mon Sep 17 00:00:00 2001 From: jpp Date: Tue, 16 Mar 2021 11:16:05 -0300 Subject: [PATCH 3/4] optimize the code and solve problems with MIMO systems converting to SISO systems from input=0 to output =0 solve problems with non stationary systems doing SteadyStateValue= nan when y_final is inf --- control/timeresp.py | 61 +++++++++++++++++++++++++-------------------- 1 file changed, 34 insertions(+), 27 deletions(-) diff --git a/control/timeresp.py b/control/timeresp.py index 7df7f34d6..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,14 +786,20 @@ 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() + InfValue = sys.dcgain().real # TODO: this could be a function step_info4data(t,y,yfinal) rise_time: float = np.NaN @@ -803,8 +810,11 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, 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]) @@ -813,29 +823,26 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, sup_margin = (1. + SettlingTimeThreshold) * InfValue inf_margin = (1. - SettlingTimeThreshold) * InfValue - if InfValue < 0.0: - # RiseTime - tr_lower_index = (np.where(yout <= RiseTimeLimits[0] * InfValue)[0])[0] - tr_upper_index = (np.where(yout <= RiseTimeLimits[1] * InfValue)[0])[0] - # SettlingTime - for i in reversed(range(T.size - 1)): - if (-yout[i] <= np.abs(inf_margin)) | (-yout[i] >= np.abs(sup_margin)): - settling_time = T[i + 1] - break - # Overshoot and Undershoot - overshoot = np.abs(100. * ((-yout).max() - np.abs(InfValue)) / np.abs(InfValue)) - undershoot = np.abs(100. * (-yout).min() / np.abs(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: - tr_lower_index = (np.where(yout >= RiseTimeLimits[0] * InfValue)[0])[0] - tr_upper_index = (np.where(yout >= RiseTimeLimits[1] * InfValue)[0])[0] - # SettlingTime - for i in reversed(range(T.size - 1)): - if (yout[i] <= inf_margin) | (yout[i] >= sup_margin): - settling_time = T[i + 1] - break - # Overshoot and Undershoot - overshoot = np.abs(100. * (yout.max() - InfValue) / InfValue) - undershoot = np.abs(100. * yout.min() / InfValue) + undershoot = 0 # RiseTime rise_time = T[tr_upper_index] - T[tr_lower_index] @@ -852,7 +859,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, 'Undershoot': undershoot, 'Peak': peak_value, 'PeakTime': peak_time, - 'SteadyStateValue': InfValue + 'SteadyStateValue': steady_state_value } def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None, From 81ae64f80222f97338af3b52d53b42a993300e3c Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Wed, 17 Mar 2021 11:04:04 +0100 Subject: [PATCH 4/4] fix comment format --- control/tests/timeresp_test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 54a349523..8705f3805 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -339,8 +339,8 @@ def test_step_info(self): [Strue[k] for k in Sktrue], rtol=rtol) - # tolerance for all parameters could be wrong for some systems ej: discrete systems time parameters tolerance - # could be +/-dt + # 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", {