Skip to content

Commit ecc3b5a

Browse files
murrayrmmurrayrm
committed
update stochastic.rst
1 parent eb70e96 commit ecc3b5a

2 files changed

Lines changed: 229 additions & 46 deletions

File tree

‎doc/statesp.rst‎

Lines changed: 5 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -162,46 +162,8 @@ calling the :func:`place` function::
162162
where `E` is the desired location of the eigenvalues and `.T` computes
163163
the transpose of a matrix.
164164

165-
Alternatively, an optimal estimator can be computed using the
166-
:func:`lqe` function. We consider a continuous time, state space
167-
system
168-
169-
.. math::
170-
171-
\frac{dx}{dt} &= Ax + Bu + Gw \\
172-
y &= Cx + Du + v
173-
174-
with unbiased process noise :math:`w` and measurement noise :math:`v`
175-
with covariances satisfying
176-
177-
.. math::
178-
179-
{\mathbb E}\{w w^T\} = QN,\qquad
180-
{\mathbb E}\{v v^T\} = RN,\qquad
181-
{\mathbb E}\{w v^T\} = NN
182-
183-
where :math:`{\mathbb E}\{\cdot\}` represents expectation.
184-
185-
The :func:`lqe` function computes the observer gain matrix :math:`L`
186-
such that the stationary (non-time-varying) Kalman filter
187-
188-
.. math::
189-
190-
\frac{d\hat x}{dt} = A \hat x + B u + L(y - C\hat x - D u),
191-
192-
produces a state estimate :math:`\hat x` that minimizes the expected
193-
squared error using the sensor measurements :math:`y`.
194-
195-
As with the :func:`lqr` function, the :func:`lqe` function can be called in several forms:
196-
197-
* `L, P, E = lqe(sys, QN, RN)`
198-
* `L, P, E = lqe(sys, QN, RN, NN)`
199-
* `L, P, E = lqe(A, G, C, QN, RN)`
200-
* `L, P, E = lqe(A, G, C, QN, RN, NN)`
201-
202-
where `sys` is an :class:`LTI` object, and `A`, `G`, `C`, `QN`, `RN`,
203-
and `NN` are 2D arrays of appropriate dimension. If `sys` is a
204-
discrete time system, the first two forms will compute the discrete
205-
time optimal controller. For the second two forms, the :func:`dlqr`
206-
function can be used. Additional arguments and details are given on
207-
the :func:`lqr` and :func:`dlqr` documentation pages.
165+
More sophisticated estimators can be constructed by modeling noise and
166+
disturbances as stochastic signals generated by a random process.
167+
Estimators constructed using these models are described in more detail
168+
in the :ref:kalman-filter section of the :ref:stochastic-systems
169+
chapter.

‎doc/stochastic.rst‎

Lines changed: 224 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,204 @@
1+
.. currentmodule:: control
2+
3+
.. _stochastic-systems:
4+
15
******************
26
Stochastic Systems
37
******************
48

5-
Optimal estimation problem setup
6-
--------------------------------
9+
The Python Control Systems Library has support for basic operations
10+
involving linear and nonlinear I/O systems with Gaussian white noise
11+
as an input.
12+
13+
14+
Stochastic signals
15+
==================
16+
17+
A stochastic signal is a representation of the output of a random
18+
process. NumPy and SciPy have a number of functions to calculate the
19+
covariance and correlation of random signals:
20+
21+
* :func:`numpy.cov` - returns the sample variance of vector random
22+
variable :math:`X \in \mathbb{R}^n` where the input argument
23+
represents samples of :math:`X`.
24+
25+
* :func:`numpy.cov` - returns the (cross-)covariance of the
26+
variables :math:`X` and :math:`Y` where the input arguments
27+
represent samples of the given random variables.
28+
29+
* :func:`scipy.signal.correlate` - the "cross-correlation" between two
30+
random (1D) sequences. If these sequences came from a random
31+
process, this is a single sample approximation of the (discrete
32+
time) correlation function. Use the function
33+
:func:`scipy.signal.correlation_lags` to compute the lag
34+
:math:`\tau` and :func:`scipy.signal.correlate` to get the (auto)
35+
correlation function :math:`r_X(\tau)`.
36+
37+
The python-control package has variants of these functions that do
38+
appropriate processing for continuous time models.
39+
40+
The :func:`white_noise` function generates a (multi-variable) white
41+
noise signal of specified intensity as either a sampled continuous time
42+
signal or a discrete time signal. A white noise signal along a 1D
43+
array of linearly spaced set of times `timepts` can be computing using
44+
45+
.. code::
46+
47+
V = ct.white_noise(timepts, Q[, dt])
48+
49+
where `Q` is a positive definite matrix providing the noise
50+
intensity and `dt` is the sampling time (or 0 for continuous time).
51+
52+
In continuous time, the white noise signal is scaled such that the
53+
integral of the covariance over a sample period is `Q`, thus
54+
approximating a white noise signal. In discrete time, the white noise
55+
signal has covariance `Q` at each point in time (without any
56+
scaling based on the sample time).
57+
58+
The python-control :func:`correlation` function computes the
59+
correlation matrix :math:`{\mathbb E}\{X^\mathsf{T}(t+\tau) X(t)\}` or the
60+
cross-correlation matrix :math:`{\mathbb E}\{X^\mathsf{T}(t+\tau) Y(t)\}`:
61+
62+
.. code::
63+
64+
tau, Rtau = ct.correlation(timepts, X[, Y])
65+
66+
where :math:`\mathbb{E}` represents expectation. The signal `X` (and
67+
`Y`, if present) represents a continuous time signal sampled at
68+
regularly spaced times `timepts`. The return value provides the
69+
correlation :math:`R_\tau` between :math:`X(t+\tau)` and :math:`X(t)`
70+
at a set of time offsets :math:`\tau` (determined based on the spacing
71+
of entries in the `timepts` vector.
772

8-
Consider a nonlinear system with discrete time dynamics of the form
73+
Note that the computation of the correlation function is based on a
74+
single time signal (or pair of time signals) and is thus a very crude
75+
approximation to the true correlation function between two random
76+
processes.
77+
78+
To compute the response of a linear (or nonlinear) system to a white
79+
noise input, use the :func:`forced_response` (or
80+
:func:`input_output_response`) function:
81+
82+
.. code::
83+
84+
a, c = 1, 1
85+
sys = ct.ss([[-a]], [[1]], [[c]], 0)
86+
timepts = np.linspace(0, 5, 1000)
87+
Q = np.array([[0.1]])
88+
V = ct.white_noise(timepts, Q)
89+
resp = ct.forced_response(sys, timepts, V)
90+
91+
The correlation function for the output can be computed using the
92+
:func:`correlation` function and compared to the analytical expression:
93+
94+
.. code::
95+
96+
tau, r_Y = ct.correlation(timepts, resp.outputs)
97+
plt.plot(tau, r_Y)
98+
plt.plot(tau, c**2 * Q.item() / (2 * a) * np.exp(-a * np.abs(tau)))
99+
100+
.. _kalman-filter:
101+
102+
Linear quadratic estimation (Kalman filter)
103+
===========================================
104+
105+
A standard application of stochastic linear systems is the computation
106+
of the linear estimator under the assumption of white Gaussian
107+
measurement and process noise. This estimator is called the linear
108+
quadratic estimator (LQE) and its gains can be computed using the
109+
:func:`lqe` function.
110+
111+
We consider a continuous time, state space system
112+
113+
.. math::
114+
115+
\frac{dx}{dt} &= Ax + Bu + Gw \\
116+
y &= Cx + Du + v
117+
118+
with unbiased process noise :math:`w` and measurement noise :math:`v`
119+
with covariances satisfying
120+
121+
.. math::
122+
123+
{\mathbb E}\{w w^T\} = QN,\qquad
124+
{\mathbb E}\{v v^T\} = RN,\qquad
125+
{\mathbb E}\{w v^T\} = NN
126+
127+
where :math:`{\mathbb E}\{\cdot\}` represents expectation.
128+
129+
The :func:`lqe` function computes the observer gain matrix :math:`L`
130+
such that the stationary (non-time-varying) Kalman filter
131+
132+
.. math::
133+
134+
\frac{d\hat x}{dt} = A \hat x + B u + L(y - C\hat x - D u),
135+
136+
produces a state estimate :math:`\hat x` that minimizes the expected
137+
squared error using the sensor measurements :math:`y`.
138+
139+
As with the :func:`lqr` function, the :func:`lqe` function can be
140+
called in several forms:
141+
142+
* `L, P, E = lqe(sys, QN, RN)`
143+
* `L, P, E = lqe(sys, QN, RN, NN)`
144+
* `L, P, E = lqe(A, G, C, QN, RN)`
145+
* `L, P, E = lqe(A, G, C, QN, RN, NN)`
146+
147+
where `sys` is an :class:`LTI` object, and `A`, `G`, `C`, `QN`, `RN`,
148+
and `NN` are 2D arrays of appropriate dimension. If `sys` is a
149+
discrete time system, the first two forms will compute the discrete
150+
time optimal controller. For the second two forms, the :func:`dlqr`
151+
function can be used. Additional arguments and details are given on
152+
the :func:`lqr` and :func:`dlqr` documentation pages.
153+
154+
The :func:`create_estimator_iosystem` function can be used to create
155+
an I/O system implementing a Kalman filter, including integration of
156+
the Riccati ODE. The command has the form
157+
158+
.. code::
159+
160+
estim = ct.create_estimator_iosystem(sys, Qv, Qw)
161+
162+
The input to the estimator is the measured outputs `Y` and the system
163+
input `U`. To run the estimator on a noisy signal, use the command
164+
165+
.. code::
166+
167+
resp = ct.input_output_response(est, timepts, [Y, U], [X0, P0])
168+
169+
If desired, the :func:`correct` parameter can be set to :func:`False`
170+
to allow prediction with no additional sensor information::
171+
172+
resp = ct.input_output_response(
173+
estim, timepts, 0, [X0, P0], param={'correct': False})
174+
175+
The :func:`create_statefbk_iosystem` function can be used to combine
176+
an estimator with a state feedback controller::
177+
178+
K, _, _ = ct.lqr(sys, Qx, Qu)
179+
estim = ct.create_estimator_iosystem(sys, Qv, Qw, P0)
180+
ctrl, clsys = ct.create_statefbk_iosystem(sys, K, estimator=estim)
181+
182+
The controller will have the same form as a full state feedback
183+
controller, but with the system state :math:`x` input replaced by the
184+
estimated state :math:`\hat x` (output of `estim`):
185+
186+
.. math::
187+
188+
u = u_\text{d} - K (\hat x - x_\text{d}).
189+
190+
The closed loop controller `clsys` includes both the state
191+
feedback and the estimator dynamics and takes as its input the desired
192+
state :math:`x_\text{d}` and input :math:`u_\text{d}`::
193+
194+
resp = ct.input_output_response(
195+
clsys, timepts, [Xd, Ud], [X0, np.zeros_like(X0), P0])
196+
197+
198+
Maximum likelihood estimation
199+
=============================
200+
201+
Consider a *nonlinear* system with discrete time dynamics of the form
9202

10203
.. math::
11204
:label: eq_fusion_nlsys-oep
@@ -104,3 +297,31 @@ In particular, we can obtain the estimated state at the end of the moving
104297
horizon window, corresponding to the current time, and we can thus
105298
implement an estimator by repeatedly solving the optimization of a window
106299
of length :math:`N` backwards in time.
300+
301+
The :mod:`optimal` module described in the :ref:`optimal-module`
302+
section implements functions for solving optimal estimation problems
303+
using maximum likelihood estimation. The
304+
:class:`optimal.OptimalEstimationProblem` class is used to define an
305+
optimal estimation problem over a finite horizon::
306+
307+
oep = opt.OptimalEstimationProblem(sys, timepts, cost[, constraints])
308+
309+
Given noisy measurements :math:`y` and control inputs :math:`u`, an
310+
estimate of the states over the time points can be computed using the
311+
:func:`optimal.OptimalEstimationProblem.compute_estimate` method::
312+
313+
estim = oep.compute_optimal(Y, U[, X0=x0, initial_guess=(xhat, v)])
314+
xhat, v, w = estim.states, estim.inputs, estim.outputs
315+
316+
For discrete time systems, the
317+
:func:`optimal.OptimalEstimationProblem.create_mhe_iosystem` method
318+
can be used to generate an input/output system that implements a
319+
moving horizon estimator.
320+
321+
Several functions are available to help set up standard optimal estimation
322+
problems:
323+
324+
.. autosummary::
325+
326+
optimal.gaussian_likelihood_cost
327+
optimal.disturbance_range_constraint

0 commit comments

Comments
 (0)