run_plot.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:654e43cc89b0518c833f2758aaeb9afc4c79079368d5fd4e53fc08cccd9c2eb4"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "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": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create the system\n",
      "sysc = LinSystem(2, 1, 1)\n",
      "\n",
      "# Pendulum parameters\n",
      "g = 9.81\n",
      "L = 1\n",
      "m = 0.1\n",
      "B = 1\n",
      "\n",
      "# Continuous time system\n",
      "sysc.A[0, 0] = 0\n",
      "sysc.A[0, 1] = 1\n",
      "sysc.A[1, 0] = -(g/L)\n",
      "sysc.A[1, 1] = -B/(m*L*L)\n",
      "\n",
      "sysc.B[0, 0] = 0\n",
      "sysc.B[1, 0] = 1/(m*L*L)\n",
      "\n",
      "sysc.C[0, 0] = 1\n",
      "sysc.C[0, 1] = 0\n",
      "\n",
      "# Timestep\n",
      "Ts = 0.1\n",
      "\n",
      "# Discretise the system for the MPC controller\n",
      "sys = c2d(sysc, Ts)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Jacobians\n",
      "def Jx(x0):\n",
      "    x1_0 = x0[0]\n",
      "    x2_0 = x0[1]\n",
      "    jac = np.array([[0, 1], [-(g/L)*np.cos(x1_0), -B/(m*L*L)]])\n",
      "    return jac\n",
      "\n",
      "def Ju(u0):\n",
      "    return np.array([[0], [1/(m*L*L)]])\n",
      "\n",
      "def V(x, u, t):\n",
      "    \"\"\"\n",
      "    V = f(x0, u0) - Jx(x0) - Ju(u0)\n",
      "    \"\"\"\n",
      "    f = nlderiv(x, t, u)\n",
      "    return f - Jx(x).dot(x) - Ju(u).dot(u)\n",
      "\n",
      "# True nonlinear system, this function calculates xdot = f(x(t), u(t), t)\n",
      "def nlderiv(x, t, u):\n",
      "    # Find xdot at (x,u) at t\n",
      "    x1 = x[0] # theta\n",
      "    x2 = x[1] # theta dot\n",
      "    xdot1 = x2\n",
      "    xdot2 = -(B*x2)/(m*L*L) - (g/L)*np.sin(x1) + u/(m*L*L)\n",
      "    re = np.array([xdot1, xdot2])\n",
      "    return re"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 27
    },
    {
     "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 = 10\n",
      "mpc.P = 1\n",
      "mpc.R = 0.001\n",
      "mpc.set_prediction_horizon(10)\n",
      "mpc.set_control_horizon(3)\n",
      "target = np.array([[np.pi], [0], [1]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Initial state and input\n",
      "x = np.array([[0.0], [0.0]])\n",
      "u = np.array([[0.0]])\n",
      "\n",
      "# Reset\n",
      "mpc.reset()\n",
      "\n",
      "# The augmented system we will use for the off-eqbm linearised system\n",
      "sysv = LinSystem(sys.order + 1, sys.numinputs, sys.numoutputs)\n",
      "\n",
      "# Pre augment\n",
      "Aa = np.vstack((np.hstack((Jx(x), V(x, u, 0))),\n",
      "                np.hstack((np.zeros([1, sys.order]), np.array([[0]])))))\n",
      "Ba = np.vstack((Ju(u), np.zeros([1, sys.numinputs])))\n",
      "xa = np.vstack((x, np.array([[1]])))\n",
      "\n",
      "# Start at what time?\n",
      "t = 0.0\n",
      "\n",
      "# Storage vectors\n",
      "tall = t\n",
      "uall = u\n",
      "xall = x\n",
      "\n",
      "# Run simulation\n",
      "while t < 10:\n",
      "    # Augment (update Jacobians around the current (x,u) pair)\n",
      "    Aa[0:sys.order, 0:sys.order] = Jx(x)\n",
      "    Aa[0:sys.order, sys.order:sys.order+1] = V(x, u, t)\n",
      "    Ba[0:sys.order, 0:sys.numinputs] = Ju(u)\n",
      "    sysv.A = Aa\n",
      "    sysv.B = Ba\n",
      "\n",
      "    # Discretise the model\n",
      "    sysv = c2d(sysv, Ts)\n",
      "    \n",
      "    # Update plant model in the MPC controller\n",
      "    mpc.set_plant(sysv)\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([[-3, np.inf, 3, np.inf]]),\n",
      "                             np.array([[-10, np.inf, 10, np.inf]]),\n",
      "                             np.array([[-pi, 1e3, pi, 1e3],\n",
      "                                       [-10, np.inf, 10, np.inf],\n",
      "                                       [1, 1e-9, 1, 1e-9]])) # augmented constrained by dynamics\n",
      "    \n",
      "    # Get control inputs\n",
      "    xa[0:sys.order] = x\n",
      "    (u, u_traj, z_traj) = mpc.controlmove(xa)\n",
      "    \n",
      "    # Simulate cont. time. plant - find new state x(k+1) at time t+Ts (second result row)\n",
      "    res = scipy.integrate.odeint(nlderiv, np.ravel(x), np.array([t, t+Ts]), args=(u,))\n",
      "    x = res[1,:].reshape(sys.order, -1)\n",
      "    \n",
      "    # Advance time\n",
      "    t += Ts\n",
      "    \n",
      "    # Store\n",
      "    xall = np.hstack((xall, x))\n",
      "    uall = np.append(uall, u)\n",
      "    tall = np.append(tall, t)\n",
      "    \n",
      "# Results\n",
      "print \"Final state is [%.6f, %.5f]\" % (x[0], x[1])\n",
      "plt.plot(tall, xall[0,:])\n",
      "plt.step(tall, uall)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Final state is [3.133626, 0.00476]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 31,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10fb3f250>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHcNJREFUeJzt3XuUVNWZ9/Hvw0XEIdogBAURiHdcaowEQY0Uii7ENxon\nONG8LoNmojFDRp3RMfqaRZPE+BpN3ogaxWswKmZGjbYRFLyU4mW8AYKCCiojoFy0QYUGufTz/rGr\n6aat7qruPlWn6tTvs9ZZdapq9zlPYfvU7ufsvY+5OyIikgyd4g5ARESio6QuIpIgSuoiIgmipC4i\nkiBK6iIiCaKkLiKSIK0mdTPb2cxeNrN5ZrbQzK7O0iZlZp+Z2dzMdmXhwhURkdZ0ae1Nd99kZqPc\nvc7MugDPm9kx7v58s6bPuvsphQtTRETykbP84u51md2dgM5AbZZmFmVQIiLSPjmTupl1MrN5wCrg\nGXdf2KyJA0eZ2RtmNt3MhhQiUBERyS2fnnq9u38T2As41sxSzZrMAQa4+2HADcDDkUcpIiJ5sbas\n/WJmvwQ2uvt1rbT5ADjC3Wubva5FZkRE2sHd8y5x5xr90tvMqjL73YETgLnN2vQ1M8vsDyN8UWSr\nu+Puid0mTpwYewz6bPp8+nzJ29qq1dEvwJ7AVDPrRPgC+Iu7P2Vm52eS9BRgHHCBmW0F6oAz2hyF\niIhEIteQxgXAt7K8PqXJ/k3ATdGHJiIibaUZpRFJpVJxh1AwSf5soM9X7pL++dqqTRdKO3QiMy/W\nuUREksLM8KgulIqISHlRUhcRSRAldRGRBFFSFxFJECV1EZEEUVIXEUkQJXURkQRRUhcRSRAldRGR\nBFFSFxFJECV1EZEEUVIXEUkQJXURkQRRUhcRSRAldRGRBFFSFxFJkFz3KBURaZU7bNkSts2bG/ez\nbVu3Nj62tG3blv35tm25tyuugKqquP9F4qWkLlLm3GHTJqirgw0bwmPTbePGxseNG0Pbhsfm25df\n7vi4eXPY//LLsN90+/LLxmTdpQt07Rq2nXYKW8Pzpu916bLj886dd3zs0iXsN7Tr3LnxecN+861b\nt8b9Tqo9KKmLxKEhEa9b17h99lnYPv98x+2LL8Lj+vXZtw0bQtL7h38I2y677Lh179742LDtvHN4\nrKoKSXHnncPWsN/wuNNOYb9bt8Zk3bDftWvY79oVLO+brUmh6R6lIh20bRvU1sKaNWH75JOwffpp\n42NtbdjWrm18BOjZMyTW3XZr3HbdtfHxa19rfOzR46uPDYm8i7pnidXWe5QqqYu0YMMGWLECPvoo\nbB9/DCtX7ritXh2S9K67wte/Dn36QO/eYdt998atV6/Gx6qqkMy7d4/7E0o5iDSpm9nOwLNAN2An\n4BF3vzxLu8nASUAdMN7d52Zpo6QuJWPr1pCw/+d/YOlS+PDDsC1bBsuXh23TJujfH/r1C9uee8Ie\nezQ+9u0btt691VOWwom8p25mu7h7nZl1AZ4HLnH355u8PxaY4O5jzexI4Hp3H57lOErqUlQbN8KS\nJWFbvDg8vv8+fPBBSNp9+sCgQTBwYNj23hsGDAhb//6hV61ascStrUk9Z//C3esyuzsBnYHaZk1O\nAaZm2r5sZlVm1tfdV+UbhEhHfPIJvPUWLFwIixbB22/DO+/AqlUweDDsvz/suy9885vw/e/DN74R\nEni3bnFHLhK9nEndzDoBc4B9gJvdfWGzJv2BZU2eLwf2ApTUJVKbN4fEPW8ezJ8ftgULwtC6gw+G\nIUPgoINgzBg44IDQC+/cOe6oRYorn556PfBNM9sNeMLMUu6ebtas+Z8GWess1dXV2/dTqRSpVKot\nsUoF2bw5JO1XX4XXXw/bO++EXvZhh4XthBPgkENCqURlEkmKdDpNOp1u98+3afSLmf0S2Oju1zV5\n7RYg7e73Z56/DYxsXn4xM6c67PfcuSe1lzWv4kglW74cXnwRXnopbAsWhJLJ0KFhO+KIkMA1YkQq\nTaQ1dTPrDWx193Vm1h04AZjUrFkNMAG438yGA+taqqf7xPAFYpPUrapk7vDee5BOw3PPhW39ejj6\naBgxAq65JiTxHj3ijlSk/OQqv+wJTM3U1TsBf3H3p8zsfAB3n+Lu081srJktATYA5xQ2ZClHH38M\nTz4Ztqefhvp6SKVg5Ei4/HI48ECVUESiEMvkI5tk23vtkkxbtoRyyowZYVu2DEaNgtGj4fjjYb/9\nlMRF8hH5kEaRfH3+eUjgNTXhcfBgGDsWbr4Zhg3TBB2RYtD/ZtIha9fCI4/AAw+E2vgxx8Cpp8Lv\nfhdGpYhIcSmpS5tt2BAS+bRpIZEffzz88Idw771hISoRiY+SuuSlvj6MVpk6NST0ESMaE/muu8Yd\nnYg0UFKXVi1fDnfeCXfdFZL3+PGhtNK3b9yRiUg2SuryFfX1MGtWuMD53HNwxhnw4INw+OEasSJS\n6pTUZbsvvgjllRtuCDM3f/YzuOceTQISKSdK6sLHH8P118Ptt4fJQLfdBt/5jnrlIuVIt2mtYO+/\nD+edF1Y43LABXnkllFmOPVYJXaRcKalXoHffhR/9KEwI6ts3PL/hhrACooiUN5VfKsjSpfCrX4UZ\nnxdeGO4EVFUVd1QiEiX11CvAJ5/ARReFlQ/79w+3dvvlL5XQRZJIST3BNm0Ky9geeGC40fLChfDr\nX4c72YtIMqn8kkDuYdbnv/97uLHEiy+G+3SKSPIpqSfM4sUwYUKYCXrLLeGWbyJSOVR+SYgvvwwX\nQUeMgBNPDDdnVkIXqTzqqSfASy/BueeG2vncuTBgQNwRiUhclNTLWF0dXHEF/Od/wuTJMG5c3BGJ\nSNxUfilTr7wSFthaswYWLFBCF5FAPfUys3Ur/OY3YQXFG2+E00+POyIRKSVK6mVk2bJwY4ru3UPt\nvF+/uCMSkVKj8kuZqKmBoUPh5JPh8ceV0EUkO/XUS9y2bWFK/z33wN/+BkcdFXdEIlLKWk3qZjYA\nuBv4OuDAre4+uVmbFPAI8H7mpQfd/TfRh1p51qwJ5RZ3eP116NMn7ohEpNTlKr9sAS5294OB4cC/\nmNlBWdo96+6HZzYl9AgsWBCWxh06FJ54QgldRPLTak/d3VcCKzP7681sEdAPWNSsqW6pEKGaGvjn\nfw53IzrzzLijEZFykveFUjMbBBwOvNzsLQeOMrM3zGy6mQ2JLrzK4g6//324N+jf/66ELiJtl9eF\nUjPrATwAXOju65u9PQcY4O51ZnYS8DCQdU3A6urqsJOG9Mg0qVSqXUEn0bZtYc3zdDpM+9dUf5HK\nlE6nSafT7f55c/fWG5h1Bf4OzHD3P+Y8oNkHwBHuXtvsdW84l00yfGLr560kGzeGC6Kffw4PPQS7\n7RZ3RCJSKswMd8+7xN1q+cXMDLgDWNhSQjezvpl2mNkwwhdFbba28lWffw4nnRQmFM2YoYQuIh2T\nq/xyNHAWMN/M5mZeuwLYG8DdpwDjgAvMbCtQB5xRoFgT55NPYMwYOPLIcOPnTpoKJiIdlLP8EtmJ\nVH7Zwccfw/HHw2mnhbVcTOOHRCSLSMsvUhgffwyjRoU6+lVXKaGLSHSU1Ivso48glYKzz4Yrr4w7\nGhFJGq39UkSrVsFxx8E558AvfhF3NCKSROqpF0ltbbh36BlnKKGLSOEoqRfBF1/A2LEwejRMnBh3\nNCKSZErqBbZ5cxjhcuihcN11uigqIoWlpF5A9fWhfv61r4Xbzymhi0ih6UJpAf3iF7B0KTz5JHTu\nHHc0IlIJlNQL5Kab4NFH4fnnwxIAIiLFoKReADNnhlmiL7wAu+8edzQiUkmU1CO2cCGcdVZYbfEb\n34g7GhGpNLpQGqFPP4XvfjeMcjnmmLijEZFKpKQekW3bwlou//iPYQkAEZE4KKlH5MorYetWuPrq\nuCMRkUqmmnoEHnwQ7rsPXnsNuuhfVERipBTUQYsXw09/Gu5a1KdP3NGISKVT+aUDNm2Cf/onqK6G\noUPjjkZEREm9Qy65BPbZB372s7gjEREJVH5ppwcfhOnTYc4crekiIqVDSb0dli8PvfNHH4Wqqrij\nERFppPJLG9XXw/jx8POfw7BhcUcjIrIjJfU2mjwZ6up09yIRKU0qv7TBm2/CVVfBf/+3xqOLSGlq\ntaduZgPM7Bkze8vM3jSzf22h3WQzW2xmb5jZ4YUJNV5bt4ayy9VXhxEvIiKlKFd/cwtwsbvPM7Me\nwOtmNsvdFzU0MLOxwL7uvp+ZHQncDAyPKsBe1/Ri7aa1Lb7fc+ee1F5WG9XpWnTttdCrF/z4xwU/\nlYhIu7Wa1N19JbAys7/ezBYB/YBFTZqdAkzNtHnZzKrMrK+7r2pvUE0Tec+de+ITvcW2Nqnw4wkX\nLYI//CEsA6DhiyJSyvKuDJvZIOBw4OVmb/UHljV5vhzYC2hTUm9LIm+q5849tyf2QvTat22Dc8+F\nX/0KBg6M9NAiIpHLK6lnSi8PABe6+/psTZo9z5qRq6urw04a0iPTpFKp7e+t3bQ270TeVNMkXohe\n+403QrducP75kR9aROQr0uk06XS63T9v7q0nUjPrCvwdmOHuf8zy/i1A2t3vzzx/GxjZvPxiZt5w\nLptkX0ng2V5rq+a9/Y722lesgMMOC/cZPfDADh1KRKRdzAx3z7vH2mpP3cwMuANYmC2hZ9QAE4D7\nzWw4sC7fenrzJNxRUffaL7wwzBxVQheRcpGr/HI0cBYw38zmZl67AtgbwN2nuPt0MxtrZkuADcA5\n+Z68vSWXYnjsMZg3D/7yl7gjERHJX67RL8+Tx6xTd58QWUQR6cgF1Lo6mDABpkyB7t0LFaGISPQS\nOy+yI6WYa68N66OfeGLUUYmIFFYsSb2hFx1FHT1qH34Y1neZMyfuSERE2i6WpF6MGaDtdemlYQVG\njUkXkXJUEas0Nvxl0OuaXq22S6fDYl3/8R/FiUtEJGqJrak31fCXQWu19W3b4KKLQj19l12KFZmI\nSLQqoqeej3vvDcn89NPjjkREpP1yziiN7ERNZpTGpaUZpxs3wgEHwLRpcPTRcUYoIrKjSGeUJk1L\nwxwnTw5DGJXQRaTcVVRSz+aTT0Id/YUX4o5ERKTjKr6m/tvfhjr6AQfEHYmISMdVbE99+zICu0HV\n13tyM6U7dl5EJF8V21OvvayWC1Y5l6x31n3Z8u3yRETKScUm9aVL4a9/1UQjEUmWik3qv/51WCu9\nT5+4IxERiU5F1tQXL4aamvAIhb/PqYhIsVRkUr/qqrBoV1VVeF7o+5yKiBRLxSX1Dz6ARx+F996L\nOxIRkehVXE39d7+Dn/60sZcuIpIkFdVTX7EijHh5552W2zS9gYdq6yJSbipqQa+LLwYz+MMfcre1\nSVayN8UWkcqhBb1asGYNTJ0Kb74ZdyQiIoVTMTX1m26CceOgX7+4IxERKZyK6KnX1cGf/gSzZ8cd\niYhIYeXsqZvZnWa2yswWtPB+ysw+M7O5me3K6MPsmLvvhhEj2rYSY8MF03zubSoiUipyXig1s+8A\n64G73f2QLO+ngH9z91NyHCeWC6XbtsGBB8Jdd8Exx7TvGLpoKiJxaeuF0pw9dXefDeRaxrBkp2HW\n1MDuu+uuRiJSGaK4UOrAUWb2hplNN7MhERwzMtddB5dcEoYyiogkXRQXSucAA9y9zsxOAh4G9s/W\nsLq6evt+KpUilUpFcPqWvfpqmHB02mkdO44W/BKRYkmn06TT6Xb/fF6Tj8xsEPBotpp6lrYfAEe4\ne22z14teUz/7bDjkELj00uiOqfq6iBRT5DX1PE7Y1ywUN8xsGOGLIvau7OrVYeGuH/847khERIon\nZ/nFzKYBI4HeZrYMmAh0BXD3KcA44AIz2wrUAWcULtz83XprmGzUS6MRRaSCJHLtly1bYPBgmD4d\nDj002mP3uqYXazeFwUCqr4tIoWntF+Bvf4N99ok+oYNuqCEipS2Ra7/ceGO4s5GISKVJXFJfuDDc\ne/TUU+OORESk+BKX1G+7Dc49F7p2jTsSEZHiS9SF0o0bYcCAMOlo8OCCngrQRVMRKbyKvlD64IMw\ndGhxEjrooqmIlJ5ElV9uvRXOOy/uKERE4pOYpN5wgfS73407EhGR+CQmqd9+O5xzTnwXSBsW/dIN\nNUQkTom4ULplC+y1F7zwAuy7b0FOkTct+CUiUSr6gl6lYMYM2H//+BO6iEjcEpHU//xnGD8+7ihE\nROJX9uWXNWtgv/3gww9h110jP3ybaey6iESp4sapT5sWRryUQkIHjV0XkXiVfflFpRcRkUZlndTn\nz4dPP4VRo+KORESkNJR1+eWee+Css6BTiX416YbVIlJsZXuhtL4eBg6EJ56AIUMiO2zBaPy6iLRH\nxYxTnz0bdt+9PBK6iEixlG355b774Ic/jDuK/KkUIyLFUJbll82boV8/mDMH9t47kkMWlUoxIpKv\niii/NNTRyzGhi4gUUlkm9XIrvYiIFEvO8ouZ3QmcDKx290NaaDMZOAmoA8a7+9wsbSIpv6xfH1Zk\nXLIEevfu8OFioaUERCRfhVgm4C7gBuDuFk44FtjX3fczsyOBm4Hh+QbQVo89BiNGlG9CBy0lICKF\nk7P84u6zgbWtNDkFmJpp+zJQZWZ9ownvq/7rv+D00wt1dBGR8hZFTb0/sKzJ8+XAXhEc9ys2bIBZ\ns+DUUwtxdBGR8hfVOPXmNYSsxfPq6urt+6lUilQq1aaTTJ8Ow4eHSUdJ0TB+XbV1EQFIp9Ok0+l2\n/3xe49TNbBDwaLYLpWZ2C5B29/szz98GRrr7qmbtOnyh9Ac/gNGj4Sc/6dBhSpLGrotINnGMU68B\nzs6cfDiwrnlCj0JdHTz+OJx2WtRHFhFJjpzlFzObBowEepvZMmAi0BXA3ae4+3QzG2tmS4ANwDmF\nCPTxx+Hb3y7vUS+t0TICIhKFslkm4MwzIZWC88+PLqZSpVKMiDRI5DIBmzbBjBnwve/FHYmISGkr\ni1Uan34aDjkE+hZs9HtpUSlGRNqrLJL6I49U1th0zTgVkfYq+aReXw81NfDcc3FHIiJS+ko+qb/6\nKvTsCfvtF3ck8VApRkTaouST+iOPVPYFUpViRKQtSn70y8MPV1Y9XUSkI0q6p754MaxbFyYdiUox\nIpJbSSf1Rx6BU06BTiX/90RxqBQjIrmUdLqsqQlJXURE8lOyywSsXQsDB8KqVdC9ewEDK1MNt8RT\nGUYk2QpxO7tYzJwJxx6rhN6ShkSuMoyINFWy5Zfp0+Hkk+OOQkSkvJRk+aW+HvbYA155BQYNKmxc\n5a6hDAMaESOSRIkov7z2GvTpo4SeD42IEZGmSjKpq/TSPhrHLiIlmdQfewyuvTbuKMqPeu0iUnIX\nSletgiVL4Oij445ERKT8lFxP/fHHYfRo6No17kjKm0oxIpWpJJP6mDFxR1H+VIoRqUwlVX6pr4cn\nn4QTT4w7kmRp6LX3uqZX3KGISIGVVE993jzo3RsGDIg7kmTR7FORylFSSX3mTPXSC0l1dpHky1l+\nMbMxZva2mS02s8uyvJ8ys8/MbG5mu7K9wcycCSec0N6fllxqL6vFJzo+0bfPQhWRZGl1mQAz6wy8\nA4wGVgCvAme6+6ImbVLAv7l7q4vk5lomYMMG6NsXVq6EHj3a9BmkHbS8gEh5iHqZgGHAEndfmjn4\n/cCpwKJm7TpcrH32WRg6VAm9WDQ6RiSZcpVf+gPLmjxfnnmtKQeOMrM3zGy6mQ1pTyCzZqn0EpeG\nWrtGyIiUv1w99XyWVZwDDHD3OjM7CXgY2D9bw+rq6u37qVSKVCq1/fnMmfDnP+dxNomceu0ipSOd\nTpNOp9v987lq6sOBancfk3l+OVDv7te08jMfAEe4e22z11usqS9fDocdBqtXQ+fO7fgUEhnV2kVK\nS9Q19deA/cxsEPAR8APgzGYn7Ausdnc3s2GEL4o2ZYKnnoLjjlNCLwXqtYuUt1aTurtvNbMJwBNA\nZ+AOd19kZudn3p8CjAMuMLOtQB1wRluDePppOP74NscuBdZQa1ePXaR8xH7nI/dwg+lZs+CAA4oS\nirSRSjIi8Sm7Ox+99x5s2wb7Z720KqVAJRmR8hF7Un/mGRg1Cky5oixoqQGR0hZ7Un/6aY1PLydN\nk3iva3opwYuUmFhr6u6w557w0ksweHBRwpACsUmGTyzO75JIJSmrmvqiRdC9uxJ6EqgsI1IaYk3q\nDfV0KX/ZyjJK7iLFF2v55fvfh9NOg7POKkoIUmQaCinScW0tv8SW1OvroU8fmD8f+jdfIkwSRwle\npH3KJqnPnx966osXF+X0UkKU4EXyVzYXSp97DkaOjOvsEidNZhIpnNiS+uzZMHZsXGeXUqFRMyLR\niqX84g79+sGLL2o4ozRqWpZpoEQvla4sauqLF4ehjMuWaXkAaZ3q71LpyqKmPns2HHusErrkpmUJ\nRNomlqT+3HMhqYu0hSY4ieQWS/lln32gpgYOPrgop5YEy1aHB/XkJTlKvqa+YkXj/Ug7dSrKqaUC\nqRYvSVHyNfXZs+GYY5TQpbBaqsU3UKKXpCp6Ulc9XYotW/LOluhByV7KX9HLL0OGwNSp8O1vF+W0\nIm2iGr2UmpKuqb/2mjNuXLgvqcovUk5aSvYNlPSlUEo6qZ93njNwIFxxRVFOKVI0SvpSKJEndTMb\nA/wR6Azc7u7XZGkzGTgJqAPGu/vcLG28qspZuDDcwk6kkuRK+qDEL9lFmtTNrDPwDjAaWAG8Cpzp\n7ouatBkLTHD3sWZ2JHC9uw/Pciw/7TTnoYfy/zDlJJ1Ok0ql4g6jIJL82aB0Pl8+iT+bXF8GpfL5\nCiXpny/qIY3DgCXuvjRz8PuBU4FFTdqcAkwFcPeXzazKzPq6+6rmBzvvvHzDKj9J/sVK8meD0vl8\n7e2ltzSSZ7tngFZuG1nufyGUyn+/UpErqfcHljV5vhw4Mo82ewFfSeonnNCOCEWkVbkScrVXUz2x\nusX3c34pRKTcvzzKRa6knu9V1Oa/EVl/rnPnPI8mIkVTrERbsC+PNEyaNGn7U59YnMEfpSpXTX04\nUO3uYzLPLwfqm14sNbNbgLS73595/jYwsnn5xcwq+19aRKSdoqypvwbsZ2aDgI+AHwBnNmtTA0wA\n7s98CazLVk9vS1AiItI+rSZ1d99qZhOAJwhDGu9w90Vmdn7m/SnuPt3MxprZEmADcE7BoxYRkayK\nNvlIREQKr+CT9c1sjJm9bWaLzeyyQp+vmMxsgJk9Y2ZvmdmbZvavccdUCGbW2czmmtmjcccStcwQ\n3AfMbJGZLcyUEBPDzC7P/H4uMLP7zKxb3DG1l5ndaWarzGxBk9d6mdksM3vXzGaaWVWcMXZEC5/v\n2szv5htm9pCZ7ZbrOAVN6pnJSzcCY4AhwJlmdlAhz1lkW4CL3f1gYDjwLwn7fA0uBBaS/2iocnI9\nMN3dDwIOZcc5GGUtcy3sJ8C33P0QQgn1jDhj6qC7CLmkqV8As9x9f+CpzPNyle3zzQQOdvfDgHeB\ny3MdpNA99e2Tl9x9C9AweSkR3H2lu8/L7K8nJIR+8UYVLTPbCxgL3M5Xh66WtUyv5zvufieEa0ju\n/lnMYUXpc0LHYxcz6wLsQpgZXpbcfTbQfMrt9smPmcfvFTWoCGX7fO4+y93rM09fJswBalWhk3q2\niUn9C3zOWGR6RYcT/uGT5P8BlwL1uRqWocHAGjO7y8zmmNltZrZL3EFFxd1rgd8DHxJGr61z9yfj\njSpyTWevrwL6xhlMgZ0LTM/VqNBJPYl/rn+FmfUAHgAuzPTYE8HM/hewOrNAW6J66RldgG8Bf3L3\nbxFGb5Xzn+87MLN9gIuAQYS/IHuY2f+ONagCytwEOZE5x8z+D7DZ3e/L1bbQSX0FMKDJ8wGE3npi\nmFlX4EHgHnd/OO54InYUcIqZfQBMA44zs7tjjilKy4Hl7v5q5vkDhCSfFEOBF939U3ffCjxE+G+a\nJKvMbA8AM9sTWB1zPJEzs/GEEmheX8iFTurbJy+Z2U6EyUs1BT5n0ZiZAXcAC939j3HHEzV3v8Ld\nB7j7YMIFtqfd/ey444qKu68ElpnZ/pmXRgNvxRhS1N4GhptZ98zv6mjCBe8kqQF+lNn/EZCojlVm\n6fNLgVPdfVM+P1PQpJ7pHTRMXloI/LXpsr0JcDRwFjAqM+RvbuY/QlIl8U/bnwP3mtkbhNEvv405\nnsi4+xvA3YTO1fzMy7fGF1HHmNk04EXgADNbZmbnAP8XOMHM3gWOyzwvS1k+37nADUAPYFYmv/wp\n53E0+UhEJDl0p1ARkQRRUhcRSRAldRGRBFFSFxFJECV1EZEEUVIXEUkQJXURkQRRUhcRSZD/D+YO\nFceC3V0rAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10fcb2490>"
       ]
      }
     ],
     "prompt_number": 31
    }
   ],
   "metadata": {}
  }
 ]
}