Skip to content

Step info improve jpp #567

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
bnavigator merged 4 commits into from
Mar 17, 2021
Merged

Step info improve jpp #567

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1263550
Added 5 test for step_info:
juanodecc Mar 12, 2021
f964170
Solve issue #337, #565 and #564
juanodecc Mar 12, 2021
018d128
optimize the code and solve problems with MIMO
juanodecc Mar 16, 2021
81ae64f
fix comment format
bnavigator Mar 17, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
105 changes: 104 additions & 1 deletion control/tests/timeresp_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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,
Expand All @@ -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]

Expand Down Expand Up @@ -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
Comment thread
bnavigator marked this conversation as resolved.
Resolved
@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
Expand Down
79 changes: 56 additions & 23 deletions control/timeresp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Comment on lines +805 to +812

@bnavigator bnavigator Mar 12, 2021

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The first use of typing in python-control?

@juanodecc juanodecc Mar 13, 2021

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, is wrong?

@bnavigator bnavigator Mar 13, 2021

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not at all! Good thing we already ditched Python 2.


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)

@bnavigator bnavigator Mar 12, 2021 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, this has a lot of duplicate code. Have you considered using a sign variable and the abs() function instead of two completely separate blocks? I think It won't turn out much more convoluted than the current InfValue < 0 part.

@juanodecc juanodecc Mar 13, 2021

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok


# 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
Expand Down