Skip to content

Step info improve jpp #583

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

Closed
wants to merge 2 commits into from
Closed

Step info improve jpp #583

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f7d3db5
new feature and bug fixes:
juanodecc Mar 18, 2021
8d5810f
modify the help according to the new characteristics
juanodecc Mar 18, 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
122 changes: 121 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 @@ -222,6 +222,15 @@ def tf2_matlab_help(self):
D = [0.06]
sys = StateSpace(A, B, C, D, 0.2)
return sys

@pytest.fixture
def mimo_matlab_help(self):
A = [[0.68, -0.34], [0.34, 0.68]];
B = [[0.18, -0.05], [0.04, 0.11]];
C = [[0, -1.53], [-1.12, -1.10]];
D = [[0, 0], [0.06, -0.37]];
sys = StateSpace(A,B,C,D,0.2);
return sys

@pytest.fixture
def tsystem(self,
Expand All @@ -233,7 +242,7 @@ def tsystem(self,
mimo_dss1, mimo_dss2, mimo_dtf1,
pole_cancellation, no_pole_cancellation, siso_tf_type1,
siso_tf_kpos, siso_tf_kneg, tf1_matlab_help,
tf2_matlab_help):
tf2_matlab_help, mimo_matlab_help):
systems = {"siso_ss1": siso_ss1,
"siso_ss2": siso_ss2,
"siso_tf1": siso_tf1,
Expand All @@ -256,6 +265,7 @@ def tsystem(self,
"siso_tf_kneg": siso_tf_kneg,
"tf1_matlab_help": tf1_matlab_help,
"tf2_matlab_help": tf2_matlab_help,
"mimo_matlab_help":mimo_matlab_help,
}
return systems[request.param]

Expand Down Expand Up @@ -341,6 +351,116 @@ def test_step_info(self):

# tolerance for all parameters could be wrong for some systems
# discrete systems time parameters tolerance could be +/-dt
def test_step_info(self):
# From matlab docs:
sys = TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2])
Strue = {
'RiseTime': 3.8456,
'SettlingTime': 27.9762,
'SettlingMin': 2.0689,
'SettlingMax': 2.6873,
'Overshoot': 7.4915,
'Undershoot': 0,
'Peak': 2.6873,
'PeakTime': 8.0530,
'SteadyStateValue': 2.50
}

S = step_info(sys)

Sk = sorted(S.keys())
Sktrue = sorted(Strue.keys())
assert Sk == Sktrue
# Very arbitrary tolerance because I don't know if the
# response from the MATLAB is really that accurate.
# maybe it is a good idea to change the Strue to match
# but I didn't do it because I don't know if it is
# accurate either...
rtol = 2e-2
np.testing.assert_allclose([S[k] for k in Sk],
[Strue[k] for k in Sktrue],
rtol=rtol)

def test_step_info_mimo(self):
# From matlab docs:
A = [[0.68, -0.34], [0.34, 0.68]];
B = [[0.18, -0.05], [0.04, 0.11]];
C = [[0, -1.53], [-1.12, -1.10]];
D = [[0, 0], [0.06, -0.37]];
sys = StateSpace(A,B,C,D,0.2);

Strue = [[{'RiseTime': 0.6000,
'SettlingTime': 3.0,
'SettlingMin': -0.5999,
'SettlingMax': -0.4689,
'Overshoot': 15.5072,
'Undershoot': 0,
'Peak': 0.5999,
'PeakTime': 1.4000,
'SteadyStateValue': -0.51935},
{'RiseTime': 0.0,
'SettlingTime': 3.6000,
'SettlingMin': -0.2797,
'SettlingMax': -0.1043,
'Overshoot': 118.9918,
'Undershoot': 0,
'Peak': 0.2797,
'PeakTime': 0.60000,
'SteadyStateValue': -0.1277}],
[{'RiseTime': 0.40000,
'SettlingTime': 2.8,
'SettlingMin': -0.6724,
'SettlingMax': -0.5188,
'Overshoot': 24.6476,
'Undershoot': 11.1224,
'Peak': 0.6724,
'PeakTime': 1.0,
'SteadyStateValue': -0.5394},
{'RiseTime': 0.0,
'SettlingTime': 3.4,
'SettlingMin': -0.435, #this value is not a result of step_info
'SettlingMax': -0.1485,
'Overshoot': 132.0170,
'Undershoot': 0,
'Peak': 0.4350,
'PeakTime': 0.2,
'SteadyStateValue': -0.18748}]]

S = step_info(sys)

S1k = sorted(S[0][0].keys())
S1ktrue = sorted(Strue[0][0].keys())
assert S1k == S1ktrue

S2k = sorted(S[0][1].keys())
S2ktrue = sorted(Strue[0][1].keys())
assert S2k == S2ktrue

S3k = sorted(S[1][0].keys())
S3ktrue = sorted(Strue[1][0].keys())
assert S3k == S3ktrue

S4k = sorted(S[1][1].keys())
S4ktrue = sorted(Strue[1][1].keys())
assert S4k == S4ktrue

rtol = 0.2 # tolerance could be better was chosen to meet the discret times parameters
np.testing.assert_allclose([S[0][0][k] for k in S1k],
[Strue[0][0][k] for k in S1ktrue],
rtol=rtol)
np.testing.assert_allclose([S[0][1][k] for k in S2k],
[Strue[0][1][k] for k in S2ktrue],
rtol=rtol)
np.testing.assert_allclose([S[1][0][k] for k in S3k],
[Strue[1][0][k] for k in S3ktrue],
rtol=rtol)
np.testing.assert_allclose([S[1][1][k] for k in S4k],
[Strue[1][1][k] for k in S4ktrue],
rtol=rtol)

# Tolerance for all parameters could be wrong for some systems
# discrete systems time parameters tolerance
# could be +/-dt
@pytest.mark.parametrize(
"tsystem, info_true, tolerance",
[("tf1_matlab_help", {
Expand Down
208 changes: 127 additions & 81 deletions control/timeresp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -746,7 +746,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,

Parameters
----------
sys : SISO dynamic system model. Dynamic systems that you can use include:
sys : SISO or MIMO dynamic system model. Dynamic systems that you can use include:
StateSpace or TransferFunction
LTI system to simulate

Expand All @@ -766,101 +766,147 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,

Returns
-------
S: a dictionary containing:
RiseTime: Time from 10% to 90% of the steady-state value.
S: a dictionary for SISO systems or a list of Ny-by-Nu of dictionaries
for MIMO systems, each dictionary containing:

RiseTime: Rise time (by default time from 10% to 90% of the steady-state value).
SettlingTime: Time to enter inside a default error of 2%
SettlingMin: Minimum value after RiseTime
SettlingMax: Maximum value after RiseTime
SettlingMin: Minimum value of `y(t)` after RiseTime
SettlingMax: Maximum value of `y(t)` after RiseTime
Overshoot: Percentage of the Peak relative to steady value
Undershoot: Percentage of undershoot
Peak: Absolute peak value
PeakTime: time of the Peak
Undershoot: Percentage of undershoot
Peak: Peak absolute value of `y(t)`
PeakTime: Time of the peak value
SteadyStateValue: Steady-state value


See Also
--------
step, lsim, initial, impulse
step_response, initial_response, impulse_response, forced_response


Examples
--------
>>> info = step_info(sys, T)
>>> from control import tf, step_info
>>> sys = tf([3.32, 0, -162.8],[1, 24.56, 186.5, 457.8, 116.2])
>>> step_info(sys)
>>>
>>> {'RiseTime': 7.70245002940221,
>>> 'SettlingTime': 14.13151114265325,
>>> 'SettlingMin': -1.3994264225116284,
>>> 'SettlingMax': -1.2612640212267976,
>>> 'Overshoot': 0,
>>> 'Undershoot': 0.6864433321178778,
>>> 'Peak': 1.3994264225116284,
>>> 'PeakTime': 24.13227287437709,
>>> 'SteadyStateValue': -1.4010327022375215}

'''

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

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

T, yout = step_response(sys, T)
T, y = step_response(sys, T)
else:
T, y = step_response(sys, T)

# Steady state value
InfValue = sys.dcgain().real

# TODO: this could be a function step_info4data(t,y,yfinal)
rise_time: float = np.NaN
settling_time: float = np.NaN
settling_min: float = np.NaN
settling_max: float = np.NaN
peak_value: float = np.Inf
peak_time: float = np.Inf
undershoot: float = np.NaN
overshoot: float = np.NaN
steady_state_value: float = np.NaN

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

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

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

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

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

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

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

return {
'RiseTime': rise_time,
'SettlingTime': settling_time,
'SettlingMin': settling_min,
'SettlingMax': settling_max,
'Overshoot': overshoot,
'Undershoot': undershoot,
'Peak': peak_value,
'PeakTime': peak_time,
'SteadyStateValue': steady_state_value
}
yfinal = sys.dcgain().real

# TODO: this could be a function step_info4data(t,y,yfinal)
# y: Step-response data, specified as:
# For SISO response data, a vector of length Ns, where Ns is the number of samples in the response data.
# For MIMO response data, an Ny-by-Nu-by-Ns array, where Ny is the number of system outputs, and Nu is the number of system inputs.
# T: Time vector corresponding to the response data in y, specified as a vector of length Ns.

# TODO: check dimentions for time vs step output and final output
if len(y.shape) == 3: # MIMO
Ny,Nu,Ns = y.shape
elif len(y.shape) == 1: # SISO
Ns = y.shape
Ny = 1
Nu = 1
else:
raise TypeError("Poner un mensaje")

S = [[[] for x in range(Ny)] for x in range(Nu)]

# Struct for each couple of inputs
for i in range(Ny): # Ny
for j in range(Nu): # Nu
rise_time: float = np.NaN
settling_time: float = np.NaN
settling_min: float = np.NaN
settling_max: float = np.NaN
peak_value: float = np.Inf
peak_time: float = np.Inf
undershoot: float = np.NaN
overshoot: float = np.NaN
steady_state_value: float = np.NaN

if len(y.shape) > 1:
yout = y[i,j,:]
InfValue = yfinal[i,j]
else:
yout = y
InfValue = yfinal

InfValue_sign = np.sign(InfValue)

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

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

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

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

y_us = (InfValue_sign*yout).min()
y_us_index = (InfValue_sign*yout).argmin()
if (InfValue_sign * yout[y_us_index]) < 0: # must have oposite sign
undershoot = np.abs(100. * np.abs(y_us) / InfValue)
else:
undershoot = 0

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

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

S[i][j] = {
'RiseTime': rise_time,
'SettlingTime': settling_time,
'SettlingMin': settling_min,
'SettlingMax': settling_max,
'Overshoot': overshoot,
'Undershoot': undershoot,
'Peak': peak_value,
'PeakTime': peak_time,
'SteadyStateValue': steady_state_value
}
if Ny > 1 or Nu > 1:
return S
else:
return S[0][0]

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