{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"import scipy.linalg"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# System Identification\n",
"\n",
"This notebook will guide you through some examples of system identification. We have a linear discrete-time system with some input $\\vec{u}[i]\\in\\mathcal{R}^{p}$ and state $\\vec{x}[i]\\in\\mathcal{R}^{n}$. At timestep $i+1$, the value of state $\\vec{x}$ evolves by\n",
"\n",
"$$\\vec{x}[i+1] = A\\vec{x}[i] + B\\vec{u}[i]$$\n",
"\n",
"where $A\\in\\mathcal{R}^{n\\times n}$ is the state transition matrix and $B\\in\\mathcal{R}^{n\\times p}$ is the input matrix. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 1: Identifying a System We Already Know\n",
"\n",
"Let's start with a simple example to show you how it works.\n",
"What we'll do is create a system whose $A$ and $B$ matrices are known, so that we know if we're right when we identify the system later. Specifically, we'll use\n",
"\\begin{align}\n",
"A &= \\begin{bmatrix} 0 & 1 \\\\ 0.3 & 0.2 \\end{bmatrix}, & B &= \\begin{bmatrix} 0 \\\\ 1 \\end{bmatrix}.\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a_known = np.array([[0, 1],\n",
" [0.3, 0.2]])\n",
"b_known = np.array([0, 1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we need to collect some data from this system. We'll use the following function to do this for us. Try following along with the code, and be sure to read the documentation."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def make_state_trace(a_matrix, b_matrix, inpt, initial_state=None):\n",
" \"\"\"\n",
" make_state_trace: for a given A and B matrices, initial condition, and input, \"run\" the system and \n",
" calculate the system state for each time step. The length of the trace will be the length of the input plus one.\n",
" \n",
" arguments:\n",
" a_matrix: A matrix of the system (numpy ndarray, n by n)\n",
" b_matrix: B matrix of the system (numpy ndarray, p by n)\n",
" inpt: input (u) (numpy ndarray, T by p)\n",
" initial_state: initial state for the system (numpy ndarray, n by 1). If no initial state is provided, \n",
" the initial state will be the origin.\n",
" \n",
" returns:\n",
" state_trace: a numpy ndarray containing the state at \n",
" each time step that was run (numpy ndarray, T+1 by n).\n",
" \"\"\"\n",
" n_states = np.shape(a_matrix)[0]\n",
" n_timesteps = np.shape(inpt)[0]\n",
" state_trace = np.zeros((n_timesteps + 1, n_states))\n",
"\n",
" state_trace[0] = initial_state if initial_state is not None else np.zeros(n_states)\n",
" for i in range(n_timesteps):\n",
" state_trace[i + 1] = a_matrix @ state_trace[i] + (b_matrix * inpt[i]).squeeze()\n",
" return state_trace"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, let's choose an input to feed into the system, and get a state trace from it.\n",
"\n",
"What kind of signal should we use as an input? In principle, it can be anything we want-- but some choices are better than others. One kind of input that generally works quite well is a *random input*, that is where each $u[i]$ is chosen at random. In this notebook, we'll be using random inputs to excite our systems. We can use the following function to generate random inputs."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def random_input(t, n_inputs, mean=0, std=1):\n",
" \"\"\"\n",
" random_input: given a time sequence t, generate a random input trace. Each entry of the trace will be sampled\n",
" from a Gaussian distribution with the given mean and standard deviation.\n",
" \"\"\"\n",
" random_trace = np.random.normal(loc=mean, scale=std, size=(len(t), n_inputs))\n",
" return random_trace"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"k = 6\n",
"t = np.arange(k)\n",
"\n",
"u_trace = random_input(t, 1)\n",
"\n",
"x_trace = make_state_trace(a_known, b_known, u_trace)\n",
"\n",
"# Let's see what the state trace looks like.\n",
"plt.plot(t, u_trace, label='input')\n",
"plt.plot(t, x_trace[:-1,0], label='state 1')\n",
"plt.plot(t, x_trace[:-1,1], label='state 2')\n",
"plt.legend(loc='best')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have a trace from our known system, let's use least squares to identify $A$ and $B$, and see if we get what we expect.\n",
"\n",
"We've written the following function to perform least-squares system identification from a trace."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def identify_system(state_trace, inpt):\n",
" \"\"\"\n",
" identify_system: given an input and the state trace it generated, determine the A and B matrices \n",
" of the system that was used to generate the state trace.\n",
" \n",
" arguments:\n",
" state_trace: the state data.\n",
" inpt: the input data:\n",
" \n",
" returns:\n",
" a_identified: the A matrix identified from the state and input data (numpy ndarray)\n",
" b_identified: the B matrix identified from the state and input data (numpy ndarray)\n",
" \"\"\"\n",
"\n",
" # set up and solve least-squares equation from state and input data\n",
" n_states = np.shape(state_trace)[1]\n",
" lsq_b = state_trace[1:, :]\n",
" lsq_a = np.concatenate((state_trace[:-1], inpt), axis=1)\n",
" ab_identified, _, _, _ = np.linalg.lstsq(lsq_a, lsq_b, rcond=None)\n",
" a_identified = ab_identified[:n_states]\n",
" b_identified = ab_identified[n_states:]\n",
" \n",
" return a_identified, b_identified\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, let's see if the system identification function works with the code you put in. Fingers crossed!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"a_identified, b_identified = identify_system(x_trace, u_trace)\n",
"print(f'Identified A matrix: \\n {np.array2string(a_identified.T, precision=4, suppress_small=True)}')\n",
"print(f'Identified B matrix: \\n {np.array2string(b_identified.T, precision=4, suppress_small=True)}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The previous cell should have printed out an $A$ and $B$ matrix that are very close to the known $A$ and $B$. It's okay if they're a tiny amount off, for example if an element that we know is 0 was identified as 1e-16 or -0 or something — that's just roundoff error.\n",
"\n",
"If the true $A$ and $B$ are close to the identified $A$ and $B$, then you have identified the system!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Variations on System Identification\n",
"\n",
"So now we know how to use state and input data to identify a model $(A,B)$ of the system that generated that data. This is big step forward. But before it becomes truly practical, we need to examine a few caveats that come up in the real world. What kinds of caveats would those be? Well, here are a few common ones:\n",
"\n",
" - Our state and input data may be corrupted with some noise.\n",
" - Instead of directly observing the state, we observe some transformed or rotated state.\n",
" - Instead of directly observing the state, we view only a single scalar observation at each time step.\n",
" - your system has too many states to be conveniently worked with, and several seem to be redundant.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generating Random Systems\n",
"\n",
"For this part of the notebook, we are going to need some new systems to test. In this case, we are going to generate *random* systems, systems whose $A$ and $B$ matrices we don't know yet. The functions provided below will allow us to do just that. Again, be sure to read the documentation. It won't be necessary for you to understand the code (some of it is out of scope), but you might want to check it out anyways. We won't be using all of the functionality of this code in this notebook, but we'll need it later."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from itertools import combinations # Huge shoutout to the combinatorial iterators in itertools, \n",
" # some of Python's lesser-known magic problem solvers.\n",
"\n",
"def make_random_system(n_states, n_inputs, z_min=0.0, z_max=1.0, must_have_real_eigenvalues=[], \n",
" must_have_complex_eigenvalues=[], ccf=False):\n",
" \"\"\"\n",
" make_random_system: make a random A and B matrix, whose eigenvalues $z$ lie in the annulus \n",
" $z_{min}<=|z|<=z_{max}$, occurring as complex conjugate pairs, plus at most one purely real. \n",
" Additionally, a number of \"must-have\" eigenvalues can be specified that the system is guaranteed to have, \n",
" which need not lie in the annulus.\n",
" \n",
" Specify ccf=True if you want the system to be given in controllable canonical form.\n",
" \n",
" arguments:\n",
" n_states: the number of states you want the sytem to have.\n",
" n_inputs: the number of inputs you want the system to have.\n",
" z_min: the smallest allowable magnitude for the random eigenvalues.\n",
" z_max: the largest allowable magnitude for the random eigenvalues.\n",
" must_have_real_eigenvalues: a list of purely eigenvalues that the random system is guaranteed to have.\n",
" must_have_complex_eigenvalues: a list of complex conjugate eigenvalues that the random system is guaranteed to have.\n",
" You only need to specify one member of the conjugate pair-- the other is implied.\n",
" ccf: whether or not you want the system matrices to be returned in controllable canonical form.\n",
" \n",
" returns:\n",
" a_matrix: A-matrix of the random system.\n",
" b_matrix: B-matrix of the random system.\n",
" \"\"\"\n",
" n_random = n_states - np.size(must_have_real_eigenvalues) - 2 * np.size(must_have_complex_eigenvalues)\n",
" assert n_random > 0, 'no random eigenvalues'\n",
" n_real = n_random % 2\n",
" n_conj_pairs = int(np.floor(n_random / 2))\n",
" # This weird probability distribution over the eigenvalue magnitudes ensures that the eigenvalues will be distributed uniformly over the annulus.\n",
" z_mags = np.sqrt(z_min ** 2 + np.random.rand(n_conj_pairs) * (z_max ** 2 - z_min ** 2))\n",
" z_angles = 2 * np.pi * np.random.rand(n_conj_pairs)\n",
" z_random_complex = z_mags * np.exp(1j * z_angles)\n",
" z_complex = np.concatenate((z_random_complex, must_have_complex_eigenvalues))\n",
" z_real = must_have_real_eigenvalues\n",
" if n_real != 0:\n",
" z_random_real = z_min + np.random.rand(n_real) * (z_max - z_min)\n",
" z_real = np.concatenate((z_real, z_random_real))\n",
" if ccf:\n",
" z = np.concatenate((z_complex, z_real))\n",
" a_matrix, b_matrix = make_ccf_system_from_eigenvalues(z)\n",
" else:\n",
" a_matrix, b_matrix = make_random_system_from_eigenvalues(z_real, z_complex)\n",
" return a_matrix, b_matrix\n",
"\n",
"\n",
"def make_random_system_from_eigenvalues(eigvals_real, eigvals_complex, n_inputs=1):\n",
" \"\"\"\n",
" Given a sequence of eigenvalues, make a random system at random that has those eigenvalues.\n",
" \n",
" arguments:\n",
" eigvals_real: sequence of real eigenvalues the system is to have.\n",
" eigvals_complex: sequence of complex eigenvalues the system is to have.\n",
" \n",
" returns:\n",
" a_matrix: A matrix of the system, in CCF\n",
" b_matrix: B matrix of the system, in CCF (i.e., all zeros except the last element, and the last element is 1)\n",
" \"\"\"\n",
" n_states = len(eigvals_real) + 2 * len(eigvals_complex)\n",
" diag_blocks = []\n",
" for z in eigvals_real:\n",
" z_block = np.array([[z]])\n",
" diag_blocks.append(z_block)\n",
" for z in eigvals_complex:\n",
" a, b = np.real(z), np.imag(z)\n",
" z_block = np.array([[a, b],\n",
" [-b, a]])\n",
" diag_blocks.append(z_block)\n",
" mcf_mat = scipy.linalg.block_diag(*diag_blocks) # modal canonical form\n",
" t_mat = np.random.rand(n_states, n_states) + np.identity(n_states) # random transformation matrix\n",
" a_matrix = np.linalg.inv(t_mat) @ mcf_mat @ t_mat\n",
" b_matrix = np.random.rand(n_states, n_inputs)\n",
" return a_matrix, b_matrix\n",
"\n",
"\n",
"def make_ccf_system_from_eigenvalues(eigvals):\n",
" \"\"\"\n",
" Given a sequence of eigenvalues, make a system in controllable canonical form that has those eigenvalues.\n",
" \n",
" arguments:\n",
" eigvals: sequence of eigenvalues\n",
" \n",
" returns:\n",
" a_matrix: A matrix of the system, in CCF\n",
" b_matrix: B matrix of the system, in CCF (i.e. all zeros except the last element, and the last element is 1)\n",
" \"\"\"\n",
" n_states = len(eigvals)\n",
" a_matrix = np.zeros((n_states, n_states))\n",
" coeffs = roots2coeffs(eigvals)\n",
" for i in range(n_states - 1):\n",
" a_matrix[i, i + 1] = 1\n",
" for i in range(n_states):\n",
" a_matrix[n_states - 1, n_states - 1 - i] = -np.real(coeffs[i])\n",
" b_matrix = np.zeros((n_states, 1))\n",
" b_matrix[n_states - 1, 0] = 1\n",
" return a_matrix, b_matrix\n",
"\n",
"\n",
"def roots2coeffs(roots):\n",
" \"\"\"\n",
" roots2coeffs: given a list of roots, find the coefficients of the monic polynomial with those roots. \n",
" The roots can be complex.\n",
" \n",
" Specifically, for given roots [r_1, r_2, ..., r_n], we have the monic polynomial\n",
" \n",
" p(s) = (s-r_1)*(s-r_2)*...*(s-r_n)\n",
" = s^n + a_{1}s^{n-1} + a_{2}s^{n-2} + ... + a_{n-1}s + a_n.\n",
" \n",
" This function returns the coefficients as list [a_{1}, a_{2}, ... a_{n}].\n",
" \n",
" arguments:\n",
" roots: the polynomial roots (list of complex numbers)\n",
" \n",
" returns:\n",
" coeffs: the coefficients, as described above (list)\n",
" \n",
" This function uses Vieta's formulas to calculate the coefficients. For additional info, check out the Wiki page:\n",
" https://en.wikipedia.org/wiki/Vieta%27s_formulas\n",
" \"\"\"\n",
" n = len(roots)\n",
" coeffs = np.zeros(n, dtype=complex)\n",
" for i in range(1, n + 1):\n",
" a = list(combinations(roots, i))\n",
" prods = np.prod(a, axis=1)\n",
" s = np.sum(prods)\n",
" coeffs[i - 1] = np.power(-1, i) * s\n",
" return coeffs\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In addition to this, it would also be nice to have a function that plots the *eigenvalues* of a given $A$ matrix. That will make it easy to visualize how well the system identification is working. The code below provides such a function. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def plot_eigenvalues(a_matrices, labels=None):\n",
" \"\"\"\n",
" plot_poles: given an list A matrix, find their eigenvalues and plot them on the complex plane. \n",
" Additionally, the unit circle is plotted as a dashed line. Additionally, each A matrix can be given a label, \n",
" which will be used to identify that matrix's eigenvalues in a legend.\n",
" \n",
" arguments:\n",
" a_matrices: the system A-matrices.\n",
" labels: list of legend entries, one for each A matrix.\n",
" \n",
" returns:\n",
" None\n",
" \"\"\"\n",
" markers=['o','x','1','8','*','+']\n",
" fig = plt.figure(figsize=(8, 8)) \n",
" for i, a_matrix in enumerate(a_matrices):\n",
" eigenvalues = np.linalg.eigvals(a_matrix)\n",
" if labels:\n",
" label = labels[i]\n",
" else:\n",
" label=''\n",
" plt.plot(np.real(eigenvalues), np.imag(eigenvalues), linestyle='None', marker=markers[i], ms=5, label=label)\n",
" t = np.linspace(0, 2 * np.pi, 1000)\n",
" plt.plot(np.cos(t), np.sin(t), 'k--', label=\"Unit Circle\")\n",
" plt.xlabel(\"Real Axis\")\n",
" plt.ylabel(\"Imaginary Axis\")\n",
" if labels:\n",
" plt.legend(loc='upper right')\n",
" plt.grid(True)\n",
" return None"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To try out our new functions, let's generate a random system and plot the eigenvalues of its $A$ matrix. Strictly speaking, there's no limit on how many states (and therefore eigenvalues) our random system can have, but you might want to keep it less than 20 for now."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a_rand_1, b_rand_1 = make_random_system(6, 1, z_min=0,z_max=0.2, must_have_real_eigenvalues=[0.3, -0.7])\n",
"a_rand_2, b_rand_2 = make_random_system(6, 1, z_min=0,z_max=0.2, must_have_real_eigenvalues=[0.31, -0.71])\n",
"plot_eigenvalues([a_rand_1, a_rand_2], labels=['Random system 1 eigenvalues', 'Random system 2 eigenvalues'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The plot should show both real systems having the two eigenvalues that we specified ($\\lambda=0.3$ and $\\lambda=-0.7$), and each one having a bunch of other random eigenvalues that all lie within $0 \\le |\\lambda| \\le 0.2$. Try running the above cell with some different parameters, and see what kind of random systems you can make!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 2: Identifying a system with using transformed data.\n",
"\n",
"The first caveat we mentioned is that the data you collect may not represent the true state of the system. Suppose that, instead of having the true state trace\n",
"$$ x[0], x[1],\\dotsc,x[\\ell],$$\n",
"you instead have a *transformed state trace*, that is\n",
"$$ Tx[0], Tx[1],\\dotsc,Tx[\\ell],$$\n",
"where $T\\in\\mathcal{R}^{n\\times n}$ is some full-rank *transformation matrix*. Suppose further that we didn't know that this transformation had occurred-- it could easily happen without our knowledge, or maybe we can't prevent it-- and we ran our system identification function on our transformed state trace and our input trace. What would we get? Well, we certainly wouldn't get the right $A$ and $B$. Instead, we would approximately identify\n",
"$$\n",
"\\begin{align}\n",
"A_{identified} &= TAT^{-1}\\\\\n",
"B_{identified} &= TB.\n",
"\\end{align}\n",
"$$\n",
"As an exercise, you should prove that this is what we want to happen (i.e. this is the right system matrices to identify if we think of these transformed states as the true states.)\n",
"\n",
"Let's look at a numerical example of this phenomenon. For now, we'll use the same $A$ and $B$ from the first example, but we'll define a transformation matrix\n",
"$$\n",
"T=\\begin{bmatrix}3 & 1 \\\\ 5 & 2 \\end{bmatrix}\n",
"$$\n",
"to apply to the state trace."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a_known = np.array([[0, 1],\n",
" [0.3, 0.2]])\n",
"b_known = np.array([0, 1])\n",
"t_mat = np.array([[3, 1],\n",
" [5, 2]])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We've defined the following simple function to apply a transformation to a state trace. This function can also supply a random transformation matrix, but we brought our own this time."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def transform_state_trace(state_trace, transformation_matrix=None):\n",
" \"\"\"\n",
" transform_state_trace: given a state trace, apply linear transformation to the states.\n",
" If no transformation is specified, one will be chosen at random.\n",
" \n",
" arguments:\n",
" state_trace: the state trace to be rotated.\n",
" transformation_matrix: if given, the matrix that specifies the transformation that is to be applied.\n",
" \n",
" returns:\n",
" transformed_state_trace: the state trace, with the transformation applied.\n",
" \"\"\"\n",
" n_states = state_trace.shape[1]\n",
" t_mat = np.random.rand(n_states, n_states) + np.identity(n_states) if transformation_matrix is None else transformation_matrix\n",
" transformed_state_trace = state_trace @ t_mat.T\n",
" return transformed_state_trace"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, let's generate a state trace, apply $T$ to it, and try to identify the system."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"k = 10\n",
"t = t = np.arange(k)\n",
"\n",
"u_trace = random_input(t, 1)\n",
"\n",
"x_trace = make_state_trace(a_known, b_known, u_trace)\n",
"tx_trace = transform_state_trace(x_trace, t_mat)\n",
"\n",
"a_identified, b_identified = identify_system(tx_trace, u_trace)\n",
"print(\"a_identified:\")\n",
"print(a_identified)\n",
"print(\"b_identified:\")\n",
"print(b_identified)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Take a look at the identified $A$ and $B$ matrix. Clearly they're wrong. Are they wrong in the way we expected?\n",
"\n",
"So, we didn't recover the supposedly underlying $A$ matrix when we use a transformed state trace. But in a way, we're not so far off: even though the identified system is wrong, the identified *eigenvalues* are actually right! We can verify this numerically:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(np.linalg.eigvals(a_known))\n",
"print(np.linalg.eigvals(a_identified))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But you should be able to explain, from your knowledge of *coordinate changes*, why this is the case. Actually, a lot more is preserved than just the eigenvalues, but those are clearly one of the most important things that are preserved."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 3: Identifying a system using noisy data.\n",
"\n",
"Similar to the above, we also often have nosiy measurements when we are observing some state $x[n]$. We will assume here that the noise is Gaussian-distributed with standard deviation $\\sigma$ and a mean $\\mu$ of zero.\n",
"\n",
"Our noisy observations will be given by\n",
"$$\\begin{align}\n",
"x_{noisy}[i] &= x[i] + w[i] \\\\ u_{noisy}[i] &= u[i] + z[i]\n",
"\\end{align}$$\n",
"where $w[i]$ and $z[i]$ are *mean-zero Gaussian noise* with some standard deviation $\\sigma$. In the presense of this noise, can we still identify $A$ and $B$?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following simple function takes a state trace and corrupts it by adding noise."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def corrupt_state_trace(state_trace, noise_std=1.0):\n",
" \"\"\"\n",
" corrupt_state_trace: corrupt a state trace by addding discrete-time Gaussian white noise to it.\n",
" \n",
" arguments:\n",
" state_trace: the state trace to be corrupted.\n",
" noise_std: standard deviation of the noise. You can think of this as the \"magnitude\" of the noise. \n",
" For example, if noise_std=1, then the noise will usually be between -3 and 3.\n",
" \n",
" returns: \n",
" corrupted_state_trace: the given state trace, plus the noise as specified.\n",
" \"\"\"\n",
" noise = np.random.normal(loc=0.0, scale=noise_std, size=np.shape(state_trace))\n",
" corrupted_state_trace = state_trace + noise\n",
" return corrupted_state_trace\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, on to the example. This time, we'll generate a random system to use as our \"known\" system."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n = 4 # number of state variables in the random system\n",
"a_known, b_known = make_random_system(n, 1, z_min=0, z_max=0.6)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, let's generate a state trace, add noise to it, and try to identify the system using the noisy data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"k = 1000\n",
"t = np.arange(k)\n",
"\n",
"sigma = 0.01 # how much noise there is\n",
"\n",
"u_trace = random_input(t, 1)\n",
"\n",
"x_trace = make_state_trace(a_known, b_known, u_trace)\n",
"noisy_x_trace = corrupt_state_trace(x_trace, noise_std=sigma)\n",
"\n",
"a_identified, b_identified = identify_system(noisy_x_trace, u_trace)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How did we do? Let's try plotting the eigenvalues to find out."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"plot_eigenvalues([a_known, a_identified], ['True eigenvalues', 'Identified eigenvalues'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the default values, it looks like the noise didn't affect our ability to identify the system too much. However, if the noise becomes too big, or if we don't take a long enough trace, or if we increase the number of state variables, our identified $A$ and $B$ can potentially be very poor estimates of the true ones."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 4: Deliberately Undermodeling a System.\n",
"\n",
"If you know that a system has $n$ states, of course the most accurate model of the system will reflect this. However, it may be the case that you have more states than are really necessary to accurately capture the dynamics of your system. In this case, it would be nice if you could discard some states. If you start off with a large amount of states (say $n=50$), it might even be necessary.\n",
"\n",
"Here we will explore what happens when we *undermodel* a system. This is when we choose to model some subset of $q