furata_nonlin.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:2444dd74fc16b0539ebda7d5b48015811f71cf6b3b1001331e4387d0d2516e6b"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Nonlinear Furata Pendulum Simulation\n",
      "\n",
      "**Jon Sowman <j.sowman@soton.ac.uk>**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from scipy import integrate\n",
      "from pprint import pprint"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Pendulum parameters\n",
      "mp = 0.024\n",
      "mr = 0.95\n",
      "Lp = 0.129\n",
      "Lr = 0.085\n",
      "Jp = 3.33e-5\n",
      "Jr = 5.72e-5\n",
      "Bp = 5e-4\n",
      "Br = 15e-5\n",
      "kt = 0.042\n",
      "km = 0.042\n",
      "Rm = 8.4\n",
      "g = 9.81\n",
      "\n",
      "# Total MoI\n",
      "Jt = Jp*mp*Lr*Lr + Jr*Jp + 0.25*Jr*mp*Lp*Lp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We have a nonlinear model that looks like $\\ddot{\\alpha} = f_1(\\alpha, \\theta, \\dot{\\alpha}, \\dot{\\theta})$ and $\\ddot{\\theta} = f_2(\\alpha, \\theta, \\dot{\\alpha}, \\dot{\\theta})$. The `nlderiv` function, given a state vector $\\mathbf{x} = [\\theta, \\alpha, \\dot{\\theta}, \\dot{\\alpha}]'$ and an input $\\mathbf{u}$, returns the derivative $\\dot{\\mathbf{x}}$, i.e. $\\mathbf{\\dot{x}} = f(\\mathbf{x}, \\mathbf{u}, t)$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def nlderiv(x, t, u):\n",
      "    # Extract information from the state vector\n",
      "    t = x[0]\n",
      "    a = x[1]\n",
      "    tdot = x[2]\n",
      "    adot = x[3]\n",
      "    \n",
      "    # Actuator model (u given to this function is motor voltage Vm)\n",
      "    u = kt*(u - km*tdot)/Rm\n",
      "    \n",
      "    # Compute theta-dot-dot (tdotdot)\n",
      "    num1 = Br*tdot \\\n",
      "            + (Lp*Lr*mp*adot**2*np.sin(a))/2 \\\n",
      "            - (Lp*Lr*mp*np.cos(a)*((mp*np.cos(a)*np.sin(a)*Lp**2*tdot**2)/4 + (g*mp*np.sin(a)*Lp)/2 - Bp*adot))/(2*((mp*Lp**2)/4 + Jp)) \\\n",
      "            + (Lp**2*adot*mp*tdot*np.cos(a)*np.sin(a))/2\n",
      "    den11 = (Lp**2*Lr**2*mp**2*(np.cos(a))**2) / ((mp*Lp**2 + 4*Jp)*((-mp*Lp**2*(np.cos(a))**2)/4 + (mp*Lp**2)/4 + mp*Lr**2 + Jr)) - 1\n",
      "    den12 = (-mp*Lp**2*(np.cos(a))**2)/4 + (mp*Lp**2)/4 + mp*Lr**2 + Jr\n",
      "    \n",
      "    # Torque term\n",
      "    den1t1 = (Lp**2*Lr**2*mp**2*(np.cos(a))**2)/(4*((mp*Lp**2)/4 + Jp)*((-mp*Lp**2*(np.cos(a))**2)/4 + (mp*Lp**2)/4 + mp*Lr**2 + Jr)) - 1\n",
      "    den1t2 = -(mp*Lp**2*(np.cos(a))**2)/4 + (mp*Lp**2)/4 + mp*Lr**2 + Jr\n",
      "    tauterm = 1/(den1t1 * den1t2)\n",
      "    \n",
      "    # Final\n",
      "    tdotdot = num1/(den11*den12) - (tauterm * u)\n",
      "    \n",
      "    # Compute alpha-dot-dot\n",
      "    num2 = Bp * adot \\\n",
      "            - (Lp*g*mp*np.sin(a))/2 \\\n",
      "            - (Lp**2*mp*tdot**2*np.cos(a)*np.sin(a))/4 \\\n",
      "            + (Lp*Lr*mp*np.cos(a)*((mp*tdot*np.cos(a)*np.sin(a)*Lp**2*adot)/2 + (Lr*mp*np.sin(a)*Lp*adot**2)/2 + Br*tdot))/(2*((-mp*Lp**2*(np.cos(a))**2)/4 + (mp*Lp**2)/4 + mp*Lr**2 + Jr))\n",
      "    den21 = (Lp**2*Lr**2*mp**2*(np.cos(a))**2)/(4*((mp*Lp**2)/4 + Jp)*((-mp*Lp**2*(np.cos(a))**2)/4 + (mp*Lp**2)/4 + mp*Lr**2 + Jr)) - 1\n",
      "    \n",
      "    # Torque term\n",
      "    den2t1 = (Lp**2*Lr**2*mp**2*(np.cos(a))**2) / (4*((mp*Lp**2)/4 + Jp) * ((-mp*Lp**2*(np.cos(a))**2)/4 + (mp*Lp**2)/4 + mp*Lr**2 + Jr)) - 1\n",
      "    den2t2 = (mp*Lp**2)/4 + Jp\n",
      "    den2t3 = (-mp*Lp**2*(np.cos(a))**2)/4 + (mp*Lp**2)/4 + mp*Lr**2 + Jr\n",
      "    \n",
      "    # Final\n",
      "    adotdot = num2/(den21 * ((mp*Lp**2)/4 + Jp)) - ((Lp*Lr*mp*np.cos(a))/(2*den2t1*den2t2*den2t3))*u\n",
      "    \n",
      "    # Return xdot\n",
      "    return np.array([tdot, adot, tdotdot, adotdot])\n",
      "\n",
      "# Verify\n",
      "nlderiv(np.array([0.5, -0.3, 1.1, 0.9]), 0, np.array([0.0]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "array([  1.1       ,   0.9       , -41.46608396, -76.48942572])"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.array([0.0, 0.0, 0.0, 0.01])\n",
      "\n",
      "# Input\n",
      "u = 0\n",
      "\n",
      "# Simulation end time in seconds\n",
      "Tend = 20\n",
      "\n",
      "x = np.float64(x) # force data type\n",
      "u = np.float64(u)\n",
      "Tend = np.float64(Tend)\n",
      "t = np.linspace(0, Tend, 1e4)\n",
      "\n",
      "# Solve\n",
      "y = integrate.odeint(nlderiv, x, t, args=(u,))\n",
      "\n",
      "# Plot \n",
      "fig = plt.figure(figsize=(12,6), dpi=100)\n",
      "plot(t, y[:,1], t, y[:,0])\n",
      "legend(['Pendulum Angle', 'Rotor Angle'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "<matplotlib.legend.Legend at 0x1096bd190>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAswAAAFwCAYAAACsMS2JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0VOW9//HPN+ESLgkk3AMJgXJRqIh3ingcLVVIBW2p\n1PYnp8fVZW0Vq57aWtoeGq0/rf3Rc6y17UGPbQFPpYhWRMFLgeClRyzloshF8ICEgNxvIQRI8vz+\nmCTmMjNkZk8yk9nv11qzMnvvZ/bzTVZW/Pjl2Xubc04AAAAAQktLdAEAAABAMiMwAwAAABEQmAEA\nAIAICMwAAABABARmAAAAIAICMwAAABBBXAKzmaWb2VozWxzm+GNmttXM1pvZBfGYEwAAAGgN8eow\n3yVpo6QmN3U2s0JJQ5xzQyV9S9Lv4jQnAAAA0OI8B2YzGyCpUNJ/SbIQQyZLmiNJzrlVkrqbWR+v\n8wIAAACtIR4d5v+Q9H1J1WGO95dUUm97l6QBcZgXAAAAaHGeArOZXSdpn3NurUJ3l+uGNtrmedwA\nAABoE9p5/PxYSZNr1ilnSMoys7nOuX+uN6ZUUl697QE1+xowM0I0AAAAWpxzLlKjtwlzLj451cyu\nlHSvc25So/2FkqY75wrNbIykR51zY0J83sWrFiCeioqKVFRUlOgygCb43UQy4/cTycrMog7MXjvM\njbmaQm6TJOfcbOfcEjMrNLNtkk5IuiXOc3p26JCUnS1ZVD86AAAA+EHcArNzbqWklTXvZzc6Nj1e\n88Tb3r1S377SrFnS976X6GoAAACQbHz/pL+//lXq0UNasCDRlSBZBQKBRJcAhMTvJpIZv59IJXFb\nw+xVotYwT58e7DA//LB05IjUvn2rlwAAAIBWkgxrmNuc//1facIEKS9P2rJF+uxnE10RAAAIxbjY\nCFGKVzPW94F5xw6poEAaMiQYngnMAAAkr2T5l3Ekv3j+D5av1zA7FwzMAwdKgwZJ27cnuiIAAAAk\nG18H5v37pU6dpMxMAjMAAABC83Vg3rdP6tMn+D4/X9q5M7H1AAAAIPn4OjAfOCD17Bl837dv8J7M\nAAAAyaqgoEDLli2L+9hkFggE9NRTTyW0Bl8H5oMHGwbmTz5JbD0AAKBtKigoUOfOnZWZmam+ffvq\nlltu0YkTJ+I+j5k1+2K2aMbGw/bt25WWlqbbb789rudt7e8jFF8H5gMHgg8tkQjMAAAgdmaml156\nScePH9eaNWu0evVqPfjgg4kuq1XNnTtXn/3sZ/XnP/9Zp0+fTnQ5ceXrwFy/w9y1a/BrWVni6gEA\nAG1fbm6uJkyYoA0bNkiS3nnnHY0dO1bZ2dkaPXq0Vq5cWTc2EAho5syZGjdunLKysnTttdfq4MGD\ndcfnzZungQMHqmfPnnrooYcazPMv//Iv+rd/+7e67eLiYuXl5YWs6WxjCwoKNGvWLI0aNUqZmZn6\n5je/qb1792rixInq1q2bvvCFL+jIkSNhv2fnnObNm6eioiL16NFDixcvbnA8LS1Ns2fP1rBhw5Sd\nna3p06fXHauurtb3vvc99erVS4MHD9bjjz+utLQ0VVdXh5zr97//vUaMGKGcnBxNmDBBO1vhIjRf\nB+b6HWaJLjMAAIhd7T2iS0pKtHTpUl144YUqLS3Vddddp5kzZ+rw4cOaNWuWpkyZ0iAUP/PMM/rj\nH/+offv26fTp05o1a5YkaePGjbr99tv13//939q9e7cOHjyoXbt21X0unsszzEzPP/+8li1bpi1b\ntuill17SxIkT9fOf/1z79u1TdXW1HnvssbCff+utt7R3714VFhbqxhtv1Jw5c5qMefnll7V69Wq9\n9957WrBggV599VVJ0hNPPKFXXnlF69ev15o1a/TCCy+ErXXRokV6+OGH9Ze//EUHDhzQFVdcoa99\n7WvN+hl44evAXL/DLAXvmEFgBgAA0XLO6YYbblB2drauuOIKBQIBzZgxQ08//bQKCws1YcIESdL4\n8eN18cUX6+WXX5YUDKq33HKLhgwZooyMDE2dOlXr1q2TJC1cuFCTJk3SuHHj1KFDB/3sZz9TWlpa\nk3mjqTGSO++8U7169VJubq6uuOIKfe5zn9P555+vjh076ktf+pLWrl0b9rNz5szRpEmTlJGRoRtv\nvFGvvPKK9u/f32DMD3/4Q2VlZSkvL09XXXWV1q9fL0lasGCB7r77buXm5qp79+6aMWNG2Fr/8z//\nUzNmzNDw4cOVlpamGTNmaN26dSopKWn2zyEWvg/MOTmfbvfsGdwHAADaHrP4vGKb27Ro0SIdPnxY\nO3bs0OOPP66MjAx9/PHHevbZZ5WdnV33evvtt/VJvQ5d375969536tRJZTXrQ3fv3q0BAwbUHevc\nubN61P+n8TjrU3uv3Zo66m9nZGTU1dXYyZMntXDhQt14442SpNGjR6ugoEB/+tOfGoyr/3127ty5\n7nx79uxpsDyk/vfc2Mcff6y77rqr7mdZ+/MoLS1t7rcZE18H5mPHpG7dPt3OyZEOHUpcPQAAIHbO\nxecVT/n5+Zo2bZoOHz5c9zp+/Lh+8IMfnPWzubm5DTqn5eXlDZZydOnSReXl5XXbn0T4Z/JoxtZq\nbvf6L3/5i44dO6bbbrtN/fr1U79+/VRSUhJyWUYoteNrReoW5+fn64knnmjw8zxx4oTGjBnTrLli\n5evAfPx48Cl/tQjMAAAgnm6++WYtXrxYr732mqqqqlRRUaHi4uIGHdFwwXTKlCl66aWX9Pbbb+v0\n6dOaOXNmgwvhRo8erSVLlujw4cP65JNP9Oijj4atI5qx0ZozZ46++c1vasOGDVq/fr3Wr1+vt99+\nW+vXr6+78LEx51zd9z116lT96le/0u7du3XkyBE98sgjYdcwf/vb39ZDDz2kjRs3SpKOHj2qZ599\nNm7fSzgEZgIzAABoIQMGDNCiRYv00EMPqXfv3srPz9cvf/nLBiG5fjisf3HeyJEj9Zvf/EZf//rX\nlZubq5ycnAZLF6ZNm6bzzz9fBQUFmjBhgm666aawQTOasWerq77S0lItX75cd999t3r37l33uvDC\nCzVhwgTNnTs37Llrz3frrbfqmmuu0ahRo3TRRRfpi1/8otLT05us15akG264Qffdd59uuukmdevW\nTeedd17dxYMtyaJZLN6SzMy1di19+kjr1kn9+gW3f/tb6f33pd/9rlXLAAAAzWBmUV3khrZp6dKl\n+s53vqMdO3Z4Ok+435ea/VGtVqfDTIcZAAAgYSoqKrRkyRJVVlaqtLRU999/v7785S8nuqwGfBuY\nKyulU6ekLl0+3UdgBgAAaF3OORUVFSknJ0cXXnihRo4cqQceeCDRZTXQLtEFJEpZWTAs11+OQ2AG\nAABoXZ06ddK7776b6DIi8m2HufFyDCkYmLkPMwAAAOojMNdDYAYAAEBjvg3MZWVNA3NWllReLlVV\nJaYmAAAAJB/fBuZQHea0NKlr1+ATAAEAAADJx4G5rCwYjhvr1k06erT16wEAAEBy8m1grqiQMjKa\n7u/WjQ4zAABAOEVFRZo2bVqiy2hVvg3Mp05JHTs23Z+VRYcZAABEp6CgQJ07d1ZmZqb69u2radOm\n6VgzO3AFBQVavnx53Gtyzmnw4MEaOXJkXM97tkdqpyJPgdnMMsxslZmtM7ONZvZwiDEBMztqZmtr\nXj/xMme8nD4tdejQdD8dZgAAEC0z00svvaTjx49r/fr1ev/99/Xggw82+7OxPvK7srIy7LE33nhD\np06d0v79+7V69eqYzh+KHx9P7ikwO+cqJF3lnBstaZSkq8xsXIihK51zF9S8mvfb08LoMAMAgJbQ\np08fXXPNNfrggw/q9r344osaOXKksrOzddVVV2nz5s2SpGnTpmnnzp2aNGmSMjMzNWvWrIjjpWBH\n+he/+IVGjRqlzMxMVVdXh6xjzpw5mjJliq6//nrNmTOnwbFAIKCZM2dq3LhxysrK0rXXXquD9e6t\nO3fuXA0cOFA9e/bUgw8+GLEL/s4772js2LHKzs7W6NGjtXLlyth+cEnM85IM51x5zdsOktIlhXpW\nXtL17sMFZi76AwAAsajtvO7atUuvvPKKLrvsMknShx9+qK9//et67LHHdODAARUWFmrSpEmqrKzU\nvHnzlJ+fX9edvvfeeyOOrzV//nwtXbpUR44cUVpa0zhXXl6u5557Tl/96lc1depUzZ8/X2fOnGkw\n5plnntEf//hH7du3T6dPn64L6xs3btQdd9yhZ555Rnv27NHRo0e1e/fukN9zaWmprrvuOs2cOVOH\nDx/WrFmzNGXKFB04cCAuP9Nk4Tkwm1mama2TtFfSCufcxkZDnKSxZrbezJaY2Qivc8ZDpMDMkgwA\nABAN55xuuOEGZWVlKT8/X5/5zGf0k58EV6H++c9/1nXXXafPf/7zSk9P17333quTJ0/qb3/7W8hz\nnW28mem73/2u+vfvr46hwoyk559/XllZWbr88st19dVXS5JefvnluuNmpltuuUVDhgxRRkaGpk6d\nqnXr1kmSFi5cqMmTJ2vs2LFq3769HnjggbDrlp9++mkVFhZqwoQJkqTx48fr4osv1pIlS2L4KSav\ndl5P4JyrljTazLpJetXMAs654npD1kjKc86Vm9lESS9IGuZ1Xq/CrWFmSQYAAG2T3R+ff9B2P41+\nja6ZadGiRbr66qv1xhtvaNKkSVq9erUuvfRS7dmzR/n5+Q3G5uXlqbS0NOS5mjM+Ly8vYj1z5szR\nl7/8ZUlSenq6brjhBs2ZM0c33HBD3Zi+ffvWve/UqZPKysokSbt379aAAQMaHOvRo0fIeT7++GM9\n++yzWrx4cd2+ysrKupCeKjwH5lrOuaNm9rKkiyUV19t/vN77pWb2WzPLcc41WbpRVFRU9z4QCCgQ\nCMSrvCZOnZK6dGm6v1s3ac+eFpsWAAC0kFiCbkv4p3/6J91555267777tGLFCuXm5ur999+vO+6c\nU0lJifr37y+p6V0nzjY+1Gfq27Vrl5YvX66///3vWrBggaTgEo2KigodOnRIOTk5EevPzc3Vli1b\n6rZPnjzZYH1zffn5+Zo2bZqeeOKJiOdMpOLiYhUXF3s6h6fAbGY9JVU6546YWSdJX5B0f6MxfSTt\nc845M7tUkoUKy1LDwNzSTp2SQv2+sCQDAAB4dffdd+s//uM/tGrVKk2dOlU///nPtXz5cl1xxRX6\n1a9+pYyMDI0dO1ZS8CLBjz76qK4re7bxZzNv3jydc845WrFiRd0+55zGjh2rP/3pT5o+fXrdvlCm\nTJmiz33uc/qf//kfXXTRRSoqKgo79uabb9Yll1yi1157TZ///Od15swZvfPOOxo6dGiDgJ9IjZuw\n999/f/jBYXhdw9xP0vKaNcyrJC12zi0zs9vM7LaaMV+R9H7NmEcl3eRxzrjgLhkAAKCl9OzZU9/4\nxjf0yCOPaNiwYXr66ad15513qlevXnr55Ze1ePFitWsX7FvOmDFDDz74oLKzs/Xv//7vZx1/NnPn\nztXtt9+u3r1717369Omjb3/725o7d27duPpdajOr2x45cqR+/etf66abblJubq4yMzPVu3fvuvXS\n9ccOGDBAixYt0kMPPaTevXsrPz9fv/zlL8PeuaOtsmS5l56Zudas5Vvfki66SLrttob7V6yQ7r9f\n8ti5BwAAceblfsWIXVlZmbKzs7Vt2zYNHDgw0eU0W7jfl5r9US1450l/jXBbOQAA4HeLFy9WeXm5\nTpw4oXvvvVejRo1qU2E53gjMjbAkAwAA+N2LL76o/v37q3///vroo480f/78RJeUUHG7S0Zbw32Y\nAQAAQnvyySf15JNPJrqMpOHbDvPZ7sPMEikAAABIPg7M4TrMHTtK6enSyZOtXxMAAACSD4E5BC78\nAwAAQC0CcwgEZgAAANTior8QCMwAACSnSI+EBlqKbwNzuIv+JAIzAADJiIeWIFFYkhECgRkAAAC1\nCMwhEJgBAABQi8AcAoEZAAAAtXwbmFnDDAAAgObwbWCmwwwAAIDm8GVgdo4OMwAAAJrHl4H59Gmp\nfXspLcx339zA/NhjUufO0s9+Ft/6AAAAkDx8G5jDdZel5gXmXbuk+++XVq6Unnwy+BUAAACpx5eB\nOdL6Zal5gfmZZ6SvfEW65BLp//0/6Qc/CC71AAAAQGohMIfQvfvZA/PzzwcDsxT8euSI9MYb8asR\nAAAAyYHAHEK3bsEAHE5ZmfTee9KVVwa309Ole++VHnkkvnUCAAAg8XwZmJu7hjncEovVq6Xzz294\njn/+Z2ndOmn9+vjWCgAAgMRql+gCEuFsHeaOHYNd45Mng3fBaOydd6QxY5p+5p57pJ//PLi+OZyS\nEmn5cungQal3b+nCC6Xhw4PzAQAAIPkQmMOo7TKHCsyrVklf/WrT/bfdJg0eLG3ZEgzB9R05In3v\ne9ILL0jXXCP17RsM3kVF0r590sUXBy8g7N8/OGdZWXD+o0elY8c+fTkntWsXvC1e+/bBoG0WvEVe\n7cvs03kbd8nrb3ORon/V/x0BAACREZjDqA3M/fo1PfbBB9KoUU33Z2VJP/5xMDgvX/7pfZ6Li6Vv\nfEP64hel7duD4+o7eFD6+9+Dr61bpfJyqWvX4LjcXOncc4PvMzOD56yslM6cCb6qqz99Offp+/qB\nqHE4inQMqY//UQIA+Nm8edF/xlyS/NfTzFxr1bJsmfR//28w1IZzySXS449Ll13WcP+pU8EwfexY\n6HXQVVXBDnKvXtKttwY7ys89J/3+99KECfH9PgAAABAdM5NzLqqWoS8v+oumw9zY1q1SQUH4iwbT\n06UXX5Ty8oIPNunUKXgxIGEZAACgbWJJRhjhAvOmTcElEpF06RJ8mAkAAADaPjrMYXTvHvpezM0J\nzAAAAEgdBOYwevaUDhxoup/ADAAA4C+eArOZZZjZKjNbZ2YbzezhMOMeM7OtZrbezC7wMmc8nO3B\nJVLwHsn79jXdT2AGAADwF0+B2TlXIekq59xoSaMkXWVm4+qPMbNCSUOcc0MlfUvS77zMGQ/N6TD3\n6iXt399wX1WV9OGHTe+xDAAAgNTleUmGc6685m0HSemSDjUaMlnSnJqxqyR1N7M+Xuf1ojmBOVSH\needOqUeP4P2QAQAA4A+eA7OZpZnZOkl7Ja1wzm1sNKS/pJJ627skDfA6rxfNDcyNO8ybN7McAwAA\nwG8831bOOVctabSZdZP0qpkFnHPFjYY1vjl0yCeUFBUV1b0PBAIKBAJeywupOWuYe/Vq2mHevJnl\nGAAAAG1JcXGxiouLPZ0jrk/6M7N/k3TSOTer3r7/lFTsnJtfs71Z0pXOub2NPttqT/qbMSO4rOJH\nPwo/pqIi+DjqU6c+fXz0bbcFH4l9xx2tUiYAAADirNWf9GdmPc2se837TpK+IGlto2EvSvrnmjFj\nJB1pHJZbW3OWZGRkBF/1H16yZYt0zjktWxsAAACSi9clGf0kzTGzNAXD9zzn3DIzu02SnHOznXNL\nzKzQzLZJOiHpFo9zetacwCxJfftKe/YEH2IisSQDAADAjzwFZufc+5IuDLF/dqPt6V7mibfmrGGW\npIICaceO4IV+hw9LJ05I/fu3dHUAAABIJjzpL4JBg6Tt24PvN24MLsewqFa8AAAAoK0jMEdQPzCv\nXStdkPBnFAIAAKC1EZgjqB+Y16yRLrqoZesCAABA8vFlYG7uGuYRI6QNG4Lv//EP6cImq7UBAACQ\n6nwZmJvbYT73XKm0VPrww+BjsUePbvnaAAAAkFwIzBG0ayeNHSvdeqs0frzUvn3L1wYAAIDk4vnR\n2G1RcwOzJP3rv0pTpkjLl7dsTQAAAEhOvgzMzV3DLEnXXiuVlbVsPQAAAEheLMkAAAAAIiAwAwAA\nABEQmAEAAIAICMwAAABABL4MzNFc9AcAAAB/82VgpsMMAACA5vJdYK6uls6cocMMAACA5vFdYK5d\njmGW6EoAAADQFvg2MAMAAADN4bvAzPplAAAARIPADAAAAERAYAYAAAAi8F1gZg0zAAAAouG7wEyH\nGQAAANEgMAMAAAAREJgBAACACHwXmFnDDAAAgGj4LjDTYQYAAEA0CMwAAABABARmAAAAIAJPgdnM\n8sxshZl9YGYbzOy7IcYEzOyoma2tef3Ey5xesYYZAAAA0Wjn8fNnJN3jnFtnZl0l/cPMXnfObWo0\nbqVzbrLHueKCDjMAAACi4anD7Jz7xDm3ruZ9maRNknJDDDUv88QTgRkAAADRiNsaZjMrkHSBpFWN\nDjlJY81svZktMbMR8ZozFgRmAAAARMPrkgxJUs1yjIWS7qrpNNe3RlKec67czCZKekHSsHjMG4vT\npwnMAAAAaD7PgdnM2kt6TtLTzrkXGh93zh2v936pmf3WzHKcc4cajy0qKqp7HwgEFAgEvJbXxKlT\nXPQHAADgF8XFxSouLvZ0DnPOxf5hM5M0R9JB59w9Ycb0kbTPOefM7FJJC5xzBSHGOS+1NNcPfyh1\n6ybNmNHiUwEAACDJmJmcc1FdX+e1w3y5pJslvWdma2v2/UhSviQ552ZL+oqk75hZpaRySTd5nNMT\n1jADAAAgGp4Cs3PuLZ3lwkHn3G8k/cbLPPFUUSFlZCS6CgAAALQVPOkPAAAAiIDADAAAAERAYAYA\nAAAiIDADAAAAEfguMHPRHwAAAKLhu8BMhxkAAADRIDADAAAAERCYAQAAgAgIzAAAAEAEvgvMFRUE\nZgAAADSf7wLzqVPcJQMAAADN58vATIcZAAAAzUVgBgAAACIgMAMAAAAR+CowV1UFX+3bJ7oSAAAA\ntBW+Csy13WWzRFcCAACAtsKXgRkAAABoLgIzAAAAEAGBGQAAAIiAwAwAAABE4KvAXFHBU/4AAAAQ\nHV8FZjrMAAAAiBaBGQAAAIiAwAwAAABEQGAGAAAAIvBVYOaiPwAAAETLV4GZDjMAAACiRWAGAAAA\nIvAUmM0sz8xWmNkHZrbBzL4bZtxjZrbVzNab2QVe5vSCwAwAAIBotfP4+TOS7nHOrTOzrpL+YWav\nO+c21Q4ws0JJQ5xzQ83sMkm/kzTG47wxITADAAAgWp46zM65T5xz62rel0naJCm30bDJkubUjFkl\nqbuZ9fEyb6y46A8AAADRitsaZjMrkHSBpFWNDvWXVFJve5ekAfGaNxp0mAEAABCtuATmmuUYCyXd\nVdNpbjKk0baLx7zRIjADAAAgWl7XMMvM2kt6TtLTzrkXQgwplZRXb3tAzb4mioqK6t4HAgEFAgGv\n5TVQUSF16xbXUwIAACCJFRcXq7i42NM5zLnYm71mZgquTz7onLsnzJhCSdOdc4VmNkbSo865Jhf9\nmZnzUktz3HGHdO650vTpLToNAAAAkpSZyTnXePVDRF47zJdLulnSe2a2tmbfjyTlS5JzbrZzbomZ\nFZrZNkknJN3icc6YVVRInTolanYAAAC0RZ4Cs3PuLTVjHbRzLil6uidPcpcMAAAARMdXT/o7eZIO\nMwAAAKJDYAYAAAAi8FVgZg0zAAAAouWrwMwaZgAAAETLd4GZDjMAAACi4avAzJIMAAAARMtXgZkl\nGQAAAIiW7wIzHWYAAABEg8AMAAAAROCbwOxccA0zSzIAAAAQDd8E5jNnpLQ0qZ2nh4EDAADAb3wT\nmFmOAQAAgFj4JjBzSzkAAADEwjeBmQ4zAAAAYuGrwMwFfwAAAIiWrwIzHWYAAABEyzeBmTXMAAAA\niIVvAjMdZgAAAMTCV4GZNcwAAACIlm8CM0syAAAAEAvfBGaWZAAAACAWBGYAAAAgAl8FZtYwAwAA\nIFq+CcysYQYAAEAsfBOYWZIBAACAWBCYAQAAgAh8E5grKljDDAAAgOj5JjDTYQYAAEAsPAdmM/u9\nme01s/fDHA+Y2VEzW1vz+onXOWNx4oTUpUsiZgYAAEBb1i4O5/iDpF9LmhthzErn3OQ4zBWz8nKp\nc+dEVgAAAIC2yHOH2Tn3pqTDZxlmXufxig4zAAAAYtEaa5idpLFmtt7MlpjZiFaYswk6zAAAAIhF\nPJZknM0aSXnOuXIzmyjpBUnDWmHeBugwAwAAIBYtHpidc8frvV9qZr81sxzn3KHGY4uKiureBwIB\nBQKBuNVx4gQdZgAAAL8pLi5WcXGxp3OYc85zIWZWIGmxc+68EMf6SNrnnHNmdqmkBc65ghDjXDxq\nCWfAAOmdd4JfAQAA4E9mJudcVNfXee4wm9kzkq6U1NPMSiT9VFJ7SXLOzZb0FUnfMbNKSeWSbvI6\nZyzoMAMAACAWcekwx0NLd5g7dpSOHQt+BQAAgD/F0mH2xZP+KiuDrw4dEl0JAAAA2hpfBOby8uAd\nMizhd4MGAABAW+OLwMz6ZQAAAMTKF4G5tsMMAAAARMsXgZkOMwAAAGLli8BMhxkAAACx8kVgpsMM\nAACAWPkiMNNhBgAAQKx8EZjpMAMAACBWvgjMdJgBAAAQK18EZjrMAAAAiJVvAjMdZgAAAMTCF4G5\nvJwOMwAAAGLji8BMhxkAAACx8kVgpsMMAACAWPkiMNNhBgAAQKx8EZjpMAMAACBWvgjMx49LmZmJ\nrgIAAABtkW8Cc9euia4CAAAAbZEvAnNZGR1mAAAAxMYXgZklGQAAAIiVbwIzSzIAAAAQC18EZpZk\nAAAAIFYpH5jPnJEqK6WMjERXAgAAgLYo5QNz7fpls0RXAgAAgLYo5QNzWRnrlwEAABC7lA/M3CED\nAAAAXhCYAQAAgAg8BWYz+72Z7TWz9yOMeczMtprZejO7wMt8seCWcgAAAPDCa4f5D5ImhDtoZoWS\nhjjnhkr6lqTfeZwvatxSDgAAAF54CszOuTclHY4wZLKkOTVjV0nqbmZ9vMwZLZZkAAAAwIuWXsPc\nX1JJve1dkga08JwNsCQDAAAAXrTGRX+N74DsWmHOOizJAAAAgBftWvj8pZLy6m0PqNkXUlFRUd37\nQCCgQCDguQCWZAAAAPhXcXGxiouLPZ3DnPPW8DWzAkmLnXPnhThWKGm6c67QzMZIetQ5NybMeZzX\nWkK56y6poEC65564nxoAAABtjJnJORfVM6A9dZjN7BlJV0rqaWYlkn4qqb0kOedmO+eWmFmhmW2T\ndELSLV5I193GAAASAUlEQVTmiwVLMgAAAOCFp8DsnPtaM8ZM9zKHVyzJAAAAgBe+eNIfd8kAAABA\nrFI+MB87JmVlJboKAAAAtFUpH5iPHpW6d090FQAAAGirfBGYu3VLdBUAAABoq1I+MB85QocZAAAA\nsUvpwFxZKZWXc9EfAAAAYpfSgbn2gr+0lP4uAQAA0JJSOkqyfhkAAABepXRgZv0yAAAAvErpwEyH\nGQAAAF6lfGCmwwwAAAAvUjowHzlChxkAAADepHRgpsMMAAAAr1I6MNNhBgAAgFcpHZjpMAMAAMCr\nlA7MdJgBAADgVUoHZjrMAAAA8CqlAzMdZgAAAHiV0oH50CEpJyfRVQAAAKAtS+nAfPCg1KNHoqsA\nAABAW0ZgBgAAACJI2cB86pR0+rTUtWuiKwEAAEBblrKB+eDB4Ppls0RXAgAAgLYspQMzyzEAAADg\nFYEZAAAAiIDADAAAAERAYAYAAAAiSNnAfOgQgRkAAADepWxgpsMMAACAePAcmM1sgpltNrOtZnZf\niOMBMztqZmtrXj/xOmdzEJgBAAAQD+28fNjM0iU9Lmm8pFJJfzezF51zmxoNXemcm+xlrmgdOEBg\nBgAAgHdeO8yXStrmnNvhnDsjab6k60OMa/XHh+zdK/Xt29qzAgAAINV4Dcz9JZXU295Vs68+J2ms\nma03syVmNsLjnM3yySdSnz6tMRMAAABSmaclGQqG4bNZIynPOVduZhMlvSBpWKiBRUVFde8DgYAC\ngUBsRblgh5nADAAA4G/FxcUqLi72dA5zrjmZN8yHzcZIKnLOTajZniGp2jn3SITPbJd0kXPuUKP9\nzkst9R05IuXnS8eOxeV0AAAASBFmJudcVMuFvS7JWC1pqJkVmFkHSV+V9GKjovqYmdW8v1TBkH6o\n6anih+4yAAAA4sXTkgznXKWZTZf0qqR0SU855zaZ2W01x2dL+oqk75hZpaRySTd5rPmsuOAPAAAA\n8eJpSUY8xXNJxoIFwdfChXE5HQAAAFJEIpZkJCU6zAAAAIiXlAzMjW8pt2n/Jo35rzF6dduriSsK\nAAAAbVJKBubSUql/vbtBf//17yunU45+tPxHiSsKAAAAbVJKBuadO6W8vOD7vWV79dbOt7TgxgXa\ndmibDpQfSGxxAAAAaFNSMjCXlATvwyxJL2x+QYVDC9W1Q1dd1v8yrdq1KrHFAQAAoE1JucDsnLRr\n16cd5hU7Vuiaz1wjSTq357nafGBzAqsDAABAW5NygXn/fqlLF6lz5+D230r+psvzLpckDe85XFsO\nbklgdQAAAGhrUi4wl5R82l0uOVqiisoKDckZIkkanD1YHx/9OIHVAQAAoK1JucBc/4K/t0ve1uX5\nl6vmydzq17Wfdh/fncDqAAAA0NakXGD++ONPL/h7e+fbdcsxJCk3M5fADAAAgKikXGDeulUaNiz4\n/u2ShoG5R+ceKjtdplOVpxJUHQAAANqalAzMQ4dKx08d14cHP9SF/S6sO5ZmaerTpY/2lO1JYIUA\nAABoS1IuMH/4YbDDvKp0lUb3Ha2O7To2ON6zc08dLD+YoOoAAADQ1qRUYD55UvrkE2ngwKbrl2tl\nd8rW4YrDCagOAAAAbVFKBeaPPpIGDZLatfv0DhmNZWdk69DJQwmoDgAAAG1RSgXmDRukESOkyupK\nrSpdFbLDnNMpR4dP0mEGAABA86RUYF67VrrgAum9ve9pQNYA9ejco8kYOswAAACIRsoF5tGjpTc/\nflPj8saFHMMaZgAAAEQjZQKzc9K6dcEO81slb+mKgVeEHJedkc2SDAAAADRbygTmHTuktDSpXz8X\n7DDnh+4wZ3XM0vHTx1u3OAAAALRZKROYV66UrrxSem/femV2zFRB94KQ47p26Kqy02WtWxwAAADa\nrJQKzIGAtHTrUk0cMjHsuGgD896yvVq5Y2UcKgQAAEBblBKB2Tlp2bJgYF6ybYkKhxaGHZvZMbPZ\ngbmqukrj543X9fOv1x/W/iFO1QIAAKAtSYnAvHq11KmT1DV3lz7Y94ECBYGwY6PpMC/bvkwd0jvo\n1Ztf1f0r71dVdVWcKgYAAEBbkRKBeeFCacoUaf6GZ/Tlc7+sjHYZYcdGE5iXbl2qL53zJV024DL1\n6tJLy7Yvi1fJAAAAaCPafGA+dUqaM0f6PzdX66m1T+kb538j4viuHbo2+y4Zb+58U1cPulqSNOXc\nKVq0eZHnegEAANC2tPnAPH++dN550oe2SFkds8LeTq5WbYfZORdxXFV1lTbu36jzep8nSbp++PVa\ntGXRWT8HAACA1NKmA/PJk9LMmdIPf3xaP1nxE/34ih/LzCJ+pl1aO3VI76CKyoqI47Yd2qa+Xfsq\ns2OmJOmcnueoY7uOem/ve3GrHwAAAMnPc2A2swlmttnMtprZfWHGPFZzfL2ZXeB1zlrf/740dqy0\n0j2ogd0GavLwyc36XHPWMX+w/wN9tvdn67bNTBOHTNTSbUs91QwAAIC2xVNgNrN0SY9LmiBphKSv\nmdm5jcYUShrinBsq6VuSfudlTil4G7lHHpFee00ad+cf9Md1f9RTk586a3e5VnMC886jOzWw28AG\n+wjMAAAA/uO1w3yppG3OuR3OuTOS5ku6vtGYyZLmSJJzbpWk7mbWJ9YJt2wJ3hFjzrMHNfbBOzTr\n7w/o1ZtfVb/Mfs0+R2aHs9+LufRYqfpn9W+w76pBV2nNnjU6WnE0ptoBAADQ9rTz+Pn+kkrqbe+S\ndFkzxgyQtDfUCSsrpfLy4Prk8nLpwAGnTdtO6t1Ne1S8plQfn1yvgf/0hnZf+ro6d/m6/vGtfyin\nU05URTfnThmlx0s1qs+oBvs6t++sy/Mu11//96+aMmLKWeepdtVKsza9TBwAAMD3vAbm5t4yovFa\niZCfsx/mSGmVUlpV8KtVSWlVSqvuoMy0XPW/Olc3DhqhKwoKdf05s6MOyrW6duiqE6dPRBxTerxp\nh1kKLstYsnVJ2MDsnNPsf8zWI28/op1Hd+oz2Z/RlHOn6OZRN2tk75Fh56t21dp9fLcOnTyk8jPl\nOl11Ws45VbvqsC+n4PGq6ipVuaqI76tdtZxzdZ+pPXc8tiMdC/XzabDdjDHJNM7L+QEAQNvkNTCX\nSsqrt52nYAc50pgBNfua+FfdqvZmSrM0Ba4M6Oqrrla6pTd7bXJzdenQpXlLMjJDBOahE/WLv/1C\nzrmQdd331/v02kev6dkbn9X5fc7X+/ve1/wN83Xt09eqT9c+mjRskgZ1H6QqV6UdR3Zoy8Et2nJg\ni7Ye2qqsjlnq1bmXOrXvpA7pHZRmaWFfJmuwnZ6WHvxq6UpPS1e6fbpd95man23tZ6Pdrp2j/vFI\nnzVZyJ+RNfr/p+aMSaZxXs4FAABa16bVm7R59WZP5zAv9xU2s3aStkj6vKTdkt6V9DXn3KZ6Ywol\nTXfOFZrZGEmPOufGhDiXa617HE/7yzRdM/gaTTt/Wsjjzjl1fqiz9n9/v7p26Nrk2JBfD9HzU5/X\n+X3Pb3DsuY3P6d7X79Wab61RdqfsBseqqqu0YscKLd++XCXHSpRu6crvlq/hPYZreM/hGtZjmLI6\nZsX3GwUAAEADZibnXFRdLU8dZudcpZlNl/SqpHRJTznnNpnZbTXHZzvnlphZoZltk3RC0i1e5oyH\nru0j3yXjcMVhdUjv0CQsS8EfcuGQQi3dtrRBYN52aJu+/fK3teTrS5qEZUlKT0vX+MHjNX7w+Ph8\nEwAAAGgVnq9Ic84tdc4Nd84Ncc49XLNvtnNudr0x02uOn++cW+N1Tq/Odlu5cMsxan1x2Be1cOPC\nurWrJ8+c1NRnp+qnV/5Ul/S/JO71AgAAIHF8eQuHswbmMBf81frC4C/oZOVJvbD5BVVVV+nWxbdq\neM/huuOSO1qiXAAAACSQbwPziTPh75Jxtg5zelq6Zl83W7cuvlUXPXGRdh3bFdWDUwAAANB2eL1L\nRpt0trtklB6PHJglaVz+OL1767vafni7riy4Uu3SfPmjBAAASHm+THnNWcPc+A4YoQzOHqzB2YPj\nWRoAAACSjG+XZEQKzLvLdp+1wwwAAAB/IDCHUHos8kV/AAAA8A8Ccwilx0uVm5nbihUBAAAgWfk2\nMIe7S8bpqtM6fPKw+nTp08pVAQAAIBn5MjB3aR/+Lhl7ju9R7y69lZ6W3spVAQAAIBn5MjBHWpJx\ntoeWAAAAwF8IzI2c7aElAAAA8BdfBuaMdhk6XXVaVdVVTY6VHCtRXlZeAqoCAABAMvJlYDazsBf+\n7Tq2S3ndCMwAAAAI8mVglsJf+Lfr2C4NyBqQgIoAAACQjHwbmMOtYyYwAwAAoD5fB+bjp4432c8a\nZgAAANTn28DcPaO7jlQcabCvsrpSe8v2ql9mvwRVBQAAgGTj28Dcs3NP7S/f32DfnuN71LNzT3VI\n75CgqgAAAJBsfBuYe3XupQPlBxrs23poq4b2GJqgigAAAJCM/BuYu/TS/hMNO8xbD27V0BwCMwAA\nAD7l28AcaknG1kNbNSRnSIIqAgAAQDLybWDu1blXyMBMhxkAAAD1+Tcwd+mlfSf2Ndi3Yd8Gjeg1\nIkEVAQAAIBn5NjAP6j5I2w9vr9s+UnFEe8v2aliPYQmsCgAAAMnGt4E5r1ue9p3Yp4rKCknSuk/W\naVSfUUpPS09wZQAAAEgmvg3M7dLaKb9bfl2X+W8lf9Nl/S9LcFUAAABINr4NzJI0rMcwbTqwSZL0\n+v++rvGDxye4IgAAACQbXwfmy/Mu15sfv6kD5Qe0Zs8aXVlwZaJLAgAAQJJpF+sHzSxH0p8lDZS0\nQ9JU59yREON2SDomqUrSGefcpbHOGW/XfOYa3fjsjcrqmKXJwyera4euiS4JAAAAScZLh/mHkl53\nzg2TtKxmOxQnKeCcuyCZwrIkXZR7kUb3Ha3H//64Hgg8kOhykKSKi4sTXQIQEr+bSGb8fiKVeAnM\nkyXNqXk/R9INEcaah3la1HNTn9P+7+/XoOxBiS4FSYo/+khW/G4imfH7iVTiJTD3cc7trXm/V1Kf\nMOOcpL+a2Wozu9XDfC3CzJRmvl7KDQAAgAgirmE2s9cl9Q1x6Mf1N5xzzsxcmNNc7pzbY2a9JL1u\nZpudc2/GVi4AAADQusy5cDn3LB8026zg2uRPzKyfpBXOuXPO8pmfSipzzv0yxLHYCgEAAACi4JyL\narlwzHfJkPSipG9IeqTm6wuNB5hZZ0npzrnjZtZF0jWS7g91smgLBwAAAFqDlw5zjqQFkvJV77Zy\nZpYr6Unn3BfNbLCk52s+0k7SfzvnHvZeNgAAANA6Yg7MAAAAgB8k/PYQZjbBzDab2VYzuy/R9QD1\nmdkOM3vPzNaa2buJrgf+ZWa/N7O9ZvZ+vX05Zva6mX1oZq+ZWfdE1gh/CvO7WWRmu2r+dq41swmJ\nrBH+ZWZ5ZrbCzD4wsw1m9t2a/VH9/UxoYDazdEmPS5ogaYSkr5nZuYmsCWgkaR+8A9/5g4J/K+tr\n7gOkgJYU6nfTSfr3mr+dFzjnXklAXYAknZF0j3NupKQxku6oyZpR/f1MdIf5UknbnHM7nHNnJM2X\ndH2CawIa44JUJFzN7TgPN9odzQOkgBYR5ndT4m8nkoBz7hPn3Lqa92WSNknqryj/fiY6MPeXVFJv\ne1fNPiBZJPWDd+B7zX2AFJAId5rZejN7iuVCSAZmViDpAkmrFOXfz0QHZq44RLK73Dl3gaSJCv4z\nzhWJLggIxQWv4OZvKpLF7yQNkjRa0h5JTZ6/ALQmM+sq6TlJdznnjtc/1py/n4kOzKWS8upt5ynY\nZQaSgnNuT83X/ZL+ouAyIiBZ7DWzvpJU8wCpfQmuB5AkOef2uRqS/kv87UQCmVl7BcPyPOdc7XND\novr7mejAvFrSUDMrMLMOkr6q4ANRgIQzs85mllnzvvbBO+9H/hTQqmofICWFeYAUkAg1AaTWl8Tf\nTiSImZmkpyRtdM49Wu9QVH8/E34fZjObKOlRSemSnuLBJkgWZjZIwa6yxIN3kGBm9oykKyX1VHC9\n3UxJixTiAVKJqhH+FOJ386eSAgoux3CStku6rd56UaDVmNk4SW9Iek+fLruYIeldRfH3M+GBGQAA\nAEhmiV6SAQAAACQ1AjMAAAAQAYEZAAAAiIDADAAAAERAYAYAAAAiIDADAAAAERCYAQAAgAgIzAAA\nAEAE/x/Fwry3X8SL3AAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1098fd3d0>"
       ]
      }
     ],
     "prompt_number": 12
    }
   ],
   "metadata": {}
  }
 ]
}