{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NOTEBOOK_HEADER-->\n",
    "*This notebook contains course material from [CBE30338](https://jckantor.github.io/CBE30338)\n",
    "by Jeffrey Kantor (jeff at nd.edu); the content is available [on Github](https://github.com/jckantor/CBE30338.git).\n",
    "The text is released under the [CC-BY-NC-ND-4.0 license](https://creativecommons.org/licenses/by-nc-nd/4.0/legalcode),\n",
    "and code is released under the [MIT license](https://opensource.org/licenses/MIT).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [Linearization](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/03.00-Linearization.ipynb) | [Contents](toc.ipynb) | [Linear Approximation of a Process Model](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/03.02-Linear-Approximation-of-a-Process-Model.ipynb) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE30338/blob/master/notebooks/03.01-Step-Response-of-a-Gravity-Drained-Tank.ipynb\"><img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open in Google Colaboratory\"></a><p><a href=\"https://raw.githubusercontent.com/jckantor/CBE30338/master/notebooks/03.01-Step-Response-of-a-Gravity-Drained-Tank.ipynb\"><img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Step Response of a Gravity Drained Tank\n",
    "\n",
    "by Jeffrey Kantor  (jeff at nd.edu). The latest version of this notebook is available at [https://github.com/jckantor/CBE30338](https://github.com/jckantor/CBE30338). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary\n",
    "\n",
    "In the example we show how to fit the step response of a nonlinear system, a gravity drained tank, to a first order linear system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Mass Balance for Tank with Constant Cross-Sectional Area\n",
    "\n",
    "For a tank with constant cross-sectional area, such as a cylindrical or rectangular tank, the liquid height is described by a differential equation\n",
    "\n",
    "$$A\\frac{dh}{dt} = q_{in}(t) - q_{out}(t)$$\n",
    "\n",
    "where $q_{out}$ is a function of liquid height. Torricelli's law tells the outlet flow from the tank is proportional to square root of the liquid height\n",
    "\n",
    "$$ q_{out}(h) = C_v\\sqrt{h} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Dividing by area we obtain a nonlinear ordinary differential equation \n",
    "\n",
    "$$ \\frac{dh}{dt} = - \\frac{C_V}{A}\\sqrt{h} + \\frac{1}{A}q_{in}(t) $$\n",
    "\n",
    "in our standard form where the LHS derivative appears with a constant coefficient of 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Nonlinear Step Response"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4XPV59vHvo92yJMursI1t2XjBNgGDwCzJC3bCHhqS\nvDSFkBTS9HJpQ5uULkDSBfK2TbolpU0baFpCCCFuU0IxDgRCYkMgGG94NzbeLVleZNmWRpK1zDzv\nH3MsD8KSRrJGZ0Zzf65rrjnL75y5fWZ8Hp3d3B0RERGAnLADiIhI+lBREBGRTioKIiLSSUVBREQ6\nqSiIiEgnFQUREemkoiAiIp1UFEREpJOKgoiIdMoLO0BfjRkzxisrK/s1bVNTE8OHDx/YQAMgXXNB\n+mZTrr5Rrr4ZirnWrFlT5+5je23o7hn1qqqq8v5atmxZv6dNpXTN5Z6+2ZSrb5Srb4ZiLmC1J7GO\n1e4jERHppKIgIiKdVBRERKSTioKIiHRSURARkU4qCiIi0klFQUREOmXcxWsiIhC/xioaczpi8feo\nO7GgP5YwPOaJ73R2Jw6PecLwGGyq64Bthzv7T7U5NV3M459/arwTH4efbuucbueJw4NHIMeC4Z7Q\nDRCLOQ7BuNNtcKegIcqCFC9XFQUReR93pz3qtLRHaW2PcrI9RmvH6ffWjuC9PUZrR4y2jtPD26Lx\n/vbOd6ctGqM9GNbZH41xpK6Fb297k/ZojI5Y/DM7ojGiMac9FiMaddpj8WGnVvKdRSCW4ufLr16V\n2vn3w0en5qf8M1QURDKcu9PcFiXS2kHjyXYaT3YQae0gErw3t0XZuLONVa3v0NQapbmtg6a2KM3B\nuJPtUZrborS0x7tbgu6zXefm5hgFuTnk5xoFebkU5Bp5QX9+bg75uTm0RqEEKC7IIy/XyMvJIS/H\nyAva5OYY+cHwU925OTnk5tDZNifHyMsxcnOMHIu/5+UauRYflxsMO91NZ7uchDY5BhYMX7/ubS65\n5JKgTbx9vG1Cd9D+1HiAnKC9cXq8dU4fH45xehzxccHknd3G6fEWtAVYvnz52X0pSVBREEkTHdEY\nx5rbOdbcRn1TG8eb2zjW3M7x5naOt7RxIuhuOBm8WjpoCIpAMn815+7cxfCCXIYX5jGsIJfiglyK\n8/MoLy5gQnkuwwpyKcrPZVjwKsrPoSg/N+GVQ2FeLoV58eGFeTkU5OV0vse748NPrdB7s3z5chYs\nuHIgFt+Aat6bS9WUkWHHCIWKgkgKdURjHIm0crihlcONrRxuPEldYxtHIvH3o02tHG1q42ikjRMt\n7d3OpyA3h/LifEYMi7/GlRYxfWweZcPyKSvKp7Qoj9KifEqK8igtzKOkKI/hBXmUFOYxvDCXNW/9\nius+vKDzL06R7qgoiPRTLOYcibRSc7yFmmMtHDjewqotrfxw/2pqT5zk4ImT1EVaz7gbprw4nzEl\nhYwpKWD2+DJGDy9g9PBCRg3PZ+TwAkYVF1BeXEB5cT4jiwsoys85qxV6Qa6pIEhSVBREenCyPcq+\n+mb21DWx92gze+ub2FffQnV9M9XHW2jriL2n/bA8mDS6iXNGDGP2OWVUjCiioqyQcaVFjCstZFxZ\nIaOHF1KQp7PBJT2pKIgAx5ra2H6okXcPR9h5JMLOI03sPBzhwImWzlMFAcqK8pgyejjnjy/lujkV\nnDtyGBNHDmNieTETyotYs+INFiy4Jrx/iMhZUlGQrNLWEePdw41srW3kndoG3jnYyDsHG6mLtHa2\nGZafy3njhlM1ZSS3jTmXaWOHUzl6OFNGF1NeXBBiepHUU1GQIas9GmPbwUY2VJ9gY81xNtacYPvB\nCG3R+C6fovwcZlaUsnDWWGZWlDKjooQZFaWMLysiJ4kzZ0SGIhUFGTKORlpZvfcYa/ce4+19x9lQ\nc5yT7fECMGJYPheeO4LPfaiSuRNGMHdCGZWjhyd12qRINlFRkIx1uOEkb+46yopd9azcfZSdR5oA\nyM815k4YwafnT2He5HIuOncEk0cV6+wbkSSoKEjGaGrt4K3dR3ltex0vr2/mwE9/DkBpYR6XTR3F\nbVWTuKxyJBdMHEFRfm7IaUUyk4qCpLXddU384p3DLN92mLd21dMWjVGYl8OMETncdfUMrjxvNHMn\njNBuIJEBoqIgacXd2VB9gpc2H+TlLYfYcTgCwPRxJdx11RSumTmOSytHsuKNX7LgmvNCTisy9Kgo\nSOhOFYKlGw7wwsaD1BxvITfHuHzqKO68fDLXzq5g0qjisGOKZAUVBQnNnrom/nddDc+tO8Duuiby\nc40PTR/Dl66dwXVzKnRNgEgIVBRkUDW1dvDCxlp+tKaalbvrMYMrpo7mnmumcePc8YwoTv394kWk\neyoKMii2HWzkqRV7efbtGiKtHUwdM5w/vXEWn7h4IuNHDAs7nogEUlYUzGwS8CRQQfxpcv/u7o90\naWPAI8DNQDNwt7uvTVUmGVwd0Rg/23KI776xh5V76inIy+GWC8dzx/zJXDplpK4bEElDqdxS6AD+\nyN3XmlkpsMbMfubuWxLa3ATMCF6XA98O3iWDRVo7WLxyH999Yw81x1uYNGoYX7l5NrdVncvI4TpO\nIJLOUlYU3L0WqA26G81sKzARSCwKtwJPevxJ1ivMrNzMxgfTSoapb2rjiTd2870393KipZ35laP4\n81vmcN2cCl1HIJIhBuWYgplVAhcDb3UZNRHYn9BfHQxTUcgg9U1tPPbaTp781V5a2qPcMLeCe645\nj4snZ+fjDEUymbmf5dO5e/sAsxLgVeCv3f3HXcYtBb7u7q8H/T8H7nf31V3aLQIWAVRUVFQtXry4\nX1kikQglJSX9mjaV0jUX9Jytqd15cXc7r+xtpzUKV4zP5dfOK2BCSeofIJOuy0y5+ka5+uZsci1c\nuHCNu1/aa0N3T9kLyAdeAu7rZvxjwB0J/duA8T3Ns6qqyvtr2bJl/Z42ldI1l/uZs51s7/DvvLbT\nL3r4Ja98YKnf+/Raf/dQQ+i50oFy9Y1y9c3Z5AJWexLr7VSefWTAfwJb3f0b3TRbAtxrZouJH2A+\n4TqekLbcnRc2HuRvXthKzfEWrp45lvtvnMXcCSPCjiYiAySVxxQ+CHwW2Ghm64JhXwYmA7j7o8AL\nxE9H3UH8lNTPpTCPnIWttQ08/PxmVuyq5/xzSnnq85fzoRljwo4lIgMslWcfvQ70eMpJsEnzhVRl\nkLN3ssP5q6VbePyN3ZQNy+evPn4Bd8yfrLOJRIYoXdEs3XplyyG+/HoL9Sd38+nLJ/OnN8zS/YhE\nhjgVBXmfY01tPPT8Zp5bd4CJJcYzv3slVVNGhR1LRAaBioK8x8+2HOLLz27kWFMbf3jtTObmVKsg\niGQRFQUBoKUtyleXbuaHK/dz/jmlPPG5y5g7YQTLl9eEHU1EBpGKgrDlQAN/sPhtdh6JcM8153Hf\ndTMpyEv9BWgikn5UFLKYu/P0yn08vGQL5cX5PPX5y/ngdJ1mKpLNVBSy1Mn2KF95dhPPrK3mmplj\n+canLmJ0SWHYsUQkZCoKWWh/fTO/8/01bD3YwBc/MoMvfmQGObruQERQUcg6q/bU8zvfX0NHNMbj\nd1/Gwlnjwo4kImlERSGLPLOmmgd/vJFzRw7jP+66lGlj0+8ukCISLhWFLODufPOVd/nnn7/LVeeN\n5tt3VjGiOD/sWCKShlQUhriOaIw/f24TP1y5n1+vOpe/+eQHyM/V6aYicmYqCkPYyfYo9z79Nq9s\nPcS9C6fzR9fPJH5HcxGRM1NRGKKaWjv4/PdW8dbueh7+2Fzuuqoy7EgikgFUFIagxpPtfO67q1i7\n7xjf/NQ8Pn7xxLAjiUiGUFEYYk60tHPX4yvZVHOCf7njEj564fiwI4lIBlFRGEIirR3c9fhKNh84\nwb/deQnXzz0n7EgikmFUFIaIlrYon39iFRtrVBBEpP90buIQ0NoR5Z6n1rByTz3f+NRF3KCCICL9\npKKQ4WIx577/Xs+r24/wtU98gFvn6aCyiPSfikKG+9qLW/nJhloevOl8bp8/Oew4IpLhVBQy2OOv\n7+Y7v9zN3VdVsujqaWHHEZEhQEUhQ/10Uy3/7ydbuGFuBX9+yxxdqSwiA0JFIQNtPnCCP/yv9cyb\nVM4jt19Mrp6FICIDREUhw9RFWln05BrKi/N57LNVFOXnhh1JRIYQXaeQQdo6YvzeU2upi7TyP/dc\nxbjSorAjicgQo6KQQb66dDMr99TzyO3z+MC5I8KOIyJDkHYfZYjn1tXw1Ip9LLp6mq5FEJGUUVHI\nADsOR3jwxxu5dMpI/uSGWWHHEZEhTEUhzbW0RfnCD9ZSlJ/Lv3z6Yj01TURSSscU0tzDz29m++FG\nvnv3ZYwfMSzsOCIyxOnPzjT20uaDLF61n9+5+jwWzBoXdhwRyQIqCmnqcONJHvzxRuZOKOO+62aG\nHUdEsoSKQhpyd+7/nw00tXbwT78xj4I8fU0iMji0tklDT6/cx7JtR3jgpvOZUVEadhwRySIpKwpm\n9riZHTazTd2MX2BmJ8xsXfD6i1RlySQ1x1v4m59s5UPTx3DXlZVhxxGRLJPKs4+eAL4FPNlDm1+6\n+y0pzJBR3J0/e3YjMYevffID5OhGdyIyyFK2peDurwH1qZr/ULRk/QGWbTvCH98wi0mjisOOIyJZ\nKOxjCleZ2QYze9HM5oacJVRHI608tGQz8yaVc/dVlWHHEZEsZe6eupmbVQJL3f2CM4wrA2LuHjGz\nm4FH3H1GN/NZBCwCqKioqFq8eHG/8kQiEUpKSvo1bSpFIhF+sCuPlbVRvnrVMCaWhl2rT0vnZaZc\nyVOuvhmKuRYuXLjG3S/ttaG7p+wFVAKbkmy7BxjTW7uqqirvr2XLlvV72lR69JlXfMr9S/3vf/pO\n2FHeJ12XmXL1jXL1zVDMBaz2JNbFof1JambnWPAMSTObT3xX1tGw8oSlIxrj+1tamVg+jC8snB52\nHBHJcik7+8jMfggsAMaYWTXwl0A+gLs/CtwG/K6ZdQAtwO1BNcsqT63YS3XEefQzsxlWoKeoiUi4\nUlYU3P2OXsZ/i/gpq1mrLtLKP/5sO3NH53DD3HPCjiMiEvrZR1nt7376Di1tUe6cXUiwJ01EJFQq\nCiHZWtvAj9ZUc/dVlUwo0dcgIulBa6OQfP3Fdygryuf3P3zGs3BFREKhohCCN3bU8er2I9y7cDoj\nivPDjiMi0klFYZDFYs7XXtzKxPJhfPbKKWHHERF5DxWFQfb8hgNsqmngT26YRVG+TkEVkfTS4ymp\nZvbJJOZx0t1fGKA8Q1p7NMY/vLyNuRPK+NhFE8KOIyLyPr1dp/Ad4Dmgp/MlrwZUFJLwzJpq9te3\n8N27L9BtsUUkLfVWFF5099/qqYGZPTWAeYasto4Y31q2g4smlbNg1tiw44iInFGPxxTc/TO9zSCZ\nNgLPrK2m+lgLX7p2hi5UE5G0ldRtLswsF/go8buedk7j7t9ITayhpa0jxrd+sYN5k8pZMFNbCSKS\nvpK999HzwElgIxBLXZyh6Zm11dQcb+GvP3GBthJEJK0lWxTOdfcLU5pkiGqPnt5KuEZbCSKS5pK9\nTuFFM7s+pUmGqBc21lJzvIV7F07XVoKIpL1ktxRWAM+aWQ7QTvwUVXf3spQlGwLcncde3cX0cSV8\n+PxxYccREelVslsK3wCuBIrdvczdS1UQevf6jjq21Daw6Oppui5BRDJCskVhP/FnLWfdk9HOxmOv\n7mJcaSG3ztPVyyKSGZLdfbQLWG5mLwKtpwbqlNTubao5wes76njgpvMpzNM9jkQkMyRbFHYHr4Lg\nJb3499d2UVKYx6cvnxx2FBGRpCVVFNz94VQHGUoOnjjJTzbW8lsfrKSsSM9LEJHM0eMxBTN7qLcZ\nJNMm2zz91l5i7vzmlZVhRxER6ZPethR+28waehhvwO3AQwOWKMO1dcR4euV+PjxrHJNGFYcdR0Sk\nT5K5dXZpEm0k8OKmWuoirXqqmohkpB6Lgo4l9N2Tb+6lcnQxV8/QLS1EJPPocZwDaFPNCdbsPcZn\nrpiii9VEJCOpKAyg77+5l6L8HH69alLYUURE+kVFYYA0nGznufU1fHzeREYU6zRUEclMPR5TMLN/\nAbq9tYW7/8GAJ8pQz68/wMn2GHfM18VqIpK5ejv7aPWgpBgCfrS6mlkVpVx47oiwo4iI9FtvZx99\nL7HfzIrdvTm1kTLP9kONrNt/nD/76Gw9M0FEMlpSxxTM7Eoz2wK8E/RfZGb/ltJkGeRHq/eTl2N8\n/OKJYUcRETkryR5o/ifgBuAogLuvB65OVahM0h6N8eO1NXxk9jjGlBSGHUdE5KwkffaRu+/vMig6\nwFky0i/eOczRpjY+dalOQxWRzJfsrbP3m9lVgJtZPvBFYGvqYmWOH62uZmxpIdfM1BXMIpL5kt1S\nuAf4AjARqAHmBf1Z7WiklWXbDvPJSyaSl6tLPkQk8yW7pWDufmdfZmxmjwO3AIfd/YIzjDfgEeBm\noBm4293X9uUzwvbCxlqiMecTOsAsIkNEsn/evmFmL5vZ582sPMlpngBu7GH8TcCM4LUI+HaS800b\nS9YfYMa4EmZV9HYjWRGRzJBUUXD3mcCfAXOBtWa21Mw+08s0rwH1PTS5FXjS41YA5WY2PsncoTtw\nvIVVe47xsYsm6NoEERky+nL20Up3vw+YT3xl/71eJunNRCDxjKbqYFhGWLrhAAC/dtGEkJOIiAwc\nc+/21kanG5mVAZ8g/pS184Bngf929zW9TFcJLO3mmMJS4Ovu/nrQ/3Pgfnd/3601zGwR8V1MVFRU\nVC1evLjXzGcSiUQoKSnp17Rd/eWvWsgx+Msrh531vAYy10BL12zK1TfK1TdDMdfChQvXuPulvTZ0\n915fwG7gm8CVybRPmK4S2NTNuMeAOxL6twHje5tnVVWV99eyZcv6PW2inYcbfcr9S/07r+0ckPkN\nVK5USNdsytU3ytU3QzEXsNqTWG8ne/bRNHd3MxvIhw4vAe41s8XA5cAJd68dwPmnzPPrazGDWy7U\nriMRGVqSPaZwRV/vfWRmPwTeBGaZWXVw5tI9ZnZP0OQFYBewg/hznn+vX/+CQebuLFlfw/zKUZwz\noijsOCIiAyrZLYVT9z5aAvF7H5lZj/c+cvc7ehnvZOAFcNsPRdh5pInPfXBq2FFERAac7n3URy9v\nPogZXD+3IuwoIiIDTvc+6qOXtxzi4knljCvVriMRGXrO5t5HGXEMYCAdON7CxpoTXD/3nLCjiIik\nRFJbCu5eB7zn3kdm9iXixxqyxs+2HALg+jnadSQiQ9PZ3NrzvgFLkSFe3nKQ6eNKmDY2/S5qEREZ\nCGdTFLLqhj8nmttZsateWwkiMqSdTVHo/f4YQ8gvth0iGnMdTxCRIa3HYwpm1siZV/4GnP1NfzLI\ny5sPUVFWyIUTR4QdRUQkZXosCu6uBwUAJ9ujvLr9CJ+8ZCI5OVm110xEsoyeIZmEt3bX09wW5SOz\ndTxBRIY2FYUkvLrtCIV5OVw5bXTYUUREUkpFIQmvbj/M5dNGU5SfG3YUEZGUUlHoRfWxZnYeaeLq\nGWPCjiIiknIqCr14bXsdAAtmjQ05iYhI6qko9OLV7YeZWD6M83QVs4hkARWFHrRHY7yx4yhXzxyL\nmU5FFZGhT0WhB2/vO06ktYNrZmrXkYhkBxWFHry6/TB5OcZV03UqqohkBxWFHry6/QiXTB5JWVF+\n2FFERAaFikI36iKtbKpp4BqddSQiWURFoRsrdh0F4IPTdX2CiGQPFYVurNh1lOEFuVwwoSzsKCIi\ng0ZFoRsrdtVz2dRR5OVqEYlI9tAa7wzqIq3sOBzh8qk660hEsouKwhm8tasegCumjQo5iYjI4FJR\nOIPO4wl6ypqIZBkVhTN4a/dRqipHka/jCSKSZbTW66Iu0sr2QxHtOhKRrKSi0MXK3aeOJ+ggs4hk\nHxWFLlbsOkpxQS4f0PEEEclCKgpdvLWrnqopI3U8QUSyktZ8Ceqb2th2qFG7jkQka6koJFiz9xgA\n86fqILOIZCcVhQRv7ztGXo7peIKIZC0VhQRv7zvOnAllFOXnhh1FRCQUKS0KZnajmW0zsx1m9sAZ\nxi8wsxNmti54/UUq8/QkGnPWVx9n3qTysCKIiIQuL1UzNrNc4F+B64BqYJWZLXH3LV2a/tLdb0lV\njmRtP9RIc1uUiyerKIhI9krllsJ8YIe773L3NmAxcGsKP++svL3vOAAXTxoZchIRkfCksihMBPYn\n9FcHw7q6ysw2mNmLZjY3hXl6tG7/MUYW5zNldHFYEUREQmfunpoZm90G3Ojuvx30fxa43N3vTWhT\nBsTcPWJmNwOPuPuMM8xrEbAIoKKiomrx4sX9yhSJRCgpKTnjuC+/3szYYTn8YVVRv+Z9NnrKFbZ0\nzaZcfaNcfTMUcy1cuHCNu1/aa0N3T8kLuBJ4KaH/QeDBXqbZA4zpqU1VVZX317Jly844/Hhzm0+5\nf6n/8yvb+z3vs9FdrnSQrtmUq2+Uq2+GYi5gtSex7k7l7qNVwAwzm2pmBcDtwJLEBmZ2jplZ0D2f\n+O6soynMdEYbqoPjCZN1PEFEslvKzj5y9w4zuxd4CcgFHnf3zWZ2TzD+UeA24HfNrANoAW4PKtqg\nenvfcczgwkm6aE1EslvKigKAu78AvNBl2KMJ3d8CvpXKDMl4e98xpo8toawoP+woIiKhyvormt2d\nt/cf1/UJIiKoKLDnaDPHm9t1PEFEBBUFNtacANBN8EREUFFga20D+bnGzIrSsKOIiIQu64vClgMN\nTB9XSkFe1i8KEREVhS21DcwZXxZ2DBGRtJDVReFIYytHGluZPV67jkREIMuLwtbaBgDmTNCWgogI\nZHlR2HKqKGj3kYgIkO1F4UADE0YUUV5cEHYUEZG0kNVFYWttg3YdiYgkyNqicLI9ys4jEe06EhFJ\nkLVFYdvBRmKug8wiIomytiicOsg8W1sKIiKdsrYobK1toKQwj0kj9UxmEZFTsrYobDnQwOzxpeTk\nWNhRRETSRlYWhVjM42ceadeRiMh7ZGVR2H+smaa2qI4niIh0kZVFYWttIwDnqyiIiLxHVhaF3XVN\nAEwbOzzkJCIi6SUri8KeuibGlBRSVpQfdhQRkbSSlUVhd10T08ZoK0FEpKusLAq76pqYqqIgIvI+\nWVcUGk62UxdpZaqOJ4iIvE/WFYU9wUFmbSmIiLxf1hWFzjOPVBRERN4n64rCriNNmMHk0brnkYhI\nV1lXFHbXNXHuyGEU5uWGHUVEJO1kZVGYOqYk7BgiImkpq4qCu+saBRGRHmRVUTjR5kRaO3TmkYhI\nN7KqKBxqckCno4qIdCerisLBphigoiAi0p3sKgrNTkFeDhPKh4UdRUQkLWVXUWiKMWVUMbl6BKeI\nyBmltCiY2Y1mts3MdpjZA2cYb2b2z8H4DWZ2SSrzHGqKadeRiEgPUlYUzCwX+FfgJmAOcIeZzenS\n7CZgRvBaBHw7VXmiMedQs+tGeCIiPUjllsJ8YIe773L3NmAxcGuXNrcCT3rcCqDczManIkzNsRai\nrnseiYj0JJVFYSKwP6G/OhjW1zYDYlddBEBXM4uI9CAv7ADJMLNFxHcvUVFRwfLly/s8j3ePRblw\nlHNw+3qW702vA82RSKRf/6bBkK7ZlKtvlKtvsjqXu6fkBVwJvJTQ/yDwYJc2jwF3JPRvA8b3NN+q\nqirvr2XLlvV72lRK11zu6ZtNufpGufpmKOYCVnsS6+5U7j5aBcwws6lmVgDcDizp0mYJ8JvBWUhX\nACfcvTaFmUREpAcp233k7h1mdi/wEpALPO7um83snmD8o8ALwM3ADqAZ+Fyq8oiISO9SekzB3V8g\nvuJPHPZoQrcDX0hlBhERSV5WXdEsIiI9U1EQEZFOKgoiItJJRUFERDqpKIiISCeLnwCUOczsCLC3\nn5OPAeoGMM5ASddckL7ZlKtvlKtvhmKuKe4+trdGGVcUzoaZrXb3S8PO0VW65oL0zaZcfaNcfZPN\nubT7SEREOqkoiIhIp2wrCv8edoBupGsuSN9sytU3ytU3WZsrq44piIhIz7JtS0FERHqQNUXBzG40\ns21mtsPMHggxxyQzW2ZmW8xss5l9MRj+kJnVmNm64HVzCNn2mNnG4PNXB8NGmdnPzOzd4H3kIGea\nlbBM1plZg5l9KYzlZWaPm9lhM9uUMKzb5WNmDwa/t21mdsMg5/p7M3vHzDaY2bNmVh4MrzSzloTl\n9mj3c05Jrm6/t5CX138lZNpjZuuC4YO5vLpbNwzubyyZhy5k+ov4rbt3AtOAAmA9MCekLOOBS4Lu\nUmA7MAd4CPjjkJfTHmBMl2F/BzwQdD8A/G3I3+NBYEoYywu4GrgE2NTb8gm+0/VAITA1+P3lDmKu\n64G8oPtvE3JVJrYLYXmd8XsLe3l1Gf+PwF+EsLy6WzcM6m8sW7YU5gM73H2Xu7cBi4Fbwwji7rXu\nvjbobgS2kqLnUg+QW4HvBd3fAz4eYpaPADvdvb8XL54Vd38NqO8yuLvlcyuw2N1b3X038WeGzB+s\nXO7+srt3BL0rgHNT8dl9zdWDUJfXKWZmwKeAH6bis3vSw7phUH9j2VIUJgL7E/qrSYMVsZlVAhcD\nbwWDfj/Y3H98sHfTBBx4xczWBM/FBqjw00/DOwhUhJDrlNt573/WsJcXdL980uk391vAiwn9U4Nd\nIa+a2f8JIc+Zvrd0WV7/Bzjk7u8mDBv05dVl3TCov7FsKQppx8xKgGeAL7l7A/Bt4ru35gG1xDdh\nB9uH3H0ecBPwBTO7OnGkx7dZQzldzeKPdP0Y8KNgUDosr/cIc/l0x8y+AnQAPwgG1QKTg+/5PuBp\nMysbxEhp9711cQfv/cNj0JfXGdYNnQbjN5YtRaEGmJTQf24wLBRmlk/8S/+Bu/8YwN0PuXvU3WPA\nd0jRpnNP3L0meD8MPBtkOGRm44Pc44HDg50rcBOw1t0PBRlDX16B7pZP6L85M7sbuAW4M1iZEOxq\nOBp0ryG+H3rmYGXq4XtLh+WVB3wS+K9TwwZ7eZ1p3cAg/8aypSisAmaY2dTgL87bgSVhBAn2Wf4n\nsNXdv5Hrr1nCAAAC4UlEQVQwfHxCs08Am7pOm+Jcw82s9FQ38QOVm4gvp7uCZncBzw1mrgTv+Qsu\n7OWVoLvlswS43cwKzWwqMANYOVihzOxG4E+Bj7l7c8LwsWaWG3RPC3LtGsRc3X1voS6vwLXAO+5e\nfWrAYC6v7tYNDPZvbDCOqqfDC7iZ+NH8ncBXQszxIeKbfxuAdcHrZuD7wMZg+BJg/CDnmkb8TIb1\nwOZTywgYDfwceBd4BRgVwjIbDhwFRiQMG/TlRbwo1QLtxPfffr6n5QN8Jfi9bQNuGuRcO4jvbz71\nG3s0aPt/g+93HbAW+LVBztXt9xbm8gqGPwHc06XtYC6v7tYNg/ob0xXNIiLSKVt2H4mISBJUFERE\npJOKgoiIdFJREBGRTioKIiLSSUVBREQ6qSjIkGZmoxNue3ywy22bf5WCz7vbzI6Y2X/0cbqvmtm1\nvbT5jeA2yUvPLqVI93SdgmQNM3sIiLj7P6TwM+4GLnX3e1M0/wXEbz19SyrmL6ItBclaZhYJ3hcE\nd8B8zsx2mdnXzexOM1tp8YcOnRe0G2tmz5jZquD1wSQ+424z+9/g4Sh7zOxeM7vPzN42sxVmNipo\n94SZ3RZ07zGzh81sbfD556dyOYgkUlEQibsIuAeYDXwWmOnu84H/AH4/aPMI8E13v4z47Q+S3UV0\nAfEbrV0G/DXQ7O4XA28Cv9nNNHXufgnxu4r+cd//OSL9kxd2AJE0scqDe9ab2U7g5WD4RmBh0H0t\nMCd+3zIAysysxN0jvcx7mccfmtJoZieA5xPmfWE305y6Q+Ya4gVFZFCoKIjEtSZ0xxL6Y5z+f5ID\nXOHuJ1Mw7+6mifbQRmTAafeRSPJe5vSuJMxsXohZRFJCRUEkeX8AXBo8SnIL8WMQIkOKTkkVGUA6\nJVUynbYURAZWC3BTXy9eS4aZ/Qbwb8CxgZ63yCnaUhARkU7aUhARkU4qCiIi0klFQUREOqkoiIhI\nJxUFERHp9P8Bsxasekr6TYkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106a6c710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "Cv  = 0.1     # Outlet valve constant [cubic meters/min/meter^1/2]\n",
    "A   = 1.0     # Tank area [meter^2]\n",
    "\n",
    "# inlet flow rate in cubic meters/min\n",
    "def qin(t):\n",
    "    return 0.15\n",
    "\n",
    "def deriv(h,t):\n",
    "    return qin(t)/A - Cv*np.sqrt(h)/A\n",
    "\n",
    "IC = [0.0]\n",
    "t = np.linspace(0,200,101)\n",
    "h = odeint(deriv,IC,t)\n",
    "\n",
    "plt.plot(t,h)\n",
    "\n",
    "plt.xlabel('Time [min]')\n",
    "plt.ylabel('Level [m]')\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear Approximation of the Step Response\n",
    "\n",
    "The step response of the gravity drained to a change in flowrate looks similar to the step response of a firat order linear system. Let's try a linear approximation\n",
    "\n",
    "$$\\tau\\frac{dx}{dt} + x = Ku$$\n",
    "\n",
    "which has a step response solution that can be written\n",
    "\n",
    "$$x(t) = x_{ss} + (x_0 - x_{ss})\\exp(-t/\\tau)$$\n",
    "\n",
    "where $x_{ss} = Ku_{ss}$. There are two parameters, $K$ and $\\tau$, which we need to estimate in order to fit the linear approximation to the nonlinear simulation results computed above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Fit for $K$\n",
    "\n",
    "The steady state gain $K$ of the linear system is given by\n",
    "\n",
    "$$ K = \\frac{x_{ss} - x(0)}{u_{ss} - u_0}$$\n",
    "\n",
    "where $u_0$ is the initial input, $u_{ss}$ is the steady-state input, and $x_0$ and $x_{ss}$ are corresponding values of the state variable. In the case of liquid level, $h\\sim x$ and $q_{in}\\sim u$, therefore an estimate of $K$ is \n",
    "\n",
    "$$ K = \\frac{h_{ss} - h_0}{q_{in,ss} - q_{in,0}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Steady-State Gain is approximately =  [ 14.98595148]\n"
     ]
    }
   ],
   "source": [
    "q0,h0 = 0,0            # initial conditions\n",
    "qss = qin(t[-1])       # final input\n",
    "hss = h[-1]            # python way to get the last element in a list\n",
    "\n",
    "K = (hss-h0)/(qss-q0)  # step change in output divided by step change in input\n",
    "print('Steady-State Gain is approximately = ', K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Fit for $\\tau$\n",
    "\n",
    "From the formula for the solution of a first-order linear equation with constant input,\n",
    "\n",
    "$$\\frac{x_{ss} - x(t)}{x_{ss} - x_0} = \\exp(-t/\\tau) \\qquad \\implies \\qquad \\tau = \\frac{-t}{\\ln\\frac{x_{ss} - x(t)}{x_{ss} - x_0}}$$\n",
    "\n",
    "We pick one point representative of the transient portion of the nonlinear response. In this case the response at $t = 25$ minutes accounts for $\\approx$60% of the ultimate response, so we choose point as a representative point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated time constant is  [ 23.92271579]\n"
     ]
    }
   ],
   "source": [
    "k = sum(t<25)        # find index in t corresponding to 25 minutes\n",
    "tk = t[k]\n",
    "hk = h[k]\n",
    "\n",
    "tau = -tk/np.log((hss-hk)/(hss-h0))\n",
    "print('Estimated time constant is ', tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Comparing the linear approximation to the nonlinear simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8FHX6wPHPs5tNIQmhhV6CFBWkSQQr4FkAC6ingBX0\nPMRez1Pvdx7e6RXPk7Ocop6IBbEX9FCxUCy0ACGh1wAhoQWSkF72+/tjJusSUxbIZpLs835lXpm2\nM8/Ozs6z8/3OfEeMMSillFIALqcDUEop1XBoUlBKKeWjSUEppZSPJgWllFI+mhSUUkr5aFJQSinl\no0nBQSKSICJGRMLs4S9EZKLTcdVGRB4Rkf8GadkLROTmY3xtVxHJExF3XcdVn0TkHBHZ6HQcjY2I\nTBeRP9bTuprsZ6RJoRYikiYi+0Qk2m/czSKyoK7XZYwZbYx5va6XeyxEZKyIJItIrogcEJHvRKQ7\ngDHmr8aYYzpw13GMaSJyfsWwMWanMSbGGFPuZFyBqhx/BWPM98aYE52IqToiMsn+ATPe6ViqY4yZ\nYoz5SzCWbb/3nn7ranCfUV3RpBAYN3C300EEQ8VZSqVxPYE3gPuBOKA78B+gURxs1bGpal/wMxE4\nCNzg0PpVPdGkEJh/Ag+ISIuqJorImSKyXERy7P9n+k1bICJ/EZEfReSwiMwTkTbVLMdXdGL/MvtB\nRJ4SkUMisl1ERvvNGycir4pIpojsFpHHK4pNRKSH/cs+y/6VP8s/dvsX6u9FJAXIr+LLOBDYboz5\n1lgOG2M+NMbstF8/VUTesvsrisBuFJFddqxTROQ0EUkRkWwRed5v3b7XVnp9Vcmp2vchIm8CXYHP\n7CKjB6sojusoInNE5KCIbBGR31aK4z0RecP+XNaKSGI1n8uLIvJUpXGfish9dv/v7c/gsIhsFJHz\nqlpOoERkhIik+w2nicgD9vbMEZF3RSTSb/ol9lldtoj8JCL9/aY9JCJb7djWicjlftMm2fvlNBHJ\nAqZWE083YDgwGRgpIu0rxypWkeIBO9Zr/abPFKtY52s7hoX28iqmGxG5XUQ2A5vtcVV+n0Sklb2u\nS+3hGPtzvcFvXY9XiutBsc70M0XkMhG5SEQ22fvEI35xDBGRxfY2zBSR50Uk3J62yJ5ttb2vja/i\nMzpZrO9vtr0vjam0Df4jIv+zt8FSEelR407gJGOMdjV0QBpwPvAR8Lg97mZggd3fCjgEXA+EAVfb\nw63t6QuArUBvIMoe/rs9LQEwQJjfvDfb/ZOAUuC3WGcqtwIZgNjTPwZeAqKBtsAy4BZ7Wk/gAiAC\niAcWAf+u9J6SgS5AVBXv+QSgCJgGnAvEVJo+FXir0nuYDkQCF9qv/cSOqxOwDxhe+bUBbINA3sf5\nNSxrEfCCHddAYD/wK784ioCL7O37N2BJNfvAMGCX37ZvCRQCHYET7Wkd/WLocTT7VhXjRwDpleZb\nZq+vFbAemGJPG2Rv36H2+5hozx9hT7/Kfp0LGA/kAx389rEy4E6sffcX+4I93x+BZXZ/KnB/pVjL\ngKftz2m4vY4T7ekzgcP2NowAngF+8Hu9Ab6231cUtX+fLgT2YO1brwAf+C1rJj9/RyviehTwYH2P\n9gNvA7FAX/sz7G7PPxg43V5ngr2N76kUZ8+qPiN7+VuAR4Bw4Ff2e/bfBlnAEHv5s4B3nD62Vbtf\nOh1AQ+/4OSmcAuRgHZz8k8L1FV8Yv9csBibZ/QuA//Obdhvwpd2fQM1JYYvf65rZ87YH2gHF+H2J\n7S/P/Grew2XAqkrv6aZa3vfpwHv2F6nI3rFj7GlT+WVS6OT32ixgvN/whxVfMI4iKQT4PqpMClgJ\nrxyI9Zv+N2CmXxzf+E3rAxRWs14BdgLD7OHfAt/Z/T2xDsrnA55j2beqGD+CXyaF6/yGnwSm2/0v\nAn+p9PqN2Em4imUnA2P99rGdAcS52e/zexhYXSnWMiDab9x7wB/t/pn4HQCBGPtz6WIPG+xEHcj3\nyR5+Dis57cZOFn7r8k8KhYDbHo611zXUb/4VwGXVvOd7gI/9hmtKCudgJSqX3/TZwFS/uP7rN+0i\nYMPR7Cv12WnxUYCMMWuAz4GHKk3qCOyoNG4H1i/kCnv8+guwvhiB8L3OGFNg98YA3bB+nWTap6vZ\nWGcNbQFEpJ2IvGMXaeQCbwGVi6x21bRiY8wSY8w4Y0w81k4/DPhDDS/Z69dfWMVwoO/ZJ8D3UZ2O\nwEFjzGG/cbV9LpFVFWMZ65v8DlbiBbgG69cexpgtWAeQqcA+O96OAcZ4NKrbh7oB91fsB/a+0AXr\n/SMiN/gVLWVj/bjx34Y17gcichZWndI79qi3gX4iMtBvtkPGmHy/4R0V66+8DmNMHlbdRJXTCez7\n9LL9PmYaY7JqCD/L/HzRQaH9v8r9UkR6i8jnIrLH3tf+ytHta7uMMd4aYj7WY0C906RwdP6E9SvR\n/8POwPpi+uuK9SsmWHZhnSm0Mca0sLvmxpi+9vS/Yv2y6WeMaQ5ch/Vr158JdGXGmOVYxWenHH/o\n5GOd9VRoX92M1P4+anoPGUArEYn1G3c8n8ts4Eq7PHwo1tmPFYQxbxtjzsbaDwzwj2Ncx7HYBTzh\ntx+0MMY0M8bMtmN9BbgD6xd1C2ANgW9DsIqjBEgWkT3AUr/xFVqK39V5WNs5w2+4S0WPiMRgFRH5\nT/ePocbvk1j1Zi9jXQhxm/hdEXScXgQ2AL3sfe0RfvmdqU4G0EVE/I+nwT4GBI0mhaNg/yp8F7jL\nb/RcoLeIXCMiYWJdstcH66wiWHFkAvOAf4lIcxFxiVUpO9yeJRbIA3JEpBPwu6NZvoicLSK/FZGK\nM4+TgDHAkjoIPxkYJtY9BXFYxRHVqe197MWq//gFY8wu4CfgbyISaVe+/gbrbOOoGWNWAQeA/wJf\nGWOyAUTkRBH5lYhEYBWzFQLe6pf0Cx47voruaK/AeQWYIiJDxRItIhfbyTAa64C73471Ro4isYtV\nmT0Oq4J5oF93J3BNpVgfE5FwETkHuAR432/aRfY+FQ78BavuprozlNq+T4/Y7+kmrAtA3pC6uS8l\nFsgF8uz9/dZK06vd17ASZQHwoIh4RGQEcCk/n101KpoUjt6fsb5sANinr5dgXb6ZBTwIXGKMORDk\nOG7AqtRah1UR9wHQwZ72GHAqVh3I/7B+5R+NbKwkkCoiecCXWBXbTx5v0MaYr7ESawpWmW5NybO2\n9/E34P/sopEHqnj91Vj1DBlY8f/JGPPNcYT/Nlbdwdt+4yKAv2MljIoK0IcBRORaEVlbyzLnYiWS\nim7q0QRkjEnCOnt9Hms/2IJVV4AxZh3wL6wy+b1AP+DHo1j8ZXZMbxhj9lR0wAyseptR9nx77HVn\nYBWrTTHGbPBbzttYZ9kHsSp0r6vh/VT7fRKRwcB9wA12sdA/sBJE5SLdY/EAVrHgYaxE+26l6VOB\n1+19bVylmEuwksBorP3gBTvGDTRCFVdTKKXUUbN/Fb9ljOlczfSZWBWy/1efcaljp2cKSimlfDQp\nKKWU8tHiI6WUUj56pqCUUsqn0TVA1aZNG5OQkOB0GEop1aisWLHigH0zao0aXVJISEggKSnJ6TCU\nUqpREZHKd4pXSYuPlFJK+WhSUEop5aNJQSmllI8mBaWUUj6aFJRSSvloUlBKKeWjSUEppZRPo7tP\nQSnVgBkD3nLwloEpt/pNOXi9lYbLwXitzr/f+Pd7reUdMfxzZ7xevF4v5d5yvOVevMaL12vwesvx\n2tO8Xi/GVPw39nSv/ehJL177v/H69+ObXu1jKzHgm9/Yb/3n/1a/scLHYM3uxb9ZIWs+ux+v/agh\n47cca9j+Awwtep1B/7MvCepHqElBqcbEGCgvhdJ8KC38uSsr+vl/WRGUFf88XF5iDVf1/4iuzNdv\nvGWYshJMeQmmrNQaLi+1DvZ2J/bBX0w5Yqxh11E9X+j4COC2u1CxpDQXNCko1QR4y6EoBwoPWf8r\nuuJcKMq1/hfnWf9L8qz+knyrvyQfSgvs4Xzr1/QxKpcwysRDKR5KJcz6b9yU4KHEuH1dsXFTatyU\nEkYZkZThogw3ZYRRZir63Xjt8eW4KK/oN1a/ETe43IjLjbjCEJcLxP7vCsPlciFua7rL5f/fhbjc\nuF1uxCW43NZ4l7hw+eYXa5zLhYgLt8tl9wvidv887BJcYs3ncom1DJcgIrjt17r8XysuXG4XbheA\nWNMqxrvE91qXWOOxp7vEenKny+1CBASxX2fP73LhEgAXuAQX2NOw57PXJ4K4BMF+LfiGAU53Bf+Q\nrUlBqWNRWgh5+6wuf7/d7YOCg5B/AAqyoPCgNVyYDcU5tS7ShMfg9URTHhZNSVg0xa4oiqQlBWEd\nyHdHkBceweHycHK9YeSUhZNbGsahMjeHSsIoMB6KCafIhFOMhyLCKTYeSuz+EsIod3loFu4hOjyM\nqHA3zcLdNPP83B8V7ibS4ybK7iI9LiI9br/ORXSYm4gwa3xEmIvwMJfvv9Vvjfe4XbhdgT7iWDUk\nmhSUqqwoB7J3QU465KZDbsbP3eE9VlfdQT48Bpq19nXeVj0pDGvOYYkl28Rw0BvFgbJm7C0JJ7Mo\ngsxiD+mFHtLzXBzKrb7oJdztokUzD3FRVtc8ykPzyDCaR3k4IdJDbGQYsZEeYiLDiI0IIyYyjOjw\nMGIiwoiOcBMdEUZEmPVrWKmaaFJQocfrhdzdcHArZG2Fg9sgewcc2mH9L6p0wBc3xHaA5h0g/kQ4\nYTjEtMMb3ZYcdyv2eJuTXhzNjsIo0vMMmTmFZOYUsWd/EQfyivFW8ciSFs08tImJoE1MOF1bRjAo\nOpzW0RG0ivbQMjqcVs3CadEsnBbNPLRsFk6kRw/oqn5oUlBNl7fcOuDvW291BzbCgU1wYAuUFf48\nnzsCWnaDFt2gyxBo0RXiukBcF4qiO7CzOIa0g0XsyCpgx8F8du4pJP1gAenZhZSUebGe9X4YgNjI\nMDrERdI+LoqT2zenXVwk7ZpH0DY2kraxEbRtHkHr6AjCw/RqcNUwaVJQTUNpIexdCxmrYE8K7Em1\nEkFZkT2DWAf7Nr0hYRi06QmtekDrHhDbkUOFZWzae5jN+/LYuj+PrRvz2brvMBk5+/B/OGHzyDC6\ntY7mpA6xXNCnHZ1bRtGpZRSdWjSjY4tIYiM9jrx9peqKJgXV+BhjFfukL4P05Va3d93PV+VEtYL2\np8BpN0O7vhB/klXsEx5NSZmXzfsOsz7zMBvW5bJhz2427NnAgbxi3+KjPG56tI1mcLeWXNmmMyfE\nR5PQOppurZvRolm4Q29aqfqhSUE1fF4v7F8P2xfBjh9h5xLrah+AiObQ6VQ4+x7oOAg6DIS4ziBC\nabmXjXsOk5KeQ+rSraTuzmHTnjxKyq0K3UiPi97tYjn3xHh6t4ulV7sYerWLpUPzSFx65YwKUZoU\nVMOUmwlbv4Mt31jJoOCANb5FV+hxHnQ7A7oMhTYngssqn8/KKyZpxyFW/rSBVTuzSdmdTVGplQDi\nojz07xzHjWcn0LdjHH07NiehdbReNqlUJZoUVMPg9UJmMmz8AjZ9YdUJAMS0g57nQffh0P0cKynY\n9uUWsTglkyXbDrJsexZb9+cD4HELfTvGcc2Qbgzs2oIBnePo2qqZXr2jVAA0KSjneL2wayms+wTW\nzYHDGSAu6HI6nD8Vep4P7U7Bvt2T/OIylm7Yy6JNB/hhywG27MsDIDYijNO6t+LKwV04LaElp3SK\nI9ITSo0fKFV3NCmo+mUM7F0DKe9C6gdwONO6JLTXBXDSo9B7JDRr5Zt9+4F8vtuwjwUb97F020FK\nyr1EhLkY0r0VVw3uzBk9WtO3Y5wWAylVRzQpqPpRcNBKBCvfhH1rwRUGvS6EUx63EkFELGC1HJmy\nK5uv1u5h3rq9vrOBnm1jmHhmN4b3bktiQks9E1AqSDQpqOAxxrpSaPkrsP4zqwXOToPhoqfglF/7\nzggqEsHnKRnMTd3D7uxC3C5haPdWXDu0K+ef3I4urZo5/GaUCg2aFFTdKyuG1Pdh6XSrwjgiDgbf\nCKfeYN0/YEs7kM8nybv5NDmD7Qfy8biFs3u24Z7ze3FBn3Z6T4BSDtCkoOpOUQ4kzYAl0yFvD7Tt\nA5f8G/qPg/BowKosnpuayfsr0lm2/SAicHr31kwZfgKj+nYgrpneEayUkzQpqONXeAiWvGh1xblw\nwrlw+XQ4YYTvyqGNew7z1pIdfLxqN3nFZXRvE82Do07k8kGd6BAX5Wj4SqmfaVJQx64oFxb/B5a8\nYCWDk8fAOfdDx4EAlJV7+XrdHl77MY1laQcJD3NxSf8OXD2kK4ndWup9A0o1QEFLCiLSBXgDaIf1\niNGXjTHPVJpHgGeAi4ACYJIxZmWwYlJ1pKwEVrwGC5+07jQ++VIY/pCvviCvuIx3lu3ktR/T2J1d\nSJdWUfzhopO5cnBnWkZrPYFSDVkwzxTKgPuNMStFJBZYISJfG2PW+c0zGuhld0OBF+3/qqHa+CV8\n9bDVJHXCOXDBn622h4CD+SXM/HE7ry/eQU5hKUMSWvHHS/pwQZ92eh+BUo1E0JKCMSYTyLT7D4vI\neqAT4J8UxgJvGGMMsEREWohIB/u1qiHJ2gpfPgybv7LaG7r2A+uOYxEO5pfw0qKtvPHTDgpLyxnZ\ntx1ThvdgUNeWTketlDpK9VKnICIJwCBgaaVJnYBdfsPp9rgjkoKITAYmA3Tt2hVVj8pK4MdnYNGT\n1p3HFz4BQ28Bt4ecwlJeXrSVmT+mUVBaztgBHbnjVz3p2TbW6aiVUsco6ElBRGKAD4F7jDG5x7IM\nY8zLwMsAiYmJVTzcUAVF+gqYc6d1B3LfK2DU3yG2HcVl5bz5/Taen7+FnMJSLunfkbvP02SgVFMQ\n1KQgIh6shDDLGPNRFbPsBrr4DXe2xyknlZXAwr/DD9Mgpj1c/S6cOApjDHNTMvnr3PXszi5kWO94\nfj/qRPp2jHM6YqVUHQnm1UcCvAqsN8Y8Xc1sc4A7ROQdrArmHK1PcNj+jfDRbyFzNQy6Hkb+FSKb\nsz4zl8c+W8uSbQc5qX0sb/1mKGf3auN0tEqpOhbMM4WzgOuBVBFJtsc9AnQFMMZMB+ZiXY66BeuS\n1BuDGI+qiTGw6k2Y+yCEN4Pxs+DkS8gvLmPa5+uY8eN2mkd5ePyyU7h6SFe9mkipJiqYVx/9ANR4\n5LCvOro9WDGoAJUUwP/uh9VvW3chX/4yxLbjm3V7efTTNWTkFHHN0K48OPJEbY9IqSZO72gOdVlb\n4d3rYN966wa04Q9yqLCcqe+s4tPkDHq3i+HDa85gcLdWtS9LKdXoaVIIZdsWwHsTraedXfch9DyP\nr9ft5ZGPUzmUX8K95/fm1hE9CA9zOR2pUqqeaFIIVctegS9+D216wzXvUBjdhT9/lMLsZbs4qX0s\nM288Ta8qUioEaVIINV4vzPuD1Yhd71FwxSusOwh3vfYDW/fnMWV4D+67oLeeHSgVojQphJKyYvh4\nCqz9CIbeirnwcd5O2s1jc9bRopmHt34zlLN66mWmSoUyTQqhoigX3r0Wti+CC/5M0Wm384cP1/Lh\nynSG947n6XEDaB0T4XSUSimHaVIIBYXZ8NYV1g1pl7/Eri5juOXFxazfk8vd5/Xi7vN64dL7DpRS\naFJo+goOwpuXwd51MO4NlkeewS3/+ZGyci8zJp3GuSe2dTpCpVQDokmhKcvPgjfHWk1XTJjFh4f7\n8vCbS+ncMor/TkzkhPgYpyNUSjUwmhSaqqIceOtyOLAZM2E209K68uy3qzmzR2tevHYwcc08Tkeo\nlGqANCk0RSUF8PZ42LuW8nGz+L/UtsxetpmrBnfmr1f0w+PWy02VUlXTpNDUlJXAe9fDrqWUXPYK\nty1twzfrd3HHuT25/8LeWI3XKqVU1TQpNCXGwKe3wZZvKL7oGSYu7cTS7Xt5bExfJp6Z4HR0SqlG\nQJNCUzL/CUh9n+Jhf+DaFb1ZufMg08YN5LJBnZyOTCnVSGhSaCpWvQWL/klJ/2sZv+5M1mRk89zV\np3Jx/w5OR6aUakQ0KTQF2xbCZ3dTljCCqzPGszYzlxeuPZUL+7Z3OjKlVCOjSaGxO5QG70/E26oH\nNxfeSXJGviYEpdQx02sTG7OSAnj3Ooy3nIc8D7FwZzFPjxvASE0ISqljpEmhsTIGPrsbs2cNL7Z+\niPe2R/C3y/sxdqBWKiuljp0mhcZq2cuQ+h6LOt3Mk9sSeHj0SUwY0tXpqJRSjZwmhcYoIxm++gM7\n2wxj0tbhTDozgcnDTnA6KqVUE6BJobEpPgwf3ERRRCvG7r6WC/t24I+X9NE7lZVSdUKTQmMz90HM\noe1MzptCQpcuPDNhEG59FoJSqo7oJamNScr7sPptZrjHs7nZAD69fjCRHrfTUSmlmhBNCo1Fbgbm\nf/ex0dOHfxWO4d0bE2kbG+l0VEqpJkaLjxoDY2DOnZSWlnBL3s387cpB9Osc53RUSqkmSJNCY7Dq\nTdjyDY8Xj2fkOWfqvQhKqaDR4qOGLnsn3i8eZpnpy/pO43h75IlOR6SUasI0KTRkxlA+526KS8v5\ni/t2/nvtYH1qmlIqqPQI05Ct/Qj3tu/4Z9lV/G78BXSIi3I6IqVUE6dnCg1VUQ7Fn/+ejd7uRJ45\nhREntnU6IqVUCNAzhQaq4MuphBVl8WqLu7n3wpOdDkcpFSI0KTRAJn0FkcmvMct7IXdceyXhYfox\nKaXqhxYfNTReLwffv5MyE4frV3+gV7tYpyNSSoUQ/QnawBxc8iatc9byUavJXDOsn9PhKKVCTNCS\ngojMEJF9IrKmmukjRCRHRJLt7tFgxdJYmOLDyLePkWJ6cMl1d+PShu6UUvUsmGcKM4FRtczzvTFm\noN39OYixNAqbPnqCluVZ7Djtj3RpHeN0OEqpEBS0pGCMWQQcDNbym5pDGVvptvG/LIoYzkUXXeZ0\nOEqpEOV0ncKZIpIiIl+ISN/qZhKRySKSJCJJ+/fvr8/46s3Wd34PBjpf9Q99PoJSyjFOJoWVQFdj\nTH/gOeCT6mY0xrxsjEk0xiTGx8fXW4D1JWXFD5ya8w2rO13DCT31ngSllHMcSwrGmFxjTJ7dPxfw\niEgbp+JxSlm5l/wv/ky+RNF//J+cDkcpFeIcSwoi0l7sBwuLyBA7liyn4nHKV19+xhllS8nsO5mo\nuNZOh6OUCnFBu3lNRGYDI4A2IpIO/AnwABhjpgNXAreKSBlQCEwwxphgxdMQHcgrpvXyJ8lxtaDX\nmAecDkcppYKXFIwxV9cy/Xng+WCtvzH46INZTGYN+8/8MxKhdy4rpZzn9NVHIWt9Rg6J2/5DTng7\n4kdMcTocpZQCNCk4Zs4nsznVtYXw4Q9AWITT4SilFKBJwRE/bjnAsMzXyQ+PJ2roRKfDUUopH00K\n9czrNXw85yPOcK8jfNjdepaglGpQNCnUs89SMrjo0CyKw1viGXKT0+EopdQRarz6SESuCGAZRfbN\nZ6oWpeVe5nw5l1fdyXjP+iOERzsdklJKHaG2S1JfAT4FamqMZxigSSEAH65I59f571EaGYtn6G+d\nDkcppX6htqTwhTGmxjIOEXmrDuNpskrKvHz07fe8416ODL0PIuOcDkkppX6hxjoFY8x1tS0gkHkU\nfLgyndH5n4IrDBk62elwlFKqSgHd0SwibuBiIMH/NcaYp4MTVtNSUuZl5rer+cSzCOn3a4ht73RI\nSilVpUCbufgMKAJSAW/wwmmaPlyZzrC8uUR5CuH0250ORymlqhVoUuhsP/dAHaXSci8vfruRDyO+\nxnQ9G+mgm1Ep1XAFep/CFyJyYVAjaaLmpmbS//BC4r37kTPucDocpZSqUaBnCkuAj0XEBZRiXaJq\njDHNgxZZE2CM4aWF23gqah4mrgfSa6TTISmlVI0CPVN4GjgDaGaMaW6MidWEULsfthxA9qymT/lG\nZMhkcOkN5Eqphi3Qo9QuYE2oPQTneL20cBs3R87HeJrBgAlOh6OUUrUKtPhoG7BARL4AiitG6iWp\n1VuzO4fVW3byWvSPyClXQlQLp0NSSqlaBZoUtttduN2pWry8aBvjIhbjKS+ERG34TinVOASUFIwx\njwU7kKZkT04R/0vNYHHcfGgxEDqd6nRISikVkBrrFERkam0LCGSeUPP20h0MYiNtC7fBab9xOhyl\nlApYbWcKN4tIbg3TBZgATK2ziBq5kjIvby/bxXMtf4Sy5nDKr50OSSmlAhZI09mxAcyjbF+syaQk\nL4uh3kUweKI+M0Ep1ajUmBS0LuHovbF4B5Oar8RVUgKnXu90OEopdVT0bqo6tGZ3Dit2HOK6yB+h\nbV9or+0cKaUaF00KdejNxTvo48mkbW4qDLwGpKYH1imlVMOjSaGO5BaV8unq3TzYfiWIG/qPczok\npZQ6ajXWKYjIc0C1TVsYY+6q84gaqc9WZ1BSWsZZ+d9Arwshpq3TISml1FGr7eqjpHqJogl4Pymd\nq1ttwVOw1yo6UkqpRqi2q49e9x8WkWbGmILghtT4bNp7mORd2fyr609gWkHvUU6HpJRSxySgOgUR\nOUNE1gEb7OEBIvJCUCNrRN5P2kULVyEnZC2EfldCmDYPpZRqnAKtaP43MBLIAjDGrAaGBSuoxqS0\n3MtHK3dzd+dNSHkx9NMKZqVU4xXw1UfGmF2VRpXXcSyN0ncb9pGVX8Kl7sUQ1xU6JzodklJKHbOA\nH7IjImcCRkQ8IvIAsD6IcTUa7yel0zOmmNZ7f4JTrtB7E5RSjVqgSWEKcDvQCdgNDLSHQ1pWXjHz\nN+7j/s4bEG+ZNn6nlGr0An3Ijhhjrg1qJI3Q3NRMyr2GYSWLoE1vaN/P6ZCUUuq4BHqm8KOIzBOR\n34hIQM+VFJEZIrJPRNZUM11E5FkR2SIiKSLS6J5EM2d1Bqe3KaZZxhLrLEGLjpRSjVxAScEY0xv4\nP6AvsFJEPheR62p52Uygpgv2RwO97G4y8GIgsTQUGdmFLE87xG1tUxEM9L3C6ZCUUuq4Hc3VR8uM\nMfcBQ4CDwOu1zL/Inq86Y4E3jGUJ0EJEOgQaj9M+T8kAYGj+AqvYKL63swEppVQdCPTmteYiMlFE\nvgB+AjJ77Xh4AAAcjklEQVSxksPx6AT4X+aabo+rav2TRSRJRJL2799/nKutG3NWZ3Beh2Ii9q7U\nswSlVJMR6JnCaqwrjv5sjOltjPm9MWZFEOM6gjHmZWNMojEmMT4+vr5WW61t+/NYszuXm9uss0b0\nGetsQEopVUcCvfroBGOMEZFmdbju3UAXv+HO9rgG77PVmYjA4MIfIf5kaN3D6ZCUUqpOBHqmcHoQ\n2j6aA9xgX4V0OpBjjMk8zmUGnTGGOat3c14XF+G7l8DJlzgdklJK1ZlAzxQq2j6aA1bbRyJSY9tH\nIjIbGAG0EZF04E+Ax379dGAucBGwBSgAbjyG+Ovdpr15bN2fz+OJm2GfF07SpKCUajoCTQoYY3bJ\nkdfh19j2kTHm6lqmGxrhXdHz1u5BBE4t+N5q66jDAKdDUkqpOqNtHx2leev2ckancCJ2LLKKjvSG\nNaVUE3I8bR/dFqygGqqM7EJSd+dwY7stUF6iRUdKqSYnoOIjY8wB4Ii2j0TkHqy6hpDx9bq9AJxR\nshiatYGupzsckVJK1a2A72iuwn11FkUjMW/dHk6OjyBmx3dw4mhwuZ0OSSml6tTxJIWQKkzPKShl\nybaDTOq0G0oOw0kXOx2SUkrVueNJCqbOomgEvtu4l3Kv4Vz3anBHQHd9GqlSqumpsU5BRA5T9cFf\ngKigRNRAzVu7l3bNI4jPXADdz4HwaKdDUkqpOldjUjDGxNZXIA1ZUWk5Czft57d9vcj6rTD0FqdD\nUkqpoAj45rVQtnT7QQpKyrkkaq01oteFzgaklFJBcjx1CiFj4cb9RIS5OOHQj9ZjN1t1dzokpZQK\nCk0KAVi4aR/DEprh3vmTniUopZo0TQq1SD9UwNb9+VzVyr6LWZOCUqoJ06RQi0WbDgAwtGwFhMdC\n1zMcjkgppYJHk0ItFm7aR6e4SJqnL4AeIyAs3OmQlFIqaDQp1KC03MuPW7K4quthJHe3Fh0ppZo8\nTQo1WLUzm7ziMkZGbbBG9PiVswEppVSQaVKowcJN+whzCT3zlkPrXhDX2emQlFIqqDQp1GDhpv0M\n6RKNZ9dPcMIIp8NRSqmg06RQjQN5xazZnctV7fdAaQH0ONfpkJRSKug0KVRjybYsAM4gBcQNCWc7\nHJFSSgWfJoVqLNmWRXS4m3b7F0OnwRAZ53RISikVdJoUqrFk20GGd/Ugmau06EgpFTI0KVThQF4x\nW/blMSZuKxgvnKBJQSkVGjQpVGHptoMADC5PhvAY6JzocERKKVU/NClUoaI+oc2+xVYFs9vjdEhK\nKVUvNClUYen2LC7sXIoc3Kb3JyilQoomhUoO5BWzaW8elzTfYo3oPszZgJRSqh5pUqhk2XarPmFg\n+RqIagnxJzsckVJK1R9NCpUs2ZZFs3A3rfYvh25ngUs3kVIqdOgRr5Kl2w5yQadSJHuH3sWslAo5\nmhT8HMwvYePew1wSt90a0e0sZwNSSql6pknBz4odhwC7PiEyDtr1dTgipZSqX5oU/KzaeYgwl9A6\nKwm6ngkut9MhKaVUvdKk4GfVzmzObleK6+BWrU9QSoUkTQq2cq9hdXo2FzffZo1I0PoEpVToCWpS\nEJFRIrJRRLaIyENVTB8hIjkikmx3jwYznpps2nuYgpJyEmUdRDSH9v2dCkUppRwTFqwFi4gb+A9w\nAZAOLBeROcaYdZVm/d4Yc0mw4gjUqp3ZAHTKWQldT9f6BKVUSArmmcIQYIsxZpsxpgR4BxgbxPUd\nl+Rdh+gRlU/4oS16KapSKmQFMyl0Anb5Dafb4yo7U0RSROQLEanyGlARmSwiSSKStH///mDEyqqd\n2VzeJt0a0KSglApRTlc0rwS6GmP6A88Bn1Q1kzHmZWNMojEmMT4+vs6DyCksZfO+PM4M3wruCOgw\noM7XoZRSjUEwk8JuoIvfcGd7nI8xJtcYk2f3zwU8ItImiDFVKSXdqk/oWbwWOp0KYeH1HYJSSjUI\nwUwKy4FeItJdRMKBCcAc/xlEpL2IiN0/xI4nK4gxVWnVzmwipITYQ2uhy5D6Xr1SSjUYQbv6yBhT\nJiJ3AF8BbmCGMWatiEyxp08HrgRuFZEyoBCYYIwxwYqpOqt2HmJUy71IQQl0GVrfq1dKqQYjaEkB\nfEVCcyuNm+7X/zzwfDBjqI0xhlW7snm8bRoUAJ31TEEpFbqcrmh2XFpWAdkFpQwwm6DVCRBT9xXZ\nSinVWIR8UkjdnQMY2uemaNGRUirkhXxSWJ+ZSw/3fjyF+7WSWSkV8kI+KazLyGVU3E5rQM8UlFIh\nTpNCZi5nRWy1GsGLP8npcJRSylEhnRT2Hy5m/+FiTixdB50TtRE8pVTIC+mksD4zlxgKaJW3RYuO\nlFKKEE8K6zJzGeDaimCg82lOh6OUUo4L7aSQkcs5UTusgU6DnQ1GKaUagJBOCuszcxkSsQNa94So\nFk6Ho5RSjgtqMxcNWVFpOVv359E7dhN0PNfpcJT6hdLSUtLT0ykqKnI6FNWIREZG0rlzZzwezzG9\nPmSTwsY9h4k3B4kp2a9FR6pBSk9PJzY2loSEBOzGhJWqkTGGrKws0tPT6d69+zEtI2SLjyoqmQHr\nGQpKNTBFRUW0bt1aE4IKmIjQunXr4zq7DNmksD4zl0RPGsYVBu37OR2OUlXShKCO1vHuMyGbFNZl\n5HJ6RBrStg94opwORymlGoSQTAper2F9Zg69yrdo0ZFSNRAR7r//ft/wU089xdSpU49pWWlpaZxy\nyikAJCUlcdddd9VFiAHbuHEjI0aMYODAgZx88slMnjy5zmOZOXMmd9xxR43zLFiwgJ9++sk3PH36\ndN544406WX9dCMmK5l2HCogv3U2U67BWMitVg4iICD766CMefvhh2rSpu8enJyYmkpiYWGfLq0p5\neTlu989N19x1113ce++9jB07FoDU1NR6i8XfggULiImJ4cwzzwRgypQp9bbuQIRkUlifeZgBYlcy\nd9QzBdXwPfbZWtZl5NbpMvt0bM6fLu1b4zxhYWFMnjyZadOm8cQTTxwxLS0tjZtuuokDBw4QHx/P\na6+9RteuXZk0aRLNmzcnKSmJPXv28OSTT3LllVce8doFCxbw1FNP8fnnnzN16lR27tzJtm3b2Llz\nJ/fcc4/vl/tbb73Fs88+S0lJCUOHDuWFF17A7XZz6623snz5cgoLC7nyyit57LHHAEhISGD8+PF8\n/fXXPPjgg0yYMMG3zszMTDp37uwb7tevX5WxbN++3RfLtGnTWLJkCV988QWdOnXis88+w+PxkJCQ\nQFJSEm3atCEpKYkHHniABQsWHPEeP/vsMx5//HFKSkpo3bo1s2bNorCwkOnTp+N2u3nrrbd47rnn\n+Pbbb4mJieGBBx4gOTmZKVOmUFBQQI8ePZgxYwYtW7ZkxIgRDB06lPnz55Odnc2rr77KOeecc3Qf\neIBCsvho+4F8Bri2YTzNtGVUpWpx++23M2vWLHJyco4Yf+eddzJx4kRSUlK49tprjyiCyczM5Icf\nfuDzzz/noYceqnUdGzZs4KuvvmLZsmU89thjlJaWsn79et59911+/PFHkpOTcbvdzJo1C4AnnniC\npKQkUlJSWLhwISkpKb5ltW7dmpUrVx6READuvfdefvWrXzF69GimTZtGdnZ2lbFs3bqV7777jjlz\n5nDddddx7rnnkpqaSlRUFP/73/8C3m5nn302S5YsYdWqVUyYMIEnn3yShIQEpkyZwr333ktycvIv\nDuw33HAD//jHP0hJSaFfv36+ZAdQVlbGsmXL+Pe//33E+LoWkmcKaQfyuSZsO9JhALhDchOoRqa2\nX/TB1Lx5c2644QaeffZZoqJ+vihj8eLFfPTRRwBcf/31PPjgg75pl112GS6Xiz59+rB3795a13Hx\nxRcTERFBREQEbdu2Ze/evXz77besWLGC006z2iUrLCykbdu2ALz33nu8/PLLlJWVkZmZybp16+jf\nvz8A48ePr3IdN954IyNHjuTLL7/k008/5aWXXmL16tW/mG/06NF4PB769etHeXk5o0aNAqwzi7S0\ntAC2mCU9PZ3x48eTmZlJSUlJrfcN5OTkkJ2dzfDhwwGYOHEiV111lW/6FVdcAcDgwYOPKo6jFZJn\nCjv353Ay27U+QakA3XPPPbz66qvk5+cHNH9ERISv3xhzVPO73W7KysowxjBx4kSSk5NJTk5m48aN\nvuKdp556im+//ZaUlBQuvvjiI67Lj46OrnY9HTt25KabbuLTTz8lLCyMNWvWVBuLy+XC4/H4LvF0\nuVyUlZUBVrGa1+sFqPaegDvvvJM77riD1NRUXnrppeO+M70irortEywhmRRcB9YTTgl0HOR0KEo1\nCq1atWLcuHG8+uqrvnFnnnkm77zzDgCzZs2q8zLu8847jw8++IB9+/YBcPDgQXbs2EFubi7R0dHE\nxcWxd+9evvjii4CW9+WXX1JaWgrAnj17yMrKolOnTscUW0JCAitWrADgww8/rHKenJwc3/Jff/11\n3/jY2FgOHz78i/nj4uJo2bIl33//PQBvvvmm76yhPoVcUsgtKqVz0SZrQJOCUgG7//77OXDggG/4\nueee47XXXqN///68+eabPPPMM3W6vj59+vD4449z4YUX0r9/fy644AIyMzMZMGAAgwYN4qSTTuKa\na67hrLPOCmh58+bN45RTTmHAgAGMHDmSf/7zn7Rv3/6YYvvTn/7E3XffTWJi4hFXOPmbOnUqV111\nFYMHDz7iyq1LL72Ujz/+mIEDB/oSQIXXX3+d3/3ud/Tv35/k5GQeffTRY4rveEggp3YNSWJioklK\nSjrm16ekZ5M8/WauifyJsEfSwRVyeVE1EuvXr+fkk092OgzVCFW174jICmNMrdfehtwRcfuBfPq6\n0iiJP0UTglJKVRJyR8Xt+3I5WXYS3mWg06EopVSDE3LXY+ZlbKSZFENHTQpKKVVZyJ0pRB6wL0Fr\n39/ZQJRSqgEKqaRgjKH14Q2USTjEn+h0OEop1eCEVFLYn1dML+82smN7gvvYHlWnlFJNWUglhe37\n8ujr2kFp/ClOh6JUoxATE/OLcU419XzgwAE8Hg/Tp0+v93X7mzNnDn//+9/rZFn//ve/KSgo8A1f\ndNFF1bbJVF9CKinsTd9KS8kjoou2jKrUsZoyZQo33HBD0JZvjPE1IeHv/fff5/TTT2f27Nl1tq5j\naS5izJgxATXyF4jKSWHu3Lm0aNGiTpZ9rELq6qOS9GQA4k7QNo9UI/PFQ7AntW6X2b4fjD76X7xT\np071NfVcXZPO5eXlPPTQQyxYsIDi4mJuv/12brnlFvLy8hg7diyHDh2itLSUxx9/nLFjx5KWlsbI\nkSMZOnQoK1asYO7cuXTr1u2I9c6ePZt//etfXHPNNaSnp/uawY6JieG3v/0t8+bNo3379rzzzjvE\nx8czYsQIBgwYwMKFCykrK2PGjBkMGTKEqVOnsnXrVrZt20bXrl157bXXuPXWW0lKSiIsLIynn36a\nc889l2nTppGamsqMGTNITU3l6quvZtmyZbz33nskJSXx/PPPM2nSJKKioli1ahX79u1jxowZvPHG\nGyxevJihQ4cyc+ZMgCqb+n722WfJyMjg3HPPpU2bNsyfP/+IJrmffvppZsyYAcDNN9/MPffcQ1pa\nGqNHj+bss8/mp59+olOnTnz66adHNFR4vELqTCFifyrluHC31+IjpepKVU06v/rqq8TFxbF8+XKW\nL1/OK6+8wvbt24mMjOTjjz9m5cqVzJ8/n/vvv9/XYN7mzZu57bbbWLt27S8Swq5du8jMzGTIkCGM\nGzeOd9991zctPz+fxMRE1q5dy/Dhw49oVrqgoIDk5GReeOEFbrrpJt/4devW8c033zB79mz+85//\nICKkpqYye/ZsJk6cSFFREXfffTdbtmzh448/5sYbb+Sll16iWbNmv3j/hw4dYvHixUybNo0xY8Zw\n7733snbtWlJTU0lOtn6IVtXU91133UXHjh2ZP38+8+fPP2KZK1as4LXXXmPp0qUsWbKEV155hVWr\nVvm20+23387atWtp0aJFtW0vHauQOlNonbeRvZ7OdAz/5QerVIN2DL/o60tVTTrPmzePlJQUPvjg\nA8BqHG7z5s107tyZRx55hEWLFuFyudi9e7evae1u3bpx+umnV7mOd999l3HjxgEwYcIEbrrpJt9j\nQl0ul6+57Ouuu84XD8DVV18NwLBhw8jNzfWV148ZM8b36/qHH37gzjvvBOCkk06iW7dubNq0if79\n+zNz5kz69+/PLbfcUm0bS5deeikiQr9+/WjXrp3v4T19+/YlLS2NgQMH1tjUd1V++OEHLr/8cl+L\nr1dccQXff/89Y8aMoXv37gwcOPAX27yuBDUpiMgo4BnADfzXGPP3StPFnn4RUABMMsasDEYs5V5D\n99ItHGyTSMdgrECpEFVVk87GGJ577jlGjhx5xLwzZ85k//79rFixwvcEs4ompWtq8nr27Nns2bPH\n95CdjIwMNm/eTK9evX4xb0VT15X7/YdrWpe/zZs3ExMTQ0ZGRrXz+De17d8EeEVT2xVNfS9fvpyW\nLVsyadKk42pGu3Iz44WFhce8rKoErfhIRNzAf4DRQB/gahHpU2m20UAvu5sMvBiseDIz0ukgBymN\n7xesVSilbCNHjuTFF1/0NVW9adMm8vPzycnJoW3btng8HubPn8+OHTtqXdamTZvIy8tj9+7dpKWl\nkZaWxsMPP+yrcPZ6vb4zkrfffpuzzz7b99qKYqYffviBuLg44uLifrH8c845x5dsNm3axM6dOznx\nxBPJycnhrrvuYtGiRWRlZfnWcbRqauq7uma0zznnHD755BMKCgrIz8/n448/DtrjNysL5pnCEGCL\nMWYbgIi8A4wF1vnNMxZ4w1iFiktEpIWIdDDGZNZ1MFlbltMZiOiqVx4pFaiCgoIjnmt83333BfS6\nm2++mbS0NE499VSMMcTHx/PJJ59w7bXXcumll9KvXz8SExM56aTaH4c7e/ZsLr/88iPG/frXv2b8\n+PE8+uijREdHs2zZMh5//HHatm17RH1DZGQkgwYNorS01FdpW9ltt93GrbfeSr9+/QgLC2PmzJlE\nRERw6623cvvtt9O7d29effVVzj33XIYNGxbQ+/fn39R3ly5djiiGmjx5MqNGjfLVLVQ49dRTmTRp\nEkOGDAGs7Tlo0KCgPnGtQtCazhaRK4FRxpib7eHrgaHGmDv85vkc+Lsx5gd7+Fvg98aYpErLmox1\nJkHXrl0HB/LrorINS7+iaOE0utw4k9bxx9aGulL1SZvODkxMTAx5eXm/GD9ixAieeuopEhNrbS26\nyTmeprMbRUWzMeZl4GWwnqdwLMs4aehIGDqy9hmVUiqEBTMp7Aa6+A13tscd7TxKKVWtqs4SABYs\nWFC/gTQRwbxPYTnQS0S6i0g4MAGYU2meOcANYjkdyAlGfYJSjVVjezKict7x7jNBO1MwxpSJyB3A\nV1iXpM4wxqwVkSn29OnAXKzLUbdgXZJ6Y7DiUaqxiYyMJCsri9atW//i0kqlqmKMISsri8jIyGNe\nRsg9o1mpxqK0tJT09PTjuqZdhZ7IyEg6d+6Mx3NkS9BNqqJZqVDk8Xjo3r2702GoEBNSbR8ppZSq\nmSYFpZRSPpoUlFJK+TS6imYR2Q8c/S3NljbAgToMp6401Lig4camcR0djevoNMW4uhlj4mubqdEl\nheMhIkmB1L7Xt4YaFzTc2DSuo6NxHZ1QjkuLj5RSSvloUlBKKeUTaknhZacDqEZDjQsabmwa19HR\nuI5OyMYVUnUKSimlahZqZwpKKaVqoElBKaWUT8gkBREZJSIbRWSLiDzkYBxdRGS+iKwTkbUicrc9\nfqqI7BaRZLu7yIHY0kQk1V5/kj2ulYh8LSKb7f8t6zmmE/22SbKI5IrIPU5sLxGZISL7RGSN37hq\nt4+IPGzvbxtFJGhPeKomrn+KyAYRSRGRj0WkhT0+QUQK/bbb9HqOq9rPzeHt9a5fTGkikmyPr8/t\nVd2xoX73MWNMk++wmu7eCpwAhAOrgT4OxdIBONXujwU2AX2AqcADDm+nNKBNpXFPAg/Z/Q8B/3D4\nc9wDdHNiewHDgFOBNbVtH/szXQ1EAN3t/c9dj3FdCITZ/f/wiyvBfz4HtleVn5vT26vS9H8Bjzqw\nvao7NtTrPhYqZwpDgC3GmG3GmBLgHWCsE4EYYzKNMSvt/sPAeqCTE7EEaCzwut3/OnCZg7GcB2w1\nxhzrHe3HxRizCDhYaXR122cs8I4xptgYsx3rmSFD6isuY8w8Y0yZPbgE66mG9aqa7VUdR7dXBbEe\nXDEOmB2MddekhmNDve5joZIUOgG7/IbTaQAHYhFJAAYBS+1Rd9qn+zPqu5jGZoBvRGSFiEy2x7Uz\nPz8Nbw/QzoG4KkzgyC+r09sLqt8+DWmfuwn4wm+4u10UslBEznEgnqo+t4ayvc4B9hpjNvuNq/ft\nVenYUK/7WKgkhQZHRGKAD4F7jDG5wItYxVsDgUysU9j6drYxZiAwGrhdRIb5TzTWOasj1zCL9UjX\nMcD79qiGsL2O4OT2qY6I/AEoA2bZozKBrvbnfB/wtog0r8eQGtznVsnVHPnDo963VxXHBp/62MdC\nJSnsBrr4DXe2xzlCRDxYH/osY8xHAMaYvcaYcmOMF3iFIJ0618QYs9v+vw/42I5hr4h0sOPuAOyr\n77hso4GVxpi9doyOby9bddvH8X1ORCYBlwDX2gcT7KKGLLt/BVY5dO/6iqmGz60hbK8w4Arg3Ypx\n9b29qjo2UM/7WKgkheVALxHpbv/inADMcSIQu8zyVWC9MeZpv/Ed/Ga7HFhT+bVBjitaRGIr+rEq\nKtdgbaeJ9mwTgU/rMy4/R/yCc3p7+alu+8wBJohIhIh0B3oBy+orKBEZBTwIjDHGFPiNjxcRt91/\ngh3XtnqMq7rPzdHtZTsf2GCMSa8YUZ/bq7pjA/W9j9VHrXpD6ICLsGrztwJ/cDCOs7FO/1KAZLu7\nCHgTSLXHzwE61HNcJ2BdybAaWFuxjYDWwLfAZuAboJUD2ywayALi/MbV+/bCSkqZQClW+e1vato+\nwB/s/W0jMLqe49qCVd5csY9Nt+f9tf35JgMrgUvrOa5qPzcnt5c9fiYwpdK89bm9qjs21Os+ps1c\nKKWU8gmV4iOllFIB0KSglFLKR5OCUkopH00KSimlfDQpKKWU8tGkoJRSykeTgmrSRKS1X7PHeyo1\n2/xTENY3SUT2i8h/j/J1fxaR82uZZ7zdTPLnxxelUtXT+xRUyBCRqUCeMeapIK5jEpBojLkjSMsf\ngdX09CXBWL5SeqagQpaI5Nn/R9gtYH4qIttE5O8icq2ILBProUM97PniReRDEVlud2cFsI5JIvKJ\n/XCUNBG5Q0TuE5FVIrJERFrZ880UkSvt/jQReUxEVtrrPymY20Epf5oUlLIMAKYAJwPXA72NMUOA\n/wJ32vM8A0wzxpyG1fxBoEVEp2A1tHYa8ARQYIwZBCwGbqjmNQeMMaditSr6wNG/HaWOTZjTASjV\nQCw3dpv1IrIVmGePTwXOtfvPB/pY7ZYB0FxEYowxebUse76xHppyWERygM/8lt2/mtdUtJC5Aiuh\nKFUvNCkoZSn26/f6DXv5+XviAk43xhQFYdnVvaa8hnmUqnNafKRU4Obxc1ESIjLQwViUCgpNCkoF\n7i4g0X6U5DqsOgilmhS9JFWpOqSXpKrGTs8UlKpbhcDoo715LRAiMh54AThU18tWqoKeKSillPLR\nMwWllFI+mhSUUkr5aFJQSinlo0lBKaWUz/8DRMvaPC/g6lsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106a6ccf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "u0 = q0\n",
    "uss = qss\n",
    "\n",
    "xss = K*(uss - u0)\n",
    "xpred = xss - xss*np.exp(-t/tau)\n",
    "\n",
    "plt.plot(t,h)\n",
    "plt.plot(t,xpred)\n",
    "plt.legend(['Nonlinear Simulation','Linear Approximation'])\n",
    "plt.xlabel('Time [min]')\n",
    "plt.ylabel('Level [m]')\n",
    "plt.title('Nonlinear Simulation vs. Linear Approximation');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [Linearization](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/03.00-Linearization.ipynb) | [Contents](toc.ipynb) | [Linear Approximation of a Process Model](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/03.02-Linear-Approximation-of-a-Process-Model.ipynb) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE30338/blob/master/notebooks/03.01-Step-Response-of-a-Gravity-Drained-Tank.ipynb\"><img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open in Google Colaboratory\"></a><p><a href=\"https://raw.githubusercontent.com/jckantor/CBE30338/master/notebooks/03.01-Step-Response-of-a-Gravity-Drained-Tank.ipynb\"><img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}