Skip to content

Commit 258c9c2

Browse files
bnavigatorbnavigator
authored
Merge pull request #567 from juanodecc/step_info_improve_jpp
Step info improve jpp
2 parents 0f951e1 + 81ae64f commit 258c9c2

File tree

2 files changed

+172
-29
lines changed

2 files changed

+172
-29
lines changed

‎control/tests/timeresp_test.py

Lines changed: 104 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -193,6 +193,35 @@ def pole_cancellation(self):
193193
def no_pole_cancellation(self):
194194
return TransferFunction([1.881e+06],
195195
[188.1, 1.881e+06])
196+
197+
@pytest.fixture
198+
def siso_tf_type1(self):
199+
# System Type 1 - Step response not stationary: G(s)=1/s(s+1)
200+
return TransferFunction(1, [1, 1, 0])
201+
202+
@pytest.fixture
203+
def siso_tf_kpos(self):
204+
# SISO under shoot response and positive final value G(s)=(-s+1)/(s²+s+1)
205+
return TransferFunction([-1, 1], [1, 1, 1])
206+
207+
@pytest.fixture
208+
def siso_tf_kneg(self):
209+
# SISO under shoot response and negative final value k=-1 G(s)=-(-s+1)/(s²+s+1)
210+
return TransferFunction([1, -1], [1, 1, 1])
211+
212+
@pytest.fixture
213+
def tf1_matlab_help(self):
214+
# example from matlab online help https://www.mathworks.com/help/control/ref/stepinfo.html
215+
return TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2])
216+
217+
@pytest.fixture
218+
def tf2_matlab_help(self):
219+
A = [[0.68, - 0.34], [0.34, 0.68]]
220+
B = [[0.18], [0.04]]
221+
C = [-1.12, - 1.10]
222+
D = [0.06]
223+
sys = StateSpace(A, B, C, D, 0.2)
224+
return sys
196225

197226
@pytest.fixture
198227
def tsystem(self,
@@ -202,7 +231,9 @@ def tsystem(self,
202231
siso_dtf0, siso_dtf1, siso_dtf2,
203232
siso_dss1, siso_dss2,
204233
mimo_dss1, mimo_dss2, mimo_dtf1,
205-
pole_cancellation, no_pole_cancellation):
234+
pole_cancellation, no_pole_cancellation, siso_tf_type1,
235+
siso_tf_kpos, siso_tf_kneg, tf1_matlab_help,
236+
tf2_matlab_help):
206237
systems = {"siso_ss1": siso_ss1,
207238
"siso_ss2": siso_ss2,
208239
"siso_tf1": siso_tf1,
@@ -220,6 +251,11 @@ def tsystem(self,
220251
"mimo_dtf1": mimo_dtf1,
221252
"pole_cancellation": pole_cancellation,
222253
"no_pole_cancellation": no_pole_cancellation,
254+
"siso_tf_type1": siso_tf_type1,
255+
"siso_tf_kpos": siso_tf_kpos,
256+
"siso_tf_kneg": siso_tf_kneg,
257+
"tf1_matlab_help": tf1_matlab_help,
258+
"tf2_matlab_help": tf2_matlab_help,
223259
}
224260
return systems[request.param]
225261

@@ -303,6 +339,73 @@ def test_step_info(self):
303339
[Strue[k] for k in Sktrue],
304340
rtol=rtol)
305341

342+
# tolerance for all parameters could be wrong for some systems
343+
# discrete systems time parameters tolerance could be +/-dt
344+
@pytest.mark.parametrize(
345+
"tsystem, info_true, tolerance",
346+
[("tf1_matlab_help", {
347+
'RiseTime': 3.8456,
348+
'SettlingTime': 27.9762,
349+
'SettlingMin': 2.0689,
350+
'SettlingMax': 2.6873,
351+
'Overshoot': 7.4915,
352+
'Undershoot': 0,
353+
'Peak': 2.6873,
354+
'PeakTime': 8.0530,
355+
'SteadyStateValue': 2.5}, 2e-2),
356+
("tf2_matlab_help", {
357+
'RiseTime': 0.4000,
358+
'SettlingTime': 2.8000,
359+
'SettlingMin': -0.6724,
360+
'SettlingMax': -0.5188,
361+
'Overshoot': 24.6476,
362+
'Undershoot': 11.1224,
363+
'Peak': 0.6724,
364+
'PeakTime': 1,
365+
'SteadyStateValue': -0.5394}, .2),
366+
("siso_tf_kpos", {
367+
'RiseTime': 1.242,
368+
'SettlingTime': 9.110,
369+
'SettlingMin': 0.950,
370+
'SettlingMax': 1.208,
371+
'Overshoot': 20.840,
372+
'Undershoot': 27.840,
373+
'Peak': 1.208,
374+
'PeakTime': 4.282,
375+
'SteadyStateValue': 1.0}, 2e-2),
376+
("siso_tf_kneg", {
377+
'RiseTime': 1.242,
378+
'SettlingTime': 9.110,
379+
'SettlingMin': -1.208,
380+
'SettlingMax': -0.950,
381+
'Overshoot': 20.840,
382+
'Undershoot': 27.840,
383+
'Peak': 1.208,
384+
'PeakTime': 4.282,
385+
'SteadyStateValue': -1.0}, 2e-2),
386+
("siso_tf_type1", {'RiseTime': np.NaN,
387+
'SettlingTime': np.NaN,
388+
'SettlingMin': np.NaN,
389+
'SettlingMax': np.NaN,
390+
'Overshoot': np.NaN,
391+
'Undershoot': np.NaN,
392+
'Peak': np.Inf,
393+
'PeakTime': np.Inf,
394+
'SteadyStateValue': np.NaN}, 2e-2)],
395+
indirect=["tsystem"])
396+
def test_step_info(self, tsystem, info_true, tolerance):
397+
""" """
398+
info = step_info(tsystem)
399+
400+
info_true_sorted = sorted(info_true.keys())
401+
info_sorted = sorted(info.keys())
402+
403+
assert info_sorted == info_true_sorted
404+
405+
np.testing.assert_allclose([info_true[k] for k in info_true_sorted],
406+
[info[k] for k in info_sorted],
407+
rtol=tolerance)
408+
306409
def test_step_pole_cancellation(self, pole_cancellation,
307410
no_pole_cancellation):
308411
# confirm that pole-zero cancellation doesn't perturb results

‎control/timeresp.py

Lines changed: 68 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -746,7 +746,8 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
746746
747747
Parameters
748748
----------
749-
sys : StateSpace or TransferFunction
749+
sys : SISO dynamic system model. Dynamic systems that you can use include:
750+
StateSpace or TransferFunction
750751
LTI system to simulate
751752
752753
T : array_like or float, optional
@@ -785,43 +786,82 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
785786
--------
786787
>>> info = step_info(sys, T)
787788
'''
788-
_, sys = _get_ss_simo(sys)
789+
790+
if not sys.issiso():
791+
sys = _mimo2siso(sys,0,0)
792+
warnings.warn(" Internal conversion from a MIMO system to a SISO system,"
793+
" the first input and the first output were used (u1 -> y1);"
794+
" it may not be the result you are looking for")
795+
789796
if T is None or np.asarray(T).size == 1:
790797
T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True)
791798

792799
T, yout = step_response(sys, T)
793800

794801
# Steady state value
795-
InfValue = sys.dcgain()
796-
797-
# RiseTime
798-
tr_lower_index = (np.where(yout >= RiseTimeLimits[0] * InfValue)[0])[0]
799-
tr_upper_index = (np.where(yout >= RiseTimeLimits[1] * InfValue)[0])[0]
800-
RiseTime = T[tr_upper_index] - T[tr_lower_index]
801-
802-
# SettlingTime
803-
sup_margin = (1. + SettlingTimeThreshold) * InfValue
804-
inf_margin = (1. - SettlingTimeThreshold) * InfValue
805-
# find Steady State looking for the first point out of specified limits
806-
for i in reversed(range(T.size)):
807-
if((yout[i] <= inf_margin) | (yout[i] >= sup_margin)):
808-
SettlingTime = T[i + 1]
809-
break
802+
InfValue = sys.dcgain().real
803+
804+
# TODO: this could be a function step_info4data(t,y,yfinal)
805+
rise_time: float = np.NaN
806+
settling_time: float = np.NaN
807+
settling_min: float = np.NaN
808+
settling_max: float = np.NaN
809+
peak_value: float = np.Inf
810+
peak_time: float = np.Inf
811+
undershoot: float = np.NaN
812+
overshoot: float = np.NaN
813+
steady_state_value: float = np.NaN
814+
815+
if not np.isnan(InfValue) and not np.isinf(InfValue):
816+
# SteadyStateValue
817+
steady_state_value = InfValue
818+
# Peak
819+
peak_index = np.abs(yout).argmax()
820+
peak_value = np.abs(yout[peak_index])
821+
peak_time = T[peak_index]
822+
823+
sup_margin = (1. + SettlingTimeThreshold) * InfValue
824+
inf_margin = (1. - SettlingTimeThreshold) * InfValue
825+
826+
# RiseTime
827+
tr_lower_index = (np.where(np.sign(InfValue.real) * (yout- RiseTimeLimits[0] * InfValue) >= 0 )[0])[0]
828+
tr_upper_index = (np.where(np.sign(InfValue.real) * yout >= np.sign(InfValue.real) * RiseTimeLimits[1] * InfValue)[0])[0]
829+
830+
# SettlingTime
831+
settling_time = T[np.where(np.abs(yout-InfValue) >= np.abs(SettlingTimeThreshold*InfValue))[0][-1]+1]
832+
# Overshoot and Undershoot
833+
y_os = (np.sign(InfValue.real)*yout).max()
834+
dy_os = np.abs(y_os) - np.abs(InfValue)
835+
if dy_os > 0:
836+
overshoot = np.abs(100. * dy_os / InfValue)
837+
else:
838+
overshoot = 0
839+
840+
y_us = (np.sign(InfValue.real)*yout).min()
841+
dy_us = np.abs(y_us)
842+
if dy_us > 0:
843+
undershoot = np.abs(100. * dy_us / InfValue)
844+
else:
845+
undershoot = 0
846+
847+
# RiseTime
848+
rise_time = T[tr_upper_index] - T[tr_lower_index]
849+
850+
settling_max = (yout[tr_upper_index:]).max()
851+
settling_min = (yout[tr_upper_index:]).min()
810852

811-
PeakIndex = np.abs(yout).argmax()
812853
return {
813-
'RiseTime': RiseTime,
814-
'SettlingTime': SettlingTime,
815-
'SettlingMin': yout[tr_upper_index:].min(),
816-
'SettlingMax': yout.max(),
817-
'Overshoot': 100. * (yout.max() - InfValue) / InfValue,
818-
'Undershoot': yout.min(), # not very confident about this
819-
'Peak': yout[PeakIndex],
820-
'PeakTime': T[PeakIndex],
821-
'SteadyStateValue': InfValue
854+
'RiseTime': rise_time,
855+
'SettlingTime': settling_time,
856+
'SettlingMin': settling_min,
857+
'SettlingMax': settling_max,
858+
'Overshoot': overshoot,
859+
'Undershoot': undershoot,
860+
'Peak': peak_value,
861+
'PeakTime': peak_time,
862+
'SteadyStateValue': steady_state_value
822863
}
823864

824-
825865
def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
826866
transpose=False, return_x=False, squeeze=None):
827867
# pylint: disable=W0622

0 commit comments

Comments
 (0)