Output and state covariance of system driven by white noise

`P = covar(sys,W) `

[P,Q] = covar(sys,W)

`covar`

calculates the stationary covariance
of the output *y* of an LTI model `sys`

driven
by Gaussian white noise inputs *w*. This function
handles both continuous- and discrete-time cases.

`P = covar(sys,W) `

returns the steady-state
output response covariance

$$P=E(y{y}^{T})$$

given the noise intensity

$$\begin{array}{cc}E(w(t)w{(\tau )}^{T})=W\delta (t-\tau )& \text{(continuoustime)}\\ E(w\left[k\right]w{\left[l\right]}^{T})=W{\delta}_{kl}& \text{(discretetime)}\end{array}$$

`[P,Q] = covar(sys,W)`

also returns the steady-state
state covariance

$$Q=E(x{x}^{T})$$

when `sys`

is a state-space model (otherwise `Q`

is
set to `[]`

).

When applied to an `N`

-dimensional LTI array `sys`

, `covar`

returns
multidimensional arrays *P*, *Q* such
that

`P(:,:,i1,...iN)`

and `Q(:,:,i1,...iN)`

are
the covariance matrices for the model `sys(:,:,i1,...iN)`

.

Compute the output response covariance of the discrete SISO system

$$\begin{array}{cc}H(z)=\frac{2z+1}{{z}^{2}+0.2z+0.5},& {T}_{s}\end{array}=0.1$$

due to Gaussian white noise of intensity `W = 5`

.
Type

sys = tf([2 1],[1 0.2 0.5],0.1); p = covar(sys,5)

These commands produce the following result.

p = 30.3167

You can compare this output of `covar`

to simulation
results.

randn('seed',0) w = sqrt(5)*randn(1,1000); % 1000 samples % Simulate response to w with LSIM: y = lsim(sys,w); % Compute covariance of y values psim = sum(y .* y)/length(w);

This yields

psim = 32.6269

The two covariance values `p`

and `psim`

do
not agree perfectly due to the finite simulation horizon.

Transfer functions and zero-pole-gain models are first converted
to state space with `ss`

.

For continuous-time state-space models

$$\begin{array}{l}\dot{x}=Ax+Bw\\ y=Cx+Dw,\end{array}$$

the steady-state state covariance *Q* is obtained
by solving the Lyapunov equation

$$AQ+Q{A}^{T}+BW{B}^{T}=0.$$

In discrete time, the state covariance *Q* solves
the discrete Lyapunov
equation

$$AQ{A}^{T}-Q+BW{B}^{T}=0.$$

In both continuous and discrete time, the output response covariance
is given by *P* = *CQC ^{T}* +

`covar`

returns `Inf`

for
the output covariance [1] Bryson, A.E. and Y.C. Ho, *Applied
Optimal Control,* Hemisphere Publishing, 1975, pp. 458-459.