DMC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementing Dynamic Matrix Control in Python\n",
    "\n",
    "#### Author: Siang Lim, February 2020\n",
    "\n",
    "** work in progress **\n",
    "\n",
    "- Contents of this notebook are derived from DMC literature:\n",
    "    - Hokanson, D. A., & Gerstle, J. G. (1992). **Dynamic Matrix Control Multivariable Controllers.** Practical Distillation Control, 248–271.\n",
    "    - Sorensen, R. C., & Cutler, C. R. (1998). **LP integrates economics into dynamic matrix control.** Hydrocarbon processing (International ed.), 77(9), 57-65.\n",
    "    - Morshedi, A. M., Cutler, C. R., & Skrovanek, T. A. (1985, June). **Optimal solution of dynamic matrix control with linear programing techniques (LDMC).** In 1985 American Control Conference (pp. 199-208). IEEE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Definitions\n",
    "\n",
    "#### Independent Variables\n",
    "\n",
    "- MV (Manipulated Variables)\n",
    "    - Variables that we can move directly, valves, feed rate etc.\n",
    "\n",
    "- FF (Feedforward Variables)\n",
    "    - Disturbances that we can't control directly - ambient temperature.\n",
    "\n",
    "#### Dependent Variables\n",
    "- Variables that change due to a change in the independent variables, e.g. temperatures, pressures\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step Response Convolution Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sequence $a_j$, where $j = 0,1,2,3 \\dots$ is the step response convolution model, which is the basis for DMC.\n",
    "\n",
    "Remembering that everything so far is in deviation variables, at is no more than a series of numbers based on\n",
    "a set length (time to steady state) and a set sampling interval.\n",
    "\n",
    "There is no set form of the model, other than that steady state is reached by the last interval. In fact, the last\n",
    "value $a_t$ is the process steady-state gain. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = -np.array([0.5,1.0,1.5,1.8,1.9,1.95,1.98,1.99,1.999,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Response')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4XVW9//H3p01LG2yBatSW0pah\nCJULlCeIDKIyKAI/aquIGBRQrHhlusogQwcoZb5c5AGRgvBjqCCjoPQnMl7A60VSBikUKFMHWqGI\njC20pd/fH+vUJiVJz0lyzjpJPq/n2U/O3tk5+9PztPl2rb32WooIzMzMitUrdwAzM+taXDjMzKwk\nLhxmZlYSFw4zMyuJC4eZmZXEhcPMzEriwmFmZiVx4TAzs5K4cJiZWUlqcgcoh0984hMxYsSI3DHM\nzLqMmTNnvh4RdcWc2y0Lx4gRI2hsbMwdw8ysy5A0t9hz3VVlZmYlceEwM7OSuHCYmVlJXDjMzKwk\nLhxmZlYSF45Vpk+HESOgV6/0dfr03InMzKpStxyOW7Lp02H8eFiyJO3PnZv2ARoa8uUyM6tCbnEA\nnHzy6qKxypIl6biZmTXjwgEwb15px83MejAXDoBhw0o7bmbWg7lwAEydCrW1zY/V1qbjZmbWjAsH\npBvg06bB8OGrjx1+uG+Mm5m1wIVjlYYGePllWL4cRo6Eu++GlStzpzIzqzouHGuqqYHJk+Fvf4Ob\nb86dxsys6rhwtOSAA2DUKJg0CT78MHcaM7Oq4sLRkt694dRTYfZsuO663GnMzKqKC0drxo2DbbdN\nBWT58txpzMyqhgtHa3r1gtNOg+efh6uvzp3GzKxqZCkckgZJukvSnMLXDVo5b5ikP0maLelpSSMq\nGnTffeFzn4MpU2DZsope2sysWuVqcfwcuCciRgL3FPZbcjVwbkRsCXwOeK1C+RIpFY25c+HXv67o\npc3MqlWuwjEGuKrw+irg62ueIGkUUBMRdwFExLsRsWTN88puzz1hl13g9NNh6dKKX97MrNrkKhyf\niohFAIWvn2zhnM2BNyXdIukxSedK6t3aG0oaL6lRUuPixYs7L6mUisbChXDppZ33vmZmXVTZCoek\nuyXNamEbU+Rb1ABfAI4Ftgc2AQ5p7eSImBYR9RFRX1dX1+H8zXzxi7D77nDmmfDee5373mZmXUzZ\nCkdE7BERW7Ww3Qa8KmkwQOFrS/cuFgCPRcSLEbEC+B2wXbnyrtWUKfDaa3DRRdkimJlVg1xdVbcD\nBxdeHwzc1sI5jwAbSFrVfNgNeLoC2Vq2446w995wzjnw9tvZYpiZ5ZarcJwF7ClpDrBnYR9J9ZIu\nB4iID0ndVPdIehIQcFmmvMlpp8Ebb8AFF2SNYWaWkyIid4ZOV19fH42NjeV583Hj4N574aWXYIMW\nHz8xM+tyJM2MiPpizvWT46U69dTUVfWf/5k7iZlZFi4cpfq3f4NvfSt1V3XmsF8zsy7ChaM9Jk9O\nDwOec07uJGZmFefC0R5bbAEHHZSG5i5alDuNmVlFuXC018SJsGJFeijQzKwHceFor003hUMPTdOQ\nzJ+fO42ZWcW4cHTEKaekr6efnjeHmVkFuXB0xLBhMH48XHEFvPhi7jRmZhXhwtFRJ50ENTXpqXIz\nsx7AhaOjBg+Gn/wErrkGnn02dxozs7Jz4egMJ5wA/funp8rNzLo5F47OUFcHRx0F118Ps2blTmNm\nVlYuHJ3l2GNhwACYNCl3EjOzsnLh6CyDBsFPfwq33AKPPpo7jZlZ2bhwdKZjjklTrU+cmDuJmVnZ\nuHB0pvXWg+OPhzvugP/939xpzMzKwoWjsx1xRLpZPmFC7iRmZmWRrXBIGiTpLklzCl9bXE5P0jmS\nnpI0W9KFklTprCX52MfgxBPh7rvhgQdypzEz63Q5Wxw/B+6JiJHAPYX9ZiTtBOwMbA1sBWwPfLGS\nIdvl8MNhyJA0l1U3XJrXzHq2nIVjDHBV4fVVwNdbOCeAfkBfYB2gD/BqRdJ1RP/+aSqSBx9MLQ8z\ns24kZ+H4VEQsAih8/eSaJ0TEX4D7gEWF7c6ImF3RlO112GFpEsQJE9zqMLNupayFQ9Ldkma1sI0p\n8uc3A7YEhgIbArtJ2rWVc8dLapTUuLga1gJfZ51UNB5+OI2yMjPrJhSZ/jcs6VngSxGxSNJg4P6I\n+Mwa5xwH9IuIKYX9icD7EdHmYt/19fXR2NhYrujFW74cttwSBg6EmTOhyu/rm1nPJWlmRNQXc27O\nrqrbgYMLrw8GbmvhnHnAFyXVSOpDujHeNbqqAPr0SVOQPPYY3Hpr7jRmZp0iZ+E4C9hT0hxgz8I+\nkuolXV445ybgBeBJ4AngiYj4fY6w7fad78AWW6SnyT/8MHcaM7MOy9ZVVU5V01W1yg03wAEHwG9+\nAwcemDuNmdlHdJWuqp7jm9+ErbdO3VYrVuROY2bWIS4cldCrV1pads4cuPba3GnMzDrEhaNS9tsP\n6utTAVm2LHcaM7N2c+GoFCkVjZdegiuvzJ3GzKzdXDgqaa+9YKed4PTT4f33c6cxM2sXF45KkmDK\nFFiwAC67LHcaM7N2ceGotN12gy99CaZOhSVLcqcxMyuZC0cOU6bAq6/CL3+ZO4mZWclcOHLYZRf4\n6lfhrLPgnXdypzEzK4kLRy5TpsA//gEXXpg7iZlZSVw4ctl++/Rsx3nnwZtv5k5jZlY0F46cTjst\nFY3zz8+dxMysaC4cOW2zDey/P1xwAbz+eu40ZmZFceHIbfJkePddOPfc3EnMzIriwpHbqFHQ0AAX\nXZSG6JqZVTkXjmowaRJ88EEanmtmVuVcOKrBZpvBwQfDJZek6UjMzKpYlsIhaX9JT0laKanVFack\n7SXpWUnPS/p5JTNW3IQJsHIlnHFG7iRmZm3K1eKYBYwDHmjtBEm9gYuBrwGjgAMljapMvAxGjIDD\nDoNLL4WhQ9PiTyNGwPTpuZOZmTVTk+OiETEbQFJbp30OeD4iXiycez0wBni67AFzGTUqtTpeeSXt\nz50L48en1w0N+XKZmTVRzfc4NgTmN9lfUDjWfZ133kePLVkCJ59c+SxmZq0oW4tD0t3Ap1v41skR\ncVsxb9HCsWjjeuOB8QDDhg0rKmPVmTevtONmZhmUrXBExB4dfIsFwEZN9ocCC9u43jRgGkB9fX2r\nBaaqDRuWuqdaOm5mViWquavqEWCkpI0l9QW+DdyeOVN5TZ0KtbXNj9XWpuNmZlUi13DcsZIWADsC\nd0i6s3B8iKQZABGxAjgCuBOYDdwQEU/lyFsxDQ0wbRoMH7762Mkn+8a4mVUVRXTNXp221NfXR2Nj\nY+4YHfPmm7DxxrDrrnBbMbeEzMzaT9LMiGj1ubqmqrmrqmdbf3049li4/Xb4619zpzEz+xcXjmp2\n1FHwiU/AxIm5k5iZ/YsLRzUbMABOOAHuvBMeeih3GjMzwIWj+v37v8OnP53msjIzqwIuHNWuthZO\nOgnuvx/uvTd3GjMzF44u4Yc/TBMfnnIKdMNRcGbWtbhwdAX9+qWi8Ze/wB//mDuNmfVwRRUOSbWS\nJki6rLA/UtK+5Y1mzRx6aHquY8IEtzrMLKtiWxxXAh+QnvSGNI/U6WVJZC3r2zcNy5050w8EmllW\nxRaOTSPiHGA5QEQspeXZa62cDjoINt88FZCVK3OnMbMeqtjCsUxSfwrTmkvalNQCsUqqqYHJk+HJ\nJ+HGG3OnMbMeqtjCMQn4I7CRpOnAPcDxZUtlrTvgAPjsZ1MB+fDD3GnMrAcqqnBExF2kNcIPAa4D\n6iPi/vLFslb16gWnnQbPPAO/+U3uNGbWAxU7qmpn4P2IuANYHzhJ0vC1/JiVy9ixMHp0anUsX547\njZn1MMV2VV0CLJG0DXAcMBe4umyprG0STJkCL74IV12VO42Z9TDFFo4VkRbuGANcGBG/AAaUL5at\n1d57ww47pALygccpmFnlFFs43pF0InAQacW+3kCf8sWytVrV6pg3Dy6/PHcaM+tBii0cB5CG3/4g\nIv4ObAicW7ZUVpw99kgrBE6dCkuX5k5jZj1EsaOq/h4R50fEg4X9eRHR7nsckvaX9JSklZJaXKpQ\n0kaS7pM0u3Du0e29Xre1qtWxaBFccknuNGbWQxQ7qmqcpDmS3pL0tqR3JL3dgevOIg3vfaCNc1YA\nP4uILYHPAz+RNKoD1+yedt01tTzOOgvefTd3GjPrAYrtqjoH2C8i1ouIgRExICIGtveiETE7Ip5d\nyzmLIuLRwut3gNmkLjJb05QpsHgxXHRR7iRm1gMUWzhejYjZZU3SBkkjgNHAw22cM15So6TGxYsX\nVypadfj852GffeCcc+Ctt3KnMbNurtjC0Sjpt5IOLHRbjZM0rq0fkHS3pFktbGNKCSjpY8DNwDER\n0Wr3WERMi4j6iKivq6sr5RLdw2mnwT//CRdckDuJmXVzNUWeNxBYAnylybEAbmntByJijw7kAkBS\nH1LRmB4RrV7LgO22g3Hj4Pzz4cgjYdCg3InMrJsqqnBExKHlDrImSQJ+DcyOiPMrff0u6dRT4dZb\n4bzz4Iwzcqcxs26q2FFVQyXdKuk1Sa9KulnS0PZeVNJYSQtIC0PdIenOwvEhkmYUTtsZ+C6wm6TH\nC9ve7b1mj7DVVvDtb8OFF8Jrr+VOY2bdVCkrAN4ODCGNbPp94Vi7RMStETE0ItaJiE9FxFcLxxdG\nxN6F1w9FhCJi64jYtrDNaPudjUmT0sOAZ5+dO4mZdVPFFo66iLgyIlYUtv8L9MA70F3AZz4D3/se\n/PKXsHBh7jRm1g0VWzhel3SQpN6F7SDgH+UMZh0wcSKsWOH7HGZWFsUWju8D3wL+Xti+WThm1Wjj\njeH734fLLkuTIJqZdaJi56qaFxH7RURdYft6RMwtdzjrgFNOSV9PPz1vDjPrdoodVbWJpN9LWlwY\nWXWbpE3KHc46YKON4Ec/giuugBdeyJ3GzLqRYruqfgPcAAwmjay6kbT2uFWzE0+EPn3SU+VmZp2k\n2MKhiLimyaiqa0lPjls1GzwYjjgCrr0Wnnkmdxoz6yaKLRz3Sfq5pBGShks6nvTg3iBJntuimh1/\nPPTvD5Mn505iZt1EKSsA/gi4D7gf+DFpVNVMoLEsyaxz1NXBMcfAb38Lf/tb7jRm1g0UO6pq4zY2\n3ySvdj/7Gay3Xnqq3Mysg4odVbW/pAGF16dIukXS6PJGs06zwQbw05/C734HM2fmTmNmXVyxXVUT\nIuIdSbsAXwWuAn5VvljW6Y45Jk21PnFi7iRm1sUVWzg+LHzdB7gkIm4D+pYnkpXFwIHpRvmMGfCX\nv+ROY2ZdWLGF4xVJl5KmHZkhaZ0SftaqxRFHwCc/CRMm5E5iZl1Ysb/8vwXcCewVEW8Cg4DjypbK\nymPdddNDgffcA/ffnzuNmXVRxY6qWgK8BuxSOLQCmFOuUFZGhx8OQ4akVkf4GU4zK12xo6omAScA\nJxYO9QGuLVcoK6N+/dIEiA89BHfdlTuNmXVBxXZVjQX2A96DtFIfMKC9Fy0M731K0kpJ9Ws5t7ek\nxyT9ob3XszX84AcwfHgqIG51mFmJii0cyyIiKMxPJWndDl53FjAOeKCIc48GZnfwetZU376pq+qR\nR+APrsdmVppiC8cNhVFV60v6IXA3cHl7LxoRsyPi2bWdJ2koaQhwu69lrfje92CzzVIBWbkydxoz\n60KKvTl+HnATcDPwGWBiRFxYzmAFFwDHA/7N1tn69ElTkDzxBNxyS+40ZtaFFP0sRkTcFRHHRcSx\nwL2SGto6X9Ldkma1sI0p5nqS9gVei4ii5siQNF5So6TGxYsXF/MjduCBsOWWqYB8+OHazzczYy2F\nQ9JASSdKukjSV5QcAbxIerajVRGxR0Rs1cJ2W5HZdgb2k/QycD2wm6RWR3JFxLSIqI+I+rq6uiIv\n0cP17g2nngpPPw3XX587jZl1EYo2RtVIug34J/AXYHdgA9JUI0dHxOMdvrh0P3BsRLQ5NbukLxXO\n27eY962vr4/GRs/2XpSVK2H0aFi6NBWQmprcicwsA0kzI6LNUa6rrK2rapOIOCQiLgUOBOqBfTta\nNCSNlbQA2JG0INSdheNDJM3oyHtbiXr1gilTYM4cuOaa3GnMrAtYW4vj0YjYrrX9auUWR4kiYIcd\n4LXX4Lnn0nBdM+tROrPFsY2ktwvbO8DWq15LervjUa0qSKnVMXcuXHFF7jRmVuXaLBwR0TsiBha2\nARFR0+T1wEqFtAr4yldg553h9NPh/fdzpzGzKuap0S1Z1ep45RW49NLcacysirlw2Gpf/jLsthuc\neSa8917uNGZWpVw4rLkpU+DVV+Hii3MnMbMq5cJhze20E2y9dVrwqVcvGDECpk/PncrMqoif9rLm\npk9PQ3JXTXw4dy6MH59eN7Q5y4yZ9RBucVhzJ5/80VFVS5ak42ZmuHDYmubNK+24mfU4LhzW3LBh\nLR/faKPK5jCzquXCYc1NnQq1tR89vvXWlc9iZlXJhcOaa2iAadPSmuRSaoHsvntaYvZyL8RoZh5V\nZS1paGg+gmrFCth3X/jxj2HjjVMhMbMeyy0OW7uaGrjhBthiC/jGN9K6HWbWY7lwWHEGDkzdVf36\nwT77pCnYzaxHcuGw4g0fDr//fZqSZMyYtGqgmfU4LhxWmu23h2uvhYcfhkMOWf2EuZn1GC4cVrpx\n4+Dss9N9jwkTcqcxswrLUjgk7S/pKUkrJbW6VKGk9SXdJOkZSbMl7VjJnNaGY4+FH/4QzjgDrrwy\ndxozq6BcLY5ZwDjggbWc9wvgjxGxBbANMLvcwaxIUpp6fc890ySI992XO5GZVUiWwhERsyPi2bbO\nkTQQ2BX4deFnlkXEm5XIZ0Xq0wduvBE23zx1Xz3zTO5EZlYB1XyPYxNgMXClpMckXS5p3dZOljRe\nUqOkxsWLF1cuZU+33npwxx3Qt28apuvP3qzbK1vhkHS3pFktbGOKfIsaYDvgkogYDbwH/Ly1kyNi\nWkTUR0R9XV1dJ/wJrGgjRsDtt8PChfD1r390WnYz61bKNuVIROzRwbdYACyIiIcL+zfRRuGwzHbY\nAa65BvbfH77//bQglJQ7lZmVQdV2VUXE34H5kj5TOLQ74Lkuqtk3vwlnngnXXQeTJuVOY2Zlkms4\n7lhJC4AdgTsk3Vk4PkTSjCanHglMl/Q3YFvgjMqntZKccAL84AcwZQpcfXXuNGZWBoqI3Bk6XX19\nfTQ2NuaO0XMtXw577QUPPgh33QVf/GLuRGa2FpJmRkSrz9U1VbVdVdaF9ekDN90Em24KY8fCc8/l\nTmRmnciFw8pjgw3SMN2amjRM9/XXcycys07iwmHls8kmcNttMH9+anl88EHuRGbWCVw4rLx23BGu\nugoeeggOOwy64T01s57GS8da+R1wADz/PJxyCmy2mYfqmnVxLhxWGSedBHPmwOTJqXg0XdPczLoU\nFw6rDAmmTYO5c9OT5cOGwRe+kDuVmbWD73FY5fTtCzffnOa2Gjs2dV+ZWZfjwmGVNWgQzChMDrDP\nPvDGG3nzmFnJXDis8jbdFH73O3j55bSOx7JluROZWQlcOCyPXXZJS87+93+nJWg9TNesy/DNccvn\nO99J9zkmTYKRI9NwXTOrei4clteECal4TJiQurAOPDB3IjNbC3dVWV4SXHYZ7LorHHoo/M//5E5k\nZmvhwmH5rbMO3HJLerZjzBh44YXcicysDS4cVh0+/vE0m+7KlWmY7j//mTuRmbXChcOqx8iRcOut\n8OKL8I1veJiuWZXKtXTs/pKekrRSUqsrTkn6j8J5syRdJ6lfJXNaBrvuCldcAffdB4cf7mG6ZlUo\nV4tjFjAOeKC1EyRtCBwF1EfEVkBv4NuViWdZHXQQTJyYnvMYNAh69UrTlEyfnjuZmZFpOG5EzAaQ\ntLZTa4D+kpYDtcDCMkezarH55tC7N7z5ZtqfOxfGj0+vPbOuWVZVe48jIl4BzgPmAYuAtyLiT3lT\nWcWcfDJ8+GHzY0uWwIkn5sljZv9StsIh6e7CvYk1tzFF/vwGwBhgY2AIsK6kg9o4f7ykRkmNixcv\n7pw/hOUzb17Lx+fPhx//GGbO9P0Ps0zKVjgiYo+I2KqF7bYi32IP4KWIWBwRy4FbgJ3auN60iKiP\niPq6urrO+CNYTsOGtXx83XXTUrT19TB6NFx0kYfumlVY1XZVkbqoPi+pVulmyO7A7MyZrFKmToXa\n2ubHamvh0kth4UL45S/TPZAjj4TBg9MN9fvuS8+BmFlZ5RqOO1bSAmBH4A5JdxaOD5E0AyAiHgZu\nAh4FnixknZYjr2XQ0JBWDBw+PE1LMnx42m9ogPXXX91d9eijcNhh8Ic/wG67pZvqZ56ZiouZlYWi\nG/YT19fXR2NjY+4YVklLl6bVBS+/PE3V3rt3egL9sMPga1+DGs/nadYWSTMjotXn6pqq5q4qs+L1\n75+6q+6/H557Do47Dv76V9hvv3S/5KSTvFStWSdx4bDuZ+TI1F01b15aabC+Hs4+Ox3fbbf0IOHS\npblTmnVZLhzWffXpk2bbvf32VESmTk0PEh50EAwZkm6sP/FE7pRmXY4Lh/UMG26YuqvmzIF774W9\n907rgGy7bWqR/OpX8NZbuVOadQkuHNaz9OoFX/5y6q5auBAuvBCWL0+jtAYPhkMOgQcf9MOFZm1w\n4bCea9Cg1F31+OPwyCPwve+lBaV23RW22ALOPRdefTUVmREjPNmiWYGH45o19d57cOONaVjvn/+c\nniHp1av5vFm1taufKTHrJjwc16y91l03dVc99BDMng0DBrQ82eLRR8PTT3/0e2Y9gJ+KMmvNFlvA\nO++0/L1//AM++9nU+thmG9huu9XbqFHQt29ls5pVkAuHWVuGDUtDeNc0eHB6NuTRR9N29dVw8cXp\ne337wlZbNS8mW2+dHlI06wZ8j8OsLdOnpwWklixZfaylexwrV8ILL6wuJI89lubSeuON9P3evWHL\nLVMRGT06fd12Wxg4sLJ/HrNWlHKPw4XDbG2mT08LS82bl1ogU6cWd2M8Iq0fsqqYrNoWLVp9zsiR\nzVsmo0fDxz/euTnMiuDC4cJh1WzRotQiWdUyefRRePnl1d8fNqx5Mdluu/TQYjEtH7N2cuFw4bCu\n5o03VheRVdtzz63+fq9eLa810to9GLMSlVI4fHPcrBoMGgS77562Vd5+O82l9dhjafhvS+bNg6FD\nUwFpbdtgg/Q8ilkncYvDrCsYMaLllsV668HYsamArNqWLWt+Tm1t24Vl6FBYZ53is/heS7fkFodZ\ndzN1asv3OC6++KOjuxYvXl1E5s9vXlSeeCJNo7KmT396dSHZaKOPFpe6utRqWXOU2dy5aR9cPHqQ\nLIVD0rnA/wGWAS8Ah0bEmy2ctxfwC6A3cHlEnFXRoGbVYtUv5bX9T79XL/jUp9K2/fYtv9f778OC\nBc0LyqoCM2sWzJjRvEAB9OuXCsr8+ennm1qyBI45JnW31da2vPXv3/mrMLrlk02WripJXwHujYgV\nks4GiIgT1jinN/AcsCewAHgEODAinl7b+7uryqwDItLN+qaFZdV2ww3tf9++fZsXktaKTDHbQw+l\nSSibFrH+/eHSS+G73+34Z1CsailenZCjS42qkjQW+GZENKxxfEdgckR8tbB/IkBEnLm293ThMCuT\n1u61DB4Mt96aWh+dubVHTU0qIm1t/fqt/Zy1bXfckVpaTVeTzDFEutiHVNeiq93j+D7w2xaObwjM\nb7K/ANihIonMrGWt3Ws591zYoZP/eUakFkVLBeULX2h9zZTjjku/zFvb3nij5eMrVnQ885IlqcXz\nox+le0KV2F566aPZlyxJLZAyFbCyFQ5JdwOfbuFbJ0fEbYVzTgZWAC0tcNDS+MFWm0eSxgPjAYYN\nG1ZyXjMrQrH3WjqDtPp/92s+Td/a8yvDh8MZZ7TveitWtF1w1twOP7zl94lIhSOiMtucOS3nmDev\nfZ9DEbJ1VUk6GDgc2D0iPtImdVeVmbWqk7pnOqS1brvhw5vPBNBFclT9ehyF0VInAPu1VDQKHgFG\nStpYUl/g28DtlcpoZlWsoSEVieHDU8tk+PDK31uYOjUVq6Zqa9PxSsqQI9dCThcBA4C7JD0u6VcA\nkoZImgEQESuAI4A7gdnADRHxVKa8ZlZtGhrS/6hXrkxfKz2aqRqKV6Yc2UdVlYO7qszMSlP1XVVm\nZtZ1uXCYmVlJXDjMzKwkLhxmZlYSFw4zMytJtxxVJWkx0N5l0T4BvN6JcboyfxbN+fNozp/Hat3h\nsxgeEXXFnNgtC0dHSGosdkhad+fPojl/Hs3581itp30W7qoyM7OSuHCYmVlJXDg+alruAFXEn0Vz\n/jya8+exWo/6LHyPw8zMSuIWh5mZlcSFo0DSXpKelfS8pJ/nzpOTpI0k3SdptqSnJB2dO1NuknpL\nekzSH3JnyU3S+pJukvRM4e/Ijrkz5STpPwr/TmZJuk5Sv9yZys2Fg/RLAbgY+BowCjhQ0qi8qbJa\nAfwsIrYEPg/8pId/HgBHk6b3N/gF8MeI2ALYhh78uUjaEDgKqI+IrYDepLWDujUXjuRzwPMR8WJE\nLAOuB8ZkzpRNRCyKiEcLr98h/WLYMG+qfCQNBfYBLs+dJTdJA4FdgV8DRMSyiHgzb6rsaoD+kmqA\nWmBh5jxl58KRbAjMb7K/gB78i7IpSSOA0cDDeZNkdQFwPLAyd5AqsAmwGLiy0HV3uaR1c4fKJSJe\nAc4D5gGLgLci4k95U5WfC0eiFo71+OFmkj4G3AwcExFv586Tg6R9gdciYmbuLFWiBtgOuCQiRgPv\nAT32nqCkDUi9ExsDQ4B1JR2UN1X5uXAkC4CNmuwPpQc0N9siqQ+paEyPiFty58loZ2A/SS+TujB3\nk3Rt3khZLQAWRMSqFuhNpELSU+0BvBQRiyNiOXALsFPmTGXnwpE8AoyUtLGkvqSbW7dnzpSNJJH6\nsGdHxPm58+QUESdGxNCIGEH6e3FvRHT7/1G2JiL+DsyX9JnCod2BpzNGym0e8HlJtYV/N7vTAwYL\n1OQOUA0iYoWkI4A7SaMiroiIpzLHymln4LvAk5IeLxw7KSJmZMxk1eNIYHrhP1kvAodmzpNNRDws\n6SbgUdJoxMfoAU+R+8lxMzN8AZWQAAACDUlEQVQribuqzMysJC4cZmZWEhcOMzMriQuHmZmVxIXD\nzMxK4sJhZmYlceEwa4Okj0t6vLD9XdIrTfb/p0zXHC2p1QkVJdVJ+mM5rm1WDD8AaNaGiPgHsC2A\npMnAuxFxXpkvexJwehuZFktaJGnniPhzmbOYfYRbHGbtJOndwtcvSfpvSTdIek7SWZIaJP1V0pOS\nNi2cVyfpZkmPFLadW3jPAcDWEfFEYf+LTVo4jxW+D/A7oKFCf1SzZlw4zDrHNqTFnv6NNF3L5hHx\nOdIaHkcWzvkF8F8RsT3wDVpe36MemNVk/1jgJxGxLfAFYGnheGNh36zi3FVl1jkeiYhFAJJeAFat\nyfAk8OXC6z2AUWkuPAAGShpQWCxrlcGk9S5W+TNwvqTpwC0RsaBw/DXSNN5mFefCYdY5PmjyemWT\n/ZWs/nfWC9gxIpbSuqXAv9asjoizJN0B7A38r6Q9IuKZwjltvY9Z2biryqxy/gQcsWpH0rYtnDMb\n2KzJOZtGxJMRcTape2qLwrc2p3mXllnFuHCYVc5RQL2kv0l6Gjh8zRMKrYn1mtwEP0bSLElPkFoY\n/69w/MvAHZUIbbYmT6tuVmUk/QfwTkS09SzHA8CYiPhn5ZKZJW5xmFWfS2h+z6QZSXXA+S4alotb\nHGZmVhK3OMzMrCQuHGZmVhIXDjMzK4kLh5mZlcSFw8zMSvL/AQWairriLUnnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "plt.plot(a, '-or')\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.ylabel(\"Response\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inputs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a lined out process with a single step in $u$ at time 0:\n",
    "\n",
    "$$ \\Delta u_0 = u_0 - u_{-1} $$\n",
    "\n",
    "then\n",
    "\n",
    "$$ y_1 - y_0 = a_1 \\Delta u_0 $$\n",
    "$$ y_2 - y_0 = a_2 \\Delta u_0 $$\n",
    "$$ y_3 - y_0 = a_3 \\Delta u_0 $$\n",
    "$$ y_4 - y_0 = a_4 \\Delta u_0 $$\n",
    "$$ \\vdots $$\n",
    "$$ y_t - y_0 = a_t \\Delta u_0 $$\n",
    "\n",
    "This is essentially the definition of the step response convolution model described before. What is interesting about this is that we can calculate what y will be at each interval for successive moves in u.\n",
    "\n",
    "For example purposes only, let us examine a series of three moves \"into the future\": one now, one at the next time interval, and one two time intervals into the future: \n",
    "\n",
    "$$ \\Delta u_0 = u_0 - u_{-1} $$\n",
    "$$ \\Delta u_1 = u_1 - u_{0} $$\n",
    "$$ \\Delta u_2 = u_2 - u_{1} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0]\n",
      " [ 1]\n",
      " [ 0]\n",
      " [ 0]\n",
      " [-1]]\n"
     ]
    }
   ],
   "source": [
    "U = np.atleast_2d(np.array([0,1,0,0,-1])).T\n",
    "print(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10e2d79e8>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEMVJREFUeJzt3X+Q3HV9x/HnCyK1/kjFSexgEk2Y\nRpsUadEDsUglNUKITugfjgNTrW2pzDCittG2MFbs0JlORYxtR4oy1uK0FhqprRmNjSLHOHaE5vAH\nCik1RiXXUDktlU6Ziozv/rEbOY5Ldu+yd3v76fMxc3P7+eyH7/v9/d7dK9/97i6bqkKS1Jbjht2A\nJGnwDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg5YNq/CKFStq7dq1wyovSSPp\nzjvv/G5Vrey1bmjhvnbtWiYmJoZVXpJGUpJv97POyzKS1CDDXZIaZLhLUoMMd0lqkOEuSQ3qGe5J\nPpTkgSRfO8L9SfLnSfYnuSvJCwff5uK5+moYH3/83Ph4Z35Qtm6FHTseP7djR2d+lGrA4hwvSXPX\nz5n7DcCWo9x/PrC++3UJcN2xtzU8p58Or3nNY4E1Pt4Zn3764Gps3gxve9tj4btjR2e8efNo1YDF\nOV6S5i79fMxekrXAJ6rqlFnu+wBwW1Xd2B3fC5xTVfcfbZtjY2O1VF/nPj4O27bBM58J998PGzbA\niScOtsbBg3DgACxfDg89BCefDGvWjF4NgAcfhH374Nxz4Y47YOdO2LRp8HUkQZI7q2qs17pBXHNf\nBRycNp7szs3W1CVJJpJMTE1NDaD0wti0qRPs990HJ500+GCHTsgeDt3lyxcmdBejBnSOz3HHwSc/\nCZdearBLS8Eg3qGaWeZmfThQVdcD10PnzH0AtRfE+HjnjP05z4GHH4Z3vnPwgbVjB3zuc3D22fD5\nz3ceKWzfPno1oHO8zjuvc7yuu65zrAx4abgGceY+CUw/J1wNHBrAdofi8DXjDRtg3brOJYbp15QH\n4fD172uu6YTvNdc8/vr4qNSAxTlekuZuEOG+C/i17qtmzgS+3+t6+1K2d28noA5fitm0qTPeu3dw\nNW65pRO2h8+it2/vjG+5ZbRqwOIcL0lz1/MJ1SQ3AucAK4DvAO8EngRQVe9PEuB9dF5R8zDwG1XV\n85nSpfyEKsA553S+33bbMLsYHR4vaXH0+4Rqz2vuVXVRj/sLeOMcepMkLTDfoSpJDTLcJalBhrsk\nNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KD\nDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchw\nl6QGGe6S1CDDXZIa1Fe4J9mS5N4k+5NcPsv9z0kynuRLSe5KsnXwrUqS+tUz3JMcD1wLnA9sBC5K\nsnHGsj8AdlbVacCFwF8MulFJUv/6OXM/A9hfVQeq6hHgJuCCGWsKWN69/VPAocG1KEmaq2V9rFkF\nHJw2ngRePGPNHwKfTvIm4KnA5oF0J0mal37O3DPLXM0YXwTcUFWrga3AXyd5wraTXJJkIsnE1NTU\n3LuVJPWln3CfBNZMG6/miZddLgZ2AlTVF4AnAytmbqiqrq+qsaoaW7ly5fw6liT11E+47wXWJ1mX\n5AQ6T5jumrHmPuDlAEk20Al3T80laUh6hntVPQpcBuwB9tF5VczdSa5Ksq277K3AG5J8BbgR+PWq\nmnnpRpK0SPp5QpWq2g3snjF35bTb9wBnDbY1SdJ8+Q5VSWqQ4S5JDTLcJalBhrskNchwl6QGGe6S\n1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkN\nMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDD\nXZIa1Fe4J9mS5N4k+5NcfoQ1r0lyT5K7k/ztYNuUJM3Fsl4LkhwPXAu8ApgE9ibZVVX3TFuzHrgC\nOKuqHkzyrIVqWJLUWz9n7mcA+6vqQFU9AtwEXDBjzRuAa6vqQYCqemCwbUqS5qKfcF8FHJw2nuzO\nTfc84HlJ/jnJ7Um2zLahJJckmUgyMTU1Nb+OJUk99RPumWWuZoyXAeuBc4CLgA8mecYT/qOq66tq\nrKrGVq5cOddeJUl96ifcJ4E108argUOzrPl4Vf2wqr4J3Esn7CVJQ9BPuO8F1idZl+QE4EJg14w1\n/whsAkiygs5lmgODbFSS1L+e4V5VjwKXAXuAfcDOqro7yVVJtnWX7QG+l+QeYBz43ar63kI1LUk6\nup4vhQSoqt3A7hlzV067XcD27pckach8h6okNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y\n3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNd\nkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqUF/hnmRLknuT\n7E9y+VHWvTpJJRkbXIuSpLnqGe5JjgeuBc4HNgIXJdk4y7qnA28G7hh0k5KkuennzP0MYH9VHaiq\nR4CbgAtmWfdHwNXA/w6wP0nSPPQT7quAg9PGk925H0tyGrCmqj4xwN4kSfPUT7hnlrn68Z3JccB7\ngbf23FBySZKJJBNTU1P9dylJmpN+wn0SWDNtvBo4NG38dOAU4LYk3wLOBHbN9qRqVV1fVWNVNbZy\n5cr5dy1JOqp+wn0vsD7JuiQnABcCuw7fWVXfr6oVVbW2qtYCtwPbqmpiQTqWJPXUM9yr6lHgMmAP\nsA/YWVV3J7kqybaFblCSNHfL+llUVbuB3TPmrjzC2nOOvS1J0rHwHaqS1CDDXZIaZLhLUoMMd0lq\nkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ\n7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEu\nSQ0y3CWpQYa7JDWor3BPsiXJvUn2J7l8lvu3J7knyV1JPpvkuYNvVZLUr57hnuR44FrgfGAjcFGS\njTOWfQkYq6pTgZuBqwfdqCSpf/2cuZ8B7K+qA1X1CHATcMH0BVU1XlUPd4e3A6sH26YkaS76CfdV\nwMFp48nu3JFcDHxqtjuSXJJkIsnE1NRU/11Kkuakn3DPLHM168LktcAY8O7Z7q+q66tqrKrGVq5c\n2X+XkqQ5WdbHmklgzbTxauDQzEVJNgNvB15WVT8YTHuSpPno58x9L7A+ybokJwAXArumL0hyGvAB\nYFtVPTD4NiVJc9Ez3KvqUeAyYA+wD9hZVXcnuSrJtu6ydwNPAz6a5MtJdh1hc5KkRdDPZRmqajew\ne8bcldNubx5wX5KkY+A7VCWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGG\nuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhL\nUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KD+gr3JFuS3Jtkf5LLZ7n/J5L8\nXff+O5KsHXSjkqT+9Qz3JMcD1wLnAxuBi5JsnLHsYuDBqvoZ4L3Auwbd6NVXw/j44+fGxzvzo1in\nFYtxvLZuhR07Hj+3Y0dnfpRqLFYd92Xp1VjMOj9WVUf9Al4C7Jk2vgK4YsaaPcBLureXAd8FcrTt\nvuhFL6q5uPXWqhUrOt9nGw/K4e2eemrVy162cHVasRjH6z3vqUo632cbj0qNxarjviy9GoOsA0xU\nj9yuqk4AH02SVwNbquq3uuPXAS+uqsumrflad81kd/yN7prvHmm7Y2NjNTExMad/iMbH4ZWvhJUr\n4f77YcMGOPHEOW2iLw8+CF/9KqxZAw8/DDt3wqZNg6/TivFxOO88OO44+NGPFubncvAgHDgAy5fD\nQw/BySd3fj6jVmOx6rgvS6/G9DrPfnYnw665BrZvn9s2ktxZVWO91vVzzT2zzM38F6GfNSS5JMlE\nkompqak+Sj/epk3wghfAfffBSSctTLBDZ7tr1nTqXHqpwd7Lpk1w7rnwgx8s3M9lzZrH/vCWL1+Y\nP7zFqLFYddyXpVdjep1Dh+ClL517sM9Jr1N7lshlmarHHvK/4x0Le6lkseq0YjGO1+GHsGefvTAP\nmRerxmLVcV+WXo1B1aHPyzL9hPsy4ACwDjgB+ArwczPWvBF4f/f2hcDOXttd6tfcF7pOKxbjeHnd\ndenVWKw6rdQYZJ2BhXtnW2wF/g34BvD27txVwLbu7ScDHwX2A/8CnNxrm3MN93e964mBceutnflB\nWqw6rViM43X++U/8A3jPezrzo1Rjseq4L0uvxiDr9BvuPZ9QXSjzeUJVkv6/G+QTqpKkEWO4S1KD\nDHdJapDhLkkNMtwlqUFDe7VMking2/P8z1fQeaNUC9yXpaeV/QD3Zak6ln15blWt7LVoaOF+LJJM\n9PNSoFHgviw9rewHuC9L1WLsi5dlJKlBhrskNWhUw/36YTcwQO7L0tPKfoD7slQt+L6M5DV3SdLR\njeqZuyTpKEYu3Ht9WPeoSLImyXiSfUnuTvKWYfd0LJIcn+RLST4x7F6ORZJnJLk5yb92fzYvGXZP\n85Xkd7q/W19LcmOSJw+7p34l+VCSB7qf8nZ47plJPpPk693vC/RxPYNzhP14d/f3664k/5DkGQtR\ne6TCvc8P6x4VjwJvraoNwJnAG0d4XwDeAuwbdhMD8GfAP1XVzwI/z4juU5JVwJuBsao6BTiezmct\njIobgC0z5i4HPltV64HPdsdL3Q08cT8+A5xSVafS+V+pX7EQhUcq3IEzgP1VdaCqHgFuAi4Yck/z\nUlX3V9UXu7f/m06IrBpuV/OTZDXwSuCDw+7lWCRZDvwS8JcAVfVIVf3XcLs6JsuAn0yyDHgKcGjI\n/fStqj4H/OeM6QuAD3dvfxj4lUVtah5m24+q+nRVPdod3g6sXojaoxbuq4CD08aTjGggTpdkLXAa\ncMdwO5m3PwV+D/jRsBs5RicDU8BfdS8xfTDJU4fd1HxU1b8D1wD3AfcD36+qTw+3q2P201V1P3RO\njoBnDbmfQfhN4FMLseFRC/e+Poh7lCR5GvD3wG9X1UPD7meukrwKeKCq7hx2LwOwDHghcF1VnQb8\nD6Px0P8JutejL6Dz8ZjPBp6a5LXD7UrTJXk7ncuzH1mI7Y9auE8C0z+XfDUj9FBzpiRPohPsH6mq\njw27n3k6C9iW5Ft0LpP9cpK/GW5L8zYJTFbV4UdQN9MJ+1G0GfhmVU1V1Q+BjwG/OOSejtV3kpwE\n0P3+wJD7mbckrwdeBfxqLdDr0Uct3PcC65OsS3ICnSeIdg25p3lJEjrXdvdV1Y5h9zNfVXVFVa2u\nqrV0fh63VtVIniFW1X8AB5M8vzv1cuCeIbZ0LO4DzkzylO7v2ssZ0SeHp9kFvL57+/XAx4fYy7wl\n2QL8Pp3PoH54oeqMVLh3n4S4DNhD5xd1Z1XdPdyu5u0s4HV0znS/3P3aOuymxJuAjyS5C/gF4I+H\n3M+8dB993Ax8Efgqnb/1kXmHZ5IbgS8Az08ymeRi4E+AVyT5OvCK7nhJO8J+vA94OvCZ7t/9+xek\ntu9QlaT2jNSZuySpP4a7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkN+j/osaWGRCJ2SgAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "I_obs = np.pad(np.cumsum(U),(0,8),'constant')\n",
    "plt.step(I_obs, '-bx', where=\"post\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "then\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "y_1 - y_0 & = a_1 \\Delta u_0 \\\\\n",
    "y_2 - y_0 & = a_2 \\Delta u_0 + a_1 \\Delta u_1 \\\\\n",
    "y_3 - y_0 & = a_3 \\Delta u_0 + a_2 \\Delta u_1 + + a_1 \\Delta u_2 \\\\\n",
    "y_4 - y_0 & = a_4 \\Delta u_0 + a_3 \\Delta u_1 + + a_2 \\Delta u_2 \\\\\n",
    "\\vdots & \\\\\n",
    "y_t - y_0 & = a_t \\Delta u_0 + a_{t-1} \\Delta u_1 + + a_{t-2} \\Delta u_2 \\\\\n",
    "\\end{align}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.5    0.     0.     0.     0.   ]\n",
      " [-1.    -0.5    0.     0.     0.   ]\n",
      " [-1.5   -1.    -0.5    0.     0.   ]\n",
      " [-1.8   -1.5   -1.    -0.5    0.   ]\n",
      " [-1.9   -1.8   -1.5   -1.    -0.5  ]\n",
      " [-1.95  -1.9   -1.8   -1.5   -1.   ]\n",
      " [-1.98  -1.95  -1.9   -1.8   -1.5  ]\n",
      " [-1.99  -1.98  -1.95  -1.9   -1.8  ]\n",
      " [-1.999 -1.99  -1.98  -1.95  -1.9  ]\n",
      " [-2.    -1.999 -1.99  -1.98  -1.95 ]]\n"
     ]
    }
   ],
   "source": [
    "# Build dynamic matrix\n",
    "t = a.shape[0]\n",
    "u = U.shape[0]\n",
    "\n",
    "A = np.tile(a, (u,1)).T\n",
    "row_indices = np.atleast_2d(np.arange(t)).T - np.arange(u)\n",
    "col_indices = np.arange(u)\n",
    "A = np.tril(A[row_indices,col_indices])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Y = np.multiply(A,U.T)\n",
    "Y_obs = np.matmul(A,U)\n",
    "U_obs = np.pad(np.cumsum(U),(0,a.shape[0]-U.shape[0]),'constant')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.   ,  0.   ,  0.   ,  0.   , -0.   ],\n",
       "       [-0.   , -0.5  ,  0.   ,  0.   , -0.   ],\n",
       "       [-0.   , -1.   , -0.   ,  0.   , -0.   ],\n",
       "       [-0.   , -1.5  , -0.   , -0.   , -0.   ],\n",
       "       [-0.   , -1.8  , -0.   , -0.   ,  0.5  ],\n",
       "       [-0.   , -1.9  , -0.   , -0.   ,  1.   ],\n",
       "       [-0.   , -1.95 , -0.   , -0.   ,  1.5  ],\n",
       "       [-0.   , -1.98 , -0.   , -0.   ,  1.8  ],\n",
       "       [-0.   , -1.99 , -0.   , -0.   ,  1.9  ],\n",
       "       [-0.   , -1.999, -0.   , -0.   ,  1.95 ]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([1, 2, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9]),\n",
       " array([1, 1, 1, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nonzero(Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 4])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nonzero(np.any(Y != 0, axis=0))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAFpCAYAAACmt+D8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VNX9//HXJxthC0lI2BPCDhEQ\nMKIigigim4BVrIq4dKG1Umv1W+vS/dvW2m/la/vFLrSuP3EBVESk4oZ1QZQAyioYQCCAsq8Bsp3f\nH2dCEggQZJK5Sd7Px2MeM3PvnTufTCDvOefee4455xAREZHgiIp0ASIiIlKewllERCRgFM4iIiIB\no3AWEREJGIWziIhIwCicRUREAkbhLCIiEjAKZxERkYBROIuIiASMwllERCRgYiL1xikpKS4jIyNS\nby8iIlKtFi1atMM5l1qZbSMWzhkZGWRnZ0fq7UVERKqVmW2o7Lan7NY2s8fMbJuZLT/BejOzv5hZ\njpktNbM+p1OsiIiIlFeZY85PAENPsn4Y0Cl0mwD87czLqln++EeYN6/8snnz/PLqNHw4TJpUftmk\nSX55XawjKL8XEZHTdcpwds69C+w6ySajgaectwBINLOW4SqwJjj3XLjmmtIgmDfPPz/33OqtY/Bg\n+K//Kg3GSZP888GD62YdQfm9iIicLqvMfM5mlgHMds51r2DdbOAPzrn3Q8/fAn7qnDvpAeWsrCxX\nm445z5sHo0ZBcjJs3QrdukFSUvXXsWkTrFsHCQmwbx+0bw9paXW3jt27YdUqGDIEPvoIpk2DQYOq\nvw4RETNb5JzLqsy24biUyipYVmHim9kEM8s2s+zt27eH4a2DY9AgH8wbN0LLlpEJZvABWBKICQmR\nCcQg1ZGUBFFR8OqrcOutCmYRqRnCcbZ2LlD2T28bYEtFGzrnpgBTwLecw/DegTFvnm8xp6dDXh78\n8peRCYJJk+Ddd+Gii+D9931r/s47624d8+bB5Zf738vf/uZ/JwpoEQm6cLScZwE3hs7aPh/Y65zb\nGob91hglxzK7dYN27XzXadljndWl5Njun/7kg/FPfyp/7Leu1RGU34uIyOmqzKVUzwIfAl3MLNfM\nvm1m3zez74c2mQOsA3KAfwI/qLJqA2rhQv+Hv6Qre9Ag/3zhwuqt4803fRCWtFDvvNM/f/PNullH\nUH4vIiKnq1InhFWF2nZCGMDFF/v7d96JZBVyLP1eRCQIqvuEMBEREQkjhbOIiEjAKJxFREQCRuEs\nIiISMApnERGRgFE4i4iIBIzCWUREJGAUziIiIgGjcBYREQkYhbOIiEjAKJxFREQCRuEsIiISMApn\nERGRgFE4i4iIBIzCWUREJGAUziIiIgGjcBYREQkYhbOIiEjAVCqczWyoma02sxwzu6eC9elmNs/M\nlpjZUjMbHv5SRURE6oZThrOZRQOPAMOATOA6M8s8ZrOfAdOcc72Ba4G/hrtQERGRuqIyLee+QI5z\nbp1zLh94Dhh9zDYOSAg9bgJsCV+JIiIidUtMJbZpDWwq8zwXOO+YbX4FvG5mPwQaAoPDUp2IiEgd\nVJmWs1WwzB3z/DrgCedcG2A48P/M7Lh9m9kEM8s2s+zt27effrUiIiJ1QGXCORdIK/O8Dcd3W38b\nmAbgnPsQiAdSjt2Rc26Kcy7LOZeVmpr69SoWERGp5SoTzguBTmbWzszi8Cd8zTpmm43ApQBm1g0f\nzmoai4iIfA2nDGfnXCEwEZgLrMKflb3CzH5jZqNCm90FfNfMPgWeBW52zh3b9S0iIiKVUJkTwnDO\nzQHmHLPsF2UerwQuDG9pIiIidZNGCBMREQkYhbOIiEjAKJxFREQCRuEsIiISMApnERGRgFE4i4iI\nBIzCWUREJGAUziIiIgGjcBYREQkYhbOIiEjAKJxFREQCRuEsIiISMApnERGRgFE4i4iIBIzCWURE\nJGAUziIiIgGjcBYREQkYhbOIiEjAVCqczWyoma02sxwzu+cE21xjZivNbIWZPRPeMkVEROqOmFNt\nYGbRwCPAZUAusNDMZjnnVpbZphNwL3Chc263mTWrqoJFRERqu8q0nPsCOc65dc65fOA5YPQx23wX\neMQ5txvAObctvGWKiIjUHZUJ59bApjLPc0PLyuoMdDazD8xsgZkNDVeBIiIidc0pu7UBq2CZq2A/\nnYCLgTbAe2bW3Tm3p9yOzCYAEwDS09NPu1gREZG6oDIt51wgrczzNsCWCrZ52TlX4JxbD6zGh3U5\nzrkpzrks51xWamrq161ZRESkVqtMOC8EOplZOzOLA64FZh2zzUxgEICZpeC7udeFs1AREZG64pTh\n7JwrBCYCc4FVwDTn3Aoz+42ZjQptNhfYaWYrgXnAT5xzO6uqaBERkdqsMseccc7NAeYcs+wXZR47\n4M7QTURERM6ARggTEREJGIWziIhIwCicRUREAkbhLCIiEjAKZxERkYBROIuIiASMwllERCRgFM4i\nIiIBo3AWEREJGIWziIhIwCicRUREAkbhLCIiEjAKZxERkYBROIuIiASMwllERCRgFM4iIiIBo3AW\nEREJGIWziIhIwFQqnM1sqJmtNrMcM7vnJNtdbWbOzLLCV6KIiEjdcspwNrNo4BFgGJAJXGdmmRVs\n1xi4Hfgo3EWKiIjUJZVpOfcFcpxz65xz+cBzwOgKtvtv4I/A4TDWJyIiUudUJpxbA5vKPM8NLTvK\nzHoDac652WGsTUREpE6qTDhbBcvc0ZVmUcD/AnedckdmE8ws28yyt2/fXvkqRURE6pDKhHMukFbm\neRtgS5nnjYHuwDtm9gVwPjCropPCnHNTnHNZzrms1NTUr1+1iIhILVaZcF4IdDKzdmYWB1wLzCpZ\n6Zzb65xLcc5lOOcygAXAKOdcdpVULCIiUsudMpydc4XARGAusAqY5pxbYWa/MbNRVV2giIhIXRNT\nmY2cc3OAOccs+8UJtr34zMsSERGpuzRCmIiISMAonEVERAJG4SwiIhIwCmcREZGAUTiLiIgEjMJZ\nREQkYBTOIiIiAaNwFhERCRiFs4iISMAonEVERAJG4SwiIhIwCmcREZGAUTiLiIgEjMJZREQkYBTO\nIiIiAaNwFhERCRiFs4iISMAonEVERAKmUuFsZkPNbLWZ5ZjZPRWsv9PMVprZUjN7y8zahr9UERGR\nuuGU4Wxm0cAjwDAgE7jOzDKP2WwJkOWc6wnMAP4Y7kJFRETqisq0nPsCOc65dc65fOA5YHTZDZxz\n85xzeaGnC4A24S1TRESk7qhMOLcGNpV5nhtadiLfBv59JkWJiIjUZTGV2MYqWOYq3NDsBiALGHiC\n9ROACQDp6emVLFFERKRuqUzLORdIK/O8DbDl2I3MbDBwPzDKOXekoh0556Y457Kcc1mpqalfp14R\nEZFarzLhvBDoZGbtzCwOuBaYVXYDM+sN/AMfzNvCX6aIiEjdccpwds4VAhOBucAqYJpzboWZ/cbM\nRoU2+x+gETDdzD4xs1kn2J2IiIicQmWOOeOcmwPMOWbZL8o8HhzmukREROosjRAmIiISMApnERGR\ngFE4i4iIBIzCWUREJGAUziIiIgGjcBYREQkYhbOIiEjAKJxFREQCRuEsIiISMApnERGRgFE4i4iI\nBIzCWUREJGAUziIiIgGjcBYREQkYhbOIiEjAKJxFREQCRuEsIiISMApnERGRgKlUOJvZUDNbbWY5\nZnZPBevrmdnzofUfmVlGuAsVERGpK04ZzmYWDTwCDAMygevMLPOYzb4N7HbOdQT+F3gw3IVW5I9/\nhHnzyi+bN88vr05BqUPKC8LvZfhwmDSp/LJJk/zy6qQ6VEfQ6whCDUGqA+fcSW/ABcDcMs/vBe49\nZpu5wAWhxzHADsBOtt9zzjnHnam333YuJcXfV/S8upS8b8+ezg0cGLk6pLwg/F4eesg5M39f0XPV\noTpUR3BqqOo6gGx3iswtuZnf/sTM7GpgqHPuO6Hn44HznHMTy2yzPLRNbuj52tA2O06036ysLJed\nnf21vlCUNW8ejBgBqamwdSt06wZJSWe829O2ezcsWwZpaZCXB9OmwaBB1V+HlDdvHlx+OURFQXFx\nZP59bNoE69ZBQgLs2wft2/t/J9VNdaiOoNcRhBrK1tGqlc+VP/0J7rzzzPdrZoucc1mV2bYyx5yt\ngmXHJnpltsHMJphZtpllb9++vTL1ndKgQdCjB2zcCC1bRiaYwb9vWpqv49ZbFcxBMWgQDBkCR45E\n7t9HWlrpH5uEhMj8sVEdqqMm1BGEGsrWsWUL9O8fnmA+badqWhPgbm3nSrsqf/7zyHYlB6UOKS8I\nv5eSbrGLLopMN53qUB01pY4g1FCVdXAa3dqVCecYYB3QDogDPgXOOmab24C/hx5fC0w71X5r4zHn\nSNch5QXh91IXjqOpDtVRW2qo6jrCGs5+fwwH1gBrgftDy34DjAo9jgemAznAx0D7U+0zHOH84IPH\n/6F9+22/vDoFpQ4pLwi/l2HDjv9P/dBDfnl1Uh2qI+h1BKGGqq7jdML5lCeEVZVwnRAmIiJSE4T7\nhDARERGpRgpnERGRgIlYt7aZbQc2hHGXKfizxMXT51GePo9S+izK0+dRnj6PUuH+LNo651Irs2HE\nwjnczCy7sn35dYE+j/L0eZTSZ1GePo/y9HmUiuRnoW5tERGRgFE4i4iIBExtCucpkS4gYPR5lKfP\no5Q+i/L0eZSnz6NUxD6LWnPMWUREpLaoTS1nERGRWqFWhLOZDTWz1WaWY2b3RLqeSDKzNDObZ2ar\nzGyFmf0o0jVFmplFm9kSM5sd6VoizcwSzWyGmX0W+jdyQaRrihQz+3Ho/8hyM3vWzOIjXVN1MrPH\nzGxbaMrfkmXJZvaGmX0euo/QPH/V7wSfx/+E/q8sNbOXzCyxuuqp8eFsZtHAI8AwIBO4zswyI1tV\nRBUCdznnugHnA7fV8c8D4EfAqkgXERB/Bl5zznUFzqaOfi5m1hq4HchyznUHovGT9tQlTwBDj1l2\nD/CWc64T8FboeV3xBMd/Hm8A3Z1zPfHzS9xbXcXU+HAG+gI5zrl1zrl84DlgdIRrihjn3Fbn3OLQ\n4/34P76tI1tV5JhZG2AE8K9I1xJpZpYADAAeBXDO5Tvn9kS2qoiKAeqbWQzQANgS4XqqlXPuXWDX\nMYtHA0+GHj8JjKnWoiKoos/DOfe6c64w9HQB0Ka66qkN4dwa2FTmeS51OIzKMrMMoDfwUWQriaiH\ngbuB4kgXEgDtge3A46Fu/n+ZWcNIFxUJzrnNwJ+AjcBWYK9z7vXIVhUIzZ1zW8F/0QeaRbieIPkW\n8O/qerPaEM5WwbI6fwq6mTUCXgDucM7ti3Q9kWBmI4FtzrlFka4lIGKAPsDfnHO9gYPUrW7Lo0LH\nUkfj56lvBTQ0sxsiW5UElZndjz9kOLW63rM2hHMukFbmeRvqWPfUscwsFh/MU51zL0a6ngi6EBhl\nZl/gD3dcYmZPR7akiMoFcp1zJT0pM/BhXRcNBtY757Y75wqAF4F+Ea4pCL4ys5YAofttEa4n4szs\nJmAkMM5V47XHtSGcFwKdzKydmcXhT+qYFeGaIsbMDH9McZVzblKk64kk59y9zrk2zrkM/L+Lt51z\ndbZ15Jz7EthkZl1Ciy4FVkawpEjaCJxvZg1C/2cupY6eHHeMWcBNocc3AS9HsJaIM7OhwE+BUc65\nvOp87xofzqGD9ROBufj/XNOccysiW1VEXQiMx7cSPwndhke6KAmMHwJTzWwp0Av4fYTriYhQ78EM\nYDGwDP+3sE6NjGVmzwIfAl3MLNfMvg38AbjMzD4HLgs9rxNO8HlMBhoDb4T+lv692urRCGEiIiLB\nUuNbziIiIrWNwllERCRgFM4iIiIBo3AWEREJGIWziIhIwCicRUREAkbhLCIiEjAKZxERkYBROIuI\niASMwllERCRgFM4iIiIBo3AWEREJGIWziIhIwCicRUREAiYmUm+ckpLiMjIyIvX2IiIi1WrRokU7\nnHOpldk2LOFsZmnAU0ALoBiY4pz788lek5GRQXZ2djjeXkREJPDMbENltw1Xy7kQuMs5t9jMGgOL\nzOwN59zKMO1fRESkzgjLMWfn3Fbn3OLQ4/3AKqB1OPYtIiJSbd5/GNa/W37Z+nf98moU9hPCzCwD\n6A18FO59i4iIVKnWfWD6zaUBvf5d/7x1n2otI6wnhJlZI+AF4A7n3L4K1k8AJgCkp6eH861FRCRo\n3n/Yh1q7AaXL1r8LmxdD/zsiU9PhvZCfBwV5UHAICg9DXENo1s2vP7ANelwDz3wTsr4Nnz4DY58o\n/zNUg7CFs5nF4oN5qnPuxYq2cc5NAaYAZGVluXC9t4iIBFBJK7Qk3EpaoWOfKN2muBiKjkBsff98\nby7k7fTBWRKgUbHQeYhfv3Q67FobWnfY3zduAZf8zK+feRt8tTz0+tA+Wp8D46b59f8YCLvXl6+z\n8zC4/jn/eO59cOAr//jD/4MBd1d7MEP4ztY24FFglXNuUjj2KSIiNUB+HhzaBU3a+Oe52bBlCeTt\n8iHbvDs8fRVceAdkPwoZF8GL3/OhWXjY3xqmwk9y/Ovn3A2rXy3/HoltofNS//iTp2HdOxBdzwd6\nbANo0b1027gG0Kh5aF3oltKldP3An0LhIf+62PoQU9+He4nvvOlb9rN/DFnf8jW3u6jaA9qcO/MG\nrJn1B94DluEvpQK4zzk350SvycrKcrqUSkSkkqqrizhvl29Z5u32oZu30y+78EdQrxEsehI+/mfp\nusLD/nX3bfXB+Np9sOARv6xeE2iQBIVHYP9W3wpt0toHeElwxtSH+CZwwQ/8azZ+BHk7StfF1od6\njaFpB7++8IhvSUdV0RhaZVv3x7b2zzCgzWyRcy6rMtuGpeXsnHsfsHDsS0REKlCZLuISxcVweA8c\n2l0arm3OhYZNIXcRLH7SLz+0u7SFe9Mr0KwrLJsO/777mB0a9LrOh3NcQ99KbtkTGiRD/WR/b6EI\nuOgu/2WhfhJEx5bWOeBu3wod+wScc/OJf870807+OcTUq9zn9XVtXlw+iNsN8M83L67W1nNYWs5f\nh1rOIiKnaf27MO1GaNMXvngPOl8OsQ19uPa/A9LPh5y3YOrV4IrLv/aGF6HjpfDZq/DKHdCgaShc\nk/zji+6EpAzYtR62f1Yaug2a+pZtVPTXq7eKWqE1UbW3nEVEJMycg72bYNPHsOkj6DQEOl0GmaNh\n0RN+m8/mhEK2KRw54Jclt/et1/rJZQI4GVI6+fVdR/jbiSS387dwCEgrtCZSOIuIBEnBYXjpez6U\n92/xy2Ib+lZtTD1Y9Qqc9z1/1vLYJ6D9wPKvT25XeuZypFV0LLzdAAVzJSicRUQi4eBOyA21ijd9\n7M9IvvJvEBvvL+XJuBDSzoO0vtDsLNg4v3yXcNeRdbqLuLZTOIuIVLXiYtiXC4mhwZeeGwefzfaP\no2Kh5dmlZyMDfOu14/ehLuI6ReEsIhJuRw7A5kWlx4tzP4biIvjpBoiOgY6DoU0WpJ0PrXqVDsBx\nMuoirlMUziIiZ8I52LPRB3GXof6a3Pn/B//5g1+f2g3OutJ3UbsiIAayboloyRJ8CmcRkdO1byus\neLH0ePH+rX75+JegwyXQY6y/rrhNFtRPjGytUiMpnEVETubgjtLu6Q6X+LOjD3zpx2BObOuHo0zr\n61vGzTL9a1I6+pvI16RwFpG651RDYRYcgtl3+kDetdavj4qFhik+nJv3gLtWlx+TWSSMFM4iUveU\nDIU55u8QEwdLp/lbp8v8+ph4+HIppHaFc27yreKWvfxlTuBP6lIwSxVSOItI3VNylvMzY0uXJWZA\nix7+sRnc+kFEShMBqKJpPUREAmbPRvjgz/6SJvDHh1ue7R/3ux3u+BQG3Re5+kTKUDiLSO2Vn+e7\nq58cBQ/3hDd+AVs+8evSz4e9uX62pE+m+mPOIgGhbm0RqZ22rYJHh8CRff6s6ovv9dMeJqYfPztS\nu4s0FKYEisJZRGqHfVth6XMQ28BPDJHSGc6+FrqNgrYXQlSZjkINhSkBp/mcRaTmKjwCq+fAkqmw\n9i0/h3HmGLjmyUhXJnIczecsInXDq3fCkqchoTX0/zH0Gld+AgmRGkrhLCI1w8EdsPR5+OQZ+MY/\noXkm9P2eH7e6/SCIio50hSJho3AWkeAqKoTPX/dnU695DYoLoVVvf5IXQMueka1PpIoonEUkeI7s\n97M7FeTBjG9BvUZw3vd9t3XzzEhXJ1LlFM4iEgyHdsPyF/zJXa4IvvcuxCfAt+f6AUOiYyNdoUi1\nUTiLSGTlLoIFj8Cq2VB0BJqdBb3H+ZG8oqJLR/ESqUMUziJS/XauhQbJUD8Jtq2EtW/7CSZ6jfNh\nbBbpCkUiSuEsItXjyH5YMdOf3LXxQ7j8AbjgB9DzGn+LqRfpCkUCI2zhbGaPASOBbc657uHar4jU\ncMVFMOt2WPGiP8GraScY/Ct/CRQolEUqEM6JL54AhoZxfyJSE7z/8PGTRiydDjO+7R9HRcPhPdBj\nLHz7DZi40A8YktCy+msVqSHC1nJ2zr1rZhnh2p+I1BCt+/hJI8b8Aw7tgg8nw5dLwaJhxENQPxGu\nnRrpKkVqFB1zFpEz024AXHAbPDMWcGBRcPb1cPE9PphF5LRV63zOZjbBzLLNLHv79u3V+dYiEk77\ntvru7LVv++c9rvHXIgP0vwuu/BsktY1cfSI1XLWGs3NuinMuyzmXlZqaWp1vLSJnqvCIP9t66lj4\n30x485ewdp5ft3s9HPgSBtwNix47/hi0iJwWdWuLSOU8MRJyP4bGreDCO/w1ySkdfRBPv7l0fuR2\nF5V/LiKnLZyXUj0LXAykmFku8Evn3KPh2r+IVKODO2DpNPhsNtzwIsTGQ/87ILoedDhmBqjNi8sH\ncbsB/vnmxQpnka/JnHMReeOsrCyXnZ0dkfcWkQoUFULOG35+5DVzobjAzwB19WOQ3D7S1YnUeGa2\nyDmXVZlt1a0tUtcVFfhJJbZ+Cs9eCw1T4bzvaQYokQhSOIvURYf2wPIZ8Mkz0KIHXPFnf73yDS/6\nrmjNACUSUQpnkbpk/Xuw6PHyM0C17OXXmUHHSyNbn4gACmeR2m/XekjK8OG7cibkvAV9bvTTMrbs\npRmgRAJI4SxSGx3ZDytfhiVTYeN8+NZcSD8fBt0PQ37nz74WkcBSOIvUJgd3wus/88FccBCadoRL\nf1l6tnWD5MjWJyKVonAWqen2bIK9m6BtP6jXGDZ8AD2ugl43QFpfdVuL1EAKZ5GaKD/PDxCy5Gk/\nQldSBty+BGLi4PZPIKpaR+YVkTBTOIvUNAsfhTd/BUf2QWK6n/3p7OtKW8gKZpEaT+EsEiTvP+yv\nNy477OWKmfDJ0zD8IT/TU5M06DLcn23dtr/CWKQWUjiLBEnrPn7SiG/8E/IPwAd/gc2hYW7XvgVZ\n34LOQ/xNRGothbNIkLQb4IP56asABxYF3a/2XdcpnSJdnYhUE/WHiUTawZ2w4G8w63b/vOOl0P5i\n/7j/nXD1owpmkTpGLWeRSCiZAeqTqbD6tdIZoPLzfDf2l0thwN2Q/Si0H6ipF0XqGIWzSHVyzp9V\nveQpmP1jaJASmgHqemh+lr8savrNpfMjt7uo/HMRqRMUziJV7dAeWP6CbyX3uRHOuRkyx0Cj5tBp\nSPkZoDYvLh/E7Qb455sXK5xF6hCFs0hVWfu2H9v6s9lQeBiaZUK9BL+uQTJ0HXH8a/rfcfyydgMU\nzCJ1jMJZJJwO7oSGTf3jt38LO9dC7/G+27pVbw2lKSKVonAWOVNHDvipGJdMhS1L4K7PoH4iXPUo\nNG6pGaBE5LQpnEW+rp1r4b2H/AheBQchuQMM/Enp+uR2katNRGo0hbPI6dizCQqPQEpHcMWwchZ0\n/wb0vgHSzlO3tYiEhcJZ5FQKDsGq2X5863X/gcxRcM1TfmCQn+So21pEwk7hLHIy8x7wo3cd2QtN\n0mHgT6HXdaXrFcwiUgUUziJl7f8Slr8Ifb/rrz+OjoUuQ6HXOMi4SDNAiUi1UDiLFObDmn/7s61z\n3gRXBC3PhowLYcB/Rbo6EamDwhbOZjYU+DMQDfzLOfeHcO1bpMrsWgf/vBQO7fKXPV14u28la6IJ\nEYmgsPTRmVk08AgwDMgErjOzzHDs+2QWPPVzlv/l15CR4bsbMzJY/pdfs+Cpn1f1W3997z/sx08u\na/27fnmQ1cS6K6p51Wx45lp/HBkgMQPOuhLGzYAfr4DBv1Iwi0jEhesAWl8gxzm3zjmXDzwHjA7T\nvk+o0Z4YWu+YwvLmceAcy5vH0XrHFBrtCXBvfes+fiKDktAomeigdZ9IVnVqNbHukprXzoM1c+Gx\nYfD8ON+FnfOm3yYqCkZOgk6XQVR0RMsVESlhzrkz34nZ1cBQ59x3Qs/HA+c55yae6DVZWVkuOzv7\nzN44I8MH8uUH2X2kIRnx2ygsMOKKHDRtCg2awsSP/bYv3Qqfv17+9YlpMOEd//j58bBhfvn1qV3h\nllf94//3Ddj6afn1rc+BcdP840cvh5055de3GwBjH/eP/94f9m31j4vy4ch+aNEd9m3xExu88iM4\nvK/863tcDcMe9I//1BmKi8qvP+cmuPQX/lKf/+1+/OdzwQ/gorvg4A545Lzj1w+828+ItHsD/POS\n49df9mt//e62VfDEyNK6Y+L9WNEDfgKX3A+52fDMN49//Zi/QufLYd07MOPbx6+/5knI6A+fvVo6\nl3FZN8zwQ14unQ6v3XP8+m+95lu52Y/7oTKP9b13YddaX1tBHmDQbZT/uVtU8HmJiABMnQr33w8b\nN0J6OvzudzBu3Bnv1swWOeeyKrNtuJqYFY28cFzqm9kEYAJAenr6mb/rxo103+BY0PUczm/3OUuL\nM/g0qiMph/bQIb8+bTt3o17JtmnnQmz98q9vkFz6uG0/aJhafn1Cy9LH7QZAUkb59WVHgOpwiZ/y\nr6zUrqWPOw3xsxOV2LwItn7i5+xtNwC6DPchW1ar3qWPu43yg16U1aKHv7doyKygoyK1m7+PqVfx\n+qYd/X1co4rXl/y89RJK15fU3bKX/5kA6idV/PrGoc+vUfOK1zds5u8TWlW8vn6Sv09Mr3h9vcb+\nPrl9xetj6/vPtusIWDYd+t87oB8KAAAdSklEQVQJg39x/HYiUnWqKOjCorgYDh+GI0f8/eHDMGMG\n/OIX/jHAhg0wYYJ/XI11h6vlfAHwK+fc5aHn9wI45x440WvC3XJevaElndtu5aWV5zEtcTCrU9oy\ncvUHTG6wAW65BTdoEBYdkG7Lki7hrG9D9qM1Z67emlh3TaxZpLaYOtUHW15e6bIGDWDKFB90RUWl\noXiyW9nwDOetoKDyP0vbtvDFF2f0cZxOyzlc4RwDrAEuBTYDC4HrnXMrTvSacITz8r/8mtY7prB5\nbkO6f/w5y/t2ovXlB9ncdALF3YYQ+8rLdHv6H6yzBowf9wDfaHyIq8cOoG3vbmf0vmekJCxKQuLY\n50FVE+uuiTWL1ATOwf79sGsX7N59/H3J46lTywdzCTOIjobCwjOrwwzi40/vVq/eydffcov/+Sp6\nr+Li45efVrnVHM6hNx0OPIy/lOox59zvTrZ9OMJ5wVM/p9GeGLpPevxol8nyO2/hQGIh59/4336j\nw4dZ+ewr/GHhdt5rnIazKPru28TYzk244qbhxCc1OaMaTtv7D/sTlcqGw/p3YfPiiufyDYqaWHdN\nrFnkZMLdRXzoUMXhemzIHnu/Z49v9Z5IXBwkJ8OXX554m/vuO/NgjY0N/3j2GRm+K/tYNbHl/HWE\npVv7NG39bB0vPvs2M3ZEs7l+Igsf/z5NRo9gx/U30fSSizCN/iQiQXWiLuK//x2GDv16IXvkyInf\nLyoKEhN9yCYl+VvJ41Pd16/vQ7MKg67KnKor/gwonE/BFRez/o33aT/9KXj+ecZc+St2JzTl6uQC\nrrr2ElpldohIXSIiFdq5EzIzYdu2039to0YnDtGTBWxCwpkPV1uFQVelAnC2dp0M57LcgQO89Nhs\npq/azYdN0jFXTP/9m/herxT6j7/CfwMUEakue/fC4sWwcCFkZ/vb+vUnf81f/nLi8I2NrZ66TyTI\nZ2tXM4Xz17Tpk8+YMf1dZuyL57b3nuH6Lxaw97rxrB99LWcPuUDd3iISXgcPwpIlpSG8cCGsWVO6\nPiMDzj0XsrLgoYcqbjkHuYtYyonEdc61Qlqvrvy4V1d+VFhE0Tut4cknmLV4Ez9vsodOrzzO2ObG\nmHGX0axDWqRLFZGa5vBhWLq0NISzs2HlytIzgFu39iE8frwP5HPOgZSU0te3bl1xF/HvTnrurdRQ\najmfwr6vdvLq068xPecAi5u0Ibq4iEEHNvHIJS2pN2qkPytRRKSsggJYsaJ81/SyZaXX1aamlraI\nS24tW558n6Au4hpO3dpVJGfBUmbMnM/GL77ir8//ClJSmHbjTzhr5MWcNahvpMsTkUgoKoLPPivf\nNf3JJ6VnQicmlgZwSSCnpYX/EiAJPIVzVSsshDfe4NATT3Fuq29woF4DMvdtYWybOMaMH0JSWiW+\nAYtIzVNcDGvXlm8RL17sjx2DPzO6T5/yreIOHRTEAiicq9Xu3C+Z9fRcpm/MZ3lCK2KLCnho+weM\nGnsxXH45xOiwvkggnaqL2Dl/jW7ZFvGiRf5savCDYPTuXb5V3LmzH/lKpAIK5whZ9Z9sZsxeyE0v\nTiZ93Ure63MJ7w8ey9gxF9DxgrMjXZ6IlKjo+tv69eHWW/1JViWBvGOHXxcbCz17lu+azsyM/GVK\nUqMonCMtPx/mzOGRWUuYlHIORVHR9NqXy9h2Ddne7Wz6bl5Jvwd+evQb+/x7H2Rp1yy+P1CDn4hU\nixONXAW+5XvWWeVP1urZ0w8dKXIGdClVpMXFwZgx3DZmDNesy2Xm1DeYvi+G+3cm0frVhTwVG89k\nl0A/55jvEpi4spjJZIPCWaRq5ebC9OknDmYz2LfPt55FIkgt52riiotZ9uZH7Lzzburt38tto+8h\ntqiQffEN+cO//8KYvA0aSECkKnz1lZ+j9/nn4b33/LLY2IqnC9SAHlKFTqflrCGvqolFRdFzyAUM\nWvkB/TYu48oV89jWuCmHY+pxx6i7+Wa/7zPjz8+St2tvpEsVqfl27PDjN196KbRqBRMn+vGpf/1r\nf9nT448f3zrWgB4SIOrWrm7p6cx3Ccw8axC3f/AsT/UZweVrPuSj9B7819YE0s8fTN/+PcgbfzP1\nB/bXkKEilbVnD8yc6VvIb7zhrz/u1MlPTfjNb0L37qXbduni7zWghwSUwrmazb/3QX+MeeYD9Nu4\njPM3LmXimHuZ3C2Keq0b0mfH2fD88/xhUz3efWEdVycX8o1rB2mmLJGK7N8Pr7wCzz0Hc+f6kzHb\ntoW77oJrr4VevU58jfG4cQpjCSyFczVb2jWLyWTT7+V9YEY/28fkzKjSs7WvuBj+/GfOf/wV1qzc\nzZ/y03noyZX03/8643qmMPSmkZopS+q2vDx49VXfQn71VT9mdevWcNttvoXct68G/ZAaTyeEBdzG\n0ExZL+yN56K12fzhw6dw113Hiitv4KzBmilL6ogjR+C113wgz5rlR+Rq3hyuvtoH8oUXnvncwyJV\nTJdS1SLpvbpyZ6+u3FFYxMG3WkHTXXzy7w+4sslIOr2smbKkFisogDff9F3WM2f6S5yaNvVd0d/8\nJgwcqNG4pNZSy7kG2r9tJ6/8v9eYvvYASxL8TFkXH8zltxe1pOU3RmimLKm5CgvhnXd8C/nFF2HX\nLmjSBK680h9DvuQSjcolNZZGCKtDchYsZfrM+by1J5rZj04kPjGBD8b/kKRRw8i8+NxIlydyasXF\n8P77PpBnzIBt2/wEEqNH+xbykCEanUtqBYVzHeQKCrA33oDHH2dowsV8lprBWfu2MDYtjtHjLyep\nTYtIlyhSyjlYsMAH8vTpsGWLP9Fx5EgfyMOH68RHqXUUznXc7twvefn/zWX6pnxWJLQirrCAO/Yu\n5QdXZvlWiGbKkkhwzk+v+Pzz/rZxo28RDxvmA3nkSN9iFqmlFM5y1Mp3FjJjdjYXvv0Cly55iy0d\nMnnyqh8y9sp+dDy/Z6TLk9rOOVi2zIfxtGmQk+O/HA4Z4gN59Gh/TFmkDlA4y/Hy8+HVV5n1wrv8\nuNUgiqKi6b0vl7HtGzJy/DASmjeNdIVSkx07N/Jtt/nrkZ9/Hlat8pc5XXqpD+Qrr4Tk5EhXLFLt\nFM5yUtvWbuLlqW8w/UvHmoQWNMrP46Ptr9Lw5vEwaJCuF5XTU9HcyCUGDvSBfNVV0KxZ9dcmEiDV\nGs5mNhb4FdAN6Oucq1TiKpwjzxUXs/T1D1n67/cZ/9QfYM8e7rrmZ7TumMbV1wwk/ewukS5Rgqyo\nCD76yB8z3rfv+PWtW/spGkUEqP5w7gYUA/8A/kvhXEMdPkzBSzP5zns7ebdxOs6iOG/fJsZ2SWT4\njcN5aul2en6WTb8Hfnq063L+vQ+WDjsqdcOBA/D6634869mz/exPJ2LmL5MSEaCaRwhzzq0KvemZ\n7koiKT6e2Ouu5cnrYMvKtbz43DxmRMXzX1sbs3fkd+jZtB63ZQznEZdAP+eY7xL8BB5kg8K5dtu0\nyQfxrFnw9tv+/IWkJH+50xVXwE9+4rc5Vnp69dcqUkvomho5TqvMDkz8TQduKy5m4Sv/odP2ZJIe\nn8I1A6K54Zu/pe+m5XzWrB1/nfmAn8Dje9+MdMkSTiWXPM2a5VvIS5b45R07+nmRR43yY1mXXJJX\nWHj8MWfNjSxyRioVzmb2JlDRKBb3O+deruybmdkEYAJAur5VB55FRdF39CAYPQgen0K/jcuY3W0A\nC9qeDc7xt/PHsmPZm1xx6BCmASNqtkOHfKv4lVf8bcsWf2Jgv37wxz/6FnKXLhXP9lQy7aLmRhYJ\nm0qFs3NucDjezDk3BZgC/phzOPYp1SQ9ndiiAg7FxnPjoleY1vMyVqW24+/nXcWoVq3guuvYdM2N\ntBnQVzNl1RRffeWnXJw1C954w7d8GzWCoUN9GA8fDikplduX5kYWCSt1a0ulzL/3QX+MeeYD9Nu4\njKFr5jNxzL3c1jwfhg9n/9TnuKz+paRPe5yxLaMYM+4yUtu3iXTZUpZzsGJFaXf1Rx/5ZenpcMst\nvrt64ECNYy0SAGfcxDGzK80sF7gAeNXM5p55WRI0S7tmMTkzin62D8zoZ/uYnBnFpmFXwtSpxOR8\nzi9aHqIhhfzuYDPO//tivnPbX1k1daY/gUgiIz/fT7v4ox9B+/bQo4fvfi4qgt/8Bj75BL74AiZP\n1gQTIgGiQUgk7HI+/JTpMz/kpUONefrZ++hMHjnjJ1Aw+kq6DazUVQRyJnbtgn//27eQX3vNX4Mc\nHw+XXea7q0eOhJYtI12lSJ2jEcIkEIryC4h+43V4/HF+nN+OlzIvpvu+LYxNj2P0+KEktm4e6RJr\nj88/L+2ufv993zJu0cIH8RVXwODB/gxqEYkYhbMEzq5NX/Ly03OZUWamrOsO5vDrKzI1U9bXUVgI\nH37ow3jWLFi92i/v2dOH8ahRkJWloVhFAuR0wln/c6VaJKe14JZ7b+LVv36XOUObMY4tNFuzHEaM\noCgjgz/f81fWfrws0mUGw9SpkJHhgzUjwz8H2L8fZsyAG2/0reIBA+Dhh/0JXf/3f/7Y8aefwm9/\nC337KphFajC1nCVy8vNh9myWPTebMRljKIqKps/eXMZ2bMTIG4fROLUOzlxU0SQSsbH+GuPVq6Gg\nwM/oNGKEbyFffjkkJESuXhGpNHVrS42zbe0mZk59g+lfOT5v3IL4giPM3PkWXcd/o+7MlLV7N2Rm\nwpdfHr8uJgbuuMN3V19wgQ4DiNRACmepsVxxMZ++/iFz5i7ip0/8kug9e/jXkFvYl3U+Y68ZSFpN\nnykrPx/WrvWt4DVr/H3JTZNIiNRqCmepHQ4dgpkzufPtzbyU3AVnUZy/byNXd01m+PhhNEgKaHeu\nc7B16/EBvGYNrF/vz6Qu0by577IuuT34IGzffvw+27b1x5RFpMZSOEuts3nlWl56bh4zdsXwRaNU\nrlr1Hx5K2Qm33ILr1y8yQ4YeOOAvYSrb+l2zxt/27y/drn596NzZ38oGcefO0KRJ+X1WdMy5QQOY\nMkXDY4rUcApnqbVcaKashNkv0/XZf7G6fgq3XvNLrkop4qrrLqVF13bhfcOiItiwoXz4ljzevLl0\nOzN/1vSx4dulC7Rpc3rHzKdO1SQSIrWQwlnqhgMH+OTpl/n90v18nJBGVHER/Q/kMrZncy6/YRiP\nLdpKz8+y6ffAT48G3fx7H2Rp1yy+f+wc1Dt3VnwcOCen/PCjiYnHh2+XLn46Rc3MJSIncTrhrFM+\npeZq1Ihe3x/HNGDDklXMmP4eL0Q35u41MKh9B3peeDm3pV/OIy6Bfs4x3yUwcUURk1dPhQ/jyreG\nd+4s3W9MDHTo4EN3xIjyYZyaWvG0iSIiYaSWs9QqRYVFrJ39Fp2nPwnPPMPF3/0HmxJbcPaWNaxJ\nTeePcx5m+JoP/cYtWlTcCm7XTpcqiUjYqVtbBCiOiua5npcx+YJr2NKkdBzvmxe9wq9mPIhLSGD1\nV/vp1Kwx0VFqDYtI1VK3tggQlZ5Gxu4tHI6N5wcfPs/TvUdwxcp3GXhgIzRpwubdeQx9+D0a1Yuh\nd3oifdKTOKdtEn3aJtGonv5riEjk6C+Q1Frz732QiSuLmTzzAfptXEb/Lz5h4ph7GZE5AIAm9WN5\n+Ju9WLRhN4s27Ob/3v6cYgd/vrYXo3u1ZtOuPD5ev4tz2ibRtmkDTMeaRaSaKJyl1lraNYvJZNPv\n5X1gRj/bx+TMKJZ2zaIf0Dg+ljG9WzOmd2sADhwp5NNNe8hs6Qc3eWf1Nn7+8goAUhrFHW1ZX9s3\nnSb1YyP1Y4lIHaBjziInUFzs+HzbARZt2E32hl0s3rCbjbvy+PSXQ2gcH8v07E2s+Wq/7wpPT6JZ\nQnykSxaRANMxZ5EwiIoyurRoTJcWjbn+vHQA9uTl0zjet5pXf7mfpxZs4J/vrQcgLbk+/Tum8sA3\nekSsZhGpHRTOIqchsUHc0cc/G5nJ3UO7snzLXhaHjlvvP1xwdP03//EhMdHGOelJnJORTK+0RHWH\ni0ilKJxFzkBcTBR90n239ncuKl3unCOzVQIfrdvF5Hk5FDs/dskt/drxiysyAdi4M4+05PrlTjT7\n+3/W0rNNE/p1SDm6bP7aHSzN3Xv8qGYiUmspnEWqgJnxyyvOAkpPNFu8YTedWzQG4Kt9hxnwP/No\n2jCOPm39iWbntE2ia4vGTHxmCZOv702/DinMX7vj6HMRqTt0QphIBOw9VMCcZVuPXsa1fsdBACZd\nczYtmsTzg6cX07NNExZt3M0PL+nEhR1SSGwQS4sm8cRGR2AGLhE5YxohTKSG2XngCIs37qFXWiKp\njetxy+MLmbd623HbzZp4IT3bJPLyJ5v52ztrSWwQS1KDOJIaxpHUIJbv9G9PUsM4cnfnsW3/EZJD\n6xLiY3SdtkiE6WxtkRqmaaN6XJbphxidv3YHn27aw3cvase07Fx+PLgTrZMasPtgPm2TGwKQUD+W\ntOQG7MnL5/NtB9h9MJ89hwq46YIMAGYsyuXhNz8/uv/oKCOpQSxv3XUxTerH8tKSXBas3UViQx/u\nyQ3iSGwQy+BuzYmKMg7lFxEbbcRUopWu4+Qi4adwFgmQo8eYx/ljzoO6Njt6zPmyzLSj2w3q0oxB\nXZqVe21xsTs6YdZVfdpwdptEdufls+tgPnvyCtidl390WNINO/OYt3obe/IKyC8qBvzJbav/eygA\nP5u5nBcW59KkfixJDWJJahhH68T6TL6+DwBzV3zJzgP5JDWIJSbKuPXpxTx4VQ+Gdm9ZI46T6wuF\nBN0Zd2ub2f8AVwD5wFrgFufcnlO9Tt3aIser7tBwzpGXX8Sug/nsP1xIZis/Otpbq75iae5edufl\nszuvgN0H84mLieKxm88F4IZ/fcT7OTvK7Sva4LZBHXn6o420bBJP7u5DxMdGUT82mvjYaHq2acIf\nrz4bgN/PWcWOA0eoHxt9dH3HZo2OjtY2d8WXFBc74kPr6sdF07RhHGnJDQDYf7iAejHRxEbb1+qu\nL/sF4tgT78p+9kGiLxQ1X7UeczazIcDbzrlCM3sQwDn301O9TuEsUnMdLihiT15BqFXuA/zVpVuY\ns/xLbr+kI0kN49iwM49D+UUcKvC3dikNuW94NwBufvxjPv/qAIcLijgcWj+wcyqP39IXgPN+/yZf\n7TtS7j1H9GzJI6GWe49fzWX/4UKioywU7lF8o0+bo/u/dsqHxMVEEx8TRf04/wVgYOdUhvVoSUFR\nMY9/sJ4tew4zPXsT/Tul8EHOTn51RSZXZ6WRX1jMp7l7iI2OIibKiIuJIjY6iqaN4kiIj6Wo2JGX\nX0hstF9eXTOa6QtF9ajKmqv1mLNz7vUyTxcAV5/pPkUk2OJjo2nRJJoWTfyQpfPX7mDB+l3cfolv\nOU++vje3XNjuhK9/IhTCJZxzFBWXNhRmfL8fB/MLOZRfxOGCYg4XFNG0UekAMHde1pmDRwo5XFDM\noVDAd2nuL1Mr2c++QwVsKxP+LZrEM6xHS/KOFPH7OZ8d3dfcFV8BsHnPYQB2HjzC2L9/eFzNPxvR\nje9c1J71Ow4weNK7R5dHGcRGR/HbMd0Zm5XGyi37+O5T2aFQNx/y0VH8ZEgX+ndKYeWWfUx6Y83R\ndbHRUcTFGDf3a0eXFo3J2baflz/ZcnRdbLT/gnD5WS2YfH1vbn16MRd0aMr7n+/guxe1Y//hQt5Y\n+RXntU8mIT6WzXsOkbPtAFEGUWZY6L5XWiLxsdF8te8wW/YcIsqs3PrOzRsREx3FzgNH2J1XcPT1\nJdu0TqxPVJSx/3ABhwqKjq6LMn/pYMkAO/mFxRQ7R5QZZ7VKYOLU0BeIjjXj0sCebZoE4nLGcB9z\n/hbwfJj3KSIBdmwL7vwOTU+7RWdmxESXtkBLuq9P5GTBHx1lPDfhghOuT6gfw4pfX857n2/npy8s\nY9TZLZn16VY6Nfcn2yU1iGPqd84jv6iYgsJiCoocBUXFnBXq8k9uWI/7h3ejoLiYgkK/rqComE6h\nLwcN4qI5r30yBUWOwtC6/CJHbOjnO1xYxJY9h46+rmT/V5zdCmhMzraDTJ6Xw7Gdmme1SqBfhxTO\nzUjiteVfAvC/ZU76+/ePLiKhZSxvrvyKX85acdzP/d7dg0hLbsALi3P542urj1u/+OeXkdwwjsc+\nWM8j89Yet/6z/x5KfFQ0D72+hifmf1FuXUyUkfP74QDc99IyZizKLbd+3L8+4oehL26dmzfipsc+\nxkLBHmVGq8T6vHnnQABum7qY+Wt3hL4U+C8GHVIbHv2d3vr0IpZv2Vv6xQHIbJVw9HyIW59exBc7\n88p9OemdlsivR3c/uv/t+49gxtEvJudmJPPjyzoD8OzHm2ib3ICbHvuY6/um88rSrRHpnahUOJvZ\nm0CLClbd75x7ObTN/UAhMPUk+5kATABIT08/7WJFJHiW5u4t98erX4cUJl/fm6W5ewPZ3WpmfJq7\nh/teWs7fbuhDvw4pDOvRkonPLCGxQRz9OqRwYccT153cMI7vDmh/wvUZKQ2ZdE2vE67vk57EnB9d\ndML1Q7u3YP0DIygqdqFg918SGsfHMn/tDhZt2M2489KZvXQr9w7rSo82TXAOMpr6LxfDerSge+sm\nOOcodlDsHMXOkdq4HgAjerSkW8sEv764ZD00rBcNwMierejcvDEu9NqS+5Lr60f0bEmHZo2gzP7L\nGt6jBe1TG/rXFfttPlq/k7+8ncPtl3SkQ7NG9E5P8q8Lvb5kvHqA89onk9wwDod/rXOO1Malk8p0\na5lAfGx0uZ+v5GcHaJ4QT0GRA0rXN4ovjbqYaCM6yvzPXQxF+C9JJXYdPEJhsSOxfixPfriB2y/p\nGJF/x2G5ztnMbgK+D1zqnMurzGt0zFlEIqUmHguticecobTuG85LP3rII8j1QtXVXN0nhA0FJgED\nnXPbK/s6hbOISOXpC0X1qMqaqzucc4B6wM7QogXOue+f6nUKZxGR2q0mfqEIytnaGr5TRESkGpxO\nOGsEfRERkYBROIuIiARMxLq1zWw7sCGMu0wBdpxyKwkHfdbVQ59z9dDnXD30OUNb51xqZTaMWDiH\nm5llV7YvX86MPuvqoc+5euhzrh76nE+PurVFREQCRuEsIiISMLUpnKdEuoA6RJ919dDnXD30OVcP\nfc6nodYccxYREaktalPLWUREpFaoFeFsZkPNbLWZ5ZjZPZGupzYyszQzm2dmq8xshZn9KNI11WZm\nFm1mS8xsdqRrqa3MLNHMZpjZZ6F/1yeeZ1LOiJn9OPR3Y7mZPWtm8ad+Vd1W48PZzKKBR4BhQCZw\nnZllRraqWqkQuMs51w04H7hNn3OV+hGwKtJF1HJ/Bl5zznUFzkafd5Uws9bA7UCWc647EA1cG9mq\ngq/GhzPQF8hxzq1zzuUDzwGjI1xTreOc2+qcWxx6vB//h6x1ZKuqncysDTAC+Feka6mtzCwBGAA8\nCuCcy3fO7YlsVbVaDFDfzGKABsCWCNcTeLUhnFsDm8o8z0WhUaXMLAPoDXwU2UpqrYeBu4HiU20o\nX1t7YDvweOjwwb/MrGGki6qNnHObgT8BG4GtwF7n3OuRrSr4akM4WwXLdAp6FTGzRsALwB3OuX2R\nrqe2MbORwDbn3KJI11LLxQB9gL8553oDBwGdr1IFzCwJ35vZDmgFNDSzGyJbVfDVhnDOBdLKPG+D\nukyqhJnF4oN5qnPuxUjXU0tdCIwysy/wh2guMbOnI1tSrZQL5DrnSnp/ZuDDWsJvMLDeObfdOVcA\nvAj0i3BNgVcbwnkh0MnM2plZHP5Eg1kRrqnWMTPDH59b5ZybFOl6aivn3L3OuTbOuQz8v+W3nXNq\nZYSZc+5LYJOZdQktuhRYGcGSarONwPlm1iD0d+RSdPLdKcVEuoAz5ZwrNLOJwFz8WYCPOedWRLis\n2uhCYDywzMw+CS27zzk3J4I1iZyJHwJTQ1/q1wG3RLieWsk595GZzQAW46/6WIJGCzsljRAmIiIS\nMLWhW1tERKRWUTiLiIgEjMJZREQkYBTOIiIiAaNwFhERCRiFs4iISMAonEVERAJG4SwiIhIw/x92\nNkUOBhyXKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6));\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.step(I_obs, '-bx', where=\"post\")\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(Y_obs,'-or')\n",
    "\n",
    "# Plot non-zero observations\n",
    "for i in np.nonzero(np.any(Y != 0, axis=0))[0]:\n",
    "    plt.plot(Y[:,i],'--x')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, the process is rarely at steady state when control is initiated. If we somehow knew the past 10 moves in $u$, we could get a series of 10 future predictions of $y$ formulated by extending the preceding logic\n",
    "and going back in time. \n",
    "\n",
    "Let us call this the vector $y_p$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "U_hist = np.array([1,0,0,0,-1,0,0,0,0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.5   -0.5   -0.5   -0.5   -0.5   -0.5   -0.5   -0.5   -0.5   -0.5  ]\n",
      " [-1.    -1.    -1.    -1.    -1.    -1.    -1.    -1.    -1.    -1.   ]\n",
      " [-1.5   -1.5   -1.5   -1.5   -1.5   -1.5   -1.5   -1.5   -1.5   -1.5  ]\n",
      " [-1.8   -1.8   -1.8   -1.8   -1.8   -1.8   -1.8   -1.8   -1.8   -1.8  ]\n",
      " [-1.9   -1.9   -1.9   -1.9   -1.9   -1.9   -1.9   -1.9   -1.9   -1.9  ]\n",
      " [-1.95  -1.95  -1.95  -1.95  -1.95  -1.95  -1.95  -1.95  -1.95  -1.95 ]\n",
      " [-1.98  -1.98  -1.98  -1.98  -1.98  -1.98  -1.98  -1.98  -1.98  -1.98 ]\n",
      " [-1.99  -1.99  -1.99  -1.99  -1.99  -1.99  -1.99  -1.99  -1.99  -1.99 ]\n",
      " [-1.999 -1.999 -1.999 -1.999 -1.999 -1.999 -1.999 -1.999 -1.999 -1.999]\n",
      " [-2.    -2.    -2.    -2.    -2.    -2.    -2.    -2.    -2.    -2.   ]]\n"
     ]
    }
   ],
   "source": [
    "A_hist = np.tile(a, (U_hist.shape[0], 1)).T\n",
    "print(A_hist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "row_shift_indices = np.atleast_2d(np.arange(a.shape[0])).T + np.flip(np.arange(A_hist.shape[0]))\n",
    "row_shift_indices = np.triu(row_shift_indices)\n",
    "col_indices = np.arange(A_hist.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-2.    -1.999 -1.99  -1.98  -1.95  -1.9   -1.8   -1.5   -1.    -0.5  ]\n",
      " [-2.    -2.    -1.999 -1.99  -1.98  -1.95  -1.9   -1.8   -1.5   -1.   ]\n",
      " [-2.    -2.    -2.    -1.999 -1.99  -1.98  -1.95  -1.9   -1.8   -1.5  ]\n",
      " [-2.    -2.    -2.    -2.    -1.999 -1.99  -1.98  -1.95  -1.9   -1.8  ]\n",
      " [-2.    -2.    -2.    -2.    -2.    -1.999 -1.99  -1.98  -1.95  -1.9  ]\n",
      " [-2.    -2.    -2.    -2.    -2.    -2.    -1.999 -1.99  -1.98  -1.95 ]\n",
      " [-2.    -2.    -2.    -2.    -2.    -2.    -2.    -1.999 -1.99  -1.98 ]\n",
      " [-2.    -2.    -2.    -2.    -2.    -2.    -2.    -2.    -1.999 -1.99 ]\n",
      " [-2.    -2.    -2.    -2.    -2.    -2.    -2.    -2.    -2.    -1.999]\n",
      " [-2.    -2.    -2.    -2.    -2.    -2.    -2.    -2.    -2.    -2.   ]]\n"
     ]
    }
   ],
   "source": [
    "M = A_hist[row_shift_indices,col_indices]\n",
    "M[np.tril_indices(M.shape[0])] = a[-1]\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_initial = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.05  -0.02  -0.01  -0.001  0.     0.     0.     0.     0.     0.   ]\n"
     ]
    }
   ],
   "source": [
    "Y_pred = np.matmul(M,U_hist) + y_initial\n",
    "print(Y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The preceding convolution model will be wrong due to unmodelled load disturbances, nonlinearity in the models, and\n",
    "modelling errors. One way to correct for these is by calculating a prediction error $b$ at the current time and then adjusting the predictions. In equation form:\n",
    "\n",
    "$$ b = y_0^{actual} - y_0^{pred} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y0_actual = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the vector y^{pred} is corrected by $b$ for predictions 1 to 10 as \n",
    "\n",
    "$$ y_i^{pc} = y_i^{pred} + b  $$\n",
    "\n",
    "Where $y^{pc}$ is the vector of corrected predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = y0_actual - Y_pred[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation becomes\n",
    "\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "y_1 - y_0 & = a_1 \\Delta u_0 + y_1^{pc} \\\\\n",
    "y_2 - y_0 & = a_2 \\Delta u_0 + a_1 \\Delta u_1 + y_2^{pc} \\\\\n",
    "y_3 - y_0 & = a_3 \\Delta u_0 + a_2 \\Delta u_1 + + a_1 \\Delta u_2 + y_3^{pc} \\\\\n",
    "y_4 - y_0 & = a_4 \\Delta u_0 + a_3 \\Delta u_1 + + a_2 \\Delta u_2 + y_4^{pc} \\\\\n",
    "\\vdots & \\\\\n",
    "y_t - y_0 & = a_t \\Delta u_0 + a_{t-1} \\Delta u_1 + + a_{t-2} \\Delta u_2 + y_t^{pc} \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "In vector form:\n",
    "\n",
    "$$ y = y^{pc} + A\\Delta u $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Control Law\n",
    "\n",
    "From the preceding text, we have a set of\n",
    "future predictions of $y$ based on past moves\n",
    "and a set of future values of $y$ with a set of\n",
    "three moves. \n",
    "\n",
    "Let us now assume that we\n",
    "have a desired target for $y$ called $y^{SP}$. (We,\n",
    "of course, need this for any control action!)\n",
    "\n",
    "One criterion we can use is to minimize the\n",
    "error squared over the time between now\n",
    "and the future (that is, to use the least squares criterion). In vector equation form:\n",
    "\n",
    "$$ \\min (y^{sp} - y)^2 $$\n",
    "\n",
    "Substituting the previous equation:\n",
    "\n",
    "$$ \\min (y^{sp} - y^{pc} + A\\Delta u)^2 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}