Skip to content

Commit 983726c

Browse files
murrayrmmurrayrm
authored
Merge pull request #714 from murrayrm/stochsys
Stochastic systems additions
2 parents cb6d9d7 + 34c2c7e commit 983726c

22 files changed

+2824
-403
lines changed

‎control/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,7 @@
6161
from .rlocus import *
6262
from .statefbk import *
6363
from .statesp import *
64+
from .stochsys import *
6465
from .timeresp import *
6566
from .xferfcn import *
6667
from .ctrlutil import *

‎control/iosys.py‎

Lines changed: 84 additions & 12 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
@@ -1656,14 +1672,14 @@ def input_output_response(
16561672
raise ValueError("ivp_method specified more than once")
16571673
solve_ivp_kwargs['method'] = kwargs.pop('solve_ivp_method')
16581674

1659-
# Make sure there were no extraneous keywords
1660-
if kwargs:
1661-
raise TypeError("unrecognized keywords: ", str(kwargs))
1662-
16631675
# Set the default method to 'RK45'
16641676
if solve_ivp_kwargs.get('method', None) is None:
16651677
solve_ivp_kwargs['method'] = 'RK45'
16661678

1679+
# Make sure there were no extraneous keywords
1680+
if kwargs:
1681+
raise TypeError("unrecognized keyword(s): ", str(kwargs))
1682+
16671683
# Sanity checking on the input
16681684
if not isinstance(sys, InputOutputSystem):
16691685
raise TypeError("System of type ", type(sys), " not valid")
@@ -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/optimal.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -776,7 +776,7 @@ def create_mpc_iosystem(self):
776776
"""Create an I/O system implementing an MPC controller"""
777777
# Check to make sure we are in discrete time
778778
if self.system.dt == 0:
779-
raise ControlNotImplemented(
779+
raise ct.ControlNotImplemented(
780780
"MPC for continuous time systems not implemented")
781781

782782
def _update(t, x, u, params={}):

‎control/sisotool.py‎

Lines changed: 11 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -215,14 +215,14 @@ def rootlocus_pid_designer(plant, gain='P', sign=+1, input_signal='r',
215215
derivative terms are given instead by Kp, Ki*dt/2*(z+1)/(z-1), and
216216
Kd/dt*(z-1)/z, respectively.
217217
218-
------> C_ff ------ d
219-
| | |
220-
r | e V V u y
221-
------->O---> C_f --->O--->O---> plant --->
222-
^- ^- |
223-
| | |
224-
| ----- C_b <-------|
225-
---------------------------------
218+
------> C_ff ------ d
219+
| | |
220+
r | e V V u y
221+
------->O---> C_f --->O--->O---> plant --->
222+
^- ^- |
223+
| | |
224+
| ----- C_b <-------|
225+
---------------------------------
226226
227227
It is also possible to move the derivative term into the feedback path
228228
`C_b` using `derivative_in_feedback_path=True`. This may be desired to
@@ -234,8 +234,8 @@ def rootlocus_pid_designer(plant, gain='P', sign=+1, input_signal='r',
234234
235235
Remark: It may be helpful to zoom in using the magnifying glass on the
236236
plot. Just ake sure to deactivate magnification mode when you are done by
237-
clicking the magnifying glass. Otherwise you will not be able to be able to choose
238-
a gain on the root locus plot.
237+
clicking the magnifying glass. Otherwise you will not be able to be able
238+
to choose a gain on the root locus plot.
239239
240240
Parameters
241241
----------
@@ -269,6 +269,7 @@ def rootlocus_pid_designer(plant, gain='P', sign=+1, input_signal='r',
269269
----------
270270
closedloop : class:`StateSpace` system
271271
The closed-loop system using initial gains.
272+
272273
"""
273274

274275
plant = _convert_to_statespace(plant)

0 commit comments

Comments
 (0)