{
"cells": [
{
"cell_type": "markdown",
"id": "03aa22e7",
"metadata": {},
"source": [
"# Stochastic Response\n",
"Richard M. Murray, 6 Feb 2022\n",
"\n",
"This notebook illustrates the implementation of random processes and stochastic response. We focus on a system of the form\n",
"$$\n",
" \\dot X = A X + F V \\qquad X \\in {\\mathbb R}^n\n",
"$$\n",
"\n",
"where $V$ is a white noise process and the system is a first order linear system."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "902af902",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pyplot as plt\n",
"import control as ct\n",
"import ctrlutil as ct_\n",
"from math import sqrt, exp"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "60192a8c",
"metadata": {},
"outputs": [],
"source": [
"# First order system\n",
"a = 1\n",
"c = 1\n",
"sys = ct.tf(c, [1, a])\n",
"\n",
"# Create the time vector that we want to use\n",
"Tf = 5\n",
"T = np.linspace(0, Tf, 1000)\n",
"dt = T[1] - T[0]\n",
"\n",
"# Create the basis for a white noise signal\n",
"# Note: use sqrt(Q/dt) for desired covariance\n",
"Q = np.array([[0.1]])\n",
"# V = np.random.normal(0, sqrt(Q[0,0]/dt), T.shape)\n",
"V = ct_.white_noise(T, Q)\n",
"\n",
"plt.plot(T, V[0])\n",
"plt.xlabel('Time [s]')\n",
"plt.ylabel('$V$');"
]
},
{
"cell_type": "markdown",
"id": "b4629e2c",
"metadata": {},
"source": [
"Note that the magnitude of the signal seems to be much larger than $Q$. This is because we have a Guassian process $\\implies$ the signal consists of a sequence of \"impulse-like\" functions that have magnitude that increases with the time step $dt$ as $1/\\sqrt{dt}$ (this gives covariance $\\mathbb{E}(V(t_1) V^T(t_2)) = Q \\delta(t_2 - t_1)$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "23319dc6",
"metadata": {},
"outputs": [],
"source": [
"# Calculate the sample properties and make sure they match\n",
"print(\"mean(V) [0.0] = \", np.mean(V))\n",
"print(\"cov(V) * dt [%0.3g] = \" % Q, np.round(np.cov(V), decimals=3) * dt)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2bdaaccf",
"metadata": {},
"outputs": [],
"source": [
"# Response of the first order system\n",
"# Scale white noise by sqrt(dt) to account for impulse\n",
"T, Y = ct.forced_response(sys, T, V)\n",
"plt.plot(T, Y)\n",
"plt.xlabel('Time [s]')\n",
"plt.ylabel('$Y$');"
]
},
{
"cell_type": "markdown",
"id": "ead0232e",
"metadata": {},
"source": [
"This is a first order system, and so we can use the calculation from the course\n",
"notes to compute the analytical correlation function and compare this to the \n",
"sampled data:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d31ce324",
"metadata": {},
"outputs": [],
"source": [
"# Compare static properties to what we expect analytically\n",
"def r(tau):\n",
" return c**2 * Q / (2 * a) * exp(-a * abs(tau))\n",
" \n",
"print(\"* mean(Y) [%0.3g] = %0.3g\" % (0, np.mean(Y)))\n",
"print(\"* cov(Y) [%0.3g] = %0.3g\" % (r(0), np.cov(Y)))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1cf5a4b1",
"metadata": {},
"outputs": [],
"source": [
"# Correlation function for the input\n",
"# Scale by dt to take time step into account\n",
"# r_V = sp.signal.correlate(V, V) * dt / Tf\n",
"# tau = sp.signal.correlation_lags(len(V), len(V)) * dt\n",
"tau, r_V = ct_.correlation(T, V)\n",
"\n",
"plt.plot(tau, r_V, 'r-')\n",
"plt.xlabel(r'$\\tau$')\n",
"plt.ylabel(r'$r_V(\\tau)$');"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "62af90a4",
"metadata": {},
"outputs": [],
"source": [
"# Correlation function for the output\n",
"# r_Y = sp.signal.correlate(Y, Y) * dt / Tf\n",
"# tau = sp.signal.correlation_lags(len(Y), len(Y)) * dt\n",
"tau, r_Y = ct_.correlation(T, Y)\n",
"plt.plot(tau, r_Y)\n",
"plt.xlabel(r'$\\tau$')\n",
"plt.ylabel(r'$r_Y(\\tau)$')\n",
"\n",
"# Compare to the analytical answer\n",
"plt.plot(tau, [r(t)[0, 0] for t in tau], 'k--');"
]
},
{
"cell_type": "markdown",
"id": "2a2785e9",
"metadata": {},
"source": [
"The analytical curve may or may not line up that well with the correlation function based on the sample. Try running the code again from the top to see how things change based on the specific random sequence chosen at the start.\n",
"\n",
"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 $R(t_1, t_2)$ across those samples."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bd5dfc75",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 5
}