docs/notebooks/research/linear_model_reduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear System Model Reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import numpy as np\n",
    "import pylab\n",
    "try:\n",
    "    import seaborn as sns  # optional; prettier graphs\n",
    "except ImportError:\n",
    "    pass\n",
    "\n",
    "from scipy.linalg import cholesky, svd\n",
    "from nengo.utils.numpy import rmse\n",
    "\n",
    "import nengolib\n",
    "from nengolib import Lowpass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem Statement\n",
    "\n",
    "Suppose we have some linear system. We'd like to reduce its order while maintaining similar characteristics. Take, for example, a lowpass filter that has a small amount of 3rd-order dynamics mixed in, resulting in a 4th-order system that consists mostly of 1st-order dynamics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "isys = Lowpass(0.05)\n",
    "noise = 0.5*Lowpass(0.2) + 0.25*Lowpass(0.007) - 0.25*Lowpass(0.003)\n",
    "p = 0.8\n",
    "sys = p*isys + (1-p)*noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exact Minimal Realizations\n",
    "\n",
    "By cancelling repeated zeros and poles from a system, we can obtain an exact version of that same system with potentially reduced order. However, for the above system, there are no poles to be cancelled, and so this does not help to reduce the order. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert nengolib.signal.pole_zero_cancel(isys/isys) == 1  # demonstration\n",
    "\n",
    "minsys = nengolib.signal.pole_zero_cancel(sys)\n",
    "assert minsys == sys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a crude way of getting around this problem, we can raise the tolerance for detecting similar poles/zeros until repeats are found. By setting the tolerance appropriately for this example, we can reduce the model to a first-order filter, with a surprisingly similar response. However, as we will soon see further down, we can do *much* better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minsys_crude = nengolib.signal.pole_zero_cancel(sys, tol=1000.0)\n",
    "assert minsys_crude.order_den == 1\n",
    "\n",
    "def test_sys(u, redsys, dt=0.001):\n",
    "    orig = sys.filt(u, dt)\n",
    "    red = redsys.filt(u, dt)\n",
    "    \n",
    "    pylab.figure()\n",
    "    pylab.title(\"(RMSE: %s)\" % rmse(orig, red))\n",
    "    pylab.plot(orig, label=\"Original\")\n",
    "    pylab.plot(red, label=\"Reduced\")\n",
    "    pylab.legend()\n",
    "    pylab.show()\n",
    "\n",
    "rng = np.random.RandomState(0)\n",
    "white = rng.normal(size=100)\n",
    "test_sys(white, minsys_crude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Balanced Realizations\n",
    "\n",
    "First we need to compute some special matrices from Lyapunov theory.\n",
    "\n",
    "The \"controllability gramian\" (a.k.a. \"reachability gramian\" for linear systems) $W_r$ satisfies: \n",
    "\n",
    "$$AW_r + W_rA' = -BB'$$\n",
    "\n",
    "The \"observability gramian\" $W_o$ satisfies:\n",
    "\n",
    "$$A'W_o + W_oA = -C'C$$\n",
    "\n",
    "See [2] for more background information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A, B, C, D = sys.ss\n",
    "\n",
    "R = nengolib.signal.control_gram(sys)\n",
    "assert np.allclose(np.dot(A, R) + np.dot(R, A.T), -np.dot(B, B.T))\n",
    "\n",
    "O = nengolib.signal.observe_gram(sys)\n",
    "assert np.allclose(np.dot(A.T, O) + np.dot(O, A), -np.dot(C.T, C))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm from [3] computes the lower cholesky factorizations of $W_r \\, ( = L_rL_r')$ and $W_o \\, ( = L_oL_o')$, and the singular value decomposition of $L_o'L_r$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LR = cholesky(R, lower=True)\n",
    "assert np.allclose(R, np.dot(LR, LR.T))\n",
    "\n",
    "LO = cholesky(O, lower=True)\n",
    "assert np.allclose(O, np.dot(LO, LO.T))\n",
    "\n",
    "M = np.dot(LO.T, LR)\n",
    "U, S, V = svd(M)\n",
    "assert np.allclose(M, np.dot(U * S, V))\n",
    "\n",
    "T = np.dot(LR, V.T) * S ** (-1. / 2)\n",
    "Tinv = (S ** (-1. / 2))[:, None] * np.dot(U.T, LO.T)\n",
    "assert np.allclose(np.dot(T, Tinv), np.eye(len(T)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives the similarity transform [4] that maps the state to the \"balanced realization\" of the same order. This system is equivalent up to a change of basis $T$ in the state-space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TA, TB, TC, TD = sys.transform(T, Tinv=Tinv).ss\n",
    "assert sys == (TA, TB, TC, TD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the reason we do this is because the singular values reflect a measure of importance for each of the states in the new realization. The order should then be reduced by removing the least important states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pylab.figure()\n",
    "pylab.plot(S)\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The short-cut to do the above procedure in `nengolib` is the function `balanced_transformation` followed by `sys.transform`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T, Tinv, S_check = nengolib.signal.balanced_transformation(sys)\n",
    "sys_check = sys.transform(T, Tinv)\n",
    "\n",
    "assert np.allclose(TA, sys_check.A)\n",
    "assert np.allclose(TB, sys_check.B)\n",
    "assert np.allclose(TC, sys_check.C)\n",
    "assert np.allclose(TD, sys_check.D)\n",
    "assert np.allclose(S, S_check)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, note that this diagonalizes the two gramiam matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "P = nengolib.signal.control_gram((TA, TB, TC, TD))\n",
    "Q = nengolib.signal.observe_gram((TA, TB, TC, TD))\n",
    "\n",
    "diag = np.diag_indices(len(P))\n",
    "offdiag = np.ones_like(P, dtype=bool)\n",
    "offdiag[diag] = False\n",
    "offdiag = np.where(offdiag)\n",
    "\n",
    "assert np.allclose(P[diag], S)\n",
    "assert np.allclose(P[offdiag], 0)\n",
    "assert np.allclose(Q[diag], S)\n",
    "assert np.allclose(Q[offdiag], 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Approximate Model Order Reduction\n",
    "\n",
    "Low singular values indicate states are less important. The method in [5] can be used to remove these states while matching the DC gain for continuous or discrete systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "redsys = nengolib.signal.modred((TA, TB, TC, TD), 0, method='dc')\n",
    "assert redsys.order_den == 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "step = np.zeros(1000)\n",
    "step[50:] = 1.0\n",
    "test_sys(step, redsys)\n",
    "test_sys(white, redsys)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, this doesn't work very well for matching the response of the system given white-noise input. If we care less about the steady-state response, then it is much more accurate to simply delete the less important states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delsys = nengolib.signal.modred((TA, TB, TC, TD), 0, method='del')\n",
    "assert delsys.order_den == 1\n",
    "test_sys(step, delsys)\n",
    "test_sys(white, delsys)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The short-cut for all of the above is to call `nengolib.signal.balred` with a desired order and method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "[1] http://www.mathworks.com/help/control/ref/minreal.html\n",
    "\n",
    "[2] https://en.wikibooks.org/wiki/Control_Systems/Controllability_and_Observability\n",
    "\n",
    "[3] Laub, A.J., M.T. Heath, C.C. Paige, and R.C. Ward, \"Computation of System Balancing Transformations and Other Applications of Simultaneous Diagonalization Algorithms,\" *IEEE® Trans. Automatic Control*, AC-32 (1987), pp. 115-122.\n",
    "\n",
    "[4] http://www.mathworks.com/help/control/ref/balreal.html\n",
    "\n",
    "[5] http://www.mathworks.com/help/control/ref/modred.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}