{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fault detection\n",
    "\n",
    "We'll consider a problem of identifying faults that have occurred in a system based on sensor measurements of system performance.\n",
    "\n",
    "# Topic references\n",
    "- [Samar, Sikandar, Dimitry Gorinevsky, and Stephen Boyd. \"Likelihood Bounds for Constrained Estimation with Uncertainty.\" Decision and Control, 2005 and 2005 European Control Conference. CDC-ECC'05. 44th IEEE Conference on. IEEE, 2005.](http://web.stanford.edu/~boyd/papers/pdf/map_bounds.pdf)\n",
    "\n",
    "# Problem statement\n",
    "\n",
    "Each of $n$ possible faults occurs independently with probability $p$.\n",
    "The vector $x \\in \\lbrace 0,1 \\rbrace^{n}$ encodes the fault occurrences, with $x_i = 1$ indicating that fault $i$ has occurred.\n",
    "System performance is measured by $m$ sensors.\n",
    "The sensor output is\n",
    "\\begin{equation}\n",
    "y = Ax + v = \\sum_{i=1}^n a_i x_i + v,\n",
    "\\end{equation}\n",
    "where $A \\in \\mathbf{R}^{m \\times n}$ is the sensing matrix with column $a_i$ being the **fault signature** of fault $i$,\n",
    "and $v \\in \\mathbf{R}^m$ is a noise vector where $v_j$ is Gaussian with mean 0 and variance $\\sigma^2$.\n",
    "\n",
    "The objective is to guess $x$ (which faults have occurred) given $y$ (sensor measurements).\n",
    "\n",
    "We are interested in the setting where $n > m$, that is, when we have more possible faults than measurements.\n",
    "In this setting, we can expect a good recovery when the vector $x$ is sparse.\n",
    "This is the subject of compressed sensing.\n",
    "\n",
    "# Solution approach\n",
    "To identify the faults, one reasonable approach is to choose $x \\in \\lbrace 0,1 \\rbrace^{n}$ to minimize the negative log-likelihood function\n",
    "\n",
    "\\begin{equation}\n",
    "\\ell(x) = \\frac{1}{2 \\sigma^2} \\|Ax-y\\|_2^2 +  \\log(1/p-1)\\mathbf{1}^T x + c.\n",
    "\\end{equation}\n",
    "\n",
    "However, this problem is nonconvex and NP-hard, due to the constraint that $x$ must be Boolean.\n",
    "\n",
    "To make this problem tractable, we can relax the Boolean constraints and instead constrain $x_i \\in [0,1]$.\n",
    "\n",
    "The optimization problem\n",
    "\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} &  \\|Ax-y\\|_2^2 + 2 \\sigma^2 \\log(1/p-1)\\mathbf{1}^T x\\\\\n",
    "\\mbox{subject to} &  0 \\leq x_i \\leq 1, \\quad i=1, \\ldots n\n",
    "\\end{array}\n",
    "\n",
    "is convex.\n",
    "We'll refer to the solution of the convex problem as the **relaxed ML** estimate.\n",
    "\n",
    "By taking the relaxed ML estimate of $x$ and rounding the entries to the nearest of 0 or 1, we recover a Boolean estimate of the fault occurrences.\n",
    "\n",
    "# Example\n",
    "\n",
    "We'll generate an example with $n = 2000$ possible faults, $m = 200$ measurements, and fault probability $p = 0.01$.\n",
    "We'll choose $\\sigma^2$ so that the signal-to-noise ratio is 5.\n",
    "That is,\n",
    "\\begin{equation}\n",
    "\\sqrt{\\frac{\\mathbf{E}\\|Ax \\|^2_2}{\\mathbf{E} \\|v\\|_2^2}} = 5.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.random.seed(1)\n",
    "\n",
    "n = 2000\n",
    "m = 200\n",
    "p = 0.01\n",
    "snr = 5\n",
    "\n",
    "sigma = np.sqrt(p*n/(snr**2))\n",
    "A = np.random.randn(m,n)\n",
    "\n",
    "x_true = (np.random.rand(n) <= p).astype(int)\n",
    "v = sigma*np.random.randn(m)\n",
    "\n",
    "y = A.dot(x_true) + v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we show $x$, $Ax$ and the noise $v$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11ae42518>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 8672 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(n),x_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11aee9630>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 52364 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(m), A.dot(x_true),range(m),v)\n",
    "plt.legend(('Ax','v'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recovery\n",
    "\n",
    "We solve the relaxed maximum likelihood problem with CVXPY and then round the result to get a Boolean solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS\n",
      "\n",
      "It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT\n",
      " 0  +7.343e+03  -3.862e+03  +5e+04  5e-01  5e-04  1e+00  1e+01    ---    ---    1  1  - |  -  - \n",
      " 1  +4.814e+02  -9.580e+02  +8e+03  1e-01  6e-05  2e-01  2e+00  0.8500  1e-02   1  2  2 |  0  0\n",
      " 2  -2.079e+02  -1.428e+03  +6e+03  1e-01  4e-05  8e-01  2e+00  0.7544  7e-01   2  2  2 |  0  0\n",
      " 3  -1.321e+02  -1.030e+03  +5e+03  8e-02  3e-05  7e-01  1e+00  0.3122  2e-01   2  2  2 |  0  0\n",
      " 4  -2.074e+02  -8.580e+02  +4e+03  6e-02  2e-05  6e-01  9e-01  0.7839  7e-01   2  2  2 |  0  0\n",
      " 5  -1.121e+02  -6.072e+02  +3e+03  5e-02  1e-05  5e-01  7e-01  0.3859  4e-01   2  3  3 |  0  0\n",
      " 6  -4.898e+01  -4.060e+02  +2e+03  3e-02  8e-06  3e-01  5e-01  0.5780  5e-01   2  2  2 |  0  0\n",
      " 7  +7.778e+01  -5.711e+01  +8e+02  1e-02  3e-06  1e-01  2e-01  0.9890  4e-01   2  3  2 |  0  0\n",
      " 8  +1.307e+02  +6.143e+01  +4e+02  6e-03  1e-06  6e-02  1e-01  0.5528  1e-01   3  3  3 |  0  0\n",
      " 9  +1.607e+02  +1.286e+02  +2e+02  3e-03  4e-07  3e-02  5e-02  0.8303  3e-01   3  3  3 |  0  0\n",
      "10  +1.741e+02  +1.557e+02  +1e+02  2e-03  2e-07  2e-02  3e-02  0.6242  3e-01   3  3  3 |  0  0\n",
      "11  +1.834e+02  +1.749e+02  +5e+01  8e-04  9e-08  8e-03  1e-02  0.8043  3e-01   3  3  3 |  0  0\n",
      "12  +1.888e+02  +1.861e+02  +2e+01  3e-04  3e-08  2e-03  4e-03  0.9175  3e-01   3  3  2 |  0  0\n",
      "13  +1.909e+02  +1.902e+02  +4e+00  7e-05  7e-09  6e-04  1e-03  0.8198  1e-01   3  3  3 |  0  0\n",
      "14  +1.914e+02  +1.912e+02  +1e+00  2e-05  2e-09  2e-04  3e-04  0.8581  2e-01   3  2  3 |  0  0\n",
      "15  +1.916e+02  +1.916e+02  +1e-01  2e-06  3e-10  2e-05  4e-05  0.9004  3e-02   3  3  3 |  0  0\n",
      "16  +1.916e+02  +1.916e+02  +4e-02  7e-07  8e-11  7e-06  1e-05  0.8174  1e-01   3  3  3 |  0  0\n",
      "17  +1.916e+02  +1.916e+02  +8e-03  1e-07  1e-11  1e-06  2e-06  0.8917  9e-02   3  2  2 |  0  0\n",
      "18  +1.916e+02  +1.916e+02  +2e-03  4e-08  4e-12  4e-07  5e-07  0.8588  2e-01   3  3  3 |  0  0\n",
      "19  +1.916e+02  +1.916e+02  +2e-04  3e-09  3e-13  3e-08  5e-08  0.9309  2e-02   3  2  2 |  0  0\n",
      "20  +1.916e+02  +1.916e+02  +2e-05  4e-10  4e-14  4e-09  6e-09  0.8768  1e-02   4  2  2 |  0  0\n",
      "21  +1.916e+02  +1.916e+02  +4e-06  6e-11  6e-15  6e-10  9e-10  0.9089  6e-02   4  2  2 |  0  0\n",
      "22  +1.916e+02  +1.916e+02  +1e-06  2e-11  2e-15  2e-10  2e-10  0.8362  1e-01   2  1  1 |  0  0\n",
      "\n",
      "OPTIMAL (within feastol=1.8e-11, reltol=5.1e-09, abstol=9.8e-07).\n",
      "Runtime: 6.538894 seconds.\n",
      "\n",
      "final objective value: 191.6347201927456\n",
      "CPU times: user 6.51 s, sys: 291 ms, total: 6.8 s\n",
      "Wall time: 7.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import cvxpy as cp\n",
    "x = cp.Variable(shape=n)\n",
    "tau = 2*cp.log(1/p - 1)*sigma**2\n",
    "obj = cp.Minimize(cp.sum_squares(A*x - y) + tau*cp.sum(x))\n",
    "const = [0 <= x, x <= 1]\n",
    "cp.Problem(obj,const).solve(verbose=True)\n",
    "print(\"final objective value: {}\".format(obj.value))\n",
    "\n",
    "# relaxed ML estimate\n",
    "x_rml = np.array(x.value).flatten()\n",
    "\n",
    "# rounded solution\n",
    "x_rnd = (x_rml >= .5).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluation\n",
    "\n",
    "We define a function for computing the estimation errors, and a function for plotting $x$, the relaxed ML estimate, and the rounded solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "\n",
    "def errors(x_true, x, threshold=.5):\n",
    "    '''Return estimation errors.\n",
    "    \n",
    "    Return the true number of faults, the number of false positives, and the number of false negatives.\n",
    "    '''\n",
    "    n = len(x_true)\n",
    "    k = sum(x_true)\n",
    "    false_pos = sum(np.logical_and(x_true < threshold, x >= threshold))\n",
    "    false_neg = sum(np.logical_and(x_true >= threshold, x < threshold))\n",
    "    return (k, false_pos, false_neg)\n",
    "\n",
    "def plotXs(x_true, x_rml, x_rnd, filename=None):\n",
    "    '''Plot true, relaxed ML, and rounded solutions.'''\n",
    "    matplotlib.rcParams.update({'font.size': 14})\n",
    "    xs = [x_true, x_rml, x_rnd]\n",
    "    titles = ['x_true', 'x_rml', 'x_rnd']\n",
    "\n",
    "    n = len(x_true)\n",
    "    k = sum(x_true)\n",
    "\n",
    "    fig, ax = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12, 3))\n",
    "\n",
    "    for i,x in enumerate(xs):\n",
    "            ax[i].plot(range(n), x)\n",
    "            ax[i].set_title(titles[i])\n",
    "            ax[i].set_ylim([0,1])\n",
    "            \n",
    "    if filename:\n",
    "        fig.savefig(filename, bbox_inches='tight')\n",
    "        \n",
    "    return errors(x_true, x_rml,.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that out of 20 actual faults, the rounded solution gives perfect recovery with 0 false negatives and 0 false positives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 0, 0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 20804 characters]\n",
      "text/plain": [
       "<Figure size 864x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotXs(x_true, x_rml, x_rnd, 'fault.pdf')"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
