From ff960c38ff162cfb70ff8b138518bc35a929abbf Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Wed, 17 Mar 2021 17:13:04 +0100 Subject: [PATCH 01/13] enhance step_info to MIMO --- control/timeresp.py | 188 +++++++++++++++++++++++--------------------- 1 file changed, 100 insertions(+), 88 deletions(-) diff --git a/control/timeresp.py b/control/timeresp.py index 9c7b7a990..7a2d1c8a4 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -746,36 +746,43 @@ 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: - StateSpace or TransferFunction - LTI system to simulate - + sys : LTI system T : array_like or float, optional Time vector, or simulation time duration if a number (time vector is autocomputed if not given, see :func:`step_response` for more detail) - T_num : int, optional Number of time steps to use in simulation if T is not provided as an array (autocomputed if not given); ignored if sys is discrete-time. - 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 + S : dict or list of list of dict + If `sys` is a SISO system, S is 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 + + If `sys` is a MIMO system, S is a 2D-List of dicts. See Also @@ -786,81 +793,86 @@ 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: 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().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 - } + ret = [[None] * sys.ninputs] * sys.noutputs + for i in range(sys.noutputs): + for j in range(sys.ninputs): + sys_siso = sys[i, j] + T, yout = step_response(sys_siso, T) + + # Steady state value + InfValue = sys_siso.dcgain() + sgnInf = np.sign(InfValue.real) + + 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: complex = np.NaN + + if not np.isnan(InfValue) and not np.isinf(InfValue): + # RiseTime + tr_lower_index = np.where( + sgnInf * (yout - RiseTimeLimits[0] * InfValue) >= 0 + )[0][0] + tr_upper_index = np.where( + sgnInf * (yout - RiseTimeLimits[1] * InfValue) >= 0 + )[0][0] + rise_time = T[tr_upper_index] - T[tr_lower_index] + + # SettlingTime + settling_th = np.abs(SettlingTimeThreshold * InfValue) + not_settled = np.where(np.abs(yout - InfValue) >= settling_th) + settling_time = T[not_settled[0][-1] + 1] + + settling_min = (yout[tr_upper_index:]).min() + settling_max = (yout[tr_upper_index:]).max() + + # Overshoot + y_os = (sgnInf * 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 + + # Undershoot + y_us = (sgnInf * yout).min() + dy_us = np.abs(y_us) + if dy_us > 0: + undershoot = np.abs(100. * dy_us / InfValue) + else: + undershoot = 0 + + # Peak + peak_index = np.abs(yout).argmax() + peak_value = np.abs(yout[peak_index]) + peak_time = T[peak_index] + + # SteadyStateValue + steady_state_value = InfValue + + ret[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 + } + + return ret[0][0] if sys.issiso() else ret def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None, transpose=False, return_x=False, squeeze=None): From b98008a47765b90661ab0838bcdd00463e2f97ba Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Wed, 17 Mar 2021 17:22:23 +0100 Subject: [PATCH 02/13] remove test function with duplicate name --- control/tests/timeresp_test.py | 46 ++++++++-------------------------- 1 file changed, 10 insertions(+), 36 deletions(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 8705f3805..bc6de9ea3 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -193,7 +193,7 @@ 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) @@ -309,36 +309,6 @@ def test_step_nostates(self, dt): t, y = step_response(sys) np.testing.assert_array_equal(y, np.ones(len(t))) - 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) - # tolerance for all parameters could be wrong for some systems # discrete systems time parameters tolerance could be +/-dt @pytest.mark.parametrize( @@ -394,17 +364,21 @@ def test_step_info(self): 'SteadyStateValue': np.NaN}, 2e-2)], indirect=["tsystem"]) def test_step_info(self, tsystem, info_true, tolerance): - """ """ + """Test step info for SISO systems""" 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) + for k in info: + np.testing.assert_allclose(info[k], info_true[k], rtol=tolerance, + err_msg=f"key {k} does not match") + + def test_step_info_mimo(self, tsystem, info_true, tolearance): + """Test step info for MIMO systems""" + # TODO: implement + pass def test_step_pole_cancellation(self, pole_cancellation, no_pole_cancellation): From ed633583d3ea2f7612e43c8d65f6604c2f812162 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Wed, 17 Mar 2021 22:01:07 +0100 Subject: [PATCH 03/13] add step_info mimo test (TransferFunction) --- control/tests/timeresp_test.py | 174 ++++++++++++++++++++------------- 1 file changed, 106 insertions(+), 68 deletions(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index bc6de9ea3..dce0a6a49 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -66,8 +66,8 @@ def siso_ss2(self, siso_ss1): T.initial = siso_ss1.yinitial - 9 T.yimpulse = np.array([86., 70.1808, 57.3753, 46.9975, 38.5766, 31.7344, 26.1668, 21.6292, 17.9245, 14.8945]) - return T + return T @pytest.fixture def siso_tf1(self): @@ -197,31 +197,98 @@ def no_pole_cancellation(self): @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]) + T = TSys(TransferFunction(1, [1, 1, 0])) + T.step_info = { + '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} + return T @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]) + T = TSys(TransferFunction([-1, 1], [1, 1, 1])) + T.step_info = { + '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} + return T @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]) + T = TSys(TransferFunction([1, -1], [1, 1, 1])) + T.step_info = { + '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} + return T @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]) + T = TSys(TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2])) + T.step_info = { + '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} + return T @pytest.fixture - def tf2_matlab_help(self): + def ss2_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 + T = TSys(StateSpace(A, B, C, D, 0.2)) + T.step_info = { + '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} + return T + + @pytest.fixture + def mimo_tf_step(self, tf1_matlab_help, + siso_tf_kpos, + siso_tf_kneg, + siso_tf_type1): + Ta = [[tf1_matlab_help, tf1_matlab_help, siso_tf_kpos], + [siso_tf_kneg, siso_tf_type1, siso_tf_type1]] + T = TSys(TransferFunction( + [[Ti.sys.num[0][0] for Ti in Tr] for Tr in Ta], + [[Ti.sys.den[0][0] for Ti in Tr] for Tr in Ta])) + T.step_info = [[Ti.step_info for Ti in Tr] for Tr in Ta] + return T @pytest.fixture def tsystem(self, @@ -233,7 +300,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): + ss2_matlab_help, mimo_tf_step): systems = {"siso_ss1": siso_ss1, "siso_ss2": siso_ss2, "siso_tf1": siso_tf1, @@ -255,7 +322,8 @@ def tsystem(self, "siso_tf_kpos": siso_tf_kpos, "siso_tf_kneg": siso_tf_kneg, "tf1_matlab_help": tf1_matlab_help, - "tf2_matlab_help": tf2_matlab_help, + "ss2_matlab_help": ss2_matlab_help, + "mimo_tf_step": mimo_tf_step, } return systems[request.param] @@ -312,73 +380,43 @@ def test_step_nostates(self, 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", { - '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)], + "tsystem, tolerance", + [("tf1_matlab_help", 2e-2), + ("ss2_matlab_help", .2), + ("siso_tf_kpos", 2e-2), + ("siso_tf_kneg", 2e-2), + ("siso_tf_type1", 2e-2)], indirect=["tsystem"]) - def test_step_info(self, tsystem, info_true, tolerance): + def test_step_info(self, tsystem, tolerance): """Test step info for SISO systems""" - info = step_info(tsystem) + info = step_info(tsystem.sys) - info_true_sorted = sorted(info_true.keys()) + info_true_sorted = sorted(tsystem.step_info.keys()) info_sorted = sorted(info.keys()) assert info_sorted == info_true_sorted for k in info: - np.testing.assert_allclose(info[k], info_true[k], rtol=tolerance, - err_msg=f"key {k} does not match") + np.testing.assert_allclose(info[k], tsystem.step_info[k], + rtol=tolerance, + err_msg=f"{k} does not match") - def test_step_info_mimo(self, tsystem, info_true, tolearance): + @pytest.mark.parametrize( + "tsystem, tolerance", + [('mimo_tf_step', 2e-2)], + indirect=["tsystem"]) + def test_step_info_mimo(self, tsystem, tolerance): """Test step info for MIMO systems""" - # TODO: implement - pass + info_dict = step_info(tsystem.sys) + from pprint import pprint + pprint(info_dict) + for i, row in enumerate(info_dict): + for j, info in enumerate(row): + for k in info: + np.testing.assert_allclose( + info[k], tsystem.step_info[i][j][k], + rtol=tolerance, + err_msg=f"{k} for input {j} to output {i} " + "does not match") def test_step_pole_cancellation(self, pole_cancellation, no_pole_cancellation): From 35d8ebaa77de7a01ee876439b4b8978944521a58 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Wed, 17 Mar 2021 22:01:48 +0100 Subject: [PATCH 04/13] fix timevector and list return population for MIMO step_info --- control/timeresp.py | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/control/timeresp.py b/control/timeresp.py index 7a2d1c8a4..bb51d3449 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -793,16 +793,17 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, -------- >>> info = step_info(sys, T) ''' - - - if T is None or np.asarray(T).size == 1: - T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True) - - ret = [[None] * sys.ninputs] * sys.noutputs + ret = [] for i in range(sys.noutputs): + retrow = [] for j in range(sys.ninputs): sys_siso = sys[i, j] - T, yout = step_response(sys_siso, T) + if T is None or np.asarray(T).size == 1: + Ti = _default_time_vector(sys_siso, N=T_num, tfinal=T, + is_step=True) + else: + Ti = T + Ti, yout = step_response(sys_siso, Ti) # Steady state value InfValue = sys_siso.dcgain() @@ -826,12 +827,12 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, tr_upper_index = np.where( sgnInf * (yout - RiseTimeLimits[1] * InfValue) >= 0 )[0][0] - rise_time = T[tr_upper_index] - T[tr_lower_index] + rise_time = Ti[tr_upper_index] - Ti[tr_lower_index] # SettlingTime settling_th = np.abs(SettlingTimeThreshold * InfValue) not_settled = np.where(np.abs(yout - InfValue) >= settling_th) - settling_time = T[not_settled[0][-1] + 1] + settling_time = Ti[not_settled[0][-1] + 1] settling_min = (yout[tr_upper_index:]).min() settling_max = (yout[tr_upper_index:]).max() @@ -855,12 +856,12 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, # Peak peak_index = np.abs(yout).argmax() peak_value = np.abs(yout[peak_index]) - peak_time = T[peak_index] + peak_time = Ti[peak_index] # SteadyStateValue steady_state_value = InfValue - ret[i][j] = { + retij = { 'RiseTime': rise_time, 'SettlingTime': settling_time, 'SettlingMin': settling_min, @@ -871,6 +872,9 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, 'PeakTime': peak_time, 'SteadyStateValue': steady_state_value } + retrow.append(retij) + + ret.append(retrow) return ret[0][0] if sys.issiso() else ret From 468ffbff098f5616559e0c9d269312be2c724cff Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 18 Mar 2021 23:44:00 +0100 Subject: [PATCH 05/13] apply isort --- control/tests/timeresp_test.py | 8 +++----- control/timeresp.py | 19 +++++++++---------- 2 files changed, 12 insertions(+), 15 deletions(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index dce0a6a49..ddbe28fe7 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -15,15 +15,13 @@ import pytest import scipy as sp - import control as ct -from control import (StateSpace, TransferFunction, c2d, isctime, isdtime, - ss2tf, tf2ss) +from control import StateSpace, TransferFunction, c2d, isctime, ss2tf, tf2ss +from control.exception import slycot_check +from control.tests.conftest import slycotonly from control.timeresp import (_default_time_vector, _ideal_tfinal_and_dt, forced_response, impulse_response, initial_response, step_info, step_response) -from control.tests.conftest import slycotonly -from control.exception import slycot_check class TSys: diff --git a/control/timeresp.py b/control/timeresp.py index bb51d3449..405f6768f 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -71,18 +71,17 @@ $Id$ """ -# Libraries that we make use of -import scipy as sp # SciPy library (used all over) -import numpy as np # NumPy library -from scipy.linalg import eig, eigvals, matrix_balance, norm -from numpy import (einsum, maximum, minimum, - atleast_1d) import warnings -from .lti import LTI # base class of StateSpace, TransferFunction -from .xferfcn import TransferFunction -from .statesp import _convert_to_statespace, _mimo2simo, _mimo2siso, ssdata -from .lti import isdtime, isctime + +import numpy as np +import scipy as sp +from numpy import atleast_1d, einsum, maximum, minimum +from scipy.linalg import eig, eigvals, matrix_balance, norm + from . import config +from .lti import (LTI, isctime, isdtime) +from .statesp import _convert_to_statespace, _mimo2simo, _mimo2siso, ssdata +from .xferfcn import TransferFunction __all__ = ['forced_response', 'step_response', 'step_info', 'initial_response', 'impulse_response'] From a1fd47e5f405e2ea0e026871a77ba34ed58adfe9 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Fri, 19 Mar 2021 00:00:51 +0100 Subject: [PATCH 06/13] step_info from MIMO step_response --- control/tests/timeresp_test.py | 206 +++++++++++++++++++++------------ control/timeresp.py | 29 ++--- 2 files changed, 149 insertions(+), 86 deletions(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index ddbe28fe7..fd22001c0 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -26,8 +26,9 @@ class TSys: """Struct of test system""" - def __init__(self, sys=None): + def __init__(self, sys=None, call_kwargs=None): self.sys = sys + self.kwargs = call_kwargs if call_kwargs else {} def __repr__(self): """Show system when debugging""" @@ -210,15 +211,16 @@ def siso_tf_type1(self): @pytest.fixture def siso_tf_kpos(self): - # SISO under shoot response and positive final value G(s)=(-s+1)/(s²+s+1) + # SISO under shoot response and positive final value + # G(s)=(-s+1)/(s²+s+1) T = TSys(TransferFunction([-1, 1], [1, 1, 1])) T.step_info = { 'RiseTime': 1.242, 'SettlingTime': 9.110, - 'SettlingMin': 0.950, + 'SettlingMin': 0.90, 'SettlingMax': 1.208, 'Overshoot': 20.840, - 'Undershoot': 27.840, + 'Undershoot': 28.0, 'Peak': 1.208, 'PeakTime': 4.282, 'SteadyStateValue': 1.0} @@ -226,23 +228,25 @@ def siso_tf_kpos(self): @pytest.fixture def siso_tf_kneg(self): - # SISO under shoot response and negative final value k=-1 G(s)=-(-s+1)/(s²+s+1) + # SISO under shoot response and negative final value + # k=-1 G(s)=-(-s+1)/(s²+s+1) T = TSys(TransferFunction([1, -1], [1, 1, 1])) T.step_info = { 'RiseTime': 1.242, 'SettlingTime': 9.110, 'SettlingMin': -1.208, - 'SettlingMax': -0.950, + 'SettlingMax': -0.90, 'Overshoot': 20.840, - 'Undershoot': 27.840, + 'Undershoot': 28.0, 'Peak': 1.208, 'PeakTime': 4.282, 'SteadyStateValue': -1.0} return T @pytest.fixture - def tf1_matlab_help(self): - # example from matlab online help https://www.mathworks.com/help/control/ref/stepinfo.html + def siso_tf_step_matlab(self): + # example from matlab online help + # https://www.mathworks.com/help/control/ref/stepinfo.html T = TSys(TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2])) T.step_info = { 'RiseTime': 3.8456, @@ -257,37 +261,82 @@ def tf1_matlab_help(self): return T @pytest.fixture - def ss2_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] + def mimo_ss_step_matlab(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]] T = TSys(StateSpace(A, B, C, D, 0.2)) - T.step_info = { - '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} + T.kwargs['step_info'] = {'T': 4.6} + T.step_info = [[{'RiseTime': 0.6000, + 'SettlingTime': 3.0000, + 'SettlingMin': -0.5999, + 'SettlingMax': -0.4689, + 'Overshoot': 15.5072, + 'Undershoot': 0., + 'Peak': 0.5999, + 'PeakTime': 1.4000, + 'SteadyStateValue': -0.5193}, + {'RiseTime': 0., + 'SettlingTime': 3.6000, + 'SettlingMin': -0.2797, + 'SettlingMax': -0.1043, + 'Overshoot': 118.9918, + 'Undershoot': 0, + 'Peak': 0.2797, + 'PeakTime': .6000, + 'SteadyStateValue': -0.1277}], + [{'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}, + {'RiseTime': 0.0000, # (*) + 'SettlingTime': 3.4000, + 'SettlingMin': -0.1034, + 'SettlingMax': -0.1485, + 'Overshoot': 132.0170, + 'Undershoot': 79.222, # 0. in MATLAB + 'Peak': 0.4350, + 'PeakTime': .2, + 'SteadyStateValue': -0.1875}]] + # (*): MATLAB gives 0.4 here, but it is unclear what + # 10% and 90% of the steady state response mean, when + # the step for this channel does not start a 0 for + # 0 initial conditions return T @pytest.fixture - def mimo_tf_step(self, tf1_matlab_help, - siso_tf_kpos, - siso_tf_kneg, - siso_tf_type1): - Ta = [[tf1_matlab_help, tf1_matlab_help, siso_tf_kpos], - [siso_tf_kneg, siso_tf_type1, siso_tf_type1]] + def siso_ss_step_matlab(self, mimo_ss_step_matlab): + T = copy(mimo_ss_step_matlab) + T.sys = T.sys[1, 0] + T.step_info = T.step_info[1][0] + return T + + @pytest.fixture + def mimo_tf_step_info(self, + siso_tf_kpos, siso_tf_kneg, + siso_tf_step_matlab): + Ta = [[siso_tf_kpos, siso_tf_kneg, siso_tf_step_matlab], + [siso_tf_step_matlab, siso_tf_kpos, siso_tf_kneg]] T = TSys(TransferFunction( [[Ti.sys.num[0][0] for Ti in Tr] for Tr in Ta], [[Ti.sys.den[0][0] for Ti in Tr] for Tr in Ta])) T.step_info = [[Ti.step_info for Ti in Tr] for Tr in Ta] + # enforce enough sample points for all channels (they have different + # characteristics) + T.kwargs['step_info'] = {'T_num': 2000} return T + @pytest.fixture def tsystem(self, request, @@ -297,8 +346,9 @@ def tsystem(self, siso_dss1, siso_dss2, mimo_dss1, mimo_dss2, mimo_dtf1, pole_cancellation, no_pole_cancellation, siso_tf_type1, - siso_tf_kpos, siso_tf_kneg, tf1_matlab_help, - ss2_matlab_help, mimo_tf_step): + siso_tf_kpos, siso_tf_kneg, + siso_tf_step_matlab, siso_ss_step_matlab, + mimo_ss_step_matlab, mimo_tf_step_info): systems = {"siso_ss1": siso_ss1, "siso_ss2": siso_ss2, "siso_tf1": siso_tf1, @@ -319,9 +369,10 @@ def tsystem(self, "siso_tf_type1": siso_tf_type1, "siso_tf_kpos": siso_tf_kpos, "siso_tf_kneg": siso_tf_kneg, - "tf1_matlab_help": tf1_matlab_help, - "ss2_matlab_help": ss2_matlab_help, - "mimo_tf_step": mimo_tf_step, + "siso_tf_step_matlab": siso_tf_step_matlab, + "siso_ss_step_matlab": siso_ss_step_matlab, + "mimo_ss_step_matlab": mimo_ss_step_matlab, + "mimo_tf_step": mimo_tf_step_info, } return systems[request.param] @@ -375,46 +426,60 @@ def test_step_nostates(self, dt): t, y = step_response(sys) np.testing.assert_array_equal(y, np.ones(len(t))) - # tolerance for all parameters could be wrong for some systems - # discrete systems time parameters tolerance could be +/-dt - @pytest.mark.parametrize( - "tsystem, tolerance", - [("tf1_matlab_help", 2e-2), - ("ss2_matlab_help", .2), - ("siso_tf_kpos", 2e-2), - ("siso_tf_kneg", 2e-2), - ("siso_tf_type1", 2e-2)], - indirect=["tsystem"]) - def test_step_info(self, tsystem, tolerance): - """Test step info for SISO systems""" - info = step_info(tsystem.sys) + def assert_step_info_match(self, sys, info, info_ref): + """Assert reasonable step_info accuracy""" - info_true_sorted = sorted(tsystem.step_info.keys()) - info_sorted = sorted(info.keys()) - assert info_sorted == info_true_sorted + if sys.isdtime(strict=True): + dt = sys.dt + else: + _, dt = _ideal_tfinal_and_dt(sys, is_step=True) - for k in info: - np.testing.assert_allclose(info[k], tsystem.step_info[k], - rtol=tolerance, + for k in ['RiseTime', 'SettlingTime', 'PeakTime']: + np.testing.assert_allclose(info[k], info_ref[k], atol=dt, err_msg=f"{k} does not match") + for k in ['Overshoot', 'Undershoot', 'Peak', 'SteadyStateValue']: + np.testing.assert_allclose(info[k], info_ref[k], rtol=5e-3, + err_msg=f"{k} does not match") + + # steep gradient right after RiseTime + absrefinf = np.abs(info_ref['SteadyStateValue']) + if info_ref['RiseTime'] > 0: + y_next_sample_max = 0.8*absrefinf/info_ref['RiseTime']*dt + else: + y_next_sample_max = 0 + for k in ['SettlingMin', 'SettlingMax']: + if (np.abs(info_ref[k]) - 0.9 * absrefinf) > y_next_sample_max: + # local min/max peak well after signal has risen + np.testing.assert_allclose(info[k], info_ref[k], rtol=1e-3) + + @pytest.mark.parametrize( + "tsystem", + ["siso_tf_step_matlab", + "siso_ss_step_matlab", + "siso_tf_kpos", + "siso_tf_kneg", + "siso_tf_type1"], + indirect=["tsystem"]) + def test_step_info(self, tsystem): + """Test step info for SISO systems""" + step_info_kwargs = tsystem.kwargs.get('step_info',{}) + info = step_info(tsystem.sys, **step_info_kwargs) + self.assert_step_info_match(tsystem.sys, info, tsystem.step_info) @pytest.mark.parametrize( - "tsystem, tolerance", - [('mimo_tf_step', 2e-2)], + "tsystem", + ['mimo_ss_step_matlab', + 'mimo_tf_step'], indirect=["tsystem"]) - def test_step_info_mimo(self, tsystem, tolerance): + def test_step_info_mimo(self, tsystem): """Test step info for MIMO systems""" - info_dict = step_info(tsystem.sys) - from pprint import pprint - pprint(info_dict) + step_info_kwargs = tsystem.kwargs.get('step_info',{}) + info_dict = step_info(tsystem.sys, **step_info_kwargs) for i, row in enumerate(info_dict): for j, info in enumerate(row): for k in info: - np.testing.assert_allclose( - info[k], tsystem.step_info[i][j][k], - rtol=tolerance, - err_msg=f"{k} for input {j} to output {i} " - "does not match") + self.assert_step_info_match(tsystem.sys, + info, tsystem.step_info[i][j]) def test_step_pole_cancellation(self, pole_cancellation, no_pole_cancellation): @@ -422,13 +487,10 @@ def test_step_pole_cancellation(self, pole_cancellation, # https://github.com/python-control/python-control/issues/440 step_info_no_cancellation = step_info(no_pole_cancellation) step_info_cancellation = step_info(pole_cancellation) - for key in step_info_no_cancellation: - if key == 'Overshoot': - # skip this test because these systems have no overshoot - # => very sensitive to parameters - continue - np.testing.assert_allclose(step_info_no_cancellation[key], - step_info_cancellation[key], rtol=1e-4) + self.assert_step_info_match(no_pole_cancellation, + step_info_no_cancellation, + step_info_cancellation) + @pytest.mark.parametrize( "tsystem, kwargs", diff --git a/control/timeresp.py b/control/timeresp.py index 405f6768f..9fbe808c3 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -792,20 +792,18 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, -------- >>> info = step_info(sys, T) ''' + 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, squeeze=False) + ret = [] for i in range(sys.noutputs): retrow = [] for j in range(sys.ninputs): - sys_siso = sys[i, j] - if T is None or np.asarray(T).size == 1: - Ti = _default_time_vector(sys_siso, N=T_num, tfinal=T, - is_step=True) - else: - Ti = T - Ti, yout = step_response(sys_siso, Ti) + yout = Yout[i, j, :] # Steady state value - InfValue = sys_siso.dcgain() + InfValue = sys.dcgain() if sys.issiso() else sys.dcgain()[i, j] sgnInf = np.sign(InfValue.real) rise_time: float = np.NaN @@ -826,12 +824,15 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, tr_upper_index = np.where( sgnInf * (yout - RiseTimeLimits[1] * InfValue) >= 0 )[0][0] - rise_time = Ti[tr_upper_index] - Ti[tr_lower_index] + rise_time = T[tr_upper_index] - T[tr_lower_index] # SettlingTime - settling_th = np.abs(SettlingTimeThreshold * InfValue) - not_settled = np.where(np.abs(yout - InfValue) >= settling_th) - settling_time = Ti[not_settled[0][-1] + 1] + settled = np.where( + np.abs(yout/InfValue -1) >= SettlingTimeThreshold)[0][-1]+1 + # MIMO systems can have unsettled channels without infinite + # InfValue + if settled < len(T): + settling_time = T[settled] settling_min = (yout[tr_upper_index:]).min() settling_max = (yout[tr_upper_index:]).max() @@ -855,10 +856,10 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, # Peak peak_index = np.abs(yout).argmax() peak_value = np.abs(yout[peak_index]) - peak_time = Ti[peak_index] + peak_time = T[peak_index] # SteadyStateValue - steady_state_value = InfValue + steady_state_value = InfValue.real retij = { 'RiseTime': rise_time, From 238482e43788c524424715ef8faf9dfe88cffca2 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Fri, 19 Mar 2021 00:37:10 +0100 Subject: [PATCH 07/13] reenable masked timevector test, and really test if we get the tfinal --- control/tests/timeresp_test.py | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index fd22001c0..23fccbef7 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -706,20 +706,23 @@ def test_step_robustness(self): (TransferFunction(1, [1, .5, 0]), 25)]) # poles at 0.5 and 0 def test_auto_generated_time_vector_tfinal(self, tfsys, tfinal): """Confirm a TF with a pole at p simulates for tfinal seconds""" - np.testing.assert_almost_equal( - _ideal_tfinal_and_dt(tfsys)[0], tfinal, decimal=4) + ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(tfsys) + np.testing.assert_allclose(ideal_tfinal, tfinal, rtol=1e-4) + T = _default_time_vector(tfsys) + np.testing.assert_allclose(T[-1], tfinal, atol=0.5*ideal_dt) @pytest.mark.parametrize("wn, zeta", [(10, 0), (100, 0), (100, .1)]) - def test_auto_generated_time_vector_dt_cont(self, wn, zeta): + def test_auto_generated_time_vector_dt_cont1(self, wn, zeta): """Confirm a TF with a natural frequency of wn rad/s gets a dt of 1/(ratio*wn)""" dtref = 0.25133 / wn tfsys = TransferFunction(1, [1, 2*zeta*wn, wn**2]) - np.testing.assert_almost_equal(_ideal_tfinal_and_dt(tfsys)[1], dtref) + np.testing.assert_almost_equal(_ideal_tfinal_and_dt(tfsys)[1], dtref, + decimal=5) - def test_auto_generated_time_vector_dt_cont(self): + def test_auto_generated_time_vector_dt_cont2(self): """A sampled tf keeps its dt""" wn = 100 zeta = .1 @@ -746,21 +749,23 @@ def test_default_timevector_long(self): def test_default_timevector_functions_c(self, fun): """Test that functions can calculate the time vector automatically""" sys = TransferFunction(1, [1, .5, 0]) + _tfinal, _dt = _ideal_tfinal_and_dt(sys) # test impose number of time steps tout, _ = fun(sys, T_num=10) assert len(tout) == 10 # test impose final time - tout, _ = fun(sys, 100) - np.testing.assert_allclose(tout[-1], 100., atol=0.5) + tout, _ = fun(sys, T=100.) + np.testing.assert_allclose(tout[-1], 100., atol=0.5*_dt) @pytest.mark.parametrize("fun", [step_response, impulse_response, initial_response]) - def test_default_timevector_functions_d(self, fun): + @pytest.mark.parametrize("dt", [0.1, 0.112]) + def test_default_timevector_functions_d(self, fun, dt): """Test that functions can calculate the time vector automatically""" - sys = TransferFunction(1, [1, .5, 0], 0.1) + sys = TransferFunction(1, [1, .5, 0], dt) # test impose number of time steps is ignored with dt given tout, _ = fun(sys, T_num=15) @@ -768,7 +773,7 @@ def test_default_timevector_functions_d(self, fun): # test impose final time tout, _ = fun(sys, 100) - np.testing.assert_allclose(tout[-1], 100., atol=0.5) + np.testing.assert_allclose(tout[-1], 100., atol=0.5*dt) @pytest.mark.parametrize("tsystem", From cc90d8d5a3024a7eae7ff5329c5665a63ba17657 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Fri, 19 Mar 2021 00:37:52 +0100 Subject: [PATCH 08/13] include tfinal in auto generated timevector --- control/tests/sisotool_test.py | 12 ++++-------- control/timeresp.py | 14 ++++++++------ 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py index 6df2493cb..14e9692c1 100644 --- a/control/tests/sisotool_test.py +++ b/control/tests/sisotool_test.py @@ -68,8 +68,8 @@ def test_sisotool(self, sys): # Check the step response before moving the point step_response_original = np.array( - [0. , 0.021 , 0.124 , 0.3146, 0.5653, 0.8385, 1.0969, 1.3095, - 1.4549, 1.5231]) + [0. , 0.0216, 0.1271, 0.3215, 0.5762, 0.8522, 1.1114, 1.3221, + 1.4633, 1.5254]) assert_array_almost_equal( ax_step.lines[0].get_data()[1][:10], step_response_original, 4) @@ -113,8 +113,8 @@ def test_sisotool(self, sys): # Check if the step response has changed step_response_moved = np.array( - [0. , 0.023 , 0.1554, 0.4401, 0.8646, 1.3722, 1.875 , 2.2709, - 2.4633, 2.3827]) + [0. , 0.0237, 0.1596, 0.4511, 0.884 , 1.3985, 1.9031, 2.2922, + 2.4676, 2.3606]) assert_array_almost_equal( ax_step.lines[0].get_data()[1][:10], step_response_moved, 4) @@ -157,7 +157,3 @@ def test_sisotool_mimo(self, sys222, sys221): # but 2 input, 1 output should with pytest.raises(ControlMIMONotImplemented): sisotool(sys221) - - - - diff --git a/control/timeresp.py b/control/timeresp.py index 9fbe808c3..aa8bfedfe 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -1303,16 +1303,18 @@ def _default_time_vector(sys, N=None, tfinal=None, is_step=True): # only need to use default_tfinal if not given; N is ignored. if tfinal is None: # for discrete time, change from ideal_tfinal if N too large/small - N = int(np.clip(ideal_tfinal/sys.dt, N_min_dt, N_max))# [N_min, N_max] - tfinal = sys.dt * N + # [N_min, N_max] + N = int(np.clip(np.ceil(ideal_tfinal/sys.dt)+1, N_min_dt, N_max)) + tfinal = sys.dt * (N-1) else: - N = int(tfinal/sys.dt) - tfinal = N * sys.dt # make tfinal an integer multiple of sys.dt + N = int(np.ceil(tfinal/sys.dt)) + 1 + tfinal = sys.dt * (N-1) # make tfinal an integer multiple of sys.dt else: if tfinal is None: # for continuous time, simulate to ideal_tfinal but limit N tfinal = ideal_tfinal if N is None: - N = int(np.clip(tfinal/ideal_dt, N_min_ct, N_max)) # N<-[N_min, N_max] + # [N_min, N_max] + N = int(np.clip(np.ceil(tfinal/ideal_dt)+1, N_min_ct, N_max)) - return np.linspace(0, tfinal, N, endpoint=False) + return np.linspace(0, tfinal, N, endpoint=True) From 01ef6668470ef4444fdc6ff274ad760d3c47b359 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Fri, 19 Mar 2021 00:58:14 +0100 Subject: [PATCH 09/13] mimo is slycot only --- control/tests/timeresp_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 23fccbef7..12ca55e46 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -469,7 +469,7 @@ def test_step_info(self, tsystem): @pytest.mark.parametrize( "tsystem", ['mimo_ss_step_matlab', - 'mimo_tf_step'], + pytest.param('mimo_tf_step', marks=slycotonly)], indirect=["tsystem"]) def test_step_info_mimo(self, tsystem): """Test step info for MIMO systems""" From 20ed3682b1977acf4b98ce22d5d116c1d10559b6 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Fri, 19 Mar 2021 12:52:55 +0100 Subject: [PATCH 10/13] Describe MIMO step_info return and add doctest example --- control/timeresp.py | 49 +++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 45 insertions(+), 4 deletions(-) diff --git a/control/timeresp.py b/control/timeresp.py index aa8bfedfe..a68b25521 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -740,7 +740,7 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None, def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9)): - ''' + """ Step response characteristics (Rise time, Settling Time, Peak and others). Parameters @@ -781,7 +781,9 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, SteadyStateValue: Steady-state value - If `sys` is a MIMO system, S is a 2D-List of dicts. + If `sys` is a MIMO system, `S` is a 2D list of dicts. To get the + step response characteristics from the j-th input to the i-th output, + access ``S[i][j]`` See Also @@ -790,8 +792,47 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, Examples -------- - >>> info = step_info(sys, T) - ''' + >>> from control import step_info, TransferFunction + >>> sys = TransferFunction([-1, 1], [1, 1, 1]) + >>> S = step_info(sys) + >>> for k in S: + ... print(f"{k}: {S[k]:3.4}") + ... + RiseTime: 1.256 + SettlingTime: 9.071 + SettlingMin: 0.9011 + SettlingMax: 1.208 + Overshoot: 20.85 + Undershoot: 27.88 + Peak: 1.208 + PeakTime: 4.187 + SteadyStateValue: 1.0 + + MIMO System: Simulate until a final time of 10. Get the step response + characteristics for the second input and specify a 5% error until the + signal is considered settled. + + >>> from numpy import sqrt + >>> from control import step_info, StateSpace + >>> sys = StateSpace([[-1., -1.], + ... [1., 0.]], + ... [[-1./sqrt(2.), 1./sqrt(2.)], + ... [0, 0]], + ... [[sqrt(2.), -sqrt(2.)]], + ... [[0, 0]]) + >>> S = step_info(sys, T=10., SettlingTimeThreshold=0.05) + >>> for k, v in S[0][1].items(): + ... print(f"{k}: {float(v):3.4}") + RiseTime: 1.212 + SettlingTime: 6.061 + SettlingMin: -1.209 + SettlingMax: -0.9184 + Overshoot: 20.87 + Undershoot: 28.02 + Peak: 1.209 + PeakTime: 4.242 + SteadyStateValue: -1.0 + """ 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, squeeze=False) From 43d73f0e08b3ad1edf7ba6541679fcb21e7194bd Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Fri, 19 Mar 2021 15:12:24 +0100 Subject: [PATCH 11/13] support time series of response data in step_info --- control/tests/timeresp_test.py | 66 +++++++++++++++++++------ control/timeresp.py | 89 ++++++++++++++++++++++------------ 2 files changed, 111 insertions(+), 44 deletions(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 12ca55e46..6ec23e27e 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -427,8 +427,7 @@ def test_step_nostates(self, dt): np.testing.assert_array_equal(y, np.ones(len(t))) def assert_step_info_match(self, sys, info, info_ref): - """Assert reasonable step_info accuracy""" - + """Assert reasonable step_info accuracy.""" if sys.isdtime(strict=True): dt = sys.dt else: @@ -460,10 +459,28 @@ def assert_step_info_match(self, sys, info, info_ref): "siso_tf_kneg", "siso_tf_type1"], indirect=["tsystem"]) - def test_step_info(self, tsystem): - """Test step info for SISO systems""" - step_info_kwargs = tsystem.kwargs.get('step_info',{}) - info = step_info(tsystem.sys, **step_info_kwargs) + @pytest.mark.parametrize( + "systype, time_2d", + [("lti", False), + ("time response data", False), + ("time response data", True), + ]) + def test_step_info(self, tsystem, systype, time_2d): + """Test step info for SISO systems.""" + step_info_kwargs = tsystem.kwargs.get('step_info', {}) + if systype == "time response data": + # simulate long enough for steady state value + tfinal = 3 * tsystem.step_info['SettlingTime'] + if np.isnan(tfinal): + pytest.skip("test system does not settle") + t, y = step_response(tsystem.sys, T=tfinal, T_num=5000) + sysdata = y + step_info_kwargs['T'] = t[np.newaxis, :] if time_2d else t + else: + sysdata = tsystem.sys + + info = step_info(sysdata, **step_info_kwargs) + self.assert_step_info_match(tsystem.sys, info, tsystem.step_info) @pytest.mark.parametrize( @@ -471,15 +488,37 @@ def test_step_info(self, tsystem): ['mimo_ss_step_matlab', pytest.param('mimo_tf_step', marks=slycotonly)], indirect=["tsystem"]) - def test_step_info_mimo(self, tsystem): - """Test step info for MIMO systems""" - step_info_kwargs = tsystem.kwargs.get('step_info',{}) - info_dict = step_info(tsystem.sys, **step_info_kwargs) + @pytest.mark.parametrize( + "systype", ["lti", "time response data"]) + def test_step_info_mimo(self, tsystem, systype): + """Test step info for MIMO systems.""" + step_info_kwargs = tsystem.kwargs.get('step_info', {}) + if systype == "time response data": + tfinal = 3 * max([S['SettlingTime'] + for Srow in tsystem.step_info for S in Srow]) + t, y = step_response(tsystem.sys, T=tfinal, T_num=5000) + sysdata = y + step_info_kwargs['T'] = t + else: + sysdata = tsystem.sys + + info_dict = step_info(sysdata, **step_info_kwargs) + for i, row in enumerate(info_dict): for j, info in enumerate(row): - for k in info: - self.assert_step_info_match(tsystem.sys, - info, tsystem.step_info[i][j]) + self.assert_step_info_match(tsystem.sys, + info, tsystem.step_info[i][j]) + + def test_step_info_invalid(self): + """Call step_info with invalid parameters.""" + with pytest.raises(ValueError, match="time series data convention"): + step_info(["not numeric data"]) + with pytest.raises(ValueError, match="time series data convention"): + step_info(np.ones((10, 15))) # invalid shape + with pytest.raises(ValueError, match="matching time vector"): + step_info(np.ones(15), T=np.linspace(0, 1, 20)) # time too long + with pytest.raises(ValueError, match="matching time vector"): + step_info(np.ones((2, 2, 15))) # no time vector def test_step_pole_cancellation(self, pole_cancellation, no_pole_cancellation): @@ -491,7 +530,6 @@ def test_step_pole_cancellation(self, pole_cancellation, step_info_no_cancellation, step_info_cancellation) - @pytest.mark.parametrize( "tsystem, kwargs", [("siso_ss2", {}), diff --git a/control/timeresp.py b/control/timeresp.py index a68b25521..53144d7d2 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -1,9 +1,7 @@ -# timeresp.py - time-domain simulation routines -# -# This file contains a collection of functions that calculate time -# responses for linear systems. +""" +timeresp.py - time-domain simulation routines. -"""The :mod:`~control.timeresp` module contains a collection of +The :mod:`~control.timeresp` module contains a collection of functions that are used to compute time-domain simulations of LTI systems. @@ -21,9 +19,7 @@ See :ref:`time-series-convention` for more information on how time series data are represented. -""" - -"""Copyright (c) 2011 by California Institute of Technology +Copyright (c) 2011 by California Institute of Technology All rights reserved. Copyright (c) 2011 by Eike Welk @@ -75,12 +71,12 @@ import numpy as np import scipy as sp -from numpy import atleast_1d, einsum, maximum, minimum +from numpy import einsum, maximum, minimum from scipy.linalg import eig, eigvals, matrix_balance, norm from . import config -from .lti import (LTI, isctime, isdtime) -from .statesp import _convert_to_statespace, _mimo2simo, _mimo2siso, ssdata +from .lti import isctime, isdtime +from .statesp import StateSpace, _convert_to_statespace, _mimo2simo, _mimo2siso from .xferfcn import TransferFunction __all__ = ['forced_response', 'step_response', 'step_info', 'initial_response', @@ -209,7 +205,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, Parameters ---------- - sys : LTI (StateSpace or TransferFunction) + sys : StateSpace or TransferFunction LTI system to simulate T : array_like, optional for discrete LTI `sys` @@ -284,9 +280,9 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, See :ref:`time-series-convention`. """ - if not isinstance(sys, LTI): - raise TypeError('Parameter ``sys``: must be a ``LTI`` object. ' - '(For example ``StateSpace`` or ``TransferFunction``)') + if not isinstance(sys, (StateSpace, TransferFunction)): + raise TypeError('Parameter ``sys``: must be a ``StateSpace`` or' + ' ``TransferFunction``)') # If return_x was not specified, figure out the default if return_x is None: @@ -738,20 +734,24 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None, squeeze=squeeze, input=input, output=output) -def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, +def step_info(sysdata, T=None, T_num=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9)): """ Step response characteristics (Rise time, Settling Time, Peak and others). Parameters ---------- - sys : LTI system + sysdata : StateSpace or TransferFunction or array_like + The system data. Either LTI system to similate (StateSpace, + TransferFunction), or a time series of step response data. T : array_like or float, optional Time vector, or simulation time duration if a number (time vector is autocomputed if not given, see :func:`step_response` for more detail) + Required, if sysdata is a time series of response data. T_num : int, optional Number of time steps to use in simulation if T is not provided as an - array (autocomputed if not given); ignored if sys is discrete-time. + array; autocomputed if not given; ignored if sysdata is a + discrete-time system or a time series or response data. SettlingTimeThreshold : float value, optional Defines the error to compute settling time (default = 0.02) RiseTimeLimits : tuple (lower_threshold, upper_theshold) @@ -760,7 +760,8 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, Returns ------- S : dict or list of list of dict - If `sys` is a SISO system, S is a dictionary containing: + If `sysdata` corresponds to a SISO system, S is a dictionary + containing: RiseTime: Time from 10% to 90% of the steady-state value. @@ -781,9 +782,9 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, SteadyStateValue: Steady-state value - If `sys` is a MIMO system, `S` is a 2D list of dicts. To get the - step response characteristics from the j-th input to the i-th output, - access ``S[i][j]`` + If `sysdata` corresponds to a MIMO system, `S` is a 2D list of dicts. + To get the step response characteristics from the j-th input to the + i-th output, access ``S[i][j]`` See Also @@ -833,18 +834,46 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, PeakTime: 4.242 SteadyStateValue: -1.0 """ - 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, squeeze=False) + if isinstance(sysdata, (StateSpace, TransferFunction)): + if T is None or np.asarray(T).size == 1: + T = _default_time_vector(sysdata, N=T_num, tfinal=T, is_step=True) + T, Yout = step_response(sysdata, T, squeeze=False) + InfValues = np.atleast_2d(sysdata.dcgain()) + retsiso = sysdata.issiso() + noutputs = sysdata.noutputs + ninputs = sysdata.ninputs + else: + # Time series of response data + errmsg = ("`sys` must be a LTI system, or time response data" + " with a shape following the python-control" + " time series data convention.") + try: + Yout = np.array(sysdata, dtype=float) + except ValueError: + raise ValueError(errmsg) + if Yout.ndim == 1 or (Yout.ndim == 2 and Yout.shape[0] == 1): + Yout = Yout[np.newaxis, np.newaxis, :] + retsiso = True + elif Yout.ndim == 3: + retsiso = False + else: + raise ValueError(errmsg) + if T is None or Yout.shape[2] != len(np.squeeze(T)): + raise ValueError("For time response data, a matching time vector" + " must be given") + T = np.squeeze(T) + noutputs = Yout.shape[0] + ninputs = Yout.shape[1] + InfValues = Yout[:, :, -1] ret = [] - for i in range(sys.noutputs): + for i in range(noutputs): retrow = [] - for j in range(sys.ninputs): + for j in range(ninputs): yout = Yout[i, j, :] # Steady state value - InfValue = sys.dcgain() if sys.issiso() else sys.dcgain()[i, j] + InfValue = InfValues[i, j] sgnInf = np.sign(InfValue.real) rise_time: float = np.NaN @@ -869,7 +898,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, # SettlingTime settled = np.where( - np.abs(yout/InfValue -1) >= SettlingTimeThreshold)[0][-1]+1 + np.abs(yout/InfValue-1) >= SettlingTimeThreshold)[0][-1]+1 # MIMO systems can have unsettled channels without infinite # InfValue if settled < len(T): @@ -917,7 +946,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, ret.append(retrow) - return ret[0][0] if sys.issiso() else ret + return ret[0][0] if retsiso else ret def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None, transpose=False, return_x=False, squeeze=None): From 6587f6922bc5c8069e22480d5ae710c8bc77c68f Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Fri, 19 Mar 2021 15:37:55 +0100 Subject: [PATCH 12/13] add yfinal parameter to step_info --- control/tests/timeresp_test.py | 35 ++++++++++++++++++++++------------ control/timeresp.py | 16 ++++++++++++---- 2 files changed, 35 insertions(+), 16 deletions(-) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 6ec23e27e..a576d0903 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -451,6 +451,15 @@ def assert_step_info_match(self, sys, info, info_ref): # local min/max peak well after signal has risen np.testing.assert_allclose(info[k], info_ref[k], rtol=1e-3) + @pytest.mark.parametrize( + "yfinal", [True, False], ids=["yfinal", "no yfinal"]) + @pytest.mark.parametrize( + "systype, time_2d", + [("ltisys", False), + ("time response", False), + ("time response", True), + ], + ids=["ltisys", "time response (n,)", "time response (1,n)"]) @pytest.mark.parametrize( "tsystem", ["siso_tf_step_matlab", @@ -459,16 +468,10 @@ def assert_step_info_match(self, sys, info, info_ref): "siso_tf_kneg", "siso_tf_type1"], indirect=["tsystem"]) - @pytest.mark.parametrize( - "systype, time_2d", - [("lti", False), - ("time response data", False), - ("time response data", True), - ]) - def test_step_info(self, tsystem, systype, time_2d): + def test_step_info(self, tsystem, systype, time_2d, yfinal): """Test step info for SISO systems.""" step_info_kwargs = tsystem.kwargs.get('step_info', {}) - if systype == "time response data": + if systype == "time response": # simulate long enough for steady state value tfinal = 3 * tsystem.step_info['SettlingTime'] if np.isnan(tfinal): @@ -478,22 +481,26 @@ def test_step_info(self, tsystem, systype, time_2d): step_info_kwargs['T'] = t[np.newaxis, :] if time_2d else t else: sysdata = tsystem.sys + if yfinal: + step_info_kwargs['yfinal'] = tsystem.step_info['SteadyStateValue'] info = step_info(sysdata, **step_info_kwargs) self.assert_step_info_match(tsystem.sys, info, tsystem.step_info) + @pytest.mark.parametrize( + "yfinal", [True, False], ids=["yfinal", "no_yfinal"]) + @pytest.mark.parametrize( + "systype", ["ltisys", "time response"]) @pytest.mark.parametrize( "tsystem", ['mimo_ss_step_matlab', pytest.param('mimo_tf_step', marks=slycotonly)], indirect=["tsystem"]) - @pytest.mark.parametrize( - "systype", ["lti", "time response data"]) - def test_step_info_mimo(self, tsystem, systype): + def test_step_info_mimo(self, tsystem, systype, yfinal): """Test step info for MIMO systems.""" step_info_kwargs = tsystem.kwargs.get('step_info', {}) - if systype == "time response data": + if systype == "time response": tfinal = 3 * max([S['SettlingTime'] for Srow in tsystem.step_info for S in Srow]) t, y = step_response(tsystem.sys, T=tfinal, T_num=5000) @@ -501,6 +508,10 @@ def test_step_info_mimo(self, tsystem, systype): step_info_kwargs['T'] = t else: sysdata = tsystem.sys + if yfinal: + step_info_kwargs['yfinal'] = [[S['SteadyStateValue'] + for S in Srow] + for Srow in tsystem.step_info] info_dict = step_info(sysdata, **step_info_kwargs) diff --git a/control/timeresp.py b/control/timeresp.py index 53144d7d2..4f407f925 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -734,8 +734,8 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None, squeeze=squeeze, input=input, output=output) -def step_info(sysdata, T=None, T_num=None, SettlingTimeThreshold=0.02, - RiseTimeLimits=(0.1, 0.9)): +def step_info(sysdata, T=None, T_num=None, yfinal=None, + SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9)): """ Step response characteristics (Rise time, Settling Time, Peak and others). @@ -752,6 +752,11 @@ def step_info(sysdata, T=None, T_num=None, SettlingTimeThreshold=0.02, Number of time steps to use in simulation if T is not provided as an array; autocomputed if not given; ignored if sysdata is a discrete-time system or a time series or response data. + yfinal: scalar or array_like, optional + Steady-state response. If not given, sysdata.dcgain() is used for + systems to simulate and the last value of the the response data is + used for a given time series of response data. Scalar for SISO, + (noutputs, ninputs) array_like for MIMO systems. SettlingTimeThreshold : float value, optional Defines the error to compute settling time (default = 0.02) RiseTimeLimits : tuple (lower_threshold, upper_theshold) @@ -838,7 +843,10 @@ def step_info(sysdata, T=None, T_num=None, SettlingTimeThreshold=0.02, if T is None or np.asarray(T).size == 1: T = _default_time_vector(sysdata, N=T_num, tfinal=T, is_step=True) T, Yout = step_response(sysdata, T, squeeze=False) - InfValues = np.atleast_2d(sysdata.dcgain()) + if yfinal: + InfValues = np.atleast_2d(yfinal) + else: + InfValues = np.atleast_2d(sysdata.dcgain()) retsiso = sysdata.issiso() noutputs = sysdata.noutputs ninputs = sysdata.ninputs @@ -864,7 +872,7 @@ def step_info(sysdata, T=None, T_num=None, SettlingTimeThreshold=0.02, T = np.squeeze(T) noutputs = Yout.shape[0] ninputs = Yout.shape[1] - InfValues = Yout[:, :, -1] + InfValues = np.atleast_2d(yfinal) if yfinal else Yout[:, :, -1] ret = [] for i in range(noutputs): From a878846400f8fbb29c2babaa02bae461ae7e40ff Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Fri, 19 Mar 2021 15:53:58 +0100 Subject: [PATCH 13/13] include yfinal in matlab.stepinfo call signature --- control/matlab/timeresp.py | 69 +++++++++++++++++++++++++------------- control/timeresp.py | 6 ++-- 2 files changed, 49 insertions(+), 26 deletions(-) diff --git a/control/matlab/timeresp.py b/control/matlab/timeresp.py index 31b761bcd..b1fa24bb0 100644 --- a/control/matlab/timeresp.py +++ b/control/matlab/timeresp.py @@ -64,38 +64,59 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False): transpose=True, return_x=return_x) return (out[1], out[0], out[2]) if return_x else (out[1], out[0]) -def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, + +def stepinfo(sysdata, T=None, yfinal=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 or number, optional + sysdata : StateSpace or TransferFunction or array_like + The system data. Either LTI system to similate (StateSpace, + TransferFunction), or a time series of step response data. + T : array_like or float, optional Time vector, or simulation time duration if a number (time vector is - autocomputed if not given) - - SettlingTimeThreshold: float value, optional + autocomputed if not given). + Required, if sysdata is a time series of response data. + yfinal : scalar or array_like, optional + Steady-state response. If not given, sysdata.dcgain() is used for + systems to simulate and the last value of the the response data is + used for a given time series of response data. Scalar for SISO, + (noutputs, ninputs) array_like for MIMO systems. + SettlingTimeThreshold : float, optional Defines the error to compute settling time (default = 0.02) - - RiseTimeLimits: tuple (lower_threshold, upper_theshold) + 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 + S : dict or list of list of dict + If `sysdata` corresponds to a SISO system, S is 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 + + If `sysdata` corresponds to a MIMO system, `S` is a 2D list of dicts. + To get the step response characteristics from the j-th input to the + i-th output, access ``S[i][j]`` See Also @@ -105,11 +126,13 @@ def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, Examples -------- >>> S = stepinfo(sys, T) - ''' + """ from ..timeresp import step_info # Call step_info with MATLAB defaults - S = step_info(sys, T, None, SettlingTimeThreshold, RiseTimeLimits) + S = step_info(sysdata, T=T, T_num=None, yfinal=yfinal, + SettlingTimeThreshold=SettlingTimeThreshold, + RiseTimeLimits=RiseTimeLimits) return S diff --git a/control/timeresp.py b/control/timeresp.py index 4f407f925..630eff03a 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -746,18 +746,18 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None, TransferFunction), or a time series of step response data. T : array_like or float, optional Time vector, or simulation time duration if a number (time vector is - autocomputed if not given, see :func:`step_response` for more detail) + autocomputed if not given, see :func:`step_response` for more detail). Required, if sysdata is a time series of response data. T_num : int, optional Number of time steps to use in simulation if T is not provided as an array; autocomputed if not given; ignored if sysdata is a discrete-time system or a time series or response data. - yfinal: scalar or array_like, optional + yfinal : scalar or array_like, optional Steady-state response. If not given, sysdata.dcgain() is used for systems to simulate and the last value of the the response data is used for a given time series of response data. Scalar for SISO, (noutputs, ninputs) array_like for MIMO systems. - SettlingTimeThreshold : float value, optional + SettlingTimeThreshold : float, 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