# Stochastic Response¶

Richard M. Murray, 6 Feb 2022

This notebook illustrates the implementation of random processes and stochastic response. We focus on a system of the form

where is a white noise process and the system is a first order linear system.

```
[83]:
```

```
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import control as ct
from math import sqrt, exp
# Fix random number seed to avoid spurious figure regeneration
np.random.seed(1)
```

We begin by defining a simple first order system

and a (scalar) white noise signal with intensity .

```
[84]:
```

```
# First order system
a = 1
c = 1
sys = ct.ss([[-a]], [[1]], [[c]], 0)
# Create the time vector that we want to use
Tf = 5
T = np.linspace(0, Tf, 1000)
dt = T[1] - T[0]
# Create the basis for a white noise signal
Q = np.array([[0.1]])
V = ct.white_noise(T, Q)
plt.plot(T, V[0])
plt.xlabel('Time [s]')
plt.ylabel('$V$');
```

Note that the magnitude of the signal seems to be much larger than . This is because we have a Guassian process the signal consists of a sequence of “impulse-like” functions that have magnitude that increases with the time step as (this gives covariance .

```
[85]:
```

```
# Calculate the sample properties and make sure they match
print("mean(V) [0.0] = ", np.mean(V))
print("cov(V) * dt [%0.3g] = " % Q, np.round(np.cov(V), decimals=3) * dt)
```

```
mean(V) [0.0] = 0.17348786109316244
cov(V) * dt [0.1] = 0.09633133133133133
```

The response of the system to white noise can be computed using the `input_output_response`

function:

```
[86]:
```

```
# Response of the first order system
T, Y = ct.input_output_response(sys, T, V)
plt.plot(T, Y)
plt.xlabel('Time [s]')
plt.ylabel('$Y$');
```

This is a first order system, and so we can compute the analytical correlation function and compare this to the sampled data:

```
[87]:
```

```
# Compare static properties to what we expect analytically
def r(tau):
return c**2 * Q / (2 * a) * exp(-a * abs(tau))
print("* mean(Y) [%0.3g] = %0.3g" % (0, np.mean(Y)))
print("* cov(Y) [%0.3g] = %0.3g" % (r(0), np.cov(Y)))
```

```
* mean(Y) [0] = 0.165
* cov(Y) [0.05] = 0.0151
```

Finally, we look at the correlation function for the input and the output:

```
[88]:
```

```
# Correlation function for the input
tau, r_V = ct.correlation(T, V)
plt.plot(tau, r_V, 'r-')
plt.xlabel(r'$\tau$')
plt.ylabel(r'$r_V(\tau)$');
```

```
[89]:
```

```
# Correlation function for the output
tau, r_Y = ct.correlation(T, Y)
plt.plot(tau, r_Y)
plt.xlabel(r'$\tau$')
plt.ylabel(r'$r_Y(\tau)$')
# Compare to the analytical answer
plt.plot(tau, [r(t)[0, 0] for t in tau], 'k--');
```

The analytical curve may or may not line up that well with the correlation function based on the sample. Try running the code again with a different random seed to see how things change based on the specific random sequence chosen at the start.

Note: the *right* way to compute the correlation function would be to run a lot of different samples of white noise filtered through the system dynamics and compute across those samples. The `correlation`

function computes the covariance between and by varying over the time range.

```
[ ]:
```

```
```