furata_lin.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:b79f74880baeceb4b13fd523fb385b7a7c32033528b6f59c9b1a78d28891c2aa"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Linearised Furata Pendulum\n",
      "\n",
      "**Jon Sowman <j.sowman@soton.ac.uk>**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import cvxopt\n",
      "from stepc.linsystem import LinSystem\n",
      "import scipy.integrate\n",
      "from stepc.c2d import c2d\n",
      "from stepc.simulator import Simulator\n",
      "from stepc.controller import *\n",
      "from pprint import pprint"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create the system (4 states, 1 ip, 2 op)\n",
      "sysc = LinSystem(4, 1, 2)\n",
      "\n",
      "# 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",
      "Dp = 0.0\n",
      "Dr = 0.0\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\n",
      "\n",
      "# Continuous time system\n",
      "sysc.A[0, 0] = 0.0\n",
      "sysc.A[0, 1] = 0.0\n",
      "sysc.A[0, 2] = Jt\n",
      "sysc.A[0, 3] = 0.0\n",
      "\n",
      "sysc.A[1, 0] = 0.0\n",
      "sysc.A[1, 1] = 0.0\n",
      "sysc.A[1, 2] = 0.0\n",
      "sysc.A[1, 3] = Jt\n",
      "\n",
      "sysc.A[2, 0] = 0.0\n",
      "sysc.A[2, 1] = 0.25*mp*mp*Lp*Lp*Lr*g\n",
      "sysc.A[2, 2] = -(Jp + 0.25*mp*Lp*Lp)*Dr\n",
      "sysc.A[2, 3] = -0.5*mp*Lp*Lr*Dp\n",
      "\n",
      "sysc.A[3, 0] = 0\n",
      "sysc.A[3, 1] = 0.5*mp*Lp*g*(Jr + mp*Lr*Lr)\n",
      "sysc.A[3, 2] = -0.5*mp*Lp*Lr*Dr\n",
      "sysc.A[3, 3] = -(Jr + mp*Lr*Lr)*Dp\n",
      "\n",
      "sysc.A = sysc.A / Jt\n",
      "\n",
      "sysc.B[0, 0] = 0.0\n",
      "sysc.B[1, 0] = 0.0\n",
      "sysc.B[2, 0] = Jp + 0.25*mp*Lp*Lp\n",
      "sysc.B[3, 0] = 0.5*mp*Lp*Lr\n",
      "sysc.B = sysc.B / Jt\n",
      "\n",
      "# Absorb actuator model\n",
      "sysc.B = sysc.B * (kt/Rm)\n",
      "sysc.A[2,2] = sysc.A[2,2] - (kt*km/Rm)*sysc.B[2, 0]\n",
      "sysc.A[3,2] = sysc.A[3,2] - (kt*km/Rm)*sysc.B[3, 0]\n",
      "\n",
      "sysc.C = np.hstack((np.eye(2), np.zeros([2, 2])))\n",
      "\n",
      "# Timestep\n",
      "Ts = 0.005\n",
      "\n",
      "# Discretise the system for the MPC controller\n",
      "sys = c2d(sysc, Ts)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# The continuous time derivative\n",
      "def nlderiv(x, t, u):\n",
      "    xdot = sysc.A.dot(x) + sysc.B.dot(u)\n",
      "    return xdot"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create a linear MPC controller and set parameters that aren't time varying\n",
      "mpc = LinearMPCController(sys)\n",
      "mpc.Q = np.diag([1e-3, 1e3, 1e-9, 1e-9])\n",
      "mpc.P = 1\n",
      "mpc.R = np.diag([1e-3])\n",
      "mpc.set_prediction_horizon(10)\n",
      "mpc.set_control_horizon(3)\n",
      "target = np.array([[0.0], [0.0], [0.0], [0.0]])\n",
      "\n",
      "# Setup MPC controller\n",
      "mpc.set_target(target)\n",
      "mpc.generate_matrices()\n",
      "# Argument order (ulim, dulim, zlim)\n",
      "# Each is a 2-col matrix of [umin, umin_weight, umax, umax_weight; ...]\n",
      "# Weights must be +ve. Use np.inf for a hard constraint.\n",
      "mpc.generate_constraints(np.array([[-12, np.inf, 12, np.inf]]), # motor voltage Vm\n",
      "                         np.array([[-1e3, np.inf, 1e3, np.inf]]), # make any adjustments you like\n",
      "                         np.array([[-pi, 1e-9, pi, 1e-9],\n",
      "                                   [-10, 1e-9, 10, 1e-9],\n",
      "                                   [1, 1e-9, 1, 1e-9],\n",
      "                                   [1, 1e-9, 1, 1e-9]]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Initial state and input\n",
      "x_angle_deg = 45\n",
      "x = np.array([[0.0], [x_angle_deg/180.0*np.pi], [0.0], [15.0]])\n",
      "u = np.array([[0.0]])\n",
      "\n",
      "# Reset\n",
      "mpc.reset()\n",
      "\n",
      "# Start at what time?\n",
      "t = 0\n",
      "\n",
      "# Storage vectors\n",
      "tall = t\n",
      "uall = u\n",
      "xall = x\n",
      "\n",
      "# Run simulation\n",
      "while t < .5:\n",
      "    # Get control inputs\n",
      "    (u, u_traj, z_traj) = mpc.controlmove(x)\n",
      "    \n",
      "    x = sys.A.dot(x) + sys.B.dot(u)\n",
      "    \n",
      "    # Store\n",
      "    xall = np.hstack((xall, x))\n",
      "    uall = np.append(uall, u)\n",
      "    tall = np.append(tall, t)\n",
      "    \n",
      "    # Advance time\n",
      "    t += Ts\n",
      "    \n",
      "# Results\n",
      "print \"Final state is [%.6f, %.5f]\" % (x[0], x[1])\n",
      "plt.plot(tall, xall[1,:] * 180.0 / np.pi)\n",
      "plt.step(tall, uall)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Final state is [-14.125520, -0.09599]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 51,
       "text": [
        "[<matplotlib.lines.Line2D at 0x115abfad0>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHRJJREFUeJzt3Xm4FNWZx/Hvy+aCKNxxAgRRxA0xbmgUY2LaNRAVNLhl\nSCC4hOBEjRoDJpnxOqPjkjDRuI8xij6K4kZwQUGlNVEjREUii6gJQYzigmjEKNs7f5y+0FzvvfS9\n1d2nu+v3eZ56uqq6lpeieetw6tQ55u6IiEhtaxc7ABERKT0lexGRFFCyFxFJASV7EZEUULIXEUkB\nJXsRkRRIlOzNbBczezFv+tDMzjSzOjObbmYLzWyamXUtVsAiItJ6Vqx29mbWDngT2A84A3jP3S83\ns7FAN3cfV5QTiYhIqxWzGucw4DV3fwMYAkzIrZ8AHFPE84iISCsVM9mfBEzMzXd396W5+aVA9yKe\nR0REWqkoyd7MOgFHA3c3/s5DPZH6ZBARiahDkY4zGHje3d/NLS81sx7u/raZ9QTeabyDmekGICLS\nBu5urd2nWNU432Z9FQ7AFGBkbn4kMLmpndxdkzsXXHBB9BgqZdK10LXQtWh5aqvEyd7MOhMezt6X\nt/pS4HAzWwgcklsWEZFIElfjuPsKYOtG65YRbgAiIlIB9AZtBchkMrFDqBi6FuvpWqyna5Fc0V6q\navWJzTzWuUVEqpWZ4REf0IqISAVTshcRSQElexGRFFCyFxFJASV7EZEUULIXEUkBJXsRkRRQshcR\nSYFi9XpZUdauhQUL4Nlnw/S3v8EXvgA9ekCvXvCtb0GfPrGjFBEpn5p6g9Yd7rkHfvYzWLMGDjgg\nTDvsAO+9B2+/Da+9FrbZf38YMwYGD4b27YsahohIybT1DdqaSfbZLPzkJyHJX3opHH5489t+8gnc\ndRdcfTVsuinccQdst13RQhERKZnUdpfgHpL7d74D55wDs2a1nOgBNt8cRo0K2x5zDHz5y3DvveWJ\nV0Qkhqou2a9aFapinn8eHnww1Me3xcyZ8O1vw5FHwhVXQLuqvwWKSK1KXTXORx/BsGGwySYwcSJ0\n6ZIsng8/DMl+wAC48kqwVl9KEZHSi1aNY2ZdzeweM5tvZvPMbH8zqzOz6Wa20MymmVnXpOfJt3o1\nnHgibLstTJ6cPNEDbLVV+N/BH/4AP/958uOJiFSSYlRYXAk87O67AnsAC4BxwHR33xl4PLdcNOec\nE5pX3nADdChi49GuXeHRR+H++8NzABGRWpGoGsfMtgJedPe+jdYvAL7u7kvNrAeQdfd+jbZpUzXO\nNdeE6ZlnQnIuhTffhK9+FX75y1BVJCJSKaLU2ZvZXsANwDxgT+B54EfAEnfvltvGgGUNy3n7tjrZ\nT58OI0bA009D374b3z6J556Do48OD3979y7tuURECtXWZJ+0EqQDMAD4obvPMrMraFRl4+5uZk1m\n9fr6+nXzmUymxXEm338fRo6EO+8sfaKH8NLV2WeHJp1PPKEXr0Qkjmw2SzabTXycpCX7HsCz7r59\nbvmrwPlAX+Bgd3/bzHoCM5JW44wYAXV1oWlkuaxZE9rsH3ww/Md/lO+8IiLNiVKyzyXzN8xsZ3df\nCBwGzM1NI4HLcp+Tk5znoYdC1c2cOUmO0nrt28Ntt8E++8Chh8JXvlLe84uIFEvidvZmtifwG6AT\n8DowCmgPTAK2BRYBJ7j78kb7FVSyX74cdt8dbr01lLBjmDQJLr441N8Xs/WPiEhr1exLVaedFhLs\nddeVIahmuIcbzYknhjd2RURiqclk39B3zfz5sOWWZQqsGS+9BEccEWKpq4sbi4ikV80le3c47DA4\n4QQYPbqMgbVgzBjo1Cl0pyAiEkPNJftp0+CMM+Dll6FjxzIG1oL33oNddw3dKe+2W+xoRCSNaqqL\n47VrYdy48FC0UhI9wNZbhyaY55wTOxIRkdapyGR/113hoWwldlUwZgy88gr88Y+xIxERKVzFVeOs\nXBmqSm68EQ45JEJgBbj2Wpg6FR54IHYkIpI2NVONM2FCGDO2UhM9wMknhzb3s2fHjkREpDAVVbJf\nuxb69w9t6mO9QFWo8eNDZ2mTJsWORETSpCZK9g8/DJ07Qwv9oVWM0aNDq5wFC2JHIiKycRWV7MeP\nh3PPrY4hAbfYAs48U4OciEh1qJhqnOefh2OPhddfr6zmli1Zvjw8X3jhBdhuu9jRiEgaVH01zvjx\noaRcLYkewkhZI0aE4RFFRCpZRZTsFy+GvfaCv/41DPxdTRYsCM8YFi8OXSmIiJRSVZfsr7oKRo2q\nvkQP0K9faEF0//2xIxERaV70kv1nn8E224Q3UnfYIUooid19dxgEvQgjh4mItKhqS/ZTpoTBSao1\n0UPohnnhQpg3L3YkIiJNS5zszWyRmc0xsxfNbGZuXZ2ZTTezhWY2zcy6Nrf/TTfBKackjSKujh3D\nn+H662NHIiLStGIMS/hXYB93X5a37nLgPXe/3MzGAt3cfVyj/XzxYmevvWDJEthss0RhRLd4Mey9\nd/js3Dl2NCJSq2JX4zQ+8RBgQm5+AnBMUzvdcksY6q/aEz3AttvCgQfCnXfGjkRE5POKkewdeMzM\n/mRmp+XWdXf3pbn5pUD3pna8+ebqr8LJd+qpoSM3EZFK06EIxzjQ3d8ys38FppvZBr3FuLubWZN1\nRStW1DNlSugqOJPJkKmGTnFaMGhQ6BFz8eJQ0hcRSSqbzZItQlO/oja9NLMLgI+B04CMu79tZj2B\nGe7er9G2fuWVzplnFu30FWH0aOjbF8aOjR2JiNSiKHX2Zra5mXXJzXcGjgD+DEwBRuY2GwlMbmr/\n4cOTnL0yDR8Ot98eOwoRkQ0lKtmb2fZAw7ujHYDb3f0SM6sDJgHbAouAE9x9eaN9WxxwvFqtXQt9\n+sBDD4X3B0REiqmtJfvob9DWorFjoV07uOSS2JGISK1Rsq8gc+bA0UeHjt3aRX9HWURqSex29pJn\n992hSxd4+unYkYiIBEr2JWCmB7UiUllUjVMiixbBvvvCW29V14AsIlLZVI1TYfr0ge23hyefjB2J\niIiSfUkNGwb33hs7ChERVeOU1Kuvwte+Bm++Ce3bx45GRGqBqnEq0E47Qffu8MwzsSMRkbRTsi8x\nVeWISCVQNU6JzZ0LgwfD3/4WmmSKiCShapwK1b8/bL45zJoVOxIRSTMl+xIzU1WOiMSnZF8Gxx0H\n99wDKai1EpEKpWRfBnvtFT5feiluHCKSXkr2ZdBQlXPPPbEjEZG0UrIvk+OPh7vvVlWOiMRRlGRv\nZu3N7EUzeyC3XGdm081soZlNM7OuxThPNdt3X/jsM3j55diRiEgaFatkfxYwD2got44Dprv7zsDj\nueVUMwsPau++O3YkIpJGiZO9mW0DfBP4DdDQ0H8IMCE3PwE4Jul5akFDqxwRkXIrRsn+V8B5wNq8\ndd3dfWlufinQvQjnqXr77w8ffxzeqhURKacOSXY2s6OAd9z9RTPLNLWNu7uZNflYsr6+ft18JpMh\nk2nyEDUjv1XObrvFjkZEqkE2myWbzSY+TqK+cczsf4DvAquBTYEtgfuALwMZd3/bzHoCM9y9X6N9\nU9E3TmPPPAPf/74e1IpI20TpG8fdf+ruvd19e+Ak4Al3/y4wBRiZ22wkMDnJeWrJwIGwfDnMnx87\nEhFJk2K3s28oql8KHG5mC4FDcssCtGunF6xEpPzUxXEEzzwDp5wC8+ap22MRaR11cVxFDjggvGA1\ne3bsSEQkLZTsIzCDf/s3uP322JGISFqoGieSefPg8MNh8WINRi4ihVM1TpXp3x++8AV46qnYkYhI\nGijZRzR8uKpyRKQ8VI0T0ZIlsMce8Pe/w6abxo5GRKqBqnGq0DbbwJ57wtSpsSMRkVqnZB+ZqnJE\npBxUjRPZ8uXQpw/85S9QVxc7GhGpdKrGqVJdu8LgwTBxYuxIRKSWKdlXgFGj4OabY0chIrVMyb4C\nHHoovPMOzJkTOxIRqVVK9hWgfXsYOVKlexEpHT2grRCvvx46SFuyBDp1ih2NiFQqPaCtcjvsELpQ\nePDB2JGISC1KlOzNbFMze87MZpvZPDO7JLe+zsymm9lCM5tmZl2LE25tO/lk+O1vY0chIrUocTWO\nmW3u7p+YWQfgD8CPgSHAe+5+uZmNBbq5+7hG+6kap5EVK8JbtXPnwhe/GDsaEalE0apx3P2T3Gwn\noD3wASHZT8itnwAck/Q8adC5M5xwgkr3IlJ8iZO9mbUzs9nAUmCGu88Furv70twmS4HuSc+TFj/4\nAfzf/8GaNbEjEZFaUoyS/Vp33wvYBjjIzA5u9L2zfiBy2Yi99w5VOA8/HDsSEaklHYp1IHf/0Mwe\nAvYBlppZD3d/28x6Au80tU99ff26+UwmQyaTKVY4VW3MGLjuOjj66NiRiEhs2WyWbDab+DiJHtCa\n2dbAandfbmabAY8CFwLfAN5398vMbBzQVQ9oC/fPf0Lv3jBrFmy/fexoRKSStPUBbdJkvzvhAWy7\n3HSbu//CzOqAScC2wCLgBHdf3mhfJfsWnHNOeLnq0ktjRyIilSRKsk9Cyb5lr7wCBx0UBiTfZJPY\n0YhIpdAbtDVml13gS1+C++6LHYmI1AIl+wp2+ulw1VWxoxCRWqBkX8GGDg2Dkc+cGTsSEal2SvYV\nrEMHOPNMuOKK2JGISLXTA9oK9+GHofnlnDmh3xwRSTc9oK1RW20F3/0uXHNN7EhEpJqpZF8FXn8d\nBg6ERYtCZ2kikl4q2dewHXaAAw+E226LHYmIVCuV7KvEk0/C6NEwbx600y1aJLVUsq9xBx0EXbrA\nlCmxIxGRaqRkXyXMYNw4uOQS0H+IRKS1lOyryLHHhqaYM2bEjkREqo2SfRVp1w7Gjg2lexGR1tAD\n2iqzcmVonXP//bDvvrGjEZFy0wPalOjUCc49V/3ci0jrqGRfhVasCF0oPPUU9OsXOxoRKacoJXsz\n621mM8xsrpm9bGZn5tbXmdl0M1toZtPMrGuS88iGOncOHaRddFHsSESkWiQdlrAH0MPdZ5vZFsDz\nwDHAKOA9d7/czMYC3TQGbXF99BHsuGN42WrXXWNHIyLlEqVk7+5vu/vs3PzHwHygFzCEMDYtuc9j\nkpxHPm/LLcM4tf/1X7EjEZFqULQ6ezPrAzwJfAlY7O7dcusNWNawnLe9SvYJffxxaJnz+ONhCEMR\nqX1RW+PkqnDuBc5y93/kf5fL6MrqJbDFFvDjH8OFF8aOREQqXYekBzCzjoREf5u7T86tXmpmPdz9\nbTPrCbzT1L719fXr5jOZDJlMJmk4qXP66TB+PLz0Euy5Z+xoRKTYstks2Ww28XGSPqA1Qp38++5+\ndt76y3PrLjOzcUBXPaAtnV/9CrJZ+N3vYkciIqXW1mqcpMn+q8BTwBzWV9WcD8wEJgHbAouAE9x9\neaN9leyL5NNPYZdd4I47Qr/3IlK7oiT7JJTsi+vWW+H66+Hpp0MPmSJSm9RdQsoNHx7erJ08eePb\nikj6qGRfQx59NLxZ+/LL0LFj7GhEpBRUsheOOAJ694abboodiYhUGpXsa8wLL8BRR8Err4RhDEWk\ntqhkLwAMGBBK+BdfHDsSEakkKtnXoLfegt13h2efhZ12ih2NiBSTSvayTs+e8JOfhI7SRERAyb5m\nnXUWLFgAjzwSOxIRqQRK9jVqk03giivgRz8K49aKSLop2dewI48MXSBfcUXsSEQkNj2grXF/+Qvs\ntx/MmhXGrRWR6qYHtNKkvn1Dn/ennw66t4qkl5J9Cpx7LixZAnffHTsSEYlF1Tgp8cwzcPzxMHcu\ndO0aOxoRaSt1cSwb9YMfhM/rr48bh4i0nZK9bNTy5eHN2ltugUMPjR2NiLRFtAe0ZvZbM1tqZn/O\nW1dnZtPNbKGZTTMzVRxUgK5d4cYb4ZRT4KOPYkcjIuVUjAe0NwODGq0bB0x3952Bx3PLUgEGDYLD\nD4fzzosdiYiUU1GqccysD/CAu++eW14AfN3dl5pZDyDr7v0a7aNqnEg++ihU59x4Y+ghU0SqR6W1\ns+/u7ktz80uB7iU6j7TBlluGAU5OPTXU44tI7etQ6hO4u5tZk0X4+vr6dfOZTIZMJlPqcCTnsMNg\n6FAYPRruvFODlItUqmw2SzabTXycUlbjZNz9bTPrCcxQNU7l+ec/Q1cK55wDo0bFjkZEClFp1ThT\ngJG5+ZHA5BKdRxLYbDOYODH0fb9wYexoRKSUEpfszWwi8HVga0L9/H8CvwMmAdsCi4AT3H15o/1U\nsq8Q11wDN98c3rLt1Cl2NCLSEr1UJW3mHurvd9wR/vd/Y0cjIi1RspdEli2DffaBX/wCjjuuvOeu\nu6yODz79AIBum3Zj2dhl5Q1ApIq0NdmXvDWOVIe6utAr5uDBsMcesPPO5Tv3B59+gF8Qbvx2oZoF\niZSCujiWdfbdFy66CIYNgxUrYkcjIsWkZC8b+P73YcCA8KlaNpHaoWQvGzCD666DBQtC/b2I1AbV\n2cvnbL45TJ4MAwdC//5w1FGxIxKRpFSylyb17g333gsnnxxGtxKR6qZkL80aOBB++UsYMgTefTd2\nNCKShKpxpEUjRsArr8DRR8MTT4QqnlLqtmm3DZpfqt29SHHopSrZKHcYORI+/BDuuw/aty/u8e1C\nW9fOvjXfiaRRpXWEJjXEDH7zm9D2/owz1CRTpBop2UtBOnUKD2yffhouvjh2NCLSWqqzl4JttRU8\n8gh87Wth/owz2n6sxv3hiEhpKdlLq/TsCY89BgcdBF26wPe+17bj5PeHIyKlp2QvrdanD0ybBgcf\nDFtsUf5eMkWk9ZTspU369YOpU2HQoPAAd9iw2BGJSEtK9oDWzAaZ2QIze9XMxpbqPBLPXnuFhP/v\n/x66RxaRylWSkr2ZtQeuBg4D3gRmmdkUd59fivNJeeU/XAVgDJwwD7iw8GPooaxIeZWqGmc/4DV3\nXwRgZncCQwEl+wQ+l2Qj6bZpt889XJ0zB77xDbjkkrY/tG3uXA1v1OptWpG2K1Wy7wW8kbe8BNi/\nROeqGkmTdVNJtlLssQfMmBES/gcfwNlnF+e4+cm9pVGsKuVGKFIqSf/tlyrZFxRVfX39uvlMJkMm\nkylROJWh1psb9usHv/89HHEEvP8+/Pd/h4e3pdA4uVfyjVAkiWw2Szab3SBftkVJ+sYxs4FAvbsP\nyi2fD6x198vytkld3zhp6efl3XfDWLYDBsC110KHIhUpNDC5SNv7xilVsu8AvAIcCvwdmAl8O/8B\nrZJ9bfvHP+DEE0M/OpMmhRewRCS5iuoIzd1XAz8EHgXmAXepJU66dOkCU6bAdtuF7hWWLIkdkUi6\nlaydvbtPdfdd3H1Hd7+kVOeRytWhQxjPdvhwOOAAmDkzdkQi6aVeL6WkzOC88+Dqq+HII+HWW2NH\nJJJO6i5BymLoUNhxx/A5ezZcfnnxHtyKyMapZC9ls9tuoSpn/vxQj//667EjEkkPJXspq7o6eOgh\nOOmkMKD5Lbdo5CuRclCyl7Jr1w7OOisMYD5+PBx2GDzwAKxZEzsykdqlAcfLKE3t7Av12WehHf6v\nfw3LlsGoUbDLLtCrV5g6doTVq8O0Zg2sXRsm91Dn37Hj+qlTpzA1LJfq7V2RmCrqpaqCTqxkL3nc\n4bnn4K67YPFiePPNMK1ZE5J6hw7hfwTt24dPWH8TWLUqTCtXrp9Wrw7bdey4fv9Cpua2b9++8GM0\nt+3G1jf1ff66xtu19F3+fLt2uvHVkrYme7WHkIpgFurwBw4szvHcw41i1ar1/ytYtSp85t8kGpYb\ntmu8Pn+58TEab9dwnobvVq6ETz7ZcLvmjpG/vpDv8s/VMJ9/jvyY3Ju/UTT12TC19H1z80m+a267\nlr5rKAAkmQo5Ri3cLJXspSaZrS/hpl3+ja+pm0JTn83NN7VdS/tt7LuVKzdcbm675qbVq0O1XiHb\ntjS1dIy1a8PvKf+m0NINorntWto//3+tTX2eemroUTYJ/VMQqXG68SXjvv6GWciNoTXzTX3XsJz/\n2bdv8j+H/vpFRFpgtr5k37Fj7GjaTsm+hJrqc11EJAYl+xKq9cFKRKR66KUqEZEUULIXEUkBJXsR\nkRRoc7I3s+PNbK6ZrTGzAY2+O9/MXjWzBWZ2RPIwRUQkiSQPaP8MHAvckL/SzPoDJwL9gV7AY2a2\ns7uvTXAuERFJoM0le3df4O4Lm/hqKDDR3Ve5+yLgNWC/tp5HRESSK0Wd/ReB/OGllxBK+CIiEkmL\n1ThmNh3o0cRXP3X3B1pxniYbm9fX16+bz2QyZDKZVhxSRKT2ZbNZstls4uMk7uLYzGYA57r7C7nl\ncQDufmlu+RHgAnd/rtF+Nd/Fsbo0FpFia2sXx8Wqxsk/8RTgJDPrZGbbAzsBM4t0HhERaYMkTS+P\nNbM3gIHAQ2Y2FcDd5wGTgHnAVOD0mi/Ci4hUOI1UVUKqxhGRYqvKkarswhoY/qUF6uVSRCqFSvYi\nIlUk9gNaERGpYEr2IiIpoGQvIpICSvYiIimgZC8ikgJK9iIiKaBkLyKSAkr2IiIpoGQvIpICSvYi\nIimgZC8ikgJK9iIiKaBkLyKSAkkGL/mFmc03s5fM7D4z2yrvu/PN7FUzW2BmRxQnVBERaaskJftp\nwG7uviewEDgfwMz6AycC/YFBwLVmpv9BtKAYgwnXCl2L9XQt1tO1SK7NSdjdp7v72tzic8A2ufmh\nwER3X+Xui4DXgP0SRVnj9ENeT9diPV2L9XQtkitWiftk4OHc/BeBJXnfLQF6Fek8IiLSBi0OS2hm\n04EeTXz1U3d/ILfNz4CV7n5HC4fSkFQiIhElGpbQzL4HnAYc6u6f5taNA3D3S3PLjwAXuPtzjfbV\nDUBEpA3aMixhm5O9mQ0CxgNfd/f38tb3B+4g1NP3Ah4DdtSAsyIi8bRYjbMRVwGdgOlmBvCsu5/u\n7vPMbBIwD1gNnK5ELyISV6JqHBERqQ4lb/9uZoNyL1e9amZjm9nm17nvXzKzvUsdUywbuxZm1s/M\nnjWzT83s3BgxlksB12J47vcwx8yeNrM9YsRZDgVci6G5a/GimT1vZofEiLMcCskXue2+bGarzexb\n5YyvnAr4XWTM7MPc7+JFM/t5iwd095JNQHtCO/s+QEdgNrBro22+CTycm98f+GMpY4o1FXgt/hXY\nF7gIODd2zJGvxQHAVrn5QSn/XXTOm98deC123LGuRd52TwAPAsNixx3xd5EBphR6zFKX7PfL/TAX\nufsq4E7CS1f5hgATADy02OlqZt1LHFcMG70W7v6uu/8JWBUjwDIq5Fo86+4f5hbzX9qrNYVcixV5\ni1sA71GbCskXAGcA9wDvljO4Miv0WhTcKqfUyb4X8EbeclMvWDW1TS3+wy7kWqRFa6/FKax/aa/W\nFHQtzOwYM5sPTAXOLFNs5bbRa2FmvQhJ77rcqlp96FjI78KBr+Sq+B7OtYRsVpLWOIUo9C+i8d2p\nFv8Ca/HP1FYFXwszO5jwhvaBpQsnqoKuhbtPBiab2deA24BdShpVHIVciyuAce7uFpoBtrq9eZUo\n5Fq8APR290/MbDAwGdi5uY1LXbJ/E+idt9ybDbtSaGqbbXLrak0h1yItCroWuYeyNwJD3P2DMsVW\nbq36Xbj774EOZvYvpQ4sgkKuxT7AnWb2V2AYoaPFIWWKr5w2ei3c/R/u/klufirQ0czqmjtgqZP9\nn4CdzKyPmXUi9IY5pdE2U4ARAGY2EFju7ktLHFcMhVyLBrVaWmmw0WthZtsC9wHfcffXIsRYLoVc\nix1ypVjMbACAu79f9khLb6PXwt37uvv27r49od5+jLs39++omhXyu+ie97vYj9CUfllzByxpNY67\nrzazHwKPEp4u3+Tu881sdO77G9z9YTP7ppm9BqwARpUyplgKuRZm1gOYBWwJrDWzs4D+7v5xtMBL\noJBrAfwn0A24Lvd7XuXuNdd7aoHXYhgwwsxWAR8DJ0ULuIQKvBapUOC1OA4YY2argU/YyO9CL1WJ\niKSABhUREUkBJXsRkRRQshcRSQElexGRFFCyFxFJASV7EZEUULIXEUkBJXsRkRT4f07QtykFLUo+\nAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x115abfb10>"
       ]
      }
     ],
     "prompt_number": 51
    }
   ],
   "metadata": {}
  }
 ]
}