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
merged 4 commits into from
Mar 17, 2021
Merged

Step info improve jpp #567

Show file tree
Hide file tree
Changes from all 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
# 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
Expand Down
96 changes: 68 additions & 28 deletions control/timeresp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -785,43 +786,82 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
--------
>>> info = step_info(sys, T)
'''
_, sys = _get_ss_simo(sys)

if not sys.issiso():
sys = _mimo2siso(sys,0,0)
warnings.warn(" Internal conversion from a MIMO system to a SISO system,"
" the first input and the first output were used (u1 -> y1);"
" it may not be the result you are looking for")

if T is None or np.asarray(T).size == 1:
T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True)

T, yout = step_response(sys, T)

# Steady state value
InfValue = sys.dcgain()

# RiseTime
tr_lower_index = (np.where(yout >= RiseTimeLimits[0] * InfValue)[0])[0]
tr_upper_index = (np.where(yout >= RiseTimeLimits[1] * InfValue)[0])[0]
RiseTime = T[tr_upper_index] - T[tr_lower_index]

# SettlingTime
sup_margin = (1. + SettlingTimeThreshold) * InfValue
inf_margin = (1. - SettlingTimeThreshold) * InfValue
# find Steady State looking for the first point out of specified limits
for i in reversed(range(T.size)):
if((yout[i] <= inf_margin) | (yout[i] >= sup_margin)):
SettlingTime = T[i + 1]
break
InfValue = sys.dcgain().real

# TODO: this could be a function step_info4data(t,y,yfinal)
rise_time: float = np.NaN
settling_time: float = np.NaN
settling_min: float = np.NaN
settling_max: float = np.NaN
peak_value: float = np.Inf
peak_time: float = np.Inf
undershoot: float = np.NaN
overshoot: float = np.NaN
Comment on lines +805 to +812
Copy link
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?

Copy link
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?

Copy link
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.

steady_state_value: float = np.NaN

if not np.isnan(InfValue) and not np.isinf(InfValue):
# SteadyStateValue
steady_state_value = InfValue
# Peak
peak_index = np.abs(yout).argmax()
peak_value = np.abs(yout[peak_index])
peak_time = T[peak_index]

sup_margin = (1. + SettlingTimeThreshold) * InfValue
inf_margin = (1. - SettlingTimeThreshold) * InfValue

# RiseTime
tr_lower_index = (np.where(np.sign(InfValue.real) * (yout- RiseTimeLimits[0] * InfValue) >= 0 )[0])[0]
tr_upper_index = (np.where(np.sign(InfValue.real) * yout >= np.sign(InfValue.real) * RiseTimeLimits[1] * InfValue)[0])[0]

# SettlingTime
settling_time = T[np.where(np.abs(yout-InfValue) >= np.abs(SettlingTimeThreshold*InfValue))[0][-1]+1]
# Overshoot and Undershoot
y_os = (np.sign(InfValue.real)*yout).max()
dy_os = np.abs(y_os) - np.abs(InfValue)
if dy_os > 0:
overshoot = np.abs(100. * dy_os / InfValue)
else:
overshoot = 0

y_us = (np.sign(InfValue.real)*yout).min()
dy_us = np.abs(y_us)
if dy_us > 0:
undershoot = np.abs(100. * dy_us / InfValue)
else:
undershoot = 0

# RiseTime
rise_time = T[tr_upper_index] - T[tr_lower_index]

settling_max = (yout[tr_upper_index:]).max()
settling_min = (yout[tr_upper_index:]).min()

PeakIndex = np.abs(yout).argmax()
return {
'RiseTime': RiseTime,
'SettlingTime': SettlingTime,
'SettlingMin': yout[tr_upper_index:].min(),
'SettlingMax': yout.max(),
'Overshoot': 100. * (yout.max() - InfValue) / InfValue,
'Undershoot': yout.min(), # not very confident about this
'Peak': yout[PeakIndex],
'PeakTime': T[PeakIndex],
'SteadyStateValue': InfValue
'RiseTime': rise_time,
'SettlingTime': settling_time,
'SettlingMin': settling_min,
'SettlingMax': settling_max,
'Overshoot': overshoot,
'Undershoot': undershoot,
'Peak': peak_value,
'PeakTime': peak_time,
'SteadyStateValue': steady_state_value
}


def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
transpose=False, return_x=False, squeeze=None):
# pylint: disable=W0622
Expand Down