Skip to content

Commit 1202257

Browse files
murrayrmmurrayrm
committed
add documentation on predict keyword + input_output_response list processing
1 parent 7821e2b commit 1202257

3 files changed

Lines changed: 136 additions & 11 deletions

File tree

‎control/iosys.py‎

Lines changed: 80 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1585,11 +1585,17 @@ def input_output_response(
15851585
T : array-like
15861586
Time steps at which the input is defined; values must be evenly spaced.
15871587
1588-
U : array-like or number, optional
1589-
Input array giving input at each time `T` (default = 0).
1590-
1591-
X0 : array-like or number, optional
1592-
Initial condition (default = 0).
1588+
U : array-like, list, or number, optional
1589+
Input array giving input at each time `T` (default = 0). If a list
1590+
is specified, each element in the list will be treated as a portion
1591+
of the input and broadcast (if necessary) to match the time vector.
1592+
1593+
X0 : array-like, list, or number, optional
1594+
Initial condition (default = 0). If a list is given, each element
1595+
in the list will be flattened and stacked into the initial
1596+
condition. If a smaller number of elements are given that the
1597+
number of states in the system, the initial condition will be padded
1598+
with zeros.
15931599
15941600
return_x : bool, optional
15951601
If True, return the state vector when assigning to a tuple (default =
@@ -1641,6 +1647,16 @@ def input_output_response(
16411647
ValueError
16421648
If time step does not match sampling time (for discrete time systems).
16431649
1650+
Notes
1651+
-----
1652+
1. If a smaller number of initial conditions are given than the number of
1653+
states in the system, the initial conditions will be padded with
1654+
zeros. This is often useful for interconnected control systems where
1655+
the process dynamics are the first system and all other components
1656+
start with zero initial condition since this can be specified as
1657+
[xsys_0, 0]. A warning is issued if the initial conditions are padded
1658+
and and the final listed initial state is not zero.
1659+
16441660
"""
16451661
#
16461662
# Process keyword arguments
@@ -1683,19 +1699,75 @@ def input_output_response(
16831699
# Use the input time points as the output time points
16841700
t_eval = T
16851701

1686-
# Check and convert the input, if needed
1687-
# TODO: improve MIMO ninputs check (choose from U)
1702+
# If we were passed a list of input, concatenate them (w/ broadcast)
1703+
if isinstance(U, (tuple, list)) and len(U) != ntimepts:
1704+
U_elements = []
1705+
for i, u in enumerate(U):
1706+
u = np.array(u) # convert everyting to an array
1707+
# Process this input
1708+
if u.ndim == 0 or (u.ndim == 1 and u.shape[0] != T.shape[0]):
1709+
# Broadcast array to the length of the time input
1710+
u = np.outer(u, np.ones_like(T))
1711+
1712+
elif (u.ndim == 1 and u.shape[0] == T.shape[0]) or \
1713+
(u.ndim == 2 and u.shape[1] == T.shape[0]):
1714+
# No processing necessary; just stack
1715+
pass
1716+
1717+
else:
1718+
raise ValueError(f"Input element {i} has inconsistent shape")
1719+
1720+
# Append this input to our list
1721+
U_elements.append(u)
1722+
1723+
# Save the newly created input vector
1724+
U = np.vstack(U_elements)
1725+
1726+
# Make sure the input has the right shape
16881727
if sys.ninputs is None or sys.ninputs == 1:
16891728
legal_shapes = [(ntimepts,), (1, ntimepts)]
16901729
else:
16911730
legal_shapes = [(sys.ninputs, ntimepts)]
1731+
16921732
U = _check_convert_array(U, legal_shapes,
16931733
'Parameter ``U``: ', squeeze=False)
1734+
1735+
# Always store the input as a 2D array
16941736
U = U.reshape(-1, ntimepts)
16951737
ninputs = U.shape[0]
16961738

1697-
# create X0 if not given, test if X0 has correct shape
1739+
# If we were passed a list of initial states, concatenate them
1740+
if isinstance(X0, (tuple, list)):
1741+
X0_list = []
1742+
for i, x0 in enumerate(X0):
1743+
x0 = np.array(x0).reshape(-1) # convert everyting to 1D array
1744+
X0_list += x0.tolist() # add elements to initial state
1745+
1746+
# Save the newly created input vector
1747+
X0 = np.array(X0_list)
1748+
1749+
# If the initial state is too short, make it longer (NB: sys.nstates
1750+
# could be None if nstates comes from size of initial condition)
1751+
if sys.nstates and isinstance(X0, np.ndarray) and X0.size < sys.nstates:
1752+
if X0[-1] != 0:
1753+
warn("initial state too short; padding with zeros")
1754+
X0 = np.hstack([X0, np.zeros(sys.nstates - X0.size)])
1755+
1756+
# Check to make sure this is not a static function
16981757
nstates = _find_size(sys.nstates, X0)
1758+
if nstates == 0:
1759+
# No states => map input to output
1760+
u = U[0] if len(U.shape) == 1 else U[:, 0]
1761+
y = np.zeros((np.shape(sys._out(T[0], X0, u))[0], len(T)))
1762+
for i in range(len(T)):
1763+
u = U[i] if len(U.shape) == 1 else U[:, i]
1764+
y[:, i] = sys._out(T[i], [], u)
1765+
return TimeResponseData(
1766+
T, y, None, U, issiso=sys.issiso(),
1767+
output_labels=sys.output_index, input_labels=sys.input_index,
1768+
transpose=transpose, return_x=return_x, squeeze=squeeze)
1769+
1770+
# create X0 if not given, test if X0 has correct shape
16991771
X0 = _check_convert_array(X0, [(nstates,), (nstates, 1)],
17001772
'Parameter ``X0``: ', squeeze=True)
17011773

‎control/stochsys.py‎

Lines changed: 15 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -323,9 +323,11 @@ def create_estimator_iosystem(
323323
324324
estim = ct.create_estimator_iosystem(sys, QN, RN)
325325
326-
where ``sys`` is the process dynamics and QN and RN are the covariance of
327-
the disturbance noise and sensor noise. The function returns the
328-
estimator ``estim`` as I/O systems.
326+
where ``sys`` is the process dynamics and QN and RN are the covariance
327+
of the disturbance noise and sensor noise. The function returns the
328+
estimator ``estim`` as I/O system with a parameter ``correct`` that can
329+
be used to turn off the correction term in the estimation (for forward
330+
predictions).
329331
330332
Parameters
331333
----------
@@ -368,6 +370,16 @@ def create_estimator_iosystem(
368370
est = ct.create_estimator_iosystem(sys, QN, RN, P0)
369371
ctrl, clsys = ct.create_statefbk_iosystem(sys, K, estimator=est)
370372
373+
The estimator can also be run on its own to process a noisy signal:
374+
375+
resp = ct.input_output_response(est, T, [Y, U], [X0, P0])
376+
377+
If desired, the ``correct`` parameter can be set to ``False`` to allow
378+
prediction with no additional sensor information:
379+
380+
resp = ct.input_output_response(
381+
est, T, 0, [X0, P0], param={'correct': False)
382+
371383
"""
372384

373385
# Make sure that we were passed an I/O system as an input

‎control/tests/iosys_test.py‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1732,3 +1732,44 @@ def test_nonuniform_timepts():
17321732
t_even, y_even = ct.input_output_response(
17331733
sys, noufpts, nonunif, t_eval=unifpts)
17341734
np.testing.assert_almost_equal(y_unif, y_even, decimal=6)
1735+
1736+
def test_input_output_broadcasting():
1737+
# Create a system, time vector, and noisy input
1738+
sys = ct.rss(6, 2, 3)
1739+
T = np.linspace(0, 10, 10)
1740+
U = np.zeros((sys.ninputs, T.size))
1741+
U[0, :] = np.sin(T)
1742+
U[1, :] = np.zeros_like(U[1, :])
1743+
U[2, :] = np.ones_like(U[2, :])
1744+
X0 = np.array([1, 2])
1745+
P0 = np.array([[3.11, 3.12], [3.21, 3.3]])
1746+
1747+
# Simulate the system with nominal input to establish baseline
1748+
resp_base = ct.input_output_response(
1749+
sys, T, U, np.hstack([X0, P0.reshape(-1)]))
1750+
1751+
# Split up the inputs into two pieces
1752+
resp_inp1 = ct.input_output_response(sys, T, [U[:1], U[1:]], [X0, P0])
1753+
np.testing.assert_equal(resp_base.states, resp_inp1.states)
1754+
1755+
# Specify two of the inputs as constants
1756+
resp_inp2 = ct.input_output_response(sys, T, [U[0], 0, 1], [X0, P0])
1757+
np.testing.assert_equal(resp_base.states, resp_inp2.states)
1758+
1759+
# Specify two of the inputs as constant vector
1760+
resp_inp3 = ct.input_output_response(sys, T, [U[0], [0, 1]], [X0, P0])
1761+
np.testing.assert_equal(resp_base.states, resp_inp3.states)
1762+
1763+
# Specify only some of the initial conditions
1764+
resp_init = ct.input_output_response(sys, T, [U[0], [0, 1]], [X0, 0])
1765+
resp_cov0 = ct.input_output_response(sys, T, U, [X0, P0 * 0])
1766+
np.testing.assert_equal(resp_cov0.states, resp_init.states)
1767+
1768+
# Specify only some of the initial conditions
1769+
with pytest.warns(UserWarning, match="initial state too short; padding"):
1770+
resp_short = ct.input_output_response(sys, T, [U[0], [0, 1]], [X0, 1])
1771+
1772+
# Make sure that inconsistent settings don't work
1773+
with pytest.raises(ValueError, match="inconsistent"):
1774+
resp_bad = ct.input_output_response(
1775+
sys, T, (U[0, :], U[:2, :-1]), [X0, P0])

0 commit comments

Comments
 (0)