Inverted pendulum on a cart.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b4d990cb",
   "metadata": {},
   "source": [
    "# Inverted pendulum on a cart\n",
    "An easy example of controlling an _inverted pendulum on a cart_ from Steve Brunton book and YT lecture [1,2]. We consider _Euler–Lagrange equations_ as our pendulum model with the full nonlinear dynamics:\n",
    "\n",
    "$\\dot{x} = v$  \n",
    "$\\dot{v} = \\frac{-m^2L^2 g cos(\\theta)sin(\\theta) + mL^2(mL\\omega^2 \\omega^2sin(\\theta)-\\delta v) +mL^2u}{mL^2(M+m(1-cos(\\theta)^2))}$  \n",
    "$\\dot{\\theta} = \\omega$  \n",
    "$\\dot{\\omega} = \\frac{(m+M)mgLsin(\\theta) - mLcos(\\theta)(mL\\omega^2sin(\\theta)-\\delta v) + mLcos(\\theta)u}{mL^2(M+m(1-cos(\\theta)^2))}$\n",
    "\\\\\n",
    "\n",
    "This dynamics can be linarized into:  \n",
    "$\n",
    "\\frac{d}{dt} \\begin{bmatrix}\n",
    "x_1\\\\\n",
    "x_2\\\\\n",
    "x_3\\\\\n",
    "x_4\n",
    "\\end{bmatrix} = \n",
    "\\begin{bmatrix}\n",
    "0 & 1 & 0 & 0 \\\\\n",
    "0 & -\\frac{\\omega}{M} & b\\frac{mg}{M} & 0\\\\\n",
    "0 & 0 & 0 & 1\\\\\n",
    "0 & -b\\frac{\\delta}{ML} & -b\\frac{(m + M)g}{ML} & 0\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "x_1\\\\\n",
    "x_2\\\\\n",
    "x_3\\\\\n",
    "x_4\n",
    "\\end{bmatrix}\n",
    "+\n",
    "\\begin{bmatrix}\n",
    "0\\\\\n",
    "\\frac{1}{M}\\\\\n",
    "0\\\\\n",
    "b\\frac{1}{ML}\n",
    "\\end{bmatrix}\n",
    "u\n",
    "$,  \n",
    "Where $m$ is pendulum mass, $M$ cart mass, $L$ pendulum length, $g$ gravitational acceleration, and $\\delta$ cart damping. This system then represents matrix system:\n",
    "\n",
    "$$\\frac{dx}{dt} = Ax + Bu$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c153eafd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import control\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "53e17cab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2AAAAHWCAYAAAARnurlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApbUlEQVR4nO3de7RXdYH//9cbj4CAyWo0EylrdfnijDPdDo1mS40DKN7R1FTMxMRJx1taWV+72c0a76Z9JcPRUBwt8X4BQXDG0eZgdjO1y3RDzTQXZhAQsn9/HOpnJXKAz9n7XB6PtVwj+Dl7v+azjnSe7s9nf0pVVQEAAKDnDWp6AAAAwEAhwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGrS1p0HlVJ+nuS5JM8nWVVVVXtPjgIAAOiPuhVga7yrqqqne2wJAABAP+cliAAAADXpboBVSeaUUh4opUzryUEAAAD9VXdfgvjOqqoeK6W8IsncUsojVVXd88IHrAmzaUkyfPjwt40ZM6bFUwEAAPqGBx544Omqqrb6698vVVWt14FKKZ9K8vuqqs5e22Pa29urRYsWrfdIAACA/qCU8sCL3bxwnS9BLKUML6Vs/qe/TzIxyQ9aPxEAAKB/685LELdOMruU8qfHX11V1R09ugoAAKAfWmeAVVX1v0neVMMWAACAfs1t6AEAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGoiwAAAAGrS7QArpWxSSnmwlHJLTw4CAADor9bnCthJSR7uqSEAAAD9XbcCrJQyOsleSS7r2TkAAAD9V3evgJ2f5MNJVvfcFAAAgP5tnQFWStk7yW+qqnpgHY+bVkpZVEpZ9NRTT7VsIAAAQH/RnStgOyfZt5Ty8yTXJBlXSpn51w+qqmp6VVXtVVW1b7XVVi2eCQAA0PetM8CqqvpoVVWjq6p6TZL3JJlfVdWUHl8GAADQz/gcMAAAgJq0rc+Dq6pakGRBjywBAADo51wBAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqIkAAwAAqMk6A6yUMrSU8j+llO+WUh4qpXy6jmEAAAD9TVs3HrMiybiqqn5fStk0yX+VUm6vqur+Ht4GAADQr6wzwKqqqpL8fs0vN13zV9WTowAAAPqj7lwBSyllkyQPJHl9kourqvpWj64CAIANsWRJct99SWdnsnBh8vjjycqVyeDByahRya67JmPHJjvtlIwc2fRaBqBuBVhVVc8neXMpZWSS2aWUHaqq+sELH1NKmZZkWpK8+tWvbvVOAABYu0WLknPOSW64oSu2li1LVq36y8c88khyzz3JsGFdUbb//smppybt7U0sZoAqXa8wXI8vKOUTSZZVVXX22h7T3t5eLVq0aGO3AQDAS1u8OJkypeuK1/LlyerV3f/aQYOSoUO7rojNnJmMHt1zOxlwSikPVFX1N3XfnbsgbrXmyldKKZslmZDkkZYvBACA7qqqZMaMZMyY5N57u654rU98JV2PX7as6+vHjOk63npenID11Z3PAdsmyd2llO8l6Uwyt6qqW3p2FgAArEVVJaeckpxwQrJ06d++1HB9rVrVdZwTTug6rgijB3XnLojfS/KWGrYAAMBL+1N8XXZZ19WrVlq2rOu4pSTnndfaY8Ma3bkCBgAAvcPllydf/WrXFauesHRpMn1618sRoQcIMAAA+obFi5MTT2z9la+/tmxZctJJyWOP9ex5GJAEGAAAfcOUKcmKFfWca/ny5PDD6zkXA4oAAwCg9+vs7PprY2+40V2rVnWdz0cr0WICDACA3u+cc7quStVp+fLk3HPrPSf9ngADAKB3W7IkufHG9f+cr421enUye3bX+aFFBBgAAL3bffclgwc3c+4hQ5L772/m3PRLAgwAgN6ts7Pn73y4NkuXdp0fWkSAAQDQuy1cWN/NN/7aqlXJggXNnJt+SYABANC7Pf74Bn/p75JcsbHnf+KJjT0C/JkAAwCgd1u5coO+bGmSvZK8P8mjG3P+uj57jAGhrekBAADwkjbgBhzLk+yX5L+TzEryfzbm/EOGbMxXw19wBQwAgN5t1Kj1evjKJO9OMi/J5UkO3tjzb7PNxh4B/kyAAQDQu+26a9LWvRdurUpyWJJbk/y/JO/d2HO3tSW77baxR4E/E2AAAPRuY8cmw4at82HPJzkyyTeTnJfk2Face/jwrvNDiwgwAAB6t512WueNOFYn+ZckVyf5fJKTW3XuFSuSHXds1dFAgAEA0MuNHJnst18y6MV/dK2SnJTksiRnJPloq847aFAyeXLX+aFFBBgAAL3faaclQ4f+zW9XST6S5MtJTk1yZivPOXRocuqprTwiCDAAAPqA9vau92L91c04zkzyb0k+sOb/lladr62t63xve1urjghJBBgAAH3FzJl/8ZlcX0ryqSRHpesKWMviK+m6+nXVVa08IiQRYAAA9BWjRycXXpgMG5aL0vXSw/ck+Wpa/EPtsGHJBRck227byqNCEgEGAEBfctRR+epOO+XEJJOTXJlkk1Yef/jwZNq0ZOrUVh4V/kyAAQDQZ8y86qocO39+Jm23XWZttlk2beXBhw1LjjkmOffcVh4V/oIAAwCgT/jGN76RI488Mrvttlu++cMfZsiXv5yMGPE3N+ZYb21tXce56KLkvPOS0tJ3k8FfEGAAAPR6t9xySw499NDstNNOuemmm7LZsGFdLxN8+OFk5527rl6t5XPC1mrQoK6v23nn5JFHvOyQWggwAAB6tblz5+bAAw/Mm9/85tx6660ZMWLE//8PR49OFixIFi5MDjmk6+6FW2yx9qtibW1d/3zo0K7HL1zY9fVuuEFNNvJ6LQAA9Jx77rkn++23X8aMGZM777wzW2yxxYs/sL09ufrqZMmS5P77k87OrrB64olkxYqu29dvs02y225dn++1447JyJH1/T8Ca5Sqqlp+0Pb29mrRokUtPy4AAAPHt771rYwfPz6jR4/OwoUL84pXvKLpSdBtpZQHqqpq/+vf9xJEAAB6nQcffDB77LFHtt5668ybN0980W8IMAAAepWHHnooEydOzOabb5558+Zl1KhRTU+ClhFgAAD0Gj/+8Y8zfvz4bLrpppk/f3622267pidBS7kJBwAAvcLPf/7zjBs3Ls8//3wWLlyY17/+9U1PgpYTYAAANO6xxx7LuHHj8vvf/z4LFizI9ttv3/Qk6BECDACARj355JPp6OjI008/nXnz5uVNb3pT05OgxwgwAAAa89vf/jbjx4/Pr371q9x5550ZO3Zs05OgRwkwAAAasWTJkkycODE//vGPc+utt+ad73xn05OgxwkwAABq99xzz2XPPffM97///dxwww3p6OhoehLUQoABAFCrZcuWZd99983//M//5Nprr82ee+7Z9CSojQADAKA2K1asyOTJk7Nw4cLMnDkzBxxwQNOToFYCDACAWvzxj3/MwQcfnDlz5uRrX/taDjvssKYnQe0GNT0AAID+7/nnn8+UKVNy00035ctf/nKmTp3a9CRohAADAKBHrV69OlOnTs21116bf/u3f8vxxx/f9CRojAADAKDHVFWV4447LldeeWXOPPPMnHbaaU1PgkYJMAAAekRVVfngBz+YSy+9NKeffnrOOOOMpidB4wQYAAA94owzzsj555+fE088MZ///OdTSml6EjROgAEA0HKf+9zn8vnPfz7Tpk3L+eefL75gDQEGAEBLnXvuuTnjjDNyxBFH5Ctf+Yr4ghcQYAAAtMxXvvKVnHrqqTnooIMyY8aMDBrkx014If9GAADQEv/+7/+e4447Lvvss09mzpyZtra2pidBryPAAADYaNdcc02OPvroTJgwIddee20GDx7c9CTolQQYAAAb5YYbbsiUKVPyzne+MzfccEOGDh3a9CTotQQYAAAb7I477sjBBx+c9vb23HLLLRk2bFjTk6BXE2AAAGyQ+fPnZ/Lkydlhhx1yxx13ZPPNN296EvR6AgwAgPV27733Zt99983rXve6zJkzJyNHjmx6EvQJAgwAgPWyaNGi7Lnnnhk1alTuuuuubLnllk1Pgj5DgAEA0G3f+973MnHixLz85S/PvHnz8spXvrLpSdCnCDAAALrlkUceyfjx4zNs2LDMnz8/r3rVq5qeBH2OAAMAYJ1++tOfpqOjI4MGDcr8+fPz2te+tulJ0Cf5eHIAAF7SL3/5y4wbNy4rVqzIggUL8sY3vrHpSdBnCTAAANbq8ccfz7hx4/Lss89m/vz52WGHHZqeBH2aAAMA4EU99dRTGT9+fJ588snMnTs3b33rW5ueBH2eAAMA4G8888wzmTBhQn7+85/n9ttvz4477tj0JOgXBBgAAH/hd7/7XfbYY488/PDDufnmm7Prrrs2PQn6DQEGAMCfLV26NHvttVcefPDBXH/99Zk4cWLTk6BfEWAAACRJ/vCHP2S//fbLf//3f+eaa67JPvvs0/Qk6HcEGAAAWblyZd797ndn/vz5ueKKK3LQQQc1PQn6JQEGADDArVq1Koceemhuu+22XHrppTniiCOangT91qCmBwAA0Jznn38+Rx55ZK6//vqcf/75mTZtWtOToF8TYAAAA9Tq1atz7LHH5uqrr84XvvCFnHTSSU1Pgn5PgAEADEBVVeXEE0/M1772tXz84x/P6aef3vQkGBAEGADAAFNVVT784Q/n4osvzmmnnZZPf/rTTU+CAUOAAQAMMJ/+9Kdz9tln57jjjsuXvvSllFKangQDhgADABhAvvjFL+bTn/50pk6dmosuukh8Qc0EGADAAHHhhRfm9NNPz6GHHprp06dn0CA/CkLd/FsHADAAfPWrX81JJ52UyZMn54orrsgmm2zS9CQYkAQYAEA/N3PmzBx77LGZNGlSZs2alU033bTpSTBgCTAAgH7suuuuy5FHHpl3vetd+eY3v5khQ4Y0PQkGNAEGANBP3XLLLTnssMOy00475cYbb8xmm23W9CQY8NYZYKWUV5VS7i6l/LCU8lApxUekAwD0cnPnzs2BBx6YN7/5zbn11lszYsSIpicBSdq68ZhVSU6tqurbpZTNkzxQSplbVdUPe3gbAAAb4J577sl+++2XMWPG5M4778wWW2zR9CRgjXVeAauq6omqqr695u+fS/Jwkm17ehgAAOvvW9/6Vvbaa69st912mTt3bl7+8pc3PQl4gfV6D1gp5TVJ3pLkWz2yBgCADfbggw9mjz32yNZbb5158+blFa94RdOTgL/S7QArpYxI8s0kJ1dV9bsX+efTSimLSimLnnrqqVZuBABgHR566KFMmDAhL3vZyzJv3ryMGjWq6UnAi+hWgJVSNk1XfF1VVdX1L/aYqqqmV1XVXlVV+1ZbbdXKjQAAvIQf/ehH6ejoyODBgzNv3rxst912TU8C1mKdN+EopZQkX0vycFVV5/b8JAAAuutnP/tZOjo6snr16tx99915/etf3/Qk4CV05wrYzkmOSDKulPKdNX/t2cO7AABYh8WLF6ejoyNLly7N3Llzs/322zc9CViHdV4Bq6rqv5KUGrYAANBNv/71r9PR0ZHf/va3ueuuu/KmN72p6UlAN3Tnc8AAAOhFnn766UyYMCGLFy/OnDlzMnbs2KYnAd0kwAAA+pAlS5Zk9913z49//OPcdttt2XnnnZueBKwHAQYA0Ec899xzmTRpUr7//e/nxhtvzLhx45qeBKwnAQYA0AcsW7Ys++yzTzo7O3Pddddl0qRJTU8CNoAAAwDo5VasWJHJkyfnnnvuyVVXXZXJkyc3PQnYQAIMAKAX++Mf/5iDDz44c+bMyYwZM3LooYc2PQnYCN35HDAAABqwatWqHH744bnpppty8cUX56ijjmp6ErCRBBgAQC+0evXqHH300bnuuuty9tln57jjjmt6EtACAgwAoJepqirHHXdcrrzyypx55pk59dRTm54EtIgAAwDoRaqqyimnnJJLL700H/3oR3PGGWc0PQloIQEGANCLnHHGGbngggty0kkn5XOf+1xKKU1PAlpIgAEA9BKf/exn8/nPfz7Tpk3LeeedJ76gHxJgAAC9wDnnnJOPf/zjOeKII/KVr3xFfEE/JcAAABp2ySWX5LTTTstBBx2UGTNmZNAgP6JBf+XfbgCABl1++eU5/vjjs+++++aqq65KW1tb05OAHiTAAAAaMmvWrBx99NGZOHFi/uM//iObbrpp05OAHibAAAAaMHv27BxxxBHZZZddMnv27AwdOrTpSUANBBgAQM1uv/32HHLIIRk7dmxuvvnmDBs2rOlJQE0EGABAjebPn58DDjggO+ywQ26//fZsvvnmTU8CaiTAAABqcu+992afffbJ61//+syZMycjR45sehJQMwEGAFCDzs7OTJo0KaNHj87cuXOz5ZZbNj0JaIAAAwDoYd/73vey++67Z8stt8y8efPyyle+sulJQEMEGABAD3r44Yczfvz4DB8+PPPmzcvo0aObngQ0SIABAPSQn/zkJ+no6MigQYMyb968vPa1r216EtAwH7UOANADfvnLX6ajoyMrV67MggUL8sY3vrHpSUAvIMAAAFrs8ccfz7hx4/Lss89m/vz52WGHHZqeBPQSAgwAoIV+85vfpKOjI08++WTmzp2bt771rU1PAnoRAQYA0CLPPPNMJkyYkF/84he54447suOOOzY9CehlBBgAQAs8++yz2X333fPII4/klltuyS677NL0JKAXEmAAABtp6dKl2WuvvfKd73wn119/fSZMmND0JKCXEmAAABvhD3/4Q/bdd9/cd999ueaaa7LPPvs0PQnoxQQYAMAGWrFiRQ488MDcfffdufLKK3PQQQc1PQno5QQYAMAGWLVqVQ499NDcfvvtufTSSzNlypSmJwF9wKCmBwAA9DXPP/983vve92b27Nm54IILMm3atKYnAX2EAAMAWA+rV6/OMccck1mzZuWss87KiSee2PQkoA8RYAAA3VRVVU488cRcfvnl+cQnPpGPfOQjTU8C+hgBBgDQDVVV5cMf/nAuvvjinHbaafnUpz7V9CSgDxJgAADd8KlPfSpnn312jj/++HzpS19KKaXpSUAfJMAAANbhrLPOyplnnpmpU6fmwgsvFF/ABhNgAAAv4YILLshHP/rRHHbYYZk+fXoGDfLjE7Dh/AkCALAW06dPz8knn5wDDjggV1xxRTbZZJOmJwF9nAADAHgRX//61/Mv//Iv2XPPPTNr1qy0tbU1PQnoBwQYAMBfue666/K+970v73rXu/KNb3wjgwcPbnoS0E8IMACAF7j55ptz2GGH5R3veEduuummbLbZZk1PAvoRAQYAsMacOXPy7ne/O295y1ty6623Zvjw4U1PAvoZAQYAkGThwoXZf//9s/322+eOO+7Iy172sqYnAf2QAAMABrz7778/e++9d17zmtdkzpw5efnLX970JKCfEmAAwID27W9/O3vssUe23nrr3HXXXXnFK17R9CSgHxNgAMCA9YMf/CATJ07MFltskfnz52fUqFFNTwL6OQEGAAxIP/rRjzJ+/PgMHjw48+fPz6tf/eqmJwEDgE8UBAAGnJ/97Gfp6OjI6tWrc/fdd+d1r3td05OAAUKAAQADyuLFizNu3LgsXbo0CxYsyPbbb9/0JGAAEWAAwIDx61//Oh0dHXnmmWcyb968/NM//VPTk4ABRoABAAPC008/nfHjx+exxx7LnXfemfb29qYnAQOQAAMA+r0lS5Zk4sSJ+elPf5pbb701O++8c9OTgAFKgAEA/dpzzz2XSZMm5Qc/+EFuvPHGjBs3rulJwAAmwACAfmvZsmXZe++909nZmeuuuy6TJk1qehIwwAkwAKBfWr58efbff//853/+Z66++upMnjy56UkAAgwA6H9WrlyZgw8+OHPnzs3ll1+e97znPU1PAkiSDGp6AABAK61atSpTpkzJzTffnEsuuSTve9/7mp4E8GcCDADoN1avXp2pU6fmuuuuyznnnJMPfOADTU8C+AsCDADoF6qqygc+8IF8/etfz2c+85l88IMfbHoSwN8QYABAn1dVVU4++eRMnz49H/vYx3LGGWc0PQngRQkwAKBPq6oqH/vYx3LhhRfm5JNPzmc/+9mmJwGslQADAPq0z372sznrrLNy7LHH5txzz00ppelJAGslwACAPuvss8/OJz7xibz3ve/NJZdcIr6AXk+AAQB90iWXXJIPfehDOfjgg/O1r30tgwb5sQbo/fxJBQD0OTNmzMjxxx+ffffdNzNnzkxbW1vTkwC6RYABAH3KrFmz8v73vz8TJ07Mtddem0033bTpSQDdJsAAgD5j9uzZOeKII7LLLrtk9uzZGTJkSNOTANaLAAMA+oTbbrsthxxySN7+9rfn5ptvzrBhw5qeBLDeBBgA0OvNmzcvBxxwQP7xH/8xt912WzbffPOmJwFsEAEGAPRq9957b/bdd9+84Q1vyJw5czJy5MimJwFsMAEGAPRanZ2dmTRpUkaPHp277rorf/d3f9f0JICNIsAAgF7pu9/9bnbfffdsueWWmTdvXrbeeuumJwFsNAEGAPQ6P/zhDzNhwoQMHz488+fPz+jRo5ueBNASAgwA6FV+8pOfZPz48dlkk00yf/78vOY1r2l6EkDL+Nh4AKDX+MUvfpGOjo6sXLkyCxcuzBve8IamJwG01DqvgJVSZpRSflNK+UEdgwCAgemxxx5LR0dHfve732Xu3Ln5h3/4h6YnAbRcd16C+O9J9ujhHQDAAPab3/wm48ePz5NPPpk77rgjb3nLW5qeBNAj1hlgVVXdk+SZGrYAAAPQM888kwkTJuQXv/hFbr311vzzP/9z05MAeoz3gAEAjXn22Wez++6759FHH83NN9+cXXbZpelJAD2qZXdBLKVMK6UsKqUseuqpp1p1WACgn/r973+fPffcM9/5znfyjW98IxMmTGh6EkCPa1mAVVU1vaqq9qqq2rfaaqtWHRYA6If+8Ic/ZL/99sv999+fWbNmZe+99256EkAtvAQRAKjVihUrcuCBB+buu+/OlVdemXe/+91NTwKoTXduQz8ryX1J/k8pZXEp5eienwUA9Ed//OMfc+ihh+b222/PpZdemilTpjQ9CaBW67wCVlXVoXUMAQD6t+effz5HHnlkZs+enQsuuCDHHHNM05MAatey94ABAKzN6tWrc8wxx2TWrFk566yzcuKJJzY9CaARAgwA6FFVVeWEE07I5Zdfnk9+8pP5yEc+0vQkgMYIMACgx1RVlQ996EO55JJL8qEPfSif/OQnm54E0CgBBgD0mE9+8pM555xz8q//+q/54he/mFJK05MAGiXAAIAe8YUvfCGf+cxncvTRR+eCCy4QXwARYABADzj//PPzsY99LIcddlguvfTSDBrkRw6ARIABAC126aWX5pRTTsmBBx6YK664IptssknTkwB6DQEGALTMlVdemQ984APZc889c/XVV6etbZ0fOQowoAgwAKAlrr322hx11FEZN25cvvnNb2bw4MFNTwLodQQYALDRbrrpphx++OF5xzvekRtvvDFDhw5tehJAr+R1AQD0W266V5c5SQ5K8pb813/dmhEjhjc9qFeqqqYXAL2BK2AAwEZYmGT/JNsnuSPJyxpdA9DbCTAAYAPdl2SvJK9JMjfJyxtdA9AXCDAAYAN8O8mkJNskmZdkq2bnAPQRAgwAWE8/SDIxych0xdc2ja4B6EsEGACwHh5NMj7JkHTF16ubnQPQx7gLIgDQTT9L0pFkdZIFSV7X6BqAvkiAAQDd8Ksk45L8IcndScY0OwegjxJgAMA6/DpdV76eSdfLDv+p2TkAfZgAAwBewtPpes/X4+n6wOX2ZucA9HECDABYiyVJJiT5aZLbkryj0TUA/YG7IAIAa7FZut7rNTvJuxreAtA/uAIGAKzFkCSzmh4B0K+4AgYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFATAQYAAFCTtqYHAEBPqaqmFwDAX3IFDAAAoCYCDAAAoCYCDAAAoCbeAwYANViyZEnuu+++dHZ2ZuHChXn88cezcuXKDB48OKNGjcquu+6asWPHZqeddsrIkSObntvreT6BvqpUPfAO5fb29mrRokUtPy4A9DWLFi3KOeeckxtuuCGDBw/OsmXLsmrVqr95XFtbW4YNG5aVK1dm//33z6mnnpr29vYGFvdunk+gryilPFBV1d/8wSPAAKAHLF68OFOmTElnZ2eWL1+e1atXd/trBw0alKFDh2bs2LGZOXNmRo8e3YNL+wbPJ9DXrC3AvAcMAFqoqqrMmDEjY8aMyb333ptly5atVywkyerVq7Ns2bLce++9GTNmTGbMmJGe+A+mfYHnE+hvBBgAtEhVVTnllFNywgknZOnSpS/60rj1sWrVqixdujQnnHBCTjnllAEXDZ5PoD8SYADQAn+KhcsuuyzLli1r6bGXLVuWyy67LB/84AdbetzezPMJ9FcCDABa4PLLL89Xv/rVLF26tEeOv3Tp0kyfPj0zZszokeP3Np5PoL9yEw4A2EiLFy/OmDFjeiwWXmjEiBF55JFHsu222/b4uZri+QT6AzfhAIAeMmXKlKxYsaKWcy1fvjyHH354LedqiucT6M8EGABshM7OznR2dm70DSK6a9WqVens7Ex/faWJ5xPo7wQYAGyEc845J8uXL6/1nMuXL8+5555b6znr4vkE+jvvAQOADbRkyZJss802tQdDkgwdOjRPPPFERo4cWfu5e4rnE+hPvAcMAFrsvvvuy+DBgxs595AhQ3L//fc3cu6e4vkEBgIBBgAbqLOzs+WfUdVdS5cuTWdnZyPn7imeT2AgEGAAsIEWLlxY280i/tqqVauyYMGCRs7dUzyfwEAgwABgAz3++OONnv+JJ55o9Pyt5vkEBgIBBgAbaOXKlY2ev67PyqqL5xMYCAQYAGygpm4Y8SdDhgxp9Pyt5vkEBgIBBgAbaNSoUY2ef5tttmn0/K3m+QQGAgEGABto1113TVtbWyPnbmtry2677dbIuXuK5xMYCAQYAGygsWPHZtiwYY2ce/jw4Rk7dmwj5+4pnk9gIBBgALCBdtppp8ZuHLFixYrsuOOOjZy7p3g+gYFAgAHABho5cmT222+/DBpU7/+cDho0KJMnT87IkSNrPW9P83wCA4EAA4CNcNppp2Xo0KG1nnPo0KE59dRTaz1nXTyfQH8nwABgI7S3t2fs2LG13Tyira0tY8eOzdve9rZazlc3zyfQ3wkwANhIM2fOrO0zpIYOHZqrrrqqlnM1xfMJ9GcCDAA20ujRo3PhhRf2+B38hg0blgsuuCDbbrttj56naZ5PoD8TYADQAkcddVSOOeaYDB8+vEeOP3z48EybNi1Tp07tkeP3Np5PoL8SYADQAqWUnHfeeXn/+9/f8is3w4YNyzHHHJNzzz23pcftzTyfQH8lwACgRf4UDRdddFFGjBix0TeSaGtry4gRI3LRRRflvPPOSymlRUv7Bs8n0B8JMABooVJKpk6dmocffjg777xzhg0btt6fazVo0KAMGzYsO++8cx555JEB/TI5zyfQ3wgwAOgBo0ePzoIFC7Jw4cIccsghGTp0aLbYYou1XsVpa2vLFltskaFDh+aQQw7JwoULs2DBAjeIWMPzCfQXpaqqlh+0vb29WrRoUcuPCwB91ZIlS3L//fens7MzCxYsyBNPPJEVK1ZkyJAh2WabbbLbbrtl7Nix2XHHHTNy5Mim5/Z6nk+gtyulPFBVVfvf/L4AAwAAaK21BVi3XoJYStmjlPJoKeUnpZTTWz8PAACg/1tngJVSNklycZJJSf4+yaGllL/v6WEAAAD9TXeugL09yU+qqvrfqqpWJrkmyX49OwsAAKD/6c4Hamyb5Fcv+PXiJP/8Ul/w6KOPZrfddtuIWQAAAP3Pxn2i4QuUUqYlmbbml79fuHDho6069gC3ZZKnmx4Ba+H7k97K9ya9le9NejPfn6213Yv9ZncC7LEkr3rBr0ev+b2/UFXV9CTTN2gaa1VKWfRid0+B3sD3J72V7016K9+b9Ga+P+vRnfeAdSZ5QynltaWUwUnek+Smnp0FAADQ/6zzClhVVatKKf+a5M4kmySZUVXVQz2+DAAAoJ/p1nvAqqq6LcltPbyFF+dlnfRmvj/prXxv0lv53qQ38/1Zg1JVVdMbAAAABoTuvAcMAACAFhBgfUgp5dRSSlVK2bLpLZAkpZR/K6U8Ukr5XilldillZNOboJSyRynl0VLKT0oppze9B5KklPKqUsrdpZQfllIeKqWc1PQmeKFSyiallAdLKbc0vaW/E2B9RCnlVUkmJvll01vgBeYm2aGqqn9K8qMkH214DwNcKWWTJBcnmZTk75McWkr5+2ZXQZJkVZJTq6r6+yQ7Jjne9ya9zElJHm56xEAgwPqO85J8OIk37dFrVFU1p6qqVWt+eX+6PicQmvT2JD+pqup/q6pameSaJPs1vAlSVdUTVVV9e83fP5euH3S3bXYVdCmljE6yV5LLmt4yEAiwPqCUsl+Sx6qq+m7TW+AlTE1ye9MjGPC2TfKrF/x6cfyQSy9TSnlNkrck+VbDU+BPzk/Xf+hf3fCOAaFbt6Gn55VS7kryyhf5R/83ycfS9fJDqN1LfW9WVXXjmsf833S9vOaqOrcB9DWllBFJvpnk5Kqqftf0Hiil7J3kN1VVPVBK2a3hOQOCAOslqqoa/2K/X0r5xySvTfLdUkrS9RKvb5dS3l5V1a9rnMgAtbbvzT8ppbwvyd5JOiqfa0HzHkvyqhf8evSa34PGlVI2TVd8XVVV1fVN74E1dk6ybyllzyRDk7yslDKzqqopDe/qt3wOWB9TSvl5kvaqqp5ueguUUvZIcm6SXauqeqrpPVBKaUvXDWE60hVenUkOq6rqoUaHMeCVrv+KekWSZ6qqOrnhOfCi1lwBO62qqr0bntKveQ8YsDG+nGTzJHNLKd8ppfy/pgcxsK25Kcy/JrkzXTc5uFZ80UvsnOSIJOPW/Hn5nTVXHIABxhUwAACAmrgCBgAAUBMBBgAAUBMBBgAAUBMBBgAAUBMBBgAAUBMBBgAAUBMBBgAAUBMBBgAAUJP/DzX/eiljCs6iAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def cart(x=0, theta=np.pi, i=-1):\n",
    "    plt.figure(figsize=(15, 8))\n",
    "    plt.xlim(-5, 5)\n",
    "    plt.ylim(0, 5)\n",
    "    plt.scatter(-.5 + x, .25, s=1000, c=\"k\")  # left wheel\n",
    "    plt.scatter(.5 + x, .25, s=1000, c=\"k\")  # right wheel\n",
    "    plt.plot([-5, 5], [0, 0], \"k-\")  # x axis\n",
    "    plt.plot([-.8 + x, .8 + x], [.65, .65], \"b-\", linewidth=30)  # cart body\n",
    "    plt.plot([0 + x, -3 * np.sin(-theta) + x], [0.7, -3 * np.cos(-theta) + .7], \"k-\")  # pendulum arm\n",
    "    plt.scatter(-3 * np.sin(-theta) + x, -3 * np.cos(-theta) + .7, s=1000, c=\"r\")  # pendulum\n",
    "    if i < 0:\n",
    "        plt.show()\n",
    "    else:\n",
    "        plt.savefig(f\"img/img_{i}.png\", dpi=120, bbox_inches=\"tight\")\n",
    "    \n",
    "\n",
    "cart(1, 3 * np.pi / 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "51d9d883",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_a_matrix(pend_mass=1, cart_mass=5, arm_length=2, g=-10, d=1, b=1):\n",
    "    return np.array([\n",
    "        [0, 1, 0, 0],\n",
    "        [0, -d / cart_mass, b * pend_mass * g / cart_mass, 0],\n",
    "        [0, 0, 0, 1],\n",
    "        [0, -b * d / (cart_mass * arm_length), -b * (pend_mass + cart_mass) * g / (cart_mass * arm_length), 0]\n",
    "    ])\n",
    "\n",
    "\n",
    "def get_b_matrix(pend_mass=1, cart_mass=5, arm_length=2, g=-10, d=1, b=1):\n",
    "    return np.array([[0, 1 / cart_mass, 0, b / (cart_mass * arm_length)]]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0e1478f8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigen values of matrix A: [0.0, -2.4311230004720983, -0.23363938083350477, 2.464762381305604]\n",
      "Rank of a control matrix C: rank=4 (if rank=4 it is controllable)\n"
     ]
    }
   ],
   "source": [
    "A = get_a_matrix()\n",
    "B = get_b_matrix()\n",
    "n = get_a_matrix().shape[0]\n",
    "\n",
    "control_matrix = control.ctrb(A, B)\n",
    "rank = np.linalg.matrix_rank(control_matrix)\n",
    "print(f\"Eigen values of matrix A: {list(np.linalg.eig(A)[0])}\")\n",
    "print(f\"Rank of a control matrix C: rank={rank} (if rank={n} it is controllable)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6dc21805",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Design LQR controller\n",
    "Q = np.eye(4)   # np.diag(np.array([1, 1, 10, 100]))\n",
    "R = .0001\n",
    "K, S, E = control.lqr(A, B, Q, R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "21d160a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pendcart(x, t, pend_mass=1, cart_mass=5, arm_length=2, g=-10, d=1, s=1):\n",
    "    Sx = np.sin(x[2])\n",
    "    Cx = np.cos(x[2])\n",
    "    D = pend_mass * (arm_length ** 2) * (pend_mass + cart_mass * (1 - Cx ** 2))\n",
    "    \n",
    "    wr = np.array([1, 0, np.pi, 0])\n",
    "    \n",
    "    u = np.dot(-K, (x - wr))[0]\n",
    "    return np.array([\n",
    "        x[1],\n",
    "        (1/D)*(-(pend_mass**2)*(arm_length**2)*g*Cx*Sx+pend_mass*(arm_length**2)*(pend_mass*arm_length*(x[3]**2)*Sx-d*x[1]))+pend_mass*(arm_length**2)*(1/D)*u,\n",
    "        x[3],\n",
    "        (1/D)*((pend_mass+cart_mass)*pend_mass*g*arm_length*Sx-pend_mass*arm_length*Cx*(pend_mass*arm_length*(x[3]**2)*Sx-d*x[1]))-pend_mass*arm_length*Cx*(1/D)*u\n",
    "    ])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f21072a",
   "metadata": {},
   "source": [
    "Initial conditions:   \n",
    "$x_0 = [-1, 0, \\pi +  0.1, 0]^T$  \n",
    "Actuation:   \n",
    "$u = -K(x - w_r)$  \n",
    "The reference position:  \n",
    "$w_r = [1, 0, \\pi, 0]^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "916ce90a",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(0, 10, 0.01)\n",
    "x0 = np.array([-1, 0, np.pi + .1, 0])\n",
    "\n",
    "sol = odeint(pendcart, x0, t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7a05d921",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2kAAAHiCAYAAAByYISUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABy+ElEQVR4nO3ddYBc1cH38d8ZWZdkJe5OCBIIBCcQtFAq0NJSo0bLQwtUnrpQ6qX+Pi2UQr2lpUgLVHALEiQJEkKIJxtbt1kbOe8f587u7GY1mdmZ3f1+2uHOXD0zczOzvzlyjbVWAAAAAIDM4Et3AQAAAAAAXQhpAAAAAJBBCGkAAAAAkEEIaQAAAACQQQhpAAAAAJBBCGkAAAAAkEEIaRhxjDHXGWP+NMzH3G6MOauPZbOMMdYYExjOMmF0MMZcaYzZb4xpNsaU9lg2w5vvT1f5DpYxZoUxpiIDyjHof5/8W85cA703xpgvGWNu6W1dY8x/jDEfSHJ5Dvl7yBjzHmPMA/0sz4h/Q2NZ4vvM5wOGGyENGckYc5kx5gXvD9S93pfsKeku12jWXxBFahhjgpJ+LOkca22BtbYmcbm1dqc3P5qGsv3OGPOt4T4uRg5jzOXGmFXDvW1vrLXfsdZ+pI9l51trf5+sYyWLtfbP1tpz4o+9ADAvnWWKM8Y8Zozp9fU0xmQbY75rjNlpjGk1xmwyxnzWGGN6bN/mfYdXG2PuMsZMHr5nQMjFyEdIQ8Yxxnxa0k8lfUfSREkzJP1S0lvSWCyk0Bj+ZXKipBxJ69NdkEQjseYOw2sM/5uF9HdJKyW9SVKhpPdJ+pikH/VY7xPW2gJJ8yQVSPrhcBZyMIbrPObfCw4GIQ0ZxRhTLOl6SVdZa++y1oastWFr7b3W2v/tY5uLjDHrjTH13q93hyUs+7wxZrcxpskYs9EYs9Kb7zPGfMEYs8UYU2OMud0YU5Kw3fuMMTu8ZV8e4nOYYoy5xxhTa4zZbIz5aMKybGPMT40xe7zbT40x2d6yFcaYCq/ZTrVXs/Wefo5zuTFmq/fctnlNZ7K84x6RsN4EY0yLMabcGFNmjLnPe61qjTFPeq/FH+XC8L3eL5+f87Y9wRjztLf+S8aYFQn7fcwY8y1vebMx5l5jTKkx5s/GmEZjzPPGmFl9lD3ebOTDxpidkh7x5n/IGLPBGFNnjLnfGDPTm2+MMT8xxlR6+37FGLPEW/Y7Y8xNxpgHvdfi8fh23vKTvLI0eNOTejyHbxpjnvK2fcAYU+YtyzHG/Mk7B+q9bSd6y4qNMbcaV8u723sdeg02fb3nxpgFkjZ6q9UbYx7p53WKN9vqr7z/McZ8osf2Lxlj3u7dX+S9RrXG/Vt4Z8J6vzPG3GiM+bcxJiTpw5LeI+lz8ffWW2+KMeZOY0yVd85dnbCPXG8/dcaY1yQd19vrkbC+NcZc7Z3D1caYG4wxvoTlvZ4LCdt+3Lhf8OuNMb8wxv2Kb4zxG2N+6O1zq6QLehy3W42x6afZWn/rJrw3HzTG7PLK+XFjzHHGmJe9cv1fP89/MJ8Fn/HO+b3GmA/2s68Peq9Vk/d6fqyfl17GmI8mrP+aMeYYb378MzE+/20J21zunXc/McbUSPqbpJskneidI/V9HKu3z6nDetvWGHOBMWatcf/Gdxljrutllx/yXq+9xpjPJhynv/exs1bI+zfRnHCzxvtcM/1/3s027rOlyRjzoKSyfl7fx40xF3v3T/aOcYH3eKUxZl3Ca7PKu/+Et3m8fJcm7O+Qz4OhnlP9HGOlpHMkXWytfdVaG7HWPivpvZKuMcbM6bmNtbZe0j8kHd3PfnONMT8y7nu3wRizyhiT6y3r7zt+u3G1eC972/3NuM/ufEn/kTQl4b2e4p0ndxj32d4o6XLTz3f2AK9Fn98Dvfx7uW4w+wS6sdZy45YxN0nnSYpICvSzznWS/uTdXyApJOlsSUFJn5O0WVKWpIWSdkma4q07S9Jc7/41kp6VNE1StqRfSbrNW7ZYUrOk07xlP/bKdFYf5ZklycbLLOkJuZq/HLkvpSpJZ3rLrveOO0FSuaSnJX3TW7bCO86PveOe7j23hb0cM19SY3yZpMmSDvfu/1LS9xPWvUbSvd7978r9cRT0bqdKMt6y7YnPUdJUSTVyv5b6vNe4RlK5t/wx77WeK6lY0muS3pB0lqSApD9I+u0Ar9kfvOeSK1dTulnSYd72X5H0tLf+uZJelDROkvHWmewt+52kpoT362eSVnnLSiTVyf3SG5D0bu9xacJz2CJ3HuV6j7/nLfuYpHsl5UnySzpWUpG37G65cybfey+fk/SxPp5rf+95/HXo9XzvuXyA8r5f0lMJ2y6WVO+9Jvly/xY+6L0OSyVVS1qc8Bo2SDrZe69zvHnfStifz3sPvib372uOpK2SzvWWf0/Sk95rPl3Sq5Iq+vl3bCU96q0/Q+7c+Yi3rM9zIWHb+7zzYYbcv7HzvGUfl/S6V4YS7xiJr+F2dT/Pr1PX50nP13sw697kvV7nSGqT+2N0gty/n0pJpx/EebFC7rPgerl/p2+S1CJpfB/7ukDu36GR+9xokXRMH+u+Q9JuuRBt5Go5ZiYsm+K915fKff7E/51d7pXpk957kuvNW9XPe9zf59QB23rP+wjv+EdK2i/prT1e79u8/R7hve9nDeJ9fEzeudXjeFfInStFGvjz7hl1fTafJveZ86d+3tv/593/kty/2e8nLPtZb6+BV+Z5PV6PpJwHB7Gvvl6z70l6vI9tdkj6aM/tJZVKekjSP/s5V37hbTNV7vP2JO+17vM7PuHf6HNy522JpA2SPp7wnCt6HOc6SWFJb/Xe51z1/519nfo+r/r8HlAv/176eu7cuPV1S3sBuHFLvMn9er9vgHUSPzS/Kun2hGU+uT9AVsj98VEpFxqCPfaxQdLKhMeTvQ/ugNwfoX9NWJYvqUODCGlyfxhGJRUmLP+upN9597dIelPCsnMlbffur/A+1PMTlt8u6au9HDNf7g/wi3t++EtaLmmnusLXC5Le6d2/XtI/lfCHQMJ229X9D9LPS/pjj3Xul/QB7/5jkr6csOxHkv6T8PjNktYN8JrNSZj3H0kf7vFetkiaKelMuT/iT5Dk67Gv3/V4vwq892C6XDh7rsf6z0i6POE5fCVh2f9I+q93/0Nyfzgf2WP7iZLaE193ufD3aB/Ptb/3vPPcGejcGkR5C+X+mJnpPf62pN949y+V9GSPff9K0tcTXsM/9PK6Joa05ZJ29ljni/KCuFxgOy9h2RUaOKSd1+O5PDzQuZCw7Sk9/p18wbv/iLw/0rzH5yi1IW1qwvIaSZcmPL5T0rUHcV6skNSaeF7IfZad0Nfr2WPf/5B0TR/L7u9rWS/rrpP0Fu/+5b28/5dr4JBWr94/p/rd1lvnp5J+0uP1XpSw/AeSbh3E+/iYegQOSad4r+kC73Gfn3dyPwT0/Gz+i/oOaSslvezd/6+kj0h61nv8uKS39/YaqPeQlpTzYKj76u018+bfooTP2x7LnpX0pYTtW+R+/LHeuTSjj+18XtmO6mVZn9/xtuvf6Ht7nBM3JTzn3kLaEwmPB/rO7vW80gDfA+rl3ws3bkO90dwRmaZGUpkZfPvtKXK/3kmSrLUxuRqDqdbazZKulfuQrTTG/NUYM8Vbdaaku73mE/VyoS0q98E7xdtHfJ8hr1ySpB5NZWb0Up5aa21Twrwdcr8OHlBe7/6UhMd13vH6Wp5Ypkvlag32GmP+ZYxZ5C1bLffluMKbN0/SPd6mN8j9CvmA1xzmCz33nWCmpHfEXyPvdTpFLtDG7U+439rL44J+9i8lvM7e8X6WcKxauV+Ep1prH5H0f3K/tlYaY242xhT1th9rbbO37RQd+HpL3d8PSdqXcL8locx/lPsj7a9e86ofGDfQx0y5X3T3JpT1V3K/pPZmoPd8qHotr3fO/UvSu7xl75b0Z+/+TEnLe7yX75E0KWFfie9Fb2bKNR1K3MeX5P7NSD3+3ejA1703PddP/PfZ67mQsH5f79vBlONQHOy/gYHOixprbSThceJz7MYYc74x5lmvuVa9XC1JX83xpssFxN72835jzLqE131Jj/0MdI5009/nVB/HX26MedS45rQN3nY9n0df58ygGWOmywX7D1hr3/Bm9/d5N0W9fzb35RlJC4xrHn20XIuB6cY1TT5eruZmsJJ5Hgx6X/2oVvfvgESTveVxV1tri+VqRcfLtVzpTZlcLVZv52Wf3/EJ6/T1WdCXxHNooO/svgzme2BI/16AnghpyDTPyP069dZBrr9H7sNSkuu7JPdHyG5Jstb+xVp7ireOlfR9b9Vdks631o5LuOVYa3dL2uvtI77PPLnmGvL2WZBw29lLeUqMMYUJ82bEy9OzvN6yPQmPx3tt6fta3slae7+19my5L8bXJf06YfHv5foIvE/SHdbaNm+bJmvtZ6y1cyRdJOnTXh8Dea9Pol1yvywnvkb51trv9Vaeg5R4zF1yTUUSj5drrX3aK/vPrbXHyjXjWyApsY9i4vtVINfsZY8OfL2l7u9H3wVzfSG/Ya1dLNf05kK5JoW75M7RsoRyFllrD+9jVwO958l0m6R3G2NOlPuj51Fv/i65JkqJr22BtfbKhG17vv+9nQ/beuyj0Fr7Jm95t383cs9zID3Xj78u/Z4LAxioHCG5Jqxxk9S3oaw7VEk5L4zrx3an3KAME6214yT9Wy7U9maXXJO4nvuZKfcZ8gm55sDj5JqsJu5noHPkAP18TvW27V/kflCa7v1xf1Mvz6Ovc2ZQvH5O/5D0U2vtfxIW9fd5t1e9fzb3ylrbItc0+BpJr1prO+Rq5T8taYu1trqvbQ/WQZwHB+shuR98Et8HGWOWy70mj/fcwFr7iqRvSersO9pDtVxT4QPOSw3wHT+Avs7PxPkDfWf3ZTDfAwP++wD6Q0hDRrHWNsg1N/yFMeatxpg8Y0zQ+4XwB71scrukC4zrjB2U9Bm5D86njTELjTFnel9ebXK/ase87W6S9G3TNTBFuTHmLd6yOyRdaIw5xRiTJddEcFD/Vqy1u+S+jL/rdV4+Um4QhniH9tskfcU7Xpn3XHt2dv+GcQOAnCoXDP7e8zjGmInGmLd4fzS0y/WhiyWs8idJb5MLan9I2O5CY8w874uuQa72ML7dfrl+Ron7eLMx5lzjBmPIMa7zeV+/hh6qmyR90RhzuFfWYmPMO7z7x3m/sgfl/nBuU/fn+6aE9+ubck2Ldsn9kbLAuEs6BIzrjL9Yrj9Tv4wxZxhjjjCuI3ijXHPYmLV2r6QHJP3IGFNk3MArc40xp/exq8G858nyb7k/aK6X9DfvV2fJPd8Fxg2IE/Rux5mEDvi96Hk+PCepybjBeHK9c2KJMSY+QMjtcu/feO8c+eQgyvu/3vrT5f6g/Zs3v89zYRBul3S1MWaaMWa8pJ61xeskvct7DZZJuqSffQ1l3aFK1nmRJdd3p0pSxBhzvlwTz77cIumzxphjjTPP+xzMl/ujskpyg1DI1aT1Z7+kad6/uwMM8DnV27aFcrUabcaY4yVd1stuv+p9Lxwu18fyb72s05/fSHrdWtvz+6TPzztr7Q65ZuPxz+ZT5Jpz9+dxucAbDy2P9Xjcm57/5oZiqOfBYAS81yF+C1prH5L0sKQ7jTGHe6/VCXKv3x+stRv72Nfv5WrdL+q5wPuc+o2kHxs3iIffGHOi993d53f8IMq/X1KpcQOS9WoQ39l9bTfU7wFgyAhpyDjW2h/J/eL4FbkvnF1yX27/6GXdjXJB5P/J/Rr3Zklv9n65zJbr5Fwt1xxiglwfGskNLnGPXLO/Jrm29Mu9fa6XdJXcr7p75QaaGMq1Vt4t13Z9j1zH4q97X2yS+zXxBUkvS3pF0hpvXtw+73h75Jqqfdxa+3ovx/DJvUZ75JqCnS6ps1bE++JZI/dH15MJ282X+yW0Wa7W8pfW2nhty3fl/misN8Z81tvHW+SatMXfh/9Vij43rLV3y9V0/tW4UbdelXS+t7hI7hf4OrmmKDVyTTfj/iLp63KvxbFy54Ssu+7YhXJf7DVync4vHOQv2ZPkAnujXHPYx+WaQEquRi1LbrCUOm+9vpoADfSeJ421tl3SXXL9MP+SML9J7g+2d8mdM/vkXuvsfnZ3q6TF3vnwD+uu1XahXPOtbXL/rm6RGzRGkr4h995sk/vj5Y8H7PFA/5SrcVgn11TzVq+8/Z0LA/m1XDPVl+Re67t6LP+q3C/2dV6Z/6K+DWXdoUrKeeG9t1fL/TFbJxds7uln/b/L9Vf8i9zgF/+QVGKtfU2uX+kzcn/cHiHpqQEO/4jc5SP2GWN6+zfV3+dUb9v+j6Trvc/kr3nPqafH5ZpsPyzph9baPi8G3Yd3SXqb6d5s/dRBfN5dJvcdUSv3WfOHA3d9QDkL1dW0sefj3lwn6ffev7l39rPeAYZ6HgzSjXI/bsZvv/XmXyxXS/9fuR/MnvHuX9FP+Trkvne/2scqn5X7d/C83Gv8fbn+x/19x/fL++68TdJW7zXtq2lsf9/Z/RnK9wAwZPGBBQCkmXHDPf/JWpuUmipjzG8k7bHWfiUZ+8tUxpjfyXUOH9XPc7QxxlhJ863rOwpghDLG/F6ub9cFgwlPAAaHmjRgFDLu+mRvl1czAQBAinxEroXGMekuCDCaENKAUcYY80255mE3WGu3pbs8AIDRyxtk6fvWXdQaQJLQ3BEAAAAAMgg1aQAAAACQQQhpAAAAAJBBAuk4aFlZmZ01a1Y6Dg0AAAAAaffiiy9WW2vLe1uWlpA2a9YsvfDCC+k4NAAAAACknTFmR1/LaO4IAAAAABmEkAYAAAAAGYSQBgAAAAAZhJAGAAAAABmEkAYAAAAAGYSQBgAAAAAZhJAGAAAAABmEkAYAAAAAGYSQBgAAAAAZhJAGAAAAABmEkAYAAAAAGYSQBgAAAAAZhJAGAAAAABmEkAYAAAAAGYSQBgAAAAAZhJAGAAAAABmEkObZF9qncDSc7mIAAAAAGOMIaZKstfr8E5/XhXdfqNV7V6e7OAAAAADGMEKa54ojr1BOIEefePgTen7f8+kuDgAAAIAxipAmyRijk6eerN+e91tNLZiqqx6+Susq16W7WAAAAADGIEJagpKcEt1y7i2akDdBVz50pdbXrE93kQAAAACMMcZaO+wHXbZsmX3hhReG/biDtS+0T5f/93I1h5v1m3N/owXjF6S7SEljrVVNW422NWzT9sbt2hfap5rWGlW1Vqm6tVot4Ra1RlrVFm1TW6RNHdEO+YxPPuNTwBeQz/iU7c9WQbBA+cF85QfzVRAsUGFWocbnjFdJTknnLfFxXjAv3U8dAAAAyBjGmBettct6XUZI692upl26/L+XKxKL6KazbtJhpYelu0hDZq3VjsYderXmVa2vXq/1Neu1qW6TmsPNnev4jE8lOSUqzy1XSW6JioJFygnkdN6CvqCstYraqGI2pqiNqi3SplA4pFA4pOZws0LhkBrbG1XXXqfWSGuvZcnx53QLbuNzxqs0p/SAMBd/nBPIGa6XCQAAABh2hLSDtL1huz764EfV1NGkn53xMy2fvDzdRRrQ3ua9enbvs3p277N6bt9zqm6tluRC0qKSRVpUskizimdpdtFszSqepYl5E+X3+ZN2/JZwi+ra61TbWqu69jrVtNaotq1WtW21qmurU217beey2tZadcQ6et1PXiBPRdlFnTV2ndMsN80L5CnLn6UsX5aC/qCy/FkK+oLK8mUpy58ln/HJyHTbpzHdH1trFbERRWNRRWxEkZi7H7VRhWPhzvuRWOSA5ZGY99gLr/EAG5/2DLbxdbqtF4vKquvfX7y88XIaGRljuubLyP2/l+Xe/Pg8n3yd8+KvRXy/nY/j2xvT/XF8HwmvYZ/7MD22847X89jx55V47MTnMpBBr2fSs7+hSNtzSdNzHuxxh3LsTH8uQzluxr9/SX5PhrJuKo6djuOm4nwY9P5GwOfIoJ/zoFfL8O+VobzHGf6cR5LyvHIdXnp4uotxAELaIdgX2qcrH7pSOxp36KsnfFVvm/+2dBepm5iN6eWql/XIzkf06K5Htb1xuyTXv2755OU6ftLxOqLsCM0dN1cBXyC9he3BWqtQOKS6tjrVtNW4ENfWFe7itXTNHc2dtXbxeaFwKC1lDvgCCvqC8hu//D6/mxq/jDHyG798xtc5jd/6nO/zd4YpWcl6/4u/NlZWif8+Ex93rustjtlY57z4Oj3nde5Tttv+4+t128YeuF63qdx28XIn7iNxn71tCwAAMJzOm3Webjj9hnQX4wCEtEPU0N6gzzz+Ga3eu1oXz79YX1z+RWX7s9NWno5oh57b95we3vmwHt35qGraahQwAR0/+XidMvUULZ+8XPPHzR+Vv4TEWWsViUXUEetQR7RD4VhYHdEOdcQ6FI6GFbXR7uv3EQ4CxgtdXuAK+AIK+AKd9xPn+Qzj7CRDbwG03/UHGewGHQAHvVqSjysl/zkn+fM72ccd0muTruecrvMwBcdO13GH8p5k/HNJ47+BdH02pfPfQNKfc4afD0OR6c95pCnKKtL0wunpLsYB+gtpmVW1kqGKs4v1q7N+pV+s+4V+/cqvta5ynb5x8jd0VPlRw1aG1kirVu1epQd3PKgnKp5QKBxSXiBPp0w9RStnrNQp005RUVbRsJUn3YwxCvqDCvqDyg/mp7s4GILEJpoAAAA4EDVpQ7Rq9yp945lvaH9ovy5deKk+dtTHVJZblpJjNXc06/GKx/XQjoe0avcqtUXbNC57nM6ccaZWzlip5ZOXp7VGDwAAAMDBobljkoXCIf1szc90+8bbleXP0rsWvkuXLrpUUwumHtJ+46MxPrXnKT21+yk9u/dZhWNhleWWaeWMlTp75tk6duKxGde3DAAAAMDQENJSZEfjDt300k3619Z/ycpq+aTlWjnT1XDNLpo9YJ+wSCyiLfVb9Er1K3ql+hWt3rtau5t3S5JmFM7Q6dNP19kzz9ZR5UfRHwoAAAAYRVIa0owxOZKekJQt18ftDmvt1/vbZrSEtLi9zXv1zy3/1D1b7tGupl2SXAfFWUWzNLVgqvKz8pXjz1E4FlZrpFVVLVWqaK7Q3ua9itiIJNfvbemEpTp5ysk6ecrJml6UeZ0bAQAAACRHqkOakZRvrW02xgQlrZJ0jbX22b62GW0hLc5aq4qmCj2771ltrN2o7Y3btad5j1rCLWqNtCrLn6XcQK7Kcss0rWCaphVO09xxc3Vk2ZGaVjhtVI/GCAAAAKBLSkd3tC7lNXsPg95tdI7fOQBjjKYXTacWDAAAAMBBS0pHJ2OM3xizTlKlpAettauTsV8AAAAAGGuSEtKstVFr7dGSpkk63hizpOc6xpgrjDEvGGNeqKqqSsZhAQAAAGDUSeqQgdbaekmPSjqvl2U3W2uXWWuXlZeXJ/OwAAAAADBqHHJIM8aUG2PGefdzJZ0t6fVD3S8AAAAAjEXJuCryZEm/N8b45ULf7dba+5KwXwAAAAAYc5IxuuPLkpYmoSwAAAAAMOYltU8aAAAAAODQENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDENIAAAAAIIMQ0gAAAAAggxDSAAAAACCDBNJdABwEa6V1f5aev0UKVUvTj5dWfEkqm5fukgEAAAA4RNSkjTTWSvddK/3zKsnGpBknSJselH51qrT5oXSXDgAAAMAhIqSNNM/8n/Ti76RTPiV99DHp4lukq56TSuZKf3uftPeldJcQAAAAwCEgpI0kNVukh6+XFl0orfy65PPevqLJ0nvvlHKKpbs/LkU60ltOAAAAAAeNkDaSPHSd5M+S3vRDyZjuywonShf+RKp8TXrqp+koHQAAAIAkIKSNFJUbpA33SCf8j6s5683C86XD3y49cYPUUDG85QMAAACQFIS0keL5W6RAjnTClf2vd/Y33OAiq346LMUCAAAAkFyHHNKMMdONMY8aY14zxqw3xlyTjIIhQbhVevnv0mEXSXkl/a87boZ09GXSmt9LDbuHp3wAAAAAkiYZNWkRSZ+x1i6WdIKkq4wxi5OwX8RtuE9qb5CWvndw65/6GTc8/zO/SG25AAAAACTdIYc0a+1ea+0a736TpA2Sph7qfpFg7R+lcTOlWacObv3xM6XFb5HW/knqCKW2bAAAAACSKql90owxsyQtlbQ6mfsd00I10vYnpSMv7RpyfzCO+4irfXv1ztSVDQAAAEDSJS2kGWMKJN0p6VprbWMvy68wxrxgjHmhqqoqWYcd/TY/5JouLjxvaNvNOFGasFh67tduIBEAAAAAI0JSQpoxJigX0P5srb2rt3WstTdba5dZa5eVl5cn47Bjwxv/lfInSJOXDm07Y6RlH5L2vexuAAAAAEaEZIzuaCTdKmmDtfbHh14kdIqGpc0PSwvOGVpTx7glF0u+oPTS35JfNgAAAAApkYyatJMlvU/SmcaYdd7tTUnYL3atdv3KFgyxqWNcXom04Fzplb9L0UhyywYAAAAgJQKHugNr7SpJJgllQU9bHpGMX5qz4uD3ceSl0uv3SVsfk+aflaySAQAAAEiRpI7uiCTb8Yw05Wgpu/Dg97HgXCmn2NWmAQAAAMh4hLRMFW6Tdr8gzTzp0PYTyJYWXiC98R8p0pGcsgEAAABIGUJaptr9ohTtkGYcYkiTpMPeLLU1SNufOPR9AQAAAEgpQlqm2vm0m8444dD3NfdMKZgvbbj30PcFAAAAIKUIaZlqx9PuYtR5JYe+r2COG8b/9X9Jseih7w8AAABAyhDSMlEsKu16XppxYvL2edibpVCVtPPZ5O0TAAAAQNIR0jJRzWapo0matix5+5x/juTPljbck7x9AgAAAEg6Qlom2rPWTacsTd4+swtd37QN90nWJm+/AAAAAJKKkJaJ9qyTgnlS2YLk7nfheVJjhVS5Ibn7BQAAAJA0hLRMtGetNOlIyedP7n7nne2mmx5I7n4BAAAAJA0hLdNEI9K+l5Pb1DGueKo0cYm06cHk7xsAAABAUhDSMk31G1K4JTUhTZLmny3tfMZd3BoAAABAxiGkZZq969w0ZSHtXMlGpS2Ppmb/AAAAAA4JIS3T7HvFDRpSOjc1+592nJRTTJNHAAAAIEMR0jLN/vVS+aLkDxoS5w9Ic1dKmx+UYrHUHAMAAADAQSOkZZrK16QJi1N7jPnnSM373QAlAAAAADIKIS2TNFdJoSppYopD2ryz3JSh+AEAAICMQ0jLJJWvuWmqa9IKyqUpxxDSAAAAgAxESMsklRvcNNUhTXJD8e9+UWqpTf2xAAAAAAwaIS2TVK6X8kqlggmpP9bcMyUbk7Y9kfpjAQAAABg0Qlom2e8NGmJM6o819Vgpq1DayvXSAAAAgExCSMsU1kpVr0sTDhue4/mD0uzTpC2PuGMDAAAAyAiEtEzRuFvqaHbXSBsuc8+Q6ndKtVuH75gAAAAA+kVIyxTVm9y0bP7wHXPumW665ZHhOyYAAACAfhHSMkXNZjctHcaQVjJHKp4hbX1s+I4JAAAAoF+EtExR/YaUVSAVThq+Yxrjmjxue0KKRobvuAAAAAD6REjLFNWbXFPH4RjZMdHcM6T2RnfNNAAAAABpR0jLFDWbh7epY9zs0yUZ+qUBAAAAGYKQlgk6WqSGXcM7aEhcXok09RiulwYAAABkCEJaJugcNGReeo4/5wyp4gWprSE9xwcAAADQiZCWCWriw+8vSM/x554p2ai07cn0HB8AAABAJ0JaJqjeLMlIpXPTc/xpx0nBfJo8AgAAABmAkJYJ6rZJRVOlYG56jh/IkmafyuAhAAAAQAYgpGWCuu3S+FnpLcOcM6Tara4sAAAAANKGkJYJMiGkzT3TTbfQ5BEAAABIJ0JauoVbpaa96Q9pZfNdk0v6pQEAAABpRUhLt7odbprukGaMNPcMaevjUiya3rIAAAAAYxghLd3ifcBKZqe1GJJcv7S2emnPunSXBAAAABizCGnpFg9p6a5Jk1xIk2GURwAAACCNCGnpVrddyiqQ8krTXRIpv1SafCT90gAAAIA0IqSlW902V4tmTLpL4sw9U9q1WmprTHdJAAAAgDGJkJZumTD8fqK5K6VYRNr+ZLpLAgAAAIxJhLR0sjbzQtr05a755eaH0l0SAAAAYEwipKVT834p0pZZIS2QJc0+3YU0a9NdGgAAAGDMIaSlU+fIjhkw/H6ieSul+p1SzeZ0lwQAAAAYcwhp6VS7zU0zqSZNciFNoskjAAAAkAaEtHSq2y7JSOOmp7sk3Y2fJZXOJ6QBAAAAaUBIS6e67VLRVCmQne6SHGjeSmn7U1K4Nd0lAQAAAMYUQlo6ZdrIjonmnSVFWqUdT6e7JAAAAMCYQkhLp4Zd0viZ6S5F72aeLPmzpc0Pp7skAAAAwJhCSEuXaFhq2isVT0t3SXqXlSfNOpl+aQAAAMAwI6SlS9NeycZcn7RMNe8sqXqjG44fAAAAwLAgpKVLQ4WbZmpNmuRCmkSTRwAAAGAYEdLSpTOkZdjw+4nKFrjy0eQRAAAAGDaEtHRp2OWmxRnc3NEYV5u25VEp3Jbu0gAAAABjAiEtXRoqpNwSKSs/3SXp38I3SeGQtH1VuksCAAAAjAmEtHRp2J3Z/dHiZp8mBfOkjf9Od0kAAACAMYGQli4NFZndHy0umCPNPVPa+B/J2nSXBgAAABj1CGnp0lAxMmrSJNfksWmPtHdduksCAAAAjHpJCWnGmN8YYyqNMa8mY3+jXluD1N4wckLagnMl43O1aQAAAABSKlk1ab+TdF6S9jX6Nex205ES0vLLpOnL6ZcGAAAADIOkhDRr7ROSapOxrzFhJFzIuqeF50v7XpHqd6W7JAAAAMCoFhiuAxljrpB0hSTNmDFjuA6bmTqvkTaSQtqbpAe/5po8Lr8i3aUBAABjjLVWMStFY1YxG7+5x9Zab746l7n5va8fs1axmFs3auPbe9t6+4nGt4n1OI63n1isa5/WumNZdZXBxufLrWulbuvG4ssS17Xq2l69r6vO56jOZUrYT+K2sVhXmbqVJ3HdeLn7WfeA55iwrnqWV9727k3r8RzVuZ4617Pd5rvt1Mv+uvahXpbFOpcduL/zDp+k7118ZIrP0OQatpBmrb1Z0s2StGzZsrE9TGDjbskXkAomprskg1c2XyqdJ238FyENAIAUicWswrGYIlGrSLTrfjgaUyRmFYnGFI5aRWLxx+5+NGY7b5GYCxYRL0REot4y22OZN43GYorG5KbeOtGoW7+v/UatW6fzGJ3Lupcl6h0/cZ1oQkCKeeGoM/wkBK3OkOWtNxYHmTZG8hkjIzeVkXy9zjPd1jXGyGd0wDzTc123ebd56vHY5zZ2+0tYZuJl8BkFOvflplLC8oR13bKu4yZuE7+vzmUmYZ2ux0rcJvE16mffS6YWD8v7lUzDFtKQoKFCKpoi+fzpLsnQLLpAeuYXUkutlFeS7tIAANCveKDpiMTUEXW3cPx+pPvjcOc8t378cTgaU3vC40jMC0xeOApHveDkBSgXrLrux8NUfL2Bto+lOYgEfEb+HrfOecbI7/emnct9CviM+0PdWyfg8yknaOQzpvsyn5sXn/qM5PcZGWPk98UDQXwd77HPW8/E13OPffF9eX+cd+4zYf1uj+PHMT229ylhvpHPO27ncUxXmbsdx3RfHg8ZPR8bmc6Q4ouHBl/3kNVz3cQQhLGLkJYOI+UaaT0d/jbpqZ9Jr98nHfP+dJcGAJChYjGrjmhM7eGY2iJRtYdjao9E1R6JqS3spu2RqNri88OxQS0bKEx1JASwcAoCj99nFPQbBX0+BfxGAb9PQZ+bBnrMD3jBJC8r4Ob5fAp226bv7YPx7f3eNp3zvfu9bNMZmrxAETggTLltfT4p4PN1Bq6Ar3uYApAZkhLSjDG3SVohqcwYUyHp69baW5Ox71GpYZc048R0l2LoJh8tjZ8lrb+bkAYAI1A4GlNrOKrWDu8Wjvb92Ju2haNq6WXdxPDV5gWpeBDriMQOqZwBn1F2wKfsoF853jTL71N20Kcsv09Bv095WW6aHXBBJivgHmcF3DqJj7vmJ6zn9ykY8Cnbm8b323377usTYgAMl6SENGvtu5OxnzEhFpUa94ysQUPijPFq034uhWqk/NJ0lwgARqVINKZQR1Sh9ohaOiJqbo+qpT3SOS/UEXHT9vjjrnVbOg4MWvFp5CCqlnKDfuVm+ZUb9Csn6FNeVkA5QZ8KsgMqzfcrO+hTTsBNswM+ZQfcetkBvxe0Epf3v8xt71PAn6wrBAHAyERzx+HWXCnFIlLR1HSX5OAc/jZp1U+k1++Vjr083aUBgIwRi1mFOiJqaovfwmpqi6jRm8bnNXcLVxEvXEW9+25++xBqovKy/MrPDig/y6+8rIDys/0qyA5oQmF2Z8DKCfqV1xm03Hq5Wb7uj4N+5Wb5uj3ODlB7BADpQEgbbo173LRoSnrLcbAmHSmVzHFNHglpAEaR9khUDa1hNbSE3bQ13BmsGnsEr/i0uT3SGcSa2yMDjj7n9xkVZAdUkB3oClfZfpUVZHfez88KKN9bXpAdUF52QAXZLjgV9JifF/QTogBgFCKkDbemvW5aODm95ThY8SaPq34qNVdJBeXpLhEAdIrFrJraImpoDau+tcNNE0KXe9zRy7ywWsPRfvcd9BsV5gRVmBNwt+ygZpTkdZ+XE0h47KZFCfdzg35GbAMADIiQNtziIW2k1qRJ0pJLpCd/JL16h3TClekuDYBRylqrUEdUtc0dqm3pUG2oXbWhsOpCHaoJdXRNW9z92pYONbaG+x3RLyfo07jcLBXnBlWc50JWcW5Q4/KC3jy3bFxuUEW5QRXlBFSQE1BRTlDZAR8BCwAwLAhpw61xj7uQdV5Zukty8CYudiM9rvsLIQ3AkLR2RFXV1K6q5nZVNbWrutndXMgKdwax2lC76kJhdUR775sV9BuNz8tSSb67LZ5SpPF5WRqf58LVuHjY8sJXPHTlBEfY9SkBAGMSIW24Ne2VCiZJvhE+ctXRl0n/+Zy0f7008fB0lwZAGrWFo17Y6ugMXonTxPuhjt6bFBblBDoD19RxOTpiapHG52epND+rWxiL3wqyA9RqAQBGLULacGvcIxWN0P5oiZZcIt3/ZVebdu63010aACkQicZU3dyhfY1t2tfQpv2NbdrX2Kb9Dd60sU1VTe1qbIv0un1xblDlhdkqK8jSEdPGqawgy3ucrfLCbJV705L8LAUZch0AgE6EtOHWtE8qX5juUhy6/FJpwbnSy7dLZ31D8nMqASNJc3tE+xpata+hvTNw7UsIX/EA1rN/V8BnNLEoRxOLsrVwUqFOnV9+QPgqK8hWaUGWsgM0LQQA4GDwl/Vwa9orzT0j3aVIjqPeLb1+n7T5IWnheekuDQBPLGZVHWrX7rpW7alv0+76Fu2pb1NFXat217dqT32rGlrDB2xXnBvUpKIcTSzO0aJJhZ33JxXlaGJRjiYV56gkL4sh3wEASDFC2nBqb5baG0fu8Ps9zT9Hyp8gvfhbQhowjMLRmPbUt2p3XasqvNC1OyGA7WloU0ePiyEX5gQ0dVyupo7L1bKZ4zVlXK6mjPPClxfCcrOo+QIAIBMQ0obTSL9GWk+BLHdB6ydukGq3SSWz010iYNRoaAlrZ22Ldta2aEdtSLu8+ztrXa1YNKEdojHShMJsTRmXqyVTi3XukkmdgWzKuFxNHZ+ropxgGp8NAAAYCkLacGrc46ajYeCQuGUfdNdMe+FW6Zxvpbs0wIgRi1ntbWzT9uqQdtS48BUPYjtqQgcMxlGan6XpJXk6ZsZ4vfXoPE0fn6dpJbmaNi5Pk4pzlBVg4A0AAEYLQtpw6qxJG8EXsu6paIp02JulNX+UVnxJyspLd4mAjFIX6tDW6pC2VYe0rbpZ26pD2lrlHrcnNEkM+o2mjc/TjJI8HT19nGaU5Gl6SZ5mlrppQTYf1wAAjBV86w+neEgbTTVpknT8FdJr/5BevUM65v3pLg0w7Noj0W7hy01dIKtr6RqgI+AzmlGSp9ll+Tp1fplmlxVoVlmeZpbma1JRjvwMyAEAAERIG16Ne6XsYikrP90lSa6ZJ0kTl0jP/FI6+r0j/0LdQB/awi6MvbG/SZsrm/XG/iZtqmzWjpqWbn3EJhXlaHZZvs4/YrLmlOVrTnm+ZpcVaNr4XK4HBgAABkRIG05Ne6TCSekuRfIZI53yKenOD0uv3ystfku6SwQcko5ITJsrm7Wpskmb9ndNt9eEOq8b5vcZzSzN04IJhbrgiMmaN6FA8yYUaFZpvvJpmggAAA4Bf0kMp8a9o6+pY9zhb5Me+64b6fGwi1xwA0aAmuZ2bdjbpA17G7Vhb6Ne29uozZXNinhpzO8zmlWap4WTCnXhkZM1f2Kh5k8s0OyyfC7WDAAAUoKQNpya9kplp6e7FKnh80unfkb6x5XSG/+VFp6f7hIB3URjVtuqm/Xa3ia9tqexM5RVNrV3rjOxKFuHTS7SGYsm6LDJRVo4sVCzyvIIYwAAYFgR0oZLLCY17Ru9NWmSdMQ7pMe+Jz3+fWnBedSmIW1iMaut1SG9srteL1c06JWKBq3f06jWcFSSG0lx3oRCnTK/TIsnF+kw71aSn5XmkgMAABDShk+oSrLR0XMh6974g9Jp/yvd8wlp/d3Skrenu0QYA6y12lHTopd3N+iVChfK1u9pVHO7u85YTtCnJVOK9a7jp2vJlGIdNrlI8yYUcF0xAACQsQhpw6XJu5D1aA5pknT0ZdLqm6SHvi4tfJMUzEl3iTDK1Ld0aO3Oeq3ZWae1O+v1ckV954WfswI+LZ5cpLcfM1VHTC3WkdPGaW55vgKMqAgAAEYQQtpwaRyl10jryeeXzv229Ie3SM/8n3TaZ9Ndor7FYlLdNmnPWqlqo9SwS2rcI3U0S+FWyfhdyMwdLxVPc7fSedKUpVLxdJpzDoNYzGpTZbPW7KzTizvqtGZnnbZWhSRJPiMtmlSkC46coiOnFeuIqcVaOKmQIe4BAMCIR0gbLvELWY/2mjRJmrPCjfD4+A/cqI+lc9Ndoi6hGmnLw9KmB920pcbNNz6pcIoL0TnF7n2yMRfWmvZJFS9IrbVd+8kvl6YdL81bKc07Sxo/Mz3PZ5RpaA1r7c46rdlZr7U767RuZ72avGaL4/OCOmbGeF18zDQtnTFOR00bx1D3AABgVOIvnOHSXCnJuD/ux4I33SBte1y652rpA/em9wLXkQ5p0/3SutvcNBaR8spcuJp1ijT5aKl8kRQYYNCI9mZX47Znjat92/6ktPFfbln5ImnJJdIRF0slc1L+lEaLyqY2Pb+tTs9tq9HqbbXauL9J1rpasoWTinTR0VN0zIzxOmbmeM0qzZOh9hIAAIwBhLTh0rxfyit1g2uMBYWTpHO/I/3zKmnVj9PT7DFULT1/q/T8r93ALQUTpROudLV7k5cOPThmF0jTjnU3SbJWqt4kbX5Q2nCf9Oi33G3qMunod0tHvFPKKUr+8xrBKupa9Ny22s7b1mrXdDEvy69jZ47Xm46YrGNnjtdR08epgFoyAAAwRvFX0HBp3u+Cy1hy9HukLY9Ij35bmnGiNOvk4Tlu9WbXH+6l26RImzT/HOm4j0pzz5T8STzljZHKF7jbiVdJ9bukV++UXvm79K/PSA98TTryHdKyD0uTj0zecUcIa90w+ImhbHd9qySpKCeg42eX6F3HT9fxs0t1+JQi+pIBAAB4jLV22A+6bNky+8ILLwz7cdPq12e6vk7vuzvdJRlebY3SzSuk1jrpIw+ltn9a5QbpiR9K6++SfEHpqEulE66SJixK3TF7Y620e430wq0utEXapGnHScdfIS1+68DNKkewnTUtempLtZ7eUqNnttSoutldKLqsIFvLZ5foeO+2cGKhfD6aLgIAgLHLGPOitXZZr8sIacPkJ0tc/6e33ZTukgy/mi3SrWdL2YXS5f9yoyQm096XpSdukDbcIwXzpeM/Ip34CalgQnKPczBaaqWX/io9f4tUu0XKnyAt+6C07EOjoma1srFNz2yt0VObXTCrqHM1ZRMKs3XS3FKdMKdUx88u0eyyfPqTAQAAJCCkpZu10rcmuP5QZ1+f7tKkR8UL0h/f1lWbWDb/0Pe5e40LZxv/LWUXScs/Jp3wP1JeyaHvO9liMdf087lfSZsekHwBV6u2/OPStGUjZjj/hpawnt1Wo6e9ULapslmSa7544txSnTyvTCfNLdPcckIZAABAf/oLafRJGw6tdVK0QyoY+TUnB23aMjfK45/eLt18hnTRz6Ulbx/6fmJRF3KevdGNHpkzTlrxJRfQcsclu9TJ4/NJ889yt5ot0nO/ltb9WXr1DnfdteUfdwOaBLLTXdJuOiIxvbCjVk+8Ua2nt1Tr1d0NilkpN+jXcbNLdPGx03Ty3DItnlIkP80XAQAAkoKatOFQ+br0y+XSxbdKR1yS7tKkV/1O6e8flHa/IM0/Vzrr69LEwwfernar9Opd0po/SPU7pKKp0vEfdYNyjNQRFNubXFPI526Wqt9wl2c41msKmcaLnm+vDumJTVV6fGOVntlao5aOqAI+o6UzxumkuWU6aW6pls4Yr6wAA30AAAAcLJo7ptvWx6U/XCR94D5p9qnpLk36RTqk1TdKj98gdTRJ05e7ERinHO0uKO3zS20NUu02qeJ5aecz0v5X3bYzTpKWXyEtunD0XM7AWmnro9LqX0lv3O+e/+K3SMd/TJp+fMqbQobaI3pmS40ef6NKT2yq0o6aFknSjJI8nb6gXKctKNeJc0sZEh8AACCJaO6Ybs373bRgYnrLkSkCWdLJ10hL3yet/ZP08u3SI9/sfd1gvmsqefY3XXPAcdOHt6zDwRh3eYC5Z7oaw+duca/Lq3e6C20v/7hrGpqkppDWWm2patZDGyr1+MYqvbCjVuGoVV6WXyfOKdWHT5mt0+aXa1ZZflKOBwAAgKGhJm04PP3/pAe+In1hpxs4AwcK1bgmf017JRtzA4EUT5PKFiT32mYjRXuz9PJfpdU3S9Ubpbwy6djLpeM+LBVNGfLuwtGYnt9eq4deq9TDr+/vrC1bNKlQpy8s1+kLynXszPHKDviT/EQAAADQG2rS0q15vxTIccEDvcsvlfJPTHcpMkd2gXTcR1yfu22Pu6aQT/5IeuqnrmnokoulhedLWX3XdjW0hPXYG5V6aEOlHttYqaa2iLICPp00t1QfPXWOVh42QZOLc4fvOQEAAGBQCGnDobnSXbOLIckxVMZIc1a4W9126flbpVfucJcdCOa5oHb4211Tyaw8VTa16f5X9+nfr+zTc9trFY1ZlRVk6fwlk7TysIk6ZV6Z8ulbBgAAkNH4a204NO076P5okaoqNdx7n5offVQdu3bJtrcrUFaqnKOOUuGZZ6rg9NNl/DRRGxPGz5LO+aZ01jeknU+7Pmvr/yG9eqcivmy9FDxKdzUv0SPRpcorn6GPnz5HKw+bqKOnjZOP4fEBAABGDELacGiulErnDmkTGw6r5pZbVP3rW2RbWpR92GHKP+EEmZxsRfbuU9P9D6jhjjsVnD5d5Vd/UkUXXsjFg8cKn0/7S5bpX+Om6oGitynY8LTO8K3VeVqnbwefk4KSchZJ4ZOlxlOkllNcTS4AAABGBELacGjeL808adCrR2prtfvaT6nluedUePbZKv/UtcqeM6fbOjYSUdPDj6jm5pu1538/p/rb/64pN/xAwUlj+ILZo1yoPaL71+/T3Wt366nN1YpZN/DHm1ZerNOOuEpTyguk6k3SG/91/dhe/pv0wq1u49L50tRj3WUOJh8tTTrC9XsDAABAxmF0x1SLdEjfKpdWfEla8fmBV6+t1Y73vk/higpN/tY3VXzRRf2ub2Mx1d95pyq/+z2ZrCxNueEHKjiVa7GNFtGY1VObq3X32t3676v71BqOanpJrt529FRddPRUzZvQT9CKRqS9L0k7Vkk7npb2rO26HISMVL7QXUi8bIFUOq9rmpU3LM8NAABgLGN0x3QKVbpp4cB90mKhkHZ99AqFd+/W9Ft+rfzjjx9wG+Pzafw73qG8Y5dp97XXatdHr9DEL35BJR/4wKGWHGm0oyakvz6/S3e+WKHKpnYV5QT01qVT9fZjpmrZzPGDa9rqD0jTjnW3k69x8xr3SnvXSXvWuenuF6VX75KU8GNN8QypZLa7Jl3xDGncDO/+dKlo6ti8JAIAAMAw4q+tVBvkhayttdp73TfUtmGDpv3i/wYV0BJlz5mtWX/7q/Z87vPa/93vKVJdo/JPf4p+aiNIRySmB1/br9ue26lVm6vl9xmdsXCCLjl2qs5YNCE51zArmuxuC8/vmhdudRfRrn7DNZes2uhGknzjga4fGeKMTyqcIhVPdYGteKpUNM1duy1+P79c8vkOvawAAABjFCEt1Zq9P3IHGLih4a671HjvvSq7+pMqPOOMgzqULzdXU3/6E+27/puq+fWvFW1u0qSvfY2gluG2V7taszte3KXq5g5NHZerz5y9QO9YNl2TinNSX4Bgrmv2OPHwA5eF26SGCqlhp1S/S2rYJdXvlBr3uOaTr/9LirZ338YX9MLgtK4w1xnovFt+GZekAAAA6AMhLdWa9rlpPzVp4f37tf8731Xe8uUq+9jHDulwxu/XpOu+Ln9hgWpuuVUmGNTEL36RoJZhrLV6clO1fvPUNj22sUp+n9HKRRP07uUzdNr8cvkzZcj8YI5UNs/demOt1FLjglzjbqlht5vG7+96zgW6WLj7dv5sr/bNq4U7oFZumpQ7niAHAADGJEJaqsVr0vL7rknb/53vykYimvzN65NyzTNjjMo/8xnF2jtU94c/ypedrfJPf5qglgFaO6K6e+1u/fapbdpU2ayygmxde9Z8veu4GcNTa5ZsxrhasfwyN3Jkb2IxKVTVPbw1Vrjw1rDbDWrSuEey0e7bBXK94Dal91q54mlSTnHKnyIAAMBwI6SlWvN+KbdECmT1vvipp9R0//0qv/YaZc2YkbTDGmM08UtflO3oUM2vb5EvL09lV16ZtP1jaCob2/T7Z7brL6t3qq4lrMWTi/SjdxylC4+anJy+ZpnM53MD5xROlKYe0/s6saj7QaNxt1crt6f7/W2PS017JRvrvl1OsTewyUxvmnArni7ljkv50wMAAEg2QlqqNe/vs6mjjcVU9eOfKDh1qko+9KGkH9oYo0lf/5pirS2q+tnPFSgv17hLLkn6cdC3nTUtuumJLbrjxQqFozGds3iiPnTybB0/u4SazUQ+f9egJtN6HYnWXVKgeX9XeGvY5frJ1e+UarZIWx6VwqHu22QXHxjeEm+EOAAAkIEIaanWvL/PQUOaHnhAbevXa/L3vitfVu81bYfK+Hya8u1vK1pTq71fv07+sjIVrliRkmOhy+v7GnXjY1t070t7FPD5dPGx0/Sx0+ZoVll+uos2cvkDXjPHqdL0XkY/tVZqqZXqd7jgVr+za6CTum3S1scIcQAAYEQgpKVa835pxokHzLaxmKr+3/8pa95cFb/5zSktggkGNfVnP9POD3xAu6/9lGb+/nfKPeqolB5zrFqzs06/fHSzHtpQqbwsvz58ymx95NQ5mlg0AvubjTTGSPml7tZbs0prpda67iEufus3xE3vaj5ZNFkqnCwVTnKXIiicJGUXMsAJAABIKkJaKlnr+tn0UpPW/Nhj6tiyRVN++MOkDBYyEH9Bvqb/6iZtf9e7tevjV2rmX/6s7NmzU37cseKlXfX60YNv6Ik3qjQuL6hPnbVAHzhppsblpaaGFAfBGCmvxN2mLD1web8hbru07Umpo+nA7YL5PcJbQoArnOQGDcovc/3nCHMAAGAQCGmp1N4oRdp67ZNWc8utCk6ZoqLzzh224gTKyjTjll9r+7sv066PXqFZt/1FgfLyYTv+aPT6vkb9+IE39MBr+zU+L6gvnr9I7z1hpvKz+ac14gwU4iSpvdldVqNpb8JtnxvcpGmfu+RA074Drx0nuevH5Zd7o2GW93I/8XGZu34dAAAYk/hLMpWa9rtpj5DWsnatWtes0cQvfUkmMLxvQdasWZr+q5u04wOXa9fHPq4Zf/iD/AX0kxqqrVXN+ulDm3Tvy3tUkBXQp89eoA+ePEuFOcF0Fw2plF0gZfdz3Tipq0YuHuJCNe4SBKEqKVTddb96kxSqdD/k9CaQ40aGzR3vgmPiNLekl/slrg+dn3MQAICRjpCWSs29h7S6226TLz9f4y5+exoKJeUeeaSm/fQn2vU/V6nik5/Q9F/9KmUDl4w2FXUt+vnDm3Tnmt3K8vt05elzdcVpc2jWiC6JNXITD+9/XWuljlCPAFfpBkBprZVa6lzga62VqjZ23Y9F+t5ndpEX3rxQl1Ps3cZ13c8d12OeN+3jUiEAAGB4EdJSqZeQFqmrU9N/79e4Sy6RLz99NVgFp5+uyd/6lvZ+8Yva89n/1dSf/HhY+saNVPsb2/SLRzfrtud2yhijD5w4S1eumKvywux0Fw0jmTFe7VyBVDLIPqLWSu1NXYGtpda7X9cV7jrv17nLFLQ1SG31UrSj/30HchMCXI9g1y3cJSzLLnTBMLtQCmTT7w4AgCQgpKVSc6WbJgwc0nDX3bIdHRr3rkvTVKgu4972VkUb6lX5ve9r79e/rsnf/CbX7uqhNtShmx7fot8/vV3RmNU7j5uuT545T5OL6S+ENDFGyilyt/EzB7+dta5pZVuD1FrvBbeGrgDXVn/gsub9UvXGrsc9Lybeky/ohbaE4HbArcgLpj3nJTwO5ruLoAMAMEYR0lIpVCX5Au4XZ0nWWtX//e/KPeYY5SxYkN6yeUovv1zRhgbV3HiTAuPGacJnP5vuImWEhtawbn1yq25dtU2t4ajeunSqrlk5XzNL6b+HEcoYNxhJMNeNOjlU8Rq8nuGuvcm7NXrT5u7zmvdLNZu75kVaB1NYF9ayCqSsfCkrzwW3+P2sAimY5z32bkFvfpY3v7f1g7nU9AEARgRCWiqFqqS8ss5fhNteeUUd27dr8kc+nOaCdVd+9dWKNTSo5pZb5cvPV9mVV6a7SGkTao/od09v181PbFVDa1gXHDFZnzp7vuZNKEx30YD0SqzB0/SD3080nBDimnoJeT1u4ZDU0eL67rVUS/Xxx81SuGXgJpzdn0T3UBfMk4I5bhrI6XHfC7SBXDc/kJswL6fH/byEdbzHDOACADgEhLRUClVJBV1D3Dfce59MMKjCc85JY6EOZIzRxK98RbFQi6p+9nPZcERln/zEmGr62BaO6s+rd+rGxzarurlDKxdN0KfPWaDDpxSnu2jA6OIPdg2skgzRsAtwHSEX2jqau0JdYsALe+vEA15HyDX/DLdI4TYXAMOt7hZpc/PCLVIsfHDlMv4Dg54/2/Xbi9/82W6wlkCO5PemSVkWv59FzSEAjFBJCWnGmPMk/UySX9It1trvJWO/I16oyl3zSJKNRNT473+r4Iwz5C8qSnPBDmR8Pk3+zrclv1/Vv/ylbCSi8k9dO+qDWkckpttf2KX/e2Sz9jW26eR5pbr5nIU6Zsb4dBcNwGD4g25Ak9xxqdl/LJoQ3lq7wlukLSHQeUEv4q3X7X5r17rRDnc/0uE1/axy9yNt3rJ2d4u29z+C51DEQ5s/6EKbP+j6Dsbvd07j97NcM/34/W7LvKmvx/r+Huv3t39foPvNH/Tu+71p/DF9EgGMbYcc0owxfkm/kHS2pApJzxtj7rHWvnao+x7xQlVSqbueUuiZZxStqVHRmy9Mc6H6Zvx+Tf7WN2UCAdXcfLNioZAmfumLo3LUx3A0prvWVOjnD2/W7vpWHTtzvH586VE6aW6ZrLWKhUKKdXRIkYhsNCobiUo2JhMMultWVteUPyaA0cvn7xqBczjFol2BLdLePcB1BrvE+fEA2Nt6Hd4t7N0SHsfCXfc7Qj3W6bE8vt1AA8gkhUkIcPEQl/DYHxjksmBCAIyHwp6B0N+1zPhdQDR+N7/btJ/5vsDg103WfACjWjJq0o6XtNlau1WSjDF/lfQWSYS0UHVnTVrDvffKV1SkgtNPT3Oh+md8Pk36xnXy5eer9re/VXjPHk394Q1pvVxAMnW0tOrB/6zWA4+skdm/X5f7Qzout0PFW+sV+0ej3mhsVLSpSYoM/ldsk5srf2GhfIWF8hcUyFdUJH9hgXwFhfIXFcpXVCx/UVHX/eIi+YuKvPUKh/2C5gBGAJ/fDXqivHSX5ECx6MBBLho5MAxGvBrCWNSbhrseR+P348sTHkcjCct63KLhhP0lzAu39rIs4XE03HtZRpI+Q11fN9PPssEsH+w+/Ckug9xUxmvOa7q27XeeDpx3wH76mjfAvge1TW/r9VVuHUQZveDe23Gl7vMS1xvSffWyv0O538dxRnkrrsFKxl+HUyXtSnhcIWl5EvY7rCqvvUQyRhN+8vfk7DDePyK/TLGWFjU99LCKL3jTiLhotDFGEz//OQWnT9P+b31b29/3Pk2/8UYFJ04ceOMMEm1uVuu6l9S2fr3aNm5U9UuvKrCnQnNsTB/31vHl5Sk4dYr8EyYqa9p0+YqL5C90gcpk58gEA5LfL+MPSMbIhjtkw2HZjrCbtrcrFgop2tykWGOTYs1NitbXK7xrl6JNTYo1NsqG++/T4svPd8ctKnZhL36/qEj+4iL5Covky8uTyc6SLydHJiu76352tkxWljuvjPG+oN0HnOnxWJIUi7lawVhUNhp1NYWxmGwkIsVrDL1lifNsNGF5NKF2MRqRjUQS7nv77byfsI9IpPv8+LrWHvii9DbPLRjku5/wAW8Gcz8+K/GL4eD2kThv6PtT//p8WfpY0Ofr2PfO7MHsK5nlGuJ733d5+z5E/69LGlGuAfi9W+L1KZNYNuv9x9ohTOVqFnsss7aP9Ye0796mA5QxcZ1u99X9+PHHnffjtaORhPfb9rjf2357lmmg+4nb9TK/r2OkS6ac+mNaX1+Mpsei3tbrWid3ySKV/vi+5BYtxYbtJ3xjzBWSrpCkGTNmDNdhB611/Zbk7jBU5ab55WpetUq2pUVFF1yQ3GOkWMlllyk4ZYp2f/oz2va2t2vKD36gglNOTnex+hSpq1PL6tVqeeFFtax5Ue2vb5Ri7ounprBUb+RPUuPS83TMimU67qQjlTVtqnzFxSnvdxdra1O0oVGxxgZFGxvd/SY3jTY2KtrYoFj8flOjwjt2qq3RPbatgxmuPIP5/a65bCAgE78fDMj43WP5/X2//kOdH5f4B2XCfZv4bWt7Wb+P7dKyj/6eY5+vS1+zD2ZfB/HaZ2K5+i1v34vSqd/XJZ0y+ZftTC1b0stlekwPdjfD/Hr18VtVv+v1ZPt84M3qLU31tl5v8/tfz31H2T6X93tc28u8ngsPWNRL+QZcZ5DL+y3PYPetAZ6/3Hs55DL3sl6vmwzyufdYLegfWRUNUnJC2m51H495mjevG2vtzZJulqRly5Zl3m8TPslGklisULWb5per+a8Py1dcrLxly5K3/2FSuGKFZt/+N+3+1Ke06yMfUelHPqyyT3xCvpycdBdNNhpV2/r1an7iSYWefFKtr7wixWIyubnKOeoo7X3zu/W3thI9GZigqVPLdfXK+br8iMny+Yb3C8qXk+Ner4kTBl65B9vRoWhjo2JtbbLt7a7mrq1dtqPH/Y4OyXq/3sbiv9jGejy2UsAv4/PLBPySPyDj93XWFJqAazLjlvm7QlUgIBOIhyq3XmLwkt/vlscfBwYRwAAAANCnZIS05yXNN8bMlgtn75J0WRL2O6yMMX03mzkYXk2azRmvpsceV+GK00ds36PsefM06/bbtf+731PNLbeq8cEHNfm665R/4onDXpZITY1CTz3lgtlTTylaVycZo9wjj1TZ//yP/MtP1N3NBbr5mZ3a39iuo+YX66dnzNNZh00c9nCWDCYrS4GysnQXAwAAAMPokFODtTZijPmEpPvlGov/xlq7/pBLNtx8RjaWxJDWXClJatm0X7GGBhWsXJm8faeBLzdXk6//horOP097r7tOOz/4IeWffpomXHONchYvTtlxbTis1pdeUvOqVQqtekpt69dL1spfWqqC005V/qmnKf/kk1Tjz9Xvn9mhP/1nh+pb9uqkuaX68TuP1klzS6nNAQAAwIiSlKoda+2/Jf07GftKF+MzfbRpPkheTVrT02tlsrNVcMopydt3GuWfeKLm/POfqv3jH1Vzy63a9vaLlX/yyRr3rktVuGKFTDB4SPu31iq8Y4dCzz6r5lWr1PLsasWamyW/X7lHHaXyqz+p/FNPU87iw2R8Pr1cUa/f3L9N9728V1FrdfZhE/XxFXO5zhkAAABGrJHZ/i4VjEnupV9C1bLBQjU9/JjyTzpJvrwMHEb5IPlyclT20Y9q/Lvepbo//Ul1f7tduz95tXzFxSo49VQVnHaqcpYcoaxZM/u9hpiNxRTes1cd27aqbcPral27Vq3r1rkmjJICUyar6PzzlX/qKco/4YTOi4BHojH9d/1+3bpqm17YUaeC7IDef+IsXX7SLM0oHT2vMwAAAMYmQponFTVp7a0liuzZq8KrrkrefjOIv7BQZVdeqdKPflTNTzyppgceUPPjj6vxPjfEqcnNVXDSJAXKy11I9ftlwx2K1jcoWl+vSGWlbFtb5/6yZs1SwYoVyj36aOUdt0xZs2d3a6q4q7ZFt7+wS397fpcqm9o1oyRPX7twsd6xbJoKcw6tBg8AAADIFIS0OJ+S2yctVKWmPbmSiahgxYrk7TcDmUBAhWeeocIzz5CNRtW+aZPaXtug9o2vK7y/UpFKd7PRqEwgIP+4ccqaPl2B8nJlzZmt7DlzlDV3rgLjD2yiGI7G9PCGSt323E49salKRtKKhRP07uNn6MxFE+QfgYOBAAAAAP0hpHmM8SX3ooWhaoV2RZWzZIkCpaVJ3HFmM36/chYtUs6iRYe0n501Lfrr8zv19xcrVNXUrklFObr6zPl653HTNXVcbpJKCwAAAGQeQlpckmvSorVVaq0IqPTCzL34c6bpiMT00Ib9uu25nXpyU7V8Rjpzkas1O31BuQL+vvu3AQAAAKMFIc1jfEmsSYvFFNrWJNnxKjj11CTtdPTaWtWsvz2/S3e8WKGaUIemFOfoU2ct0DuPm6bJxdSaAQAAYGwhpMX5kji6Y2utmvdmyZebpdwjj0zSTkeXtnBU96/fp9ue26lnt9bK7zNa6dWanbagnL5mAAAAGLMIaXEmeaM72uZKhfbmKP+o+TIBXuJEe+pb9ftntutvz+9SfUtYM0ry9L/nLtQ7jp2mCUU56S4eAAAAkHYkCI/xmaSNwN+x4WVFWv3KP35pcnY4CqzdWadbV23Tf17dJ2utzlsySe9ZPlMnzimVj1ozAAAAoBMhLc5npCQ1d2x+5jlJGvP90ay1WrW5Wj97aJNe2FGnwpyAPnzKbL3/xJmaNp6LTgMAAAC9IaR5XE1acqrSQmteU1ZRWME5i5Oyv5HGWqtnttToxw++oRd21GlKcY6ue/NiXbJsugqyOeUAAACA/vAXc5wxSRnd0XZ0qGVjhcbN7JByD7w482i3cV+Trr9vvZ7aXKNJRTn65luX6J3Lpik74E930QAAAIARgZDmMT5fUkZ3bH11vWxHVHkzsiXf2Akm9S0d+smDb+hPq3eqIDugr124WJctn6Gc4Nh5DQAAAIBkIKTF+ZJTk9by3GpJUt6ckkPf2Qhx//p9+vLdr6g21KH3LJ+pT5+9QOPzs9JdLAAAAGBEIqR5kjW6Y2j1amWXBxUom3DoO8twDS1hXXfvet29drcWTy7S7z54vJZMLU53sQAAAIARjZAW5/Mdck1arKNDrWvXadz8mJQ/ukPa2p11uurPa1TZ1K5rVs7XVWfMU1bAl+5iAQAAACMeIc1jvOaO1loZc3DX7Wp7+WXZtjbll7ZK+eVJLmFmsNbqT6t36vp712tiUY7u+p+TdOS0cekuFgAAADBqENLi4hdUjsUk/8ENdhFavVoyRnkl9VJ+WfLKliE6IjF96e5XdMeLFTpjYbl+cunRGpdH3zMAAAAgmQhpHuPzmupFowcd0lpWP6fs+XPlz9o96mrSmtsjuvJPL+rJTdW6euV8Xbtyvny+g6txBAAAANA3OhHFeYHDxg5uHP5Ye7ta161T/pHz3YxRFNKqmtr1rpuf0dNbavSDS47Up89eQEADAAAAUoSaNE+3mrSD0LruJdmODuUtmiZt1qgJafsb23Tpr57R/sZ2/fr9x+rMRRPTXSQAAABgVKMmLe4Qa9JaVq+WfD7lzfaGoB8FfdKqmtp12a+fVVVTu/70keUENAAAAGAYUJPmOdSatJbnn1fO4sXy2yY3Y4TXpNWGOvTeW1ZrT32bfv+h43XszPHpLhIAAAAwJlCTFncINWm2o0OtL7+svGOPlUJVUjBPyi5IdgmHTag9og/85jltrwnp1g8s0/GzS9JdJAAAAGDMIKR5DqUmre2112Tb25V7zDEupI3gpo7RmNU1f12r9XsadON7j9FJ80bucwEAAABGIkJanBfSDqYmrWXNWklS3jFLvZA2cps6fuffG/TQhkp946LD6YMGAAAApAEhzWPiQ8ofRE1a69o1Cs6YoUB5+YgOaX96doduXbVNl580S+87cVa6iwMAAACMSYS0OJ+7gPVQa9KstWp5cY3yli51M0LVI7K54wvba/X1e9brjIXl+uqFi9NdHAAAAGDMIqR5DrYmLbxjh6K1ta4/mrUjsiatprldn/jLWk0bn6ufvXup/FyoGgAAAEgbhuCP83t90qJDq0nr7I927DFSW70Ui4yokBaLWX3q9pdU29Khu648SUU5wXQXCQAAABjTqEnzGOO9FLGh1aS1rl0jX3GxsubMcU0dpREV0m58fIueeKNKX3/zYi2ZWpzu4gAAAABjHiEt7mBr0l5co7yjj3ZD+DdXupkjpE/amp11+tEDG3XRUVN02fEz0l0cAAAAACKkdYn3wxpCTVqkrk4dW7e6/miS648mSfkTkly45GvtiOqzt7+kycW5+vbblsgY+qEBAAAAmYA+aR7j90Z3HEJNWuvadZK8/mhSQkjL/OaOP7j/dW2tDukvH1muQvqhAQAAABmDmrQ4M/SatNY1L0rBoHKWLHEz4n3S8kqTXLjkenpLtX771HZdftIsnTRvZDTNBAAAAMYKQpqnqyZt8CGtZc1a5S5eLF9OjpsRqpJySyR/5lZQhtoj+twdL2t2Wb4+f96idBcHAAAAQA+EtLjOPmmDa+4Ya29X2yuvKPfYY7tmjoBrpP3kwTdUUdeqGy45UrlZ/nQXBwAAAEAPhDTPUGvS2tavlw2HlXfM0q6ZoeqMDmmv7m7Qb5/ersuWz9CyWSXpLg4AAACAXhDS4jqvkza4mrTWNWskSblLE0NaVcYOvx+NWX357lc0Pi+oz59LM0cAAAAgUxHSPKbzOmmDq0lrWbNWWbNmKVCaMEhIBjd3/PPqHXqpokFfvXCxivMYzREAAADIVIS0OC+kaRAhzVqr1jVruq6PJkmRDqmtXirIvGukVTa16Yb/btQp88p00VFT0l0cAAAAAP0gpHmMidekRQZct2PbNkXr67v3R2vxht/PwOaON/x3o9oiUV3/lsO5aDUAAACQ4Qhpcd7AIYoMXJPW2R8tsSYtQy9k/UpFg+5YU6EPnjxbc8oL0l0cAAAAAAMgpHmM39UwDaYmrWXtWvnHjVPW7NldMzMwpFlrdf1961WSl6VPnDkv3cUBAAAAMAiEtDifV5MWG0xN2lrlLl3avelgKN7cMXNC2r9e2avnt9fps+cuVFEOg4UAAAAAIwEhzWN88YFD+q9Ji9TVqWPbtu5D70sJNWmZ0SetLRzVd//9ug6bXKR3Lpue7uIAAAAAGCRCWlz8YtYD9ElrXbtOkroPGiK5kObPkrKLUlG6Ifv1E1u1u75VX7twsfw+BgsBAAAARgpCmqerJm2gkLZGCgaVs2RJ9wWhatfUMQNGT6xqateNj2/ReYdP0olzSwfeAAAAAEDGIKTFxWvSBuiT1rJ2rXIWHyZfTk73Bc2VGdMf7f89skkdkZg+f/6idBcFAAAAwBAR0jyD6ZNmOzrU9sqryjt66YELQ1UZEdJ21rToL6t36tLjpmt2WX66iwMAAABgiAhpcfGatH6aO7Zt2CDb3t79+mhx8eaOafajBzcq4De6euX8dBcFAAAAwEEgpHkG0yetZc1aSVLu0qO7L7DWq0lL78iO6/c06J/r9uhDJ8/WxKKcgTcAAAAAkHEIaXFeSOuvJq117VoFp01TcMKE7gvam6Roe9pr0m64f6OKc4P62Olz01oOAAAAAAePkOYx/oC700dIs9aqZe0a5fYcel9KuEZa+kLas1tr9NjGKv3PirkqzuXC1QAAAMBIRUiL8/dfkxauqFC0qlp5PS9iLbn+aFLaQpq1Vj/47+uaVJSjD5w0Ky1lAAAAAJAchDSP8fc/umPrWq8/Wq+DhsRr0tLTJ+3B1/Zrzc56XXvWfOUE/WkpAwAAAIDkIKTFDdAnrWXNGvkKCpQ9b96BC0OVbpqGmrRozOqG+zdqTnm+Ljl22rAfHwAAAEByEdI8A/VJa127TrlHHSXj76WmKo3NHe9aU6FNlc3633MWKuDn7QQAAABGOv6qj4v3SYsdGNKiTU1qf+MN5fbWH01yzR1ziqVAVipLeIC2cFQ/fWiTjppWrPOWTBrWYwMAAABIjUMKacaYdxhj1htjYsaYZckqVFr4vItZRw4Maa0vvSxZq7zeRnaUvGukDX8t2p+e3aHd9a36/HmLZIwZ9uMDAAAASL5DrUl7VdLbJT2RhLKklTE+GX9Mtr3tgGWta9ZIPp9yjjyq941D1cMe0prawvrFo5t16vwynTQvvRfRBgAAAJA8gUPZ2Fq7QdKoqMVpj1r5AlaxltYDlrWsXaPshQvlL8jvfeNQlVQ2P8Ul7O7XT2xVXUtYnzt30bAeFwAAAEBq0SdNUixmdcuq7fIFrMKhlm7LbDis1pdeVt7So/vewTA3d6xqatctq7bpgiMn64hpxcN2XAAAAACpN2BNmjHmIUm9jUrxZWvtPwd7IGPMFZKukKQZM2YMuoDDweczOmleuXxBq+0VNUosXduGDbItLco77rjeN45GpJbaYQ1p/++RTeqIxPTZcxYO2zEBAAAADI8BQ5q19qxkHMhae7OkmyVp2bJlNhn7TKalM0q0PWBVV9OoulCHxue7kRpbnn9ekpS3rI9xUVprJdlhC2k7akL6y+qduvS46Zpd1kfzSwAAAAAjFs0d44xPvkBMwXCH/vDMjs7ZLc89r6zZsxUo7yOENQ/vhax/9MAbCvp9umbl8PaBAwAAADA8DnUI/rcZYyoknSjpX8aY+5NTrDQwRr6A1XhfVL9/ZrtaO6Ky0ahaXnyx76aOkuuPJkn5qR9h8dXdDbrnpT360CmzNKEoJ+XHAwAAADD8DnV0x7sl3Z2ksqSdL2BVbKKqDXXo7rW79baCJsWam5V3XD+XgAtVu2n+hJSX7wf3b9S4vKA+dvrclB8LAAAAQHrQ3DHO+OQLWgXCHVoytUi3rtqq0PMvSNIANWnx5o6prUl7eku1nnijSletmKeinGBKjwUAAAAgfQhpccbIF4gp1tqmj5wyR1uqQtr16CoFp09XcFJvg1t6QlWSLyDljk9Z0ay1+v5/N2pycY7ed+LMlB0HAAAAQPoR0uKMT76Ale0I6/xFZZpcmKXYS2v7r0WTpGbvGmkpvKD3fS/v1Uu76vWpsxYoJ+hP2XEAAAAApB8hrZMbOESSAuF2fXymlNcWUsP8Jf1vluILWbeFo/ref17X4slFuvjYaSk7DgAAAIDMQEiLMz4ZL6TFWlq0snGrJOm2WD9NHSXXJy2FIe2WJ7dqd32rvnrhYvl9qautAwAAAJAZCGlxXp80yYW06LNPq27yTN22tVX7Gtr63i5ULRWkZmTH/Y1t+uVjW3Te4ZN04tzSlBwDAAAAQGYhpMV5fdIkKVJZqZY1azRh5QpZK930+Jbet7HWXcw6RSM7/uC/GxWJWn3pTYelZP8AAAAAMg8hrVNXn7SmBx6UIhFNPudMvf2YqbrtuZ2qbOylNq29SYq2p+QaaWt31unONRX60CmzNaM0L+n7BwAAAJCZCGlxxqfscRHJZ1T3l7/IV1CgvKVLddUZ8xSJWf3qia0HbhOqctMk90mLRGP60t2valJRjj5x5ryk7hsAAABAZiOkxRmjQHZMgfHjJElF558nEwxqZmm+3nr0VP159Q5VNbV33yYe0gqSG9J+9/R2bdjbqOsuWqyC7EBS9w0AAAAgsxHS4rzrnJVe9mZJUtknPtm56BNnzlM4avXzhzd136a50k2TWJO2u75VP37wDa1cNEHnHj7AyJIAAAAARh1CWicX0sZfuEILX35JwYld/cxml+Xrvctn6M+rd2jjvqauTTqbOyanT5q1Vl+482VZK33jLYfLpPAC2QAAAAAyEyEtzriXwkjyZWUdsPjasxaoIDugb/3rNVnrBhjpCmnJGd3xT8/u0JObqvWlCw7TtPEMFgIAAACMRYS0uHitlY31unh8fpauOWuBntxUrUc3es0cQ1VS7njJHzzkw2+tata3/71Bpy0o13uXzzjk/QEAAAAYmQhpneJNC22fa7zvhJmaU56vr/1zvULtEe8aaYfeH60tHNU1f12n7IBfN1xyJM0cAQAAgDGMkBbnNXeU7TukZQV8+v7FR2p3fau+95/XpVD1IfdHs9bqK/94Va/sbtANlxypiUU5h7Q/AAAAACMbIS1ugOaOccfNKtEHT5qtPz67Q6G6vYfcH+2Pz+7QHS9W6OqV83UOozkCAAAAYx4hLS5ek9ZPc8e4z523UIsnFynaVKmmQMlBH/LhDft1/b2vaeWiCbp25fyD3g8AAACA0YOQ1mlwNWmSlBP066Z3L1GRQrpzY5tqQx1DPtpTm6t15Z/XaPGUIv30XUfL56MfGgAAAABCWpdB9ElLNCO7RZK0pSVP7//N6iEFtfte3qMP/vZ5zS7N1+8/eLwKcw59dEgAAAAAowMhLa6zT9rgQppCbhj+d55+jN7Y36y3/GKVNuxt7HeT9khU37rvNX3iL2t11PRi3f6xEzU+/8BrsgEAAAAYuwhpcUPokybJjewo6YiF8/W3K05QWzimN/+/VfrWfa+poq6l+6rtEd3xYoXO/vETumXVNr3/xJn644eXqziPGjQAAAAA3QXSXYDMMfg+aZLcNdIkKb9MS0vG67/XnKob7t+oW5/apltWbdP8CQWaVJyjxtawXt/XpPZITIdPKdIfP3y8Tp1/6NdWAwAAADA6EdLihtzcscpNveuklRZk63sXH6mrzpine17ao5d21Wt/Y5uKcoN6z/KZOnvxRC2fXcIAIQAAAAD6RUiLG+R10jqFqqRArpSV32329JI8XXXGvCQXDgAAAMBYQZ+0uCH3SauSCsq7wh0AAAAAJAEhrdNB9EnLp28ZAAAAgOQipMUNuU9adWd/NAAAAABIFkJaXOfFrAfbJ61Syi9LXXkAAAAAjEmEtE7xvmWDqEmLxVxNWgE1aQAAAACSi5AW11mTNoiQ1lon2Sh90gAAAAAkHSEtbihD8HdeI42QBgAAACC5CGlxZggvRajSTQlpAAAAAJKMkNZpCDVpzV5Io08aAAAAgCQjpMUNZQj+zpA2MXXlAQAAADAmEdLihtInrXmf5AtKueNTWyYAAAAAYw4hLa6zT9oga9IKJnYFOwAAAABIEkJapyHUpDXtkwpp6ggAAAAg+QhpcUO5Tlq8Jg0AAAAAkoyQFjekPmn7CWkAAAAAUoKQFjfYPmnRsNRSTUgDAAAAkBKEtE6DrEkLVbkpfdIAAAAApAAhLW6w10lr2uem1KQBAAAASAFCWlznwCED1KR1Xsh6UmrLAwAAAGBMIqR1GuQ1z5rjNWkTUlcUAAAAAGMWIS1usEPwd9akEdIAAAAAJB8hLW6wQ/A37ZNyx0uB7NSXCQAAAMCYQ0iLi4e0gYbgb95PfzQAAAAAKUNI6zTImrTm/TR1BAAAAJAyhLS4QfdJ2y8VUpMGAAAAIDUIaXGD6ZNmrdRETRoAAACA1CGkxcVr0vrrk9bWIEXb6ZMGAAAAIGUIaZ0GUZPWvN9NCyamvjgAAAAAxqRAuguQMQbTJy0e0goJaQAAAMBwC4fDqqioUFtbW7qLMmg5OTmaNm2agsHgoLchpMUNpk9aEzVpAAAAQLpUVFSosLBQs2bNkum8hFbmstaqpqZGFRUVmj179qC3o7lj3GD6pNHcEQAAAEibtrY2lZaWjoiAJknGGJWWlg655o+Q1ilek9ZfSNsn+bOlnOLhKRIAAACAbkZKQIs7mPIS0uLMIEJa037XH22EnRgAAAAARg5CWtxgmjs27ZUKpwxLcQAAAACMTYcU0owxNxhjXjfGvGyMudsYMy5J5Rp+gxk4pHGPVERIAwAAAMai559/XkceeaTa2toUCoV0+OGH69VXX036cQ51dMcHJX3RWhsxxnxf0hclff7Qi5Umxi/FIr0vs9aFtIXnD2+ZAAAAABzgG/eu12t7GpO6z8VTivT1Nx/e5/LjjjtOF110kb7yla+otbVV733ve7VkyZKklkE6xJBmrX0g4eGzki45tOKkmT8oRcO9L2urlyKtUuHkYS0SAAAAgMzxta99Tccdd5xycnL085//PCXHSOZ10j4k6W99LTTGXCHpCkmaMWNGEg+bRP6svmvSGve6Kc0dAQAAgLTrr8YrlWpqatTc3KxwOKy2tjbl5+cn/RgD9kkzxjxkjHm1l9tbEtb5sqSIpD/3tR9r7c3W2mXW2mXl5eXJKX2y+QJStKP3ZY173JSQBgAAAIxZH/vYx/TNb35T73nPe/T5z6emp9eANWnW2rP6W26MuVzShZJWWtvf+PUjgD+r7+aOTV5Io7kjAAAAMCb94Q9/UDAY1GWXXaZoNKqTTjpJjzzyiM4888ykHueQmjsaY86T9DlJp1trW5JTpDTqr09avLkjIQ0AAAAYk97//vfr/e9/vyTJ7/dr9erVKTnOoV4n7f8kFUp60BizzhhzUxLKlD7+oBTrpyYtv1wKZA1vmQAAAACMKYc6uuO8ZBUkI/iC/fdJoxYNAAAAQIodak3a6OLPkqL9jO7IoCEAAAAAUoyQlsjfz+iOTXsIaQAAAABSjpCWyJ/Ve5+0cJvUUiMVEtIAAAAApBYhLZGvj9Edm+IXsqZPGgAAAIDUIqQl6msI/oZdblo8bXjLAwAAAGDMIaQl8vcxumNDhZsWTx/e8gAAAAAYcwhpifxZUqyX0R3rqUkDAAAAxrovfOEL+sUvftH5+LrrrtMPf/jDpB/nkK6TNur4+hjdsWGnVDBJCmQPf5kAAAAAHOg/X5D2vZLcfU46Qjr/e30uvvTSS3XttdfqqquukiTdfvvtuv/++5NbBhHSuvNn9d4nrX6XNI6mjgAAAMBYtnTpUlVWVmrPnj2qqqrS+PHjNX168nMCIS1RfwOHTD5q+MsDAAAAoHf91Hil0jve8Q7dcccd2rdvny699NKUHIOQlsgfPPA6abGYGzhk0YXpKRMAAACAjHHppZfqox/9qKqrq/X444+n5BgMHJLI18vojqFKN2/cjPSUCQAAAEDGOPzww9XU1KSpU6dq8uTUXEeZmrREvfVJY/h9AAAAAAleeSXJA5b0QE1aIn/gwJBWv9NNGTgEAAAAwDAgpCXyZx3Y3LGBa6QBAAAAGD6EtES+oCQrxaJd8+p3SdnFUk5x2ooFAAAAYOwgpCXyB900sTatfgeDhgAAAAAYNoS0RJ0hLaFfWu1WqXROesoDAAAAYMwhpCXyZ7lpPKRFI1LddqmEkAYAAABgeBDSEvm8KxLEL2jdsFOKRaSSuekrEwAAAIAxhZCWqLMmzeuTVrvVTUsJaQAAAACGByEtUc8+aTVeSKO5IwAAAABPNBrVNddco8MPP1xHHHGEtm7dmtT9E9IS9QxptVukrAKpYGL6ygQAAAAgo3z3u9/VnDlztH79el199dX65S9/mdT9B5K6t5HO54W0eJ+02q1SyWzJmPSVCQAAAEDGCIVCuvvuu/Xiiy9KkmbPnq1//etfST0GIS1Rzz5pNVukSUvSVx4AAAAAvfr+c9/X67WvJ3Wfi0oW6fPHf77fdR566CHt2rVLRx99tCSptrZWZ511VlLLQXPHRH4vs0Yj7la/g5EdAQAAAHRat26drr/+eq1bt07r1q3TOeec0xnYkoWatESJNWnx4fcZ2REAAADIOAPVeKVKXV2dZs+eLUmKRCJ64IEH9OUvfzmpx6AmLVEgx00jbVKlV3VatiB95QEAAACQURYsWKBnn31WkvSTn/xEF1xwQWdoSxZCWqL8cjdt3i9VvubuTzgsfeUBAAAAkFHe/e53a82aNZo3b55efvll/fjHP076MWjumKhwkps27XMhbdwMKbswvWUCAAAAkDHGjx/fWZOWKtSkJQrmSjnjXEjb/5o04fB0lwgAAADAGENI66lwklSzWap+g+H3AQAAAAw7mjv2VDhJ2vqouz/jhPSWBQAAAMCYQ01aT4WT3dT4pGnHp7csAAAAAMYcQlpPeaVuOuUYKacovWUBAAAAMObQ3LGnJRdLzZXSmV9Jd0kAAAAAjEGEtJ6mHiNd/Ot0lwIAAADAGEVzRwAAAADIIIQ0AAAAABiCe+65RxdffHG3eTfeeKM++clPJmX/NHcEAAAAMOLs+8531L7h9aTuM/uwRZr0pS8NuN6Xv/xl3Xbbbd3mzZ07V3feeWdSykFNGgAAAAAM0ksvvaRYLKYlS5Zox44duvHGGyVJ4XBYxpikHIOaNAAAAAAjzmBqvFJh3bp1OvbYYyVJDz74oDZt2iRJeu2113TUUUcl5RjUpAEAAADAIMViMTU3Nysajequu+5SU1OTWltb9bvf/U6XXXZZUo5BSAMAAACAQXrTm96krVu36uijj9bHP/5xrV+/XsuWLdMVV1yhY445JinHoLkjAAAAAAzSxIkTtW7dus7HF110UdKPQU0aAAAAAGQQQhoAAAAAZBBCGgAAAABkEEIaAAAAgBHDWpvuIgzJwZSXkAYAAABgRMjJyVFNTc2ICWrWWtXU1CgnJ2dI2zG6IwAAAIARYdq0aaqoqFBVVVW6izJoOTk5mjZt2pC2IaQBAAAAGBGCwaBmz56d7mKkHM0dAQAAACCDENIAAAAAIIMQ0gAAAAAgg5h0jIxijKmStGPYDzywMknV6S4ERi3OL6Qa5xhSifMLqcT5hVTLxHNsprW2vLcFaQlpmcoY84K1dlm6y4HRifMLqcY5hlTi/EIqcX4h1UbaOUZzRwAAAADIIIQ0AAAAAMgghLTubk53ATCqcX4h1TjHkEqcX0glzi+k2og6x+iTBgAAAAAZhJo0AAAAAMgghDRJxpjzjDEbjTGbjTFfSHd5MDIZY6YbYx41xrxmjFlvjLnGm19ijHnQGLPJm4735htjzM+98+5lY8wx6X0GGAmMMX5jzFpjzH3e49nGmNXeefQ3Y0yWNz/be7zZWz4rrQVHxjPGjDPG3GGMed0Ys8EYcyKfX0gmY8ynvO/HV40xtxljcvgMw8EyxvzGGFNpjHk1Yd6QP7OMMR/w1t9kjPlAOp5Lb8Z8SDPG+CX9QtL5khZLercxZnF6S4URKiLpM9baxZJOkHSVdy59QdLD1tr5kh72HkvunJvv3a6QdOPwFxkj0DWSNiQ8/r6kn1hr50mqk/Rhb/6HJdV583/irQf052eS/mutXSTpKLnzjM8vJIUxZqqkqyUts9YukeSX9C7xGYaD9ztJ5/WYN6TPLGNMiaSvS1ou6XhJX48Hu3Qb8yFN7g3ZbK3daq3tkPRXSW9Jc5kwAllr91pr13j3m+T+wJkqdz793lvt95Le6t1/i6Q/WOdZSeOMMZOHt9QYSYwx0yRdIOkW77GRdKakO7xVep5f8fPuDkkrvfWBAxhjiiWdJulWSbLWdlhr68XnF5IrICnXGBOQlCdpr/gMw0Gy1j4hqbbH7KF+Zp0r6UFrba21tk7Sgzow+KUFIc39Eb0r4XGFNw84aF6zjKWSVkuaaK3d6y3aJ2mid59zD0P1U0mfkxTzHpdKqrfWRrzHiedQ5/nlLW/w1gd6M1tSlaTfes1pbzHG5IvPLySJtXa3pB9K2ikXzhokvSg+w5BcQ/3MytjPMkIakGTGmAJJd0q61lrbmLjMuuFUGVIVQ2aMuVBSpbX2xXSXBaNSQNIxkm601i6VFFJXMyFJfH7h0HhNyN4i94PAFEn5ypAaC4xOI/0zi5Am7ZY0PeHxNG8eMGTGmKBcQPuztfYub/b+eDMgb1rpzefcw1CcLOkiY8x2uWbZZ8r1IRrnNR2Sup9DneeXt7xYUs1wFhgjSoWkCmvtau/xHXKhjc8vJMtZkrZZa6ustWFJd8l9rvEZhmQa6mdWxn6WEdKk5yXN90YXypLrxHpPmsuEEchrK3+rpA3W2h8nLLpHUny0oA9I+mfC/Pd7Iw6dIKkhoYoe6MZa+0Vr7TRr7Sy5z6lHrLXvkfSopEu81XqeX/Hz7hJv/RH7iyJSy1q7T9IuY8xCb9ZKSa+Jzy8kz05JJxhj8rzvy/g5xmcYkmmon1n3SzrHGDPeq+09x5uXdlzMWpIx5k1yfT38kn5jrf12ekuEkcgYc4qkJyW9oq4+Q1+S65d2u6QZknZIeqe1ttb7kvo/ueYeLZI+aK19YdgLjhHHGLNC0mettRcaY+bI1ayVSFor6b3W2nZjTI6kP8r1jayV9C5r7dY0FRkjgDHmaLlBabIkbZX0Qbkfc/n8QlIYY74h6VK50ZDXSvqIXP8fPsMwZMaY2yStkFQmab/cKI3/0BA/s4wxH5L7e02Svm2t/e0wPo0+EdIAAAAAIIPQ3BEAAAAAMgghDQAAAAAyCCENAAAAADIIIQ0AAAAAMgghDQAAAAAyCCENAAAAADIIIQ0AAAAAMgghDQAAAAAyyP8HAgRhH8a9GzIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15, 8))\n",
    "plt.title(\"Closed-loop system response of inverted pendulum on a cart stabilized with an LQR controller\")\n",
    "plt.plot(sol, label=[\"x\", \"v\", r\"$\\theta$\", r\"$\\omega$\"])\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "034f9025",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for i in range(len(t)):\n",
    "    cart(sol[i][0], sol[i][2], i)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0a120197",
   "metadata": {},
   "source": [
    "![pendulum](https://user-images.githubusercontent.com/45606539/119007993-a6800d00-b989-11eb-9bbd-1009e790448a.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c616232",
   "metadata": {},
   "source": [
    "## References\n",
    "[1] Brunton, S., & Kutz, J. (2019). _Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control._ Cambridge: Cambridge University Press. doi:10.1017/9781108380690  \n",
    "[2] Brunton, S, (2017). _Linear Quadratic Regulator (LQR) Control for the Inverted Pendulum on a Cart [Control Bootcamp]_, https://www.youtube.com/watch?v=1_UobILf3cc  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}