Skip to content

Commit 43d73f0

Browse files
bnavigatorbnavigator
committed
support time series of response data in step_info
1 parent 20ed368 commit 43d73f0

File tree

2 files changed

+111
-44
lines changed

2 files changed

+111
-44
lines changed

‎control/tests/timeresp_test.py

Lines changed: 52 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -427,8 +427,7 @@ def test_step_nostates(self, dt):
427427
np.testing.assert_array_equal(y, np.ones(len(t)))
428428

429429
def assert_step_info_match(self, sys, info, info_ref):
430-
"""Assert reasonable step_info accuracy"""
431-
430+
"""Assert reasonable step_info accuracy."""
432431
if sys.isdtime(strict=True):
433432
dt = sys.dt
434433
else:
@@ -460,26 +459,66 @@ def assert_step_info_match(self, sys, info, info_ref):
460459
"siso_tf_kneg",
461460
"siso_tf_type1"],
462461
indirect=["tsystem"])
463-
def test_step_info(self, tsystem):
464-
"""Test step info for SISO systems"""
465-
step_info_kwargs = tsystem.kwargs.get('step_info',{})
466-
info = step_info(tsystem.sys, **step_info_kwargs)
462+
@pytest.mark.parametrize(
463+
"systype, time_2d",
464+
[("lti", False),
465+
("time response data", False),
466+
("time response data", True),
467+
])
468+
def test_step_info(self, tsystem, systype, time_2d):
469+
"""Test step info for SISO systems."""
470+
step_info_kwargs = tsystem.kwargs.get('step_info', {})
471+
if systype == "time response data":
472+
# simulate long enough for steady state value
473+
tfinal = 3 * tsystem.step_info['SettlingTime']
474+
if np.isnan(tfinal):
475+
pytest.skip("test system does not settle")
476+
t, y = step_response(tsystem.sys, T=tfinal, T_num=5000)
477+
sysdata = y
478+
step_info_kwargs['T'] = t[np.newaxis, :] if time_2d else t
479+
else:
480+
sysdata = tsystem.sys
481+
482+
info = step_info(sysdata, **step_info_kwargs)
483+
467484
self.assert_step_info_match(tsystem.sys, info, tsystem.step_info)
468485

469486
@pytest.mark.parametrize(
470487
"tsystem",
471488
['mimo_ss_step_matlab',
472489
pytest.param('mimo_tf_step', marks=slycotonly)],
473490
indirect=["tsystem"])
474-
def test_step_info_mimo(self, tsystem):
475-
"""Test step info for MIMO systems"""
476-
step_info_kwargs = tsystem.kwargs.get('step_info',{})
477-
info_dict = step_info(tsystem.sys, **step_info_kwargs)
491+
@pytest.mark.parametrize(
492+
"systype", ["lti", "time response data"])
493+
def test_step_info_mimo(self, tsystem, systype):
494+
"""Test step info for MIMO systems."""
495+
step_info_kwargs = tsystem.kwargs.get('step_info', {})
496+
if systype == "time response data":
497+
tfinal = 3 * max([S['SettlingTime']
498+
for Srow in tsystem.step_info for S in Srow])
499+
t, y = step_response(tsystem.sys, T=tfinal, T_num=5000)
500+
sysdata = y
501+
step_info_kwargs['T'] = t
502+
else:
503+
sysdata = tsystem.sys
504+
505+
info_dict = step_info(sysdata, **step_info_kwargs)
506+
478507
for i, row in enumerate(info_dict):
479508
for j, info in enumerate(row):
480-
for k in info:
481-
self.assert_step_info_match(tsystem.sys,
482-
info, tsystem.step_info[i][j])
509+
self.assert_step_info_match(tsystem.sys,
510+
info, tsystem.step_info[i][j])
511+
512+
def test_step_info_invalid(self):
513+
"""Call step_info with invalid parameters."""
514+
with pytest.raises(ValueError, match="time series data convention"):
515+
step_info(["not numeric data"])
516+
with pytest.raises(ValueError, match="time series data convention"):
517+
step_info(np.ones((10, 15))) # invalid shape
518+
with pytest.raises(ValueError, match="matching time vector"):
519+
step_info(np.ones(15), T=np.linspace(0, 1, 20)) # time too long
520+
with pytest.raises(ValueError, match="matching time vector"):
521+
step_info(np.ones((2, 2, 15))) # no time vector
483522

484523
def test_step_pole_cancellation(self, pole_cancellation,
485524
no_pole_cancellation):
@@ -491,7 +530,6 @@ def test_step_pole_cancellation(self, pole_cancellation,
491530
step_info_no_cancellation,
492531
step_info_cancellation)
493532

494-
495533
@pytest.mark.parametrize(
496534
"tsystem, kwargs",
497535
[("siso_ss2", {}),

‎control/timeresp.py

Lines changed: 59 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,7 @@
1-
# timeresp.py - time-domain simulation routines
2-
#
3-
# This file contains a collection of functions that calculate time
4-
# responses for linear systems.
1+
"""
2+
timeresp.py - time-domain simulation routines.
53
6-
"""The :mod:`~control.timeresp` module contains a collection of
4+
The :mod:`~control.timeresp` module contains a collection of
75
functions that are used to compute time-domain simulations of LTI
86
systems.
97
@@ -21,9 +19,7 @@
2119
See :ref:`time-series-convention` for more information on how time
2220
series data are represented.
2321
24-
"""
25-
26-
"""Copyright (c) 2011 by California Institute of Technology
22+
Copyright (c) 2011 by California Institute of Technology
2723
All rights reserved.
2824
2925
Copyright (c) 2011 by Eike Welk
@@ -75,12 +71,12 @@
7571

7672
import numpy as np
7773
import scipy as sp
78-
from numpy import atleast_1d, einsum, maximum, minimum
74+
from numpy import einsum, maximum, minimum
7975
from scipy.linalg import eig, eigvals, matrix_balance, norm
8076

8177
from . import config
82-
from .lti import (LTI, isctime, isdtime)
83-
from .statesp import _convert_to_statespace, _mimo2simo, _mimo2siso, ssdata
78+
from .lti import isctime, isdtime
79+
from .statesp import StateSpace, _convert_to_statespace, _mimo2simo, _mimo2siso
8480
from .xferfcn import TransferFunction
8581

8682
__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,
209205
210206
Parameters
211207
----------
212-
sys : LTI (StateSpace or TransferFunction)
208+
sys : StateSpace or TransferFunction
213209
LTI system to simulate
214210
215211
T : array_like, optional for discrete LTI `sys`
@@ -284,9 +280,9 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
284280
See :ref:`time-series-convention`.
285281
286282
"""
287-
if not isinstance(sys, LTI):
288-
raise TypeError('Parameter ``sys``: must be a ``LTI`` object. '
289-
'(For example ``StateSpace`` or ``TransferFunction``)')
283+
if not isinstance(sys, (StateSpace, TransferFunction)):
284+
raise TypeError('Parameter ``sys``: must be a ``StateSpace`` or'
285+
' ``TransferFunction``)')
290286

291287
# If return_x was not specified, figure out the default
292288
if return_x is None:
@@ -738,20 +734,24 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
738734
squeeze=squeeze, input=input, output=output)
739735

740736

741-
def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
737+
def step_info(sysdata, T=None, T_num=None, SettlingTimeThreshold=0.02,
742738
RiseTimeLimits=(0.1, 0.9)):
743739
"""
744740
Step response characteristics (Rise time, Settling Time, Peak and others).
745741
746742
Parameters
747743
----------
748-
sys : LTI system
744+
sysdata : StateSpace or TransferFunction or array_like
745+
The system data. Either LTI system to similate (StateSpace,
746+
TransferFunction), or a time series of step response data.
749747
T : array_like or float, optional
750748
Time vector, or simulation time duration if a number (time vector is
751749
autocomputed if not given, see :func:`step_response` for more detail)
750+
Required, if sysdata is a time series of response data.
752751
T_num : int, optional
753752
Number of time steps to use in simulation if T is not provided as an
754-
array (autocomputed if not given); ignored if sys is discrete-time.
753+
array; autocomputed if not given; ignored if sysdata is a
754+
discrete-time system or a time series or response data.
755755
SettlingTimeThreshold : float value, optional
756756
Defines the error to compute settling time (default = 0.02)
757757
RiseTimeLimits : tuple (lower_threshold, upper_theshold)
@@ -760,7 +760,8 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
760760
Returns
761761
-------
762762
S : dict or list of list of dict
763-
If `sys` is a SISO system, S is a dictionary containing:
763+
If `sysdata` corresponds to a SISO system, S is a dictionary
764+
containing:
764765
765766
RiseTime:
766767
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,
781782
SteadyStateValue:
782783
Steady-state value
783784
784-
If `sys` is a MIMO system, `S` is a 2D list of dicts. To get the
785-
step response characteristics from the j-th input to the i-th output,
786-
access ``S[i][j]``
785+
If `sysdata` corresponds to a MIMO system, `S` is a 2D list of dicts.
786+
To get the step response characteristics from the j-th input to the
787+
i-th output, access ``S[i][j]``
787788
788789
789790
See Also
@@ -833,18 +834,46 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
833834
PeakTime: 4.242
834835
SteadyStateValue: -1.0
835836
"""
836-
if T is None or np.asarray(T).size == 1:
837-
T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True)
838-
T, Yout = step_response(sys, T, squeeze=False)
837+
if isinstance(sysdata, (StateSpace, TransferFunction)):
838+
if T is None or np.asarray(T).size == 1:
839+
T = _default_time_vector(sysdata, N=T_num, tfinal=T, is_step=True)
840+
T, Yout = step_response(sysdata, T, squeeze=False)
841+
InfValues = np.atleast_2d(sysdata.dcgain())
842+
retsiso = sysdata.issiso()
843+
noutputs = sysdata.noutputs
844+
ninputs = sysdata.ninputs
845+
else:
846+
# Time series of response data
847+
errmsg = ("`sys` must be a LTI system, or time response data"
848+
" with a shape following the python-control"
849+
" time series data convention.")
850+
try:
851+
Yout = np.array(sysdata, dtype=float)
852+
except ValueError:
853+
raise ValueError(errmsg)
854+
if Yout.ndim == 1 or (Yout.ndim == 2 and Yout.shape[0] == 1):
855+
Yout = Yout[np.newaxis, np.newaxis, :]
856+
retsiso = True
857+
elif Yout.ndim == 3:
858+
retsiso = False
859+
else:
860+
raise ValueError(errmsg)
861+
if T is None or Yout.shape[2] != len(np.squeeze(T)):
862+
raise ValueError("For time response data, a matching time vector"
863+
" must be given")
864+
T = np.squeeze(T)
865+
noutputs = Yout.shape[0]
866+
ninputs = Yout.shape[1]
867+
InfValues = Yout[:, :, -1]
839868

840869
ret = []
841-
for i in range(sys.noutputs):
870+
for i in range(noutputs):
842871
retrow = []
843-
for j in range(sys.ninputs):
872+
for j in range(ninputs):
844873
yout = Yout[i, j, :]
845874

846875
# Steady state value
847-
InfValue = sys.dcgain() if sys.issiso() else sys.dcgain()[i, j]
876+
InfValue = InfValues[i, j]
848877
sgnInf = np.sign(InfValue.real)
849878

850879
rise_time: float = np.NaN
@@ -869,7 +898,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
869898

870899
# SettlingTime
871900
settled = np.where(
872-
np.abs(yout/InfValue -1) >= SettlingTimeThreshold)[0][-1]+1
901+
np.abs(yout/InfValue-1) >= SettlingTimeThreshold)[0][-1]+1
873902
# MIMO systems can have unsettled channels without infinite
874903
# InfValue
875904
if settled < len(T):
@@ -917,7 +946,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
917946

918947
ret.append(retrow)
919948

920-
return ret[0][0] if sys.issiso() else ret
949+
return ret[0][0] if retsiso else ret
921950

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

0 commit comments

Comments
 (0)