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.
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 .
# First order system a = 1 c = 1 sys = ct.ss([[-a]], [], [[c]], 0) # Create the time vector that we want to use Tf = 5 T = np.linspace(0, Tf, 1000) dt = T - T # Create the basis for a white noise signal Q = np.array([[0.1]]) V = ct.white_noise(T, Q) plt.plot(T, V) 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 .
# 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
# 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:
# 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.165 * cov(Y) [0.05] = 0.0151
Finally, we look at the correlation function for the input and the output:
# 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)$');
# 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.