linear_SMPC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "id": "XfJCAjmGzuSF"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import casadi as ca\n",
    "import chaospy\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "U1lEeVzDGb4j"
   },
   "outputs": [],
   "source": [
    "normal = chaospy.Normal(mu=0, sigma=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "id": "NniqSRfkzx7C",
    "outputId": "7bf00b3e-201b-410e-d942-73bce72ead7c"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZP0lEQVR4nO3de5BU5ZnH8d/TjSgE5CIDymUAE6JiJV4YwWtiKmuC5EKylagkrhpiiKmYS23VliapJFWbpGqz2c0mbjRIvERzI6nEJJRLouZqthQFjDcgIhAFBGUQQRC5TJ9n/+juYXbs7vOeoaffnub7qaLonj7T/XCY/s3p57znfc3dBQAY+HKxCwAA1AeBDgAtgkAHgBZBoANAiyDQAaBFDIr1wmPGjPEpU6bEenkAGJBWrly53d3bKj0WLdCnTJmiFStWxHp5ABiQzOzZao/RcgGAFkGgA0CLINABoEUQ6ADQIgh0AGgRBDoAtAgCHQBaRLRx6EAz2L3voNa+sFsvvLxfO145oJdeOaCCuwYPymlwPqfjhg1W++ihmjR6qNqGHS0zi10yUBWBjiPK1l2v6k9PdeovT3fqied2adOOV4O/d/TrBmvG5FE6a8oonfv6MTp1/LEEPJoKgY6W9+Ke/frlX5/Tz1du1t+e3y1JGj/iGJ0xeZQuO6tdJx8/XBNGDdHooYM1YuhRGpTL6WAh0f6uRJ2792vTjr3auGOvnnxul1Y8+5LuW/2CJGnKcUP1rjefoPedPkHTxg2P+U8EJEkWa8Wijo4O59J/9KfHNu3Uovs36N7Vz+tgwXX6pJGa86bjdeFJYzVt7LA+H11v271Pf1izTXc/vlUPrN+uxKULpo3R/POm6q1vbFMux1E7+o+ZrXT3joqPEehoNSue2aEb/rBO96/t1IghR+mDMybqkrMm6Y39cBS9fc9+/XT5Jt354DN64eX9esPYYfqXd56kd0wfRzsG/YJAxxFh0469+tr/rNFvVz2v4143WFdfcKIuP7tdw485qt9f+2Ah0dInturbv39aGzpf0ZntI/W5OaforCmj+/21cWQ57EA3s9mSvi0pL+kWd/+3Xo+PkPRDSe0q9uX/w91vr/WcBDrqZd/Bgm7603rd/Of1ypnpk297veafP1VDBzf+FFFXIdHPV27Wt373tJ5/eZ8u6Zioz885RSOHDm54LWhNhxXoZpaXtFbSRZI2S1ouaZ67r+6xzecljXD368ysTdJTko539wPVnpdARz08vnmn/vlnj2ndtj1672nj9bk5J+uEEUNil6VXDxR0wx+e1qL7N2jkkKP0pfdM19zTJ8QuCy2gVqCHXFg0U9I6d99QCujFkub22sYlDbdi03CYpB2Sug6jZqCmg4VE37xvrd5/0wPas69Ld8yfqRvmndEUYS5JQwbndd3sk3X3p87XxNFD9ZnFj+qzi/+q3fsOxi4NLSzkM+kESZt63N8saVavbb4jaYmkLZKGS7rU3ZPeT2RmCyQtkKT29va+1Avo+V37dO2PH9GKZ1/SP545QV9+z6kaMaT/++R9ccoJx+quT5yrG/+4Tt/63Vo9snGnbph3hk6fNDJ2aWhBIUfolU7V9+7TvFPSo5LGSzpd0nfM7NjXfJP7InfvcPeOtraKKygBNT2wbrve/d9/0eqtL+uGeWfom5ec3rRhXpbPmT799mn62cfPUSFxfeC7D+jHD22MXRZaUEigb5Y0qcf9iSoeiff0EUl3edE6SX+XdHJ9SgQkd9et//t3XX7rQxo5dLCWXHue3nva+NhlZdIxZbSWfvoCnfuGMfr8L5/QF3/1pA4WXvNBFuizkEBfLmmamU01s8GSLlOxvdLTRklvlyQzGyfpJEkb6lkojlxdhURfXrJKX7l7tS6aPk6//uR5esPYgXll5oihR+n2q87Sx99yon6w7FldfstD2rm36tgBIJPUQHf3LknXSrpH0hpJP3P3VWZ2jZldU9rsK5LONbMnJP1e0nXuvr2/isaRY8/+Ln3szhW688Fn9fG3nKjvfniGXnf0wJ6xIp8zfW7OKfqvS0/TXzfu1AcXPqitu8LnlAGq4cIiNK2dew/oytuX68nndulf556qD8+aHLukuntg/XYtuHOljj1mkO6YP5M5YZDqcIctAg237eV9uvTmZVqz9WUtvHxGS4a5JJ37+jH66cfP1sHE9YGFD+qxTTtjl4QBjEBH09n80l5dcvOD2vTSXt1+1Vm6aPq42CX1q1PHj9BdnzhXxw4ZpMtveUiPbHwpdkkYoAh0NJUtO1/VvO8t04uvHNAPPjpL571hTOySGmLS6KH66YJzNHrYYF1x68Na8cyO2CVhACLQ0TSe37VP8763TDtfOagffnSWZkweFbukhho/coh+uuAcjR1+tK647WGtfJZQRzYEOprCtt379KFblmn77v36/vyZOu0IvZLy+BHHaPGCszXu2GN01e3LtWrLrtglYQAh0BHdrlcP6opbH9bWnfv0/fkzj7gj897GHnuMfnj1LA0/epCuuPVhre/cE7skDBAEOqLad7Cgj925Qus792jRFTOYP7xkwsgh+uHVs2QmXX7LQ3puJ+PUkY5ARzRdhUSf/slftfyZHfrPS07XBdOY36enE9uG6c75s7Rnf5euvO1hrihFKgId0dx8/wbdu/oFffnd0wfcvCyNMn38sfreFR3a+OJefe6uJ2KXgyZHoCOa9Z17NGHkEF113tTYpTS1s088Tm87uU0bOl+JXQqaHIGOaNylQXkWUg4xKJdTEmmaDgwcBDqiKSSunBHoIcykAoGOFAQ6okncRZ6HyZmJPEcaAh3RuIsj9EA5Ey0XpCLQEU3irhx5HiRnRqAjFYGOaOihh8vlTAmr1SEFgY5oElouwWi5IASBjmjcXTl+AoPQckEI3k6IpthD5wg9hJkpIc+RgkBHNAVaLsHyOSkh0ZGCQEc0ziiXYLRcEIJARzS0XMLlaLkgAIGOaJKElksoY5QLAhDoiCZhlEuwPJf+IwBvJ0RDyyVcLmcq0HNBCgId0XBhUThaLghBoCMaZlsMx2yLCEGgI5rEpTzjFoPkGbaIAAQ6okmYnCtYjgUuEIBARzRMnxvOSi0XJ9RRA4GOaBIvBhXSlT/JkOeohUBHNFz6H668n+ijoxYCHdEUEuekaKBcaT/RR0ctBDqiKQ5bJNBD0HJBCAId0bBIdDhaLghBoCMaRrmEK//i4+p/1EKgI5qCu/IcoQfp7qGT6KiBQEc0ScKwxVDlTzKMQ0ctBDqiYdhiOFouCBEU6GY228yeMrN1ZnZ9lW0uNLNHzWyVmf25vmWiFTHbYjhOiiLEoLQNzCwv6UZJF0naLGm5mS1x99U9thkp6SZJs919o5mN7ad60UIK7t29YdRW3k8sFI1aQo7QZ0pa5+4b3P2ApMWS5vba5kOS7nL3jZLk7tvqWyZaES2XcLRcECIk0CdI2tTj/ubS13p6o6RRZvYnM1tpZldUeiIzW2BmK8xsRWdnZ98qRsug5RKOlgtChAR6pXdc75+qQZJmSHqXpHdK+qKZvfE13+S+yN073L2jra0tc7FoLYxDD2fdR+gEOqpL7aGreEQ+qcf9iZK2VNhmu7u/IukVM7tf0mmS1talSrSkJKGHHirPpf8IEHKEvlzSNDObamaDJV0maUmvbX4t6QIzG2RmQyXNkrSmvqWi1dByCZcrvVO5sAi1pB6hu3uXmV0r6R5JeUm3ufsqM7um9PhCd19jZr+V9LikRNIt7v5kfxaOgY+WS7gcLRcECGm5yN2XSlra62sLe93/hqRv1K80tLpioJPoIYxRLgjAlaKIhhWLwh3qoZPoqI5ARzRJ4srzExik3JpigQvUwtsJ0dByCdfdckkiF4KmRqAjGlou4biwCCEIdERR7gUzyiUMS9AhBIGOKMrjqVngIkyeRaIRgEBHFOXhd1wpGsZouSAAgY4oysHEAXqYHMMWEYBARxTlXGKUSximz0UIAh1RlHvB9NDDMJcLQhDoiIKWSzbM5YIQBDqi8NIFMrRcwjBsESEIdESRMA49Ey4sQggCHVGUgylPogfpXiSaPEcNBDqiKHT30An0EN09dBIdNRDoiIJhi9nQckEIAh1R0EPPhnHoCEGgIwou/c+GYYsIQaAjinIvmJZLmPKFRfTQUQuBjihouWRDywUhCHREkXBSNBNOiiIEgY4ouPQ/G6OHjgAEOqIo94K5sChMnkBHAAIdUdByySbHItEIQKAjCk6KZsOKRQhBoCOKhEv/MymP1yfPUQuBjijKrQMWuAhT3k8sEo1aCHRE0d1y4ScwCMMWEYK3E6Kg5ZKNcWERAhDoiIJRLtmUj9CdI3TUQKAjioRFojMpj9dnkWjUQqAjikOTc0UuZICg5YIQBDqiKAcTPfQwtFwQgkBHFM6FRZkwHzpCEOiIonyEzlwuYfIsEo0ABDqiYJHobMq7iZOiqIVARxTM5ZJNueVCDx21EOiI4lAPnUQPwYpFCBEU6GY228yeMrN1ZnZ9je3OMrOCmX2gfiWiFZXnciHQw3DpP0KkBrqZ5SXdKOliSdMlzTOz6VW2+7qke+pdJFpPgblcMjEzmbFINGoLeTvNlLTO3Te4+wFJiyXNrbDdpyT9QtK2OtaHFkXLJbucGS0X1BQS6BMkbepxf3Ppa93MbIKk90taWOuJzGyBma0wsxWdnZ1Za0ULYS6X7HJGywW1hQR6pXdc75+qb0m6zt0LtZ7I3Re5e4e7d7S1tQWWiFbEKJfsjCN0pBgUsM1mSZN63J8oaUuvbTokLS6NKR4jaY6Zdbn7r+pRJFpPeTx1jkQPljfjCB01hQT6cknTzGyqpOckXSbpQz03cPep5dtm9n1JdxPmqMVpuWSW46QoUqQGurt3mdm1Ko5eyUu6zd1Xmdk1pcdr9s2BSmi5ZMdJUaQJOUKXuy+VtLTX1yoGubtfdfhlodVxUjQ746QoUjAKGFEk9NAzy+fooaM2Ah1R0HLJLsdJUaQg0BEFLZfsGLaINAQ6oki6p8+NXMgAkjNmW0RtBDqicBaJziyfs+5JzYBKCHRE0X1hEYEeLGfWPakZUAmBjijooWfHsEWkIdARRXcPnZ/AYDkzkeeohbcToigHEz30cIxDRxoCHVEUmA89MzMWiUZtBDqiYNhidrRckIZARxTMtpgdC1wgDYGOKLrnciHPg3HpP9IQ6Iii3EPPk+jBcmYqcGERaiDQEUX53J7RcgmWy3HpP2oj0BGFu9NuyYiWC9IQ6IgiceeEaEbMtog0BDqiKCQsbpFVnlEuSEGgIwpaLtnRckEaAh1R0HLJLmdMn4vaCHREkTgXFWXFbItIQ6AjikJCyyUrJudCGgIdUbg7J0UzyjHKBSkIdERByyU7Wi5IQ6AjioRRLplxhI40BDqi4Ag9u3zOuPQfNRHoiCJJGLaYVY4FLpCCQEcUtFyy49J/pCHQEUXizLSYVc6YbRG1EeiIojhsMXYVAwuX/iMNbylEUXBXniP0THI5o4eOmgh0RMEol+xYJBppCHREkbiLPM+GRaKRhkBHFM5si5lxYRHSEOiIopA4C0RnVFwkmkRHdQQ6omDYYnYMW0QaAh1RsGJRdrRckIZARxSMcskul+OkKGoLCnQzm21mT5nZOjO7vsLjHzazx0t/HjCz0+pfKlpJIWE+9Ky4sAhpUgPdzPKSbpR0saTpkuaZ2fRem/1d0lvd/c2SviJpUb0LRWthLpfsaLkgTcgR+kxJ69x9g7sfkLRY0tyeG7j7A+7+UunuMkkT61smWo3TcsmMcehIExLoEyRt6nF/c+lr1XxU0m8Opyi0Po7QszMzJRyio4ZBAdtUettV/Kkys7epGOjnV3l8gaQFktTe3h5YIlpRwoVFmRUXuIhdBZpZyBH6ZkmTetyfKGlL743M7M2SbpE0191frPRE7r7I3TvcvaOtra0v9aJFJAktl6xyVpzUDKgmJNCXS5pmZlPNbLCkyyQt6bmBmbVLukvSP7n72vqXiVaTMH1uZoxyQZrUlou7d5nZtZLukZSXdJu7rzKza0qPL5T0JUnHSbqpdPVfl7t39F/ZGOhouWTHikVIE9JDl7svlbS019cW9rh9taSr61saWhkXFmWXz3HpP2rjQy+iYJRLdkzOhTQEOqKg5ZIdLRekIdARRZIw22JW5U80tF1QDYGOKGi5ZFf+RMNROqoh0BFF4ixwkVV5f9FHRzUEOqJglEt25d3FWHRUQ6AjChaJzq78C5A8RzUEOqJgtsXschyhIwWBjihYJDq78i9A5nNBNQQ6oqDlkl13yyWJXAiaFoGOKGi5ZEfLBWkIdETBOPTsymuwEuiohkBHFIxDz44Li5CGQEcUBS79z+xQoJPoqIxARxROyyUzeuhIQ6AjCmZbzI6WC9IQ6IiCS/+z6z4pSqKjCgIdUSQJR+hZ0XJBGgIdUTBsMTtaLkhDoCOKxA+1EBCG2RaRhkBHFFz6n92h2RYJdFRGoCOKxF15Ej2TQwtcRC4ETYtARxSMcsmOk6JIQ6AjCk6KZmdcKYoUBDoazt3lzqX/WbFiEdIQ6Gi48rA7JufKJl96t7JINKoh0NFw5ZYBeZ4NLRekIdDRcOVAouWSDRcWIQ2BjoYrH2AyyiWb8icaxqGjGgIdDVfuAef56cukPG6fHjqq4S2FhjvUQ+cIPQuj5YIUBDoarhxI9NCzoeWCNAQ6Gs4Z5dInhxaJjlwImhaBjoZjHHrfsKYo0hDoaLjyST1aLtmUf/8VCHRUQaCj4Wi59A3T5yINgY6GSxiH3ifdLRemz0UVBDoajkv/+4YVi5CGQEfDlXvoHKFnk89xUhS1BQW6mc02s6fMbJ2ZXV/hcTOzG0qPP25mZ9a/VLQKLv3vG+ZyQZrUQDezvKQbJV0sabqkeWY2vddmF0uaVvqzQNJ361wnWkh3y4XPh5mwYhHSDArYZqakde6+QZLMbLGkuZJW99hmrqQ7vXj6fZmZjTSzE9x9a70L/vPaTn317tXpG6JpHSgtiskRejblYZ5fvXuNvv27pyNXg8Nx6VmTdPUFJ9b9eUMCfYKkTT3ub5Y0K2CbCZL+X6Cb2QIVj+DV3t6etVZJ0rCjB2nauGF9+l40jzPbR+nsE4+LXcaAMvm4ofrwrHa9tPdA7FJwmMYMO7pfnjck0CsdRvX+zBeyjdx9kaRFktTR0dGnz40zJo/SjMkz+vKtwIB2VD6nr73/TbHLQBML6WJuljSpx/2Jkrb0YRsAQD8KCfTlkqaZ2VQzGyzpMklLem2zRNIVpdEuZ0va1R/9cwBAdaktF3fvMrNrJd0jKS/pNndfZWbXlB5fKGmppDmS1knaK+kj/VcyAKCSkB663H2piqHd82sLe9x2SZ+sb2kAgCwYCQwALYJAB4AWQaADQIsg0AGgRVisyfLNrFPSs3389jGSttexnHpp1rqk5q2NurKhrmxasa7J7t5W6YFogX44zGyFu3fErqO3Zq1Lat7aqCsb6srmSKuLlgsAtAgCHQBaxEAN9EWxC6iiWeuSmrc26sqGurI5ouoakD10AMBrDdQjdABALwQ6ALSIARHoZvYNM/tbaQHqX5rZyCrb1VzMuh/q+qCZrTKzxMyqDkEys2fM7Akze9TMVjRRXY3eX6PN7D4ze7r096gq2zVkfzXr4ucBdV1oZrtK++dRM/tSg+q6zcy2mdmTVR6Ptb/S6oq1vyaZ2R/NbE3p/fiZCtvUd5+5e9P/kfQOSYNKt78u6esVtslLWi/pREmDJT0maXo/13WKpJMk/UlSR43tnpE0poH7K7WuSPvr3yVdX7p9faX/x0btr5B/v4pTQv9GxRW5zpb0UAP+70LqulDS3Y36eerxum+RdKakJ6s83vD9FVhXrP11gqQzS7eHS1rb3z9jA+II3d3vdfeu0t1lKq6I1Fv3YtbufkBSeTHr/qxrjbs/1Z+v0ReBdTV8f5We/47S7Tskva+fX6+WkH9/9+Ln7r5M0kgzO6EJ6orC3e+XtKPGJjH2V0hdUbj7Vnd/pHR7t6Q1Kq613FNd99mACPRe5qv4G623agtVNwOXdK+ZrSwtlN0MYuyvcV5ayar099gq2zVif4X8+2Pso9DXPMfMHjOz35jZqf1cU6hmfg9G3V9mNkXSGZIe6vVQXfdZ0AIXjWBmv5N0fIWHvuDuvy5t8wVJXZJ+VOkpKnztsMdkhtQV4Dx332JmYyXdZ2Z/Kx1VxKyr4fsrw9PUfX9VULfFz+ss5DUfUXE+jz1mNkfSryRN6+e6QsTYXyGi7i8zGybpF5I+6+4v9364wrf0eZ81TaC7+z/UetzMrpT0bklv91LzqZd+Wag6ra7A59hS+nubmf1SxY/VhxVQdair4fvLzF4wsxPcfWvpY+W2Ks9R9/1VQbMufp76mj1Dwd2XmtlNZjbG3WNPQtWUi8XH3F9mdpSKYf4jd7+rwiZ13WcDouViZrMlXSfpve6+t8pmIYtZN5yZvc7Mhpdvq3iCt+LZ+AaLsb+WSLqydPtKSa/5JNHA/dWsi5+n1mVmx5uZlW7PVPF9/GI/1xWiKReLj7W/Sq95q6Q17v7NKpvVd581+sxvX/6ouPj0JkmPlv4sLH19vKSlPbabo+KZ5PUqth76u673q/gbdr+kFyTd07suFUcrPFb6s6pZ6oq0v46T9HtJT5f+Hh1zf1X690u6RtI1pdsm6cbS40+oxkimBtd1bWnfPKbiIIFzG1TXTyRtlXSw9PP10SbZX2l1xdpf56vYPnm8R3bN6c99xqX/ANAiBkTLBQCQjkAHgBZBoANAiyDQAaBFEOgA0CIIdABoEQQ6ALSI/wNLeFCqFWTwsgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "normal_trunc2 = chaospy.Trunc(normal, lower=-0.5, upper=0.5)\n",
    "q_loc = np.linspace(-2, 2, 200)\n",
    "pyplot.plot(q_loc, normal_trunc2.pdf(q_loc))\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "id": "bHndPsbHH8ws"
   },
   "outputs": [],
   "source": [
    "normal = chaospy.Normal(mu=0, sigma=0.43)\n",
    "dis1 = chaospy.Trunc(normal, lower=-0.5, upper=0.5)\n",
    "dis2 = chaospy.Trunc(normal, lower=-0.5, upper=0.5)\n",
    "\n",
    "joint = chaospy.J(dis1, dis2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "n7es7M7CE1Ou",
    "outputId": "7fe210ad-46d8-4c02-e767-ee8144f994ea"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.17583247],\n",
       "       [ 0.32460843]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dis1.sample(2).reshape(-1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "gLb9cElRJkIw",
    "outputId": "4e4567a9-8d43-4f52-81d0-dba386aca87a"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.006944444444444444"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/12**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "0VK0hyvgHwVb"
   },
   "outputs": [],
   "source": [
    "cov_w = chaospy.Cov(joint)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "id": "FdStIcaZEjTs"
   },
   "outputs": [],
   "source": [
    "# joint.sample(size=1,rule = 'random')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "id": "6zghYPJaKcWN"
   },
   "outputs": [],
   "source": [
    "def dlqr(A,B,Q,R):\n",
    "    \"\"\"\n",
    "    Solve the discrete time lqr controller.\n",
    "    x[k+1] = A x[k] + B u[k]\n",
    "    cost = sum x[k].T*Q*x[k] + u[k].T*R*u[k]\n",
    "    \"\"\"\n",
    "    # first, solve the ricatti equation\n",
    "    P = np.matrix(scipy.linalg.solve_discrete_are(A, B, Q, R))\n",
    "    # compute the LQR gain\n",
    "    K = np.matrix(scipy.linalg.inv(B.T*P*B+R)*(B.T*P*A))\n",
    "    return -K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "id": "UnVT3m8P57Fv"
   },
   "outputs": [],
   "source": [
    "# K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "id": "M94PGvxjKXdO"
   },
   "outputs": [],
   "source": [
    "N_pred = 6\n",
    "nu = 1\n",
    "nx = 2\n",
    "\n",
    "A = np.array([[1.6,1.1],[-0.7,1.2]])\n",
    "B = np.array([[1],[1]])\n",
    "D = np.eye(2)\n",
    "E_ww = 1/12**2 * np.eye(2)\n",
    "Q = np.eye(2)\n",
    "R = np.eye(1)\n",
    "\n",
    "p = 0.8\n",
    "\n",
    "F = np.array([[1, 0.136]])\n",
    "G = np.zeros([1, nu])\n",
    "# E = np.hstack([np.eye(nu), np.zeros(nu, N_pred - 2)])\n",
    "K = np.array(dlqr(A,B,Q,R))\n",
    "\n",
    "Phi = A + B @ K\n",
    "\n",
    "F_tilde = F + G @ K\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "E_temp1 = np.eye(nu)\n",
    "E_temp2 = np.zeros([nu, (N_pred-2) * nu])\n",
    "E = np.hstack([E_temp1,E_temp2])\n",
    "\n",
    "M_temp1 = np.eye(nu * (N_pred - 2))\n",
    "M_temp2 = np.vstack([M_temp1, np.zeros([nu, M_temp1.shape[1]])])\n",
    "\n",
    "M = np.hstack([np.zeros([M_temp2.shape[0],nu]),M_temp2])\n",
    "\n",
    "F_bar = np.hstack([F_tilde, G@E])\n",
    "Psi = np.vstack([np.hstack([Phi, B @ E]), np.hstack([np.zeros([np.shape(M)[0],np.shape(Phi)[1]]), M])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "JV_7SEz8CFVr",
    "outputId": "f0eee5a8-85c7-4704-ebcd-4a1e16d18493"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 5)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RN0oGLiZxPjW"
   },
   "source": [
    "System"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "id": "ku_eM8hTKanW"
   },
   "outputs": [],
   "source": [
    "A_DM = ca.DM(A)\n",
    "B_DM = ca.DM(B)\n",
    "D_DM = ca.DM(D)\n",
    "E_ww_DM = ca.DM(E_ww)\n",
    "Q_DM = ca.DM(Q)\n",
    "R_DM = ca.DM(R)\n",
    "\n",
    "Phi_DM = ca.DM(Phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NIjuH-6hxRtk"
   },
   "source": [
    "Offline-calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BWK5fULJxobf"
   },
   "source": [
    "Beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "YzjHvtmpy4nl"
   },
   "source": [
    "gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "id": "RKSjrjRXKhZM"
   },
   "outputs": [],
   "source": [
    "N_gamma = 10\n",
    "\n",
    "gamma_list = []\n",
    "\n",
    "kappa = np.sqrt(p/(1-p))\n",
    "\n",
    "g = F.T\n",
    "\n",
    "\n",
    "P1 = D @ cov_w @ D.T\n",
    "Pi = P1\n",
    "\n",
    "for i in range(N_gamma):\n",
    "    gamma_i = kappa * np.sqrt(g.T @ Pi @ g)\n",
    "    gamma_list += [gamma_i]\n",
    "    \n",
    "    Pi_p1 = Phi @ Pi @ Phi.T + D @ cov_w @ D.T\n",
    "    Pi = Pi_p1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "BCOistE9MwGU",
    "outputId": "c650f385-865d-480d-b49a-2441b97cce90"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([[0.53152455]]),\n",
       " array([[0.55929024]]),\n",
       " array([[0.56016487]]),\n",
       " array([[0.56027406]]),\n",
       " array([[0.56040883]]),\n",
       " array([[0.5604517]]),\n",
       " array([[0.56045857]]),\n",
       " array([[0.56045908]]),\n",
       " array([[0.56045908]]),\n",
       " array([[0.56045909]])]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gamma_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "whmKXa85zjvQ",
    "outputId": "8eb98610-65e1-4b2c-c316-84b07bb6b5ef"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "******************************************************************************\n",
      "This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
      " Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
      "         For more information visit https://github.com/coin-or/Ipopt\n",
      "******************************************************************************\n",
      "\n"
     ]
    }
   ],
   "source": [
    "solver_opt = {}\n",
    "solver_opt['print_time'] = False\n",
    "solver_opt['ipopt'] = {\n",
    "    'max_iter': 300,\n",
    "    'print_level': 1,\n",
    "    'acceptable_tol': 1e-6,\n",
    "    'acceptable_obj_change_tol': 1e-6}\n",
    "\n",
    "i_max = 10\n",
    "\n",
    "w_sym = ca.SX.sym(\"w\",2)\n",
    "a_list = []\n",
    "\n",
    "nl = {}\n",
    "\n",
    "lbx = [-0.5] * 2\n",
    "ubx = [0.5] * 2\n",
    "\n",
    "nl['lbx'] = lbx\n",
    "nl['ubx'] = ubx\n",
    "\n",
    "for i in range(1,i_max):\n",
    "    obj = - F_tilde @ np.linalg.matrix_power(Phi, i) @ D @ w_sym\n",
    "\n",
    "    nlp_prob = {\n",
    "    'f': obj,\n",
    "    'x': w_sym\n",
    "    }\n",
    "\n",
    "    solver = ca.nlpsol(\"solver\", \"ipopt\", nlp_prob, solver_opt)\n",
    "\n",
    "    sol = solver(**nl)\n",
    "\n",
    "\n",
    "    a_list += [-sol['f'].full()]\n",
    "\n",
    "\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "psNzKsdy4CIb",
    "outputId": "68f3aa11-af76-4864-8201-f9acb8c7f193"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([[0.19778557]]),\n",
       " array([[0.04023669]]),\n",
       " array([[0.01341676]]),\n",
       " array([[0.01168454]]),\n",
       " array([[0.00716333]]),\n",
       " array([[0.00303603]]),\n",
       " array([[0.00089288]]),\n",
       " array([[9.8153436e-05]]),\n",
       " array([[0.00011195]])]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "id": "_j8RZrj34U4g"
   },
   "outputs": [],
   "source": [
    "beta_list = []\n",
    "gamma_1 = gamma_list[0]\n",
    "beta_temp = gamma_1\n",
    "for i in range(1,i_max):\n",
    "    beta_temp = beta_temp + a_list[i-1]\n",
    "    beta_list += [beta_temp]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "VGbP3cdJ8NN1",
    "outputId": "99b4a565-d814-4c26-c3ba-041ea4ba0e6a"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([[0.72931012]]),\n",
       " array([[0.7695468]]),\n",
       " array([[0.78296356]]),\n",
       " array([[0.7946481]]),\n",
       " array([[0.80181143]]),\n",
       " array([[0.80484746]]),\n",
       " array([[0.80574034]]),\n",
       " array([[0.80583849]]),\n",
       " array([[0.80595044]])]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "olc9kz6k8zVj"
   },
   "source": [
    "Objective function\n",
    "\n",
    "$$\n",
    "J\\left(x_{k}, \\mathbf{c}_{k}\\right)=\\sum_{i=0}^{\\infty} \\mathbb{E}_{k}\\left(\\left\\|x_{i \\mid k}\\right\\|_{Q}^{2}+\\left\\|u_{i \\mid k}\\right\\|_{R}^{2}-l_{s s}\\right)=x_{k}^{T} W_{x} x_{k}+\\mathbf{c}_{k}^{T} W_{c} \\mathbf{c}_{k}+w_{1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sdWdQi2AWiGr"
   },
   "source": [
    "$$\n",
    "W_{x}-\\mathbb{E}\\left(\\left(A_{k}+B_{k} K\\right)^{T} W_{x}\\left(A_{k}+B_{k} K\\right)\\right)=Q+K^{T} R K\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eizsnoLwWn3E"
   },
   "source": [
    "$$\n",
    "A^{H} X A-X-\\left(A^{H} X B\\right)\\left(R+B^{H} X B\\right)^{-1}\\left(B^{H} X A\\right)+Q=0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "id": "NAtJgvf2MxSn"
   },
   "outputs": [],
   "source": [
    "Wx = scipy.linalg.solve_discrete_are(A+B@K, np.zeros([2,2]), Q + K.T @ R @ K, np.diag([0.01,0.01]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "L3KIJfuKX1g6",
    "outputId": "420fa185-07f5-43fb-9c87-81c130e7527c"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10.39419575,  0.80911714],\n",
       "       [ 0.80911714,  2.18978722]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Wx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "id": "GhqJHFYvX8Bi"
   },
   "outputs": [],
   "source": [
    "nc = N_pred - 1\n",
    "\n",
    "Wc = R + B.T @ Wx @ B\n",
    "for i in range(nc-1):\n",
    "    Wc = scipy.linalg.block_diag(Wc, R + B.T @ Wx @ B)\n",
    "\n",
    "\n",
    "\n",
    "# temp = scipy.linalg.block_diag(Wx,Wx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "1Ry9gDwnZU1I",
    "outputId": "75933b0e-be51-4f25-d8ac-7cad778f5aa4"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[15.20221724,  0.        ,  0.        ,  0.        ,  0.        ],\n",
       "       [ 0.        , 15.20221724,  0.        ,  0.        ,  0.        ],\n",
       "       [ 0.        ,  0.        , 15.20221724,  0.        ,  0.        ],\n",
       "       [ 0.        ,  0.        ,  0.        , 15.20221724,  0.        ],\n",
       "       [ 0.        ,  0.        ,  0.        ,  0.        , 15.20221724]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Wc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "u6sSjatQ-i9D",
    "outputId": "1f34ac27-a7d6-450d-c8e1-77a28cecc096"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SX([a_0, a_1, b_0, b_1])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = ca.SX.sym(\"a\",2)\n",
    "b = ca.SX.sym(\"b\",2)\n",
    "ca.vertcat(a,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "id": "2mEOOjsGZ_iB"
   },
   "outputs": [],
   "source": [
    "# Not necesssary to solve the Lyapunov eq."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "3sKLbiwzAys9",
    "outputId": "7fc5d3da-07ae-4240-84bb-48bfa3aed598"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 7)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(F_bar @ np.linalg.matrix_power(Psi, i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "id": "ZJemTjZlBJu2"
   },
   "outputs": [],
   "source": [
    "# ca.vertcat( x0_para, ck_var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Tc-Mi3wzBezn",
    "outputId": "87b02a19-0f13-4829-8487-60e1f69d2b51"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7, 7)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Psi.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1kgo2m4YBev3"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "id": "qgFQRGa6xnVs"
   },
   "outputs": [],
   "source": [
    "x0 = np.array([[-5],[60]])\n",
    "\n",
    "\n",
    "v = 5 # i_max\n",
    "\n",
    "\n",
    "# u_ub = [1]\n",
    "# u_nom_ub = [u_ub-Z_proj_max]\n",
    "\n",
    "# u_lb = [-1]\n",
    "# u_nom_lb = [-1-Z_proj_min]\n",
    "\n",
    "# x_ub = [ca.inf, 2]\n",
    "# x_lb = [-ca.inf, -ca.inf]\n",
    "\n",
    "# x_nom_ub = [ca.inf, 2 - Z.V[:,1].max()]\n",
    "# x_nom_lb = [-ca.inf, -ca.inf]\n",
    "\n",
    "# xi_var = ca.SX.sym('x', nx, N_pred)\n",
    "x0_para = ca.SX.sym('x0', nx, 1)\n",
    "ck_var = ca.SX.sym('ck', N_pred - 1,nu)\n",
    "\n",
    "obj = 0\n",
    "g = []\n",
    "for i in range(v):\n",
    "    g += [F_bar @ np.linalg.matrix_power(Psi, i) @ ca.vertcat( x0_para, ck_var) ]\n",
    "\n",
    "# obj += x0_para.T @ Wx @ x0_para + ck_var.T @ Wc @ ck_var\n",
    "obj += ck_var.T @ Wc @ ck_var\n",
    "# g += [xi_var[:, -1] - (A_para @ xi_var[:, -2] + B_para @ ui_var[:,  -1])] # Terminal state\n",
    "\n",
    "\n",
    "dec_list = [ca.reshape(ck_var, -1, 1)]\n",
    "para_list = [x0_para]\n",
    "\n",
    "nlp_prob = {\n",
    "    'f': obj,\n",
    "    'x': ca.vertcat(*dec_list),\n",
    "    'g': ca.vertcat(*g),\n",
    "    'p': ca.vertcat(*para_list)\n",
    "}\n",
    "\n",
    "solver_opt = {}\n",
    "solver_opt['print_time'] = False\n",
    "solver_opt['ipopt'] = {\n",
    "    'max_iter': 300,\n",
    "    'print_level': 1,\n",
    "    'acceptable_tol': 1e-6,\n",
    "    'acceptable_obj_change_tol': 1e-6}\n",
    "\n",
    "solver = ca.nlpsol(\"solver\", \"ipopt\", nlp_prob, solver_opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6n6iLbWNACos"
   },
   "outputs": [],
   "source": [
    "N_sim = 20\n",
    "\n",
    "nl = {}\n",
    "\n",
    "\n",
    "\n",
    "lbg = []\n",
    "ubg = []\n",
    "\n",
    "one_vec = np.ones([np.shape(F_bar)[0],1])\n",
    "for i in range(v):\n",
    "    ubg += [one_vec - beta_list[i]]\n",
    "\n",
    "lbg = [-ca.inf] * v\n",
    "\n",
    "nl['lbg'] = lbg\n",
    "nl['ubg'] = ubg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "id": "6n6iLbWNACos"
   },
   "outputs": [],
   "source": [
    "N_loop = 30\n",
    "\n",
    "x_diff_traj_list = []\n",
    "for i in range(N_loop):\n",
    "\n",
    "    x0 = np.array([[-5],[60]])\n",
    "    x_traj_list = []\n",
    "    x_traj_list += [x0]\n",
    "    xk = x0\n",
    "\n",
    "    for i in range(N_sim):  \n",
    "        nl['p'] = xk\n",
    "        sol = solver(**nl)\n",
    "        ck = sol['x'][0].full()\n",
    "        # x0_sol = sol['x'][:nx].full()\n",
    "        uk = ck + K @ xk\n",
    "        wk = dis1.sample(2).reshape(-1,1)\n",
    "        xk_next = A @ xk + B @ uk + D @ wk\n",
    "        x_traj_list += [xk_next]\n",
    "\n",
    "        xk = xk_next\n",
    "        \n",
    "        \n",
    "    x_diff_traj_list += [x_traj_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABwJElEQVR4nO3ddXhUZ/r/8fczmpm4u7tAAkGCJ1hxqQuUKqUthSqlTnVpu9Vtty5UKXWjuLtLkBhxd5nI2Pn9Abvf7m+3SugQeF7XlWsyk5Nz7vNcw4eTe855jlAUBUmSJOncp3J0AZIkSdJfQwa+JEnSeUIGviRJ0nlCBr4kSdJ5Qga+JEnSeUIGviRJ0nmiWwJfCOEhhPhCCHFcCHFMCDFICOElhFgthMg79ejZHduSJEmS/pzuOsJ/CVihKEoCkAocAxYCaxVFiQXWnnouSZIkOYg43QuvhBBuwEEgSvnZyoQQOUCmoiiVQohAYIOiKPGntTFJkiTpT9N0wzqigFrgPSFEKrAXmA/4K4pSCXAq9P3+1y8LIWYDswFUBrd0Jy8f/DQ2dHp9N5T217Lb7ahUPfdjEVm/Y/Xk+nty7dDz68/Nza1TFMX3t5brjiP8fsAOYIiiKDuFEC8BLcBtiqJ4/Gy5RkVRfrWP7xIepXhf8Q9cjV3c69fCzDnXnFZtf7UNGzaQmZnp6DL+NFm/Y/Xk+nty7dDz6xdC7FUUpd9vLdcd/6WVAWWKouw89fwLoC9QfaqVw6nHmt9aUYBOTVRiE61mJx4r9eLxh5/rhvIkSZIk6IbAVxSlCigVQvyrPz8KOAp8B8w69dos4NvfWpdapWLpxLHE9mqky6DnHXM8t977D+x2++mWKUmSdN7rrqbVbcDHQohDQBrwFLAYGCOEyAPGnHr+m/x9fHk/cyRxSfXY/A38KKK47L4PqW1o6aZSJUmSzk/dEviKohxQFKWfoii9FUWZpihKo6Io9YqijFIUJfbUY8PvXV+IfyBvDhxKbEQdllg3dgtvLnpmFRt3H+2OciVJks5LZ+3H0lHBYbySmk6sdzXmvj5UqrXc+mU+z7//naNLkyRJ6pHO2sAHSIqK4fn4RKL1lbQNCkRt6OCV4ypmL3pH9vUlSZL+oLM68AH6JibzdEQ4UbZqagZFEeJVw6rOAKbe+z5NLW2OLk+SJKnHOOsDH2Bwal8eC/AhvKOW/PTe9AnLIVvtx6THf2B3doGjy5MkSeoRekTgA4waOIgHPZwIaatnV8JwRiftpFHtwrUfHuaNL9c7ujxJkqSzXo8JfIBJI0Zyt14hoK2R5SFTmZi+Ci+h4eldJuY884ns60uSJP2KHhX4AJeOG8d8pQM/Uwufes9gdMZXZKjsrGhwZ8J9H9As+/qSJEn/U48LfICrp07hZlM9XqZ23nKbQ9qQpVyl7iRH+DDmse84cLzQ0SVKkiSddXpk4APMvvwSrm+owMPUySvGeXgP/ox7dU20a9y58t0DfPDDVkeXKEmSdFbpsYEPMG/WlcyqLMKt3crLhjvoyviKB/XF+Kn0LNrcyPwXPpd9fUmSpFN6dOAD3D37aq4qzMfYofCi093UZ6xkrnE/Q4WWb6uNTHzgQ1paTY4uU5IkyeF6fOAD3HfbNVyRcxynDsFz+nuo7b+NMa6ruQE9x+3ejHr0G7Jzix1dpiRJkkOdE4EP8MDt13DpkaNoOrU8p19IR3o2Xp6f8bRwol3tziVv7eHT5dscXaYkSZLDnDOBr1KpeOTua7n0wBFEl55ndfehSyuiyOtd3hRq/FROPLCxnrtf/kr29SVJOi+dM4EPJ0N/0b1Xc8nuoyhdRp7RPoBvahVrvd7iddHOEKHhiwo9Ux/+iFZTu6PLlSRJ+kudU4EPoFFrePi+y7l41zGsZmee1jxIRFo9n3p/wMOqEq5Hz2GrNyMXfc2x/BJHlytJkvSXOecCH8BJ68QDC6Zz6bZcuqxuLFY/TFxqIx96fMM03S6exoBJ5cb0N3bz+codji5XkiTpL3FOBj6Aq8GVBfeM59Kt+bRbPVisfoTEPk186bqZAOcfeFu44qvSsmBdLfe9+rXs60uSdM47ZwMfwMvViztuz+LSLSdotXnxlGoRCWkNbHE6RonXR7wjnBgEfFqqY/ojH8u+viRJ57RzOvABArz8uXXuQC7dUkyzzZfFYhFxfeooUpXzVeASnhF2rlPUHLR4MWrR1+QUlDq6ZEmSpDPinA98gHD/MG6c3YtLt5RSr/jzlLKI2LQa1JZGnoz4hJtUVfwNJ1pULkx7fSdfr9nl6JIlSZK63XkR+ABxIbHMnBXFpVsqqCWQvymLiOlTQaypjdmxPzBKdYB3hDseQs2dq6t58LVvZF9fkqRzSrcEvhCiSAhxWAhxQAix59RrXkKI1UKIvFOPnt2xrdORGt2LSy7155KtVVSKYJ6yPUp0Whmjmtu5KHYrMZrlvKfypp9i5aNiLRc/+gkmU4ejy5YkSeoW3XmEn6UoSpqiKP1OPV8IrFUUJRZYe+q5w2Uk92fKZBcu2VJLhSqMp6yLiOhdymUtnUyOycGo/4DnVO7MtNvZ1+XJyEe/IvdEmaPLliRJOm1nsqUzFVhy6vslwLQzuK0/JCt9OGPHCC7aVkepJoLF1kVE9Crm2lYL0yKqaHJ5mzkqDU8qKhpxZuprO/hurezrS5LUswlFUU5/JUIUAo2AAryhKMqbQogmRVE8frZMo6Io/9XWEULMBmYD+Pr6pi9btuy06/m99hdlk1+WyNeDvIgy53Of/hHKjkXwhpOdxQ1e9G6YSp7dhzvsjTSodExwq+PCQRGoVP/7/8m2tjZcXFz+svq7m6zfsXpy/T25duj59WdlZe39WXflF3VX4AcpilIhhPADVgO3Ad/9nsD/ufj4eCUnJ+e06/kjPvphGdsO+vB1hiexXbksdFpERU4Eb6jN3GCOYnp5HxqUFO62VpKtcaGfUxNL7rkQZ2fDf61rw4YNZGZm/qX1dydZv2P15Pp7cu3Q8+sXQvyuwO+Wlo6iKBWnHmuAr4EBQLUQIvBUMYFATXdsq7vNmHQpfePLmbKzmTynOJ7teIjg+EJuVrS8p8nh73EF+KrW8k9NIBfaOtjT6cGoR78ir6jc0aVLkiT9Iacd+EIIZyGE67++B8YC2cB3wKxTi80Cvj3dbZ0pN1w8k95hJ5i0q4VjxiT+bnqQoNhCblIb2Na+g1tSBG7qj7hT7c8DNhv1ODPl1e38sH6Po0uXJEn63brjCN8f2CKEOAjsAn5UFGUFsBgYI4TIA8acen7WunXGdST5HGfi7haOuKTwfOv9BEUXMtvJSFX9SqalRWPUvMB4tZHXccJJsTFvRQWPv/W9PF9fkqQe4bQDX1GUE4qipJ76SlYU5clTr9crijJKUZTYU48Np1/umXX79bOJdz7MuL2tHHJL5cWW+wiILOIGF2fUFV8xMn0wGu1zJKo6+UjtQ5S1mXcKVFz++KfyfH1Jks56582Vtr+HEIK7b76ZePUBxuxrY797H15qupeA8GKu83TGt+hzhvcbicnwNj6qEt7RhDLSUs+uDg9GP/oV1fWtjt4FSZKkXyQD//+jUqlYMP8WEiz7GXnAxF7Pfvyj/h78Qkq51tdITN5XXJCaQZ77WpzVu3hMG8nNljbqMLJoRxc/bpB9fUmSzk4y8P8HtUrNPXffRJJpHyMOt7PLZwCv1t6Jb1ApVwfqSMlbzqz4MFb6V+Gi/oKrtEE8Y9egVSzc9lMFT77zg+zrS5J01pGB/wt0Gh13LbiBlPr9DD3SwQ6/QbxWdQc+gWVcFaIhNXcTj4TYeCUyFHf1P+ivceFdlRfBljreyhNc+cRSTB1djt4NSZKkf5OB/ysMOgPzF84kteoAg451sC1wCG+Uz8fbv5wrIiA95wAfuB3n7sQJuIvHCVTZeU8bTrqllB3t7oxd9AX5xRWO3g1JkiRABv5vcjO4ccvCS0kvPczAnE42hwznrdJ5ePpWckmshX65+WzUrOaK9Ntw4RGcVbW8oEnkQlsF1YozU17Zxk8b9zl6NyRJkmTg/x7eLt5ct2ASA08cpV9eJxvDRvBu0a24e1VzUUIXg3IrKTR/yNhBj6GxPY9Bnc2d6gTuVDpQ2y3curyMv737o+zrS5LkUDLwf6cgzyCuuDuLIbk59C3oZF3kSJYUzMHds5bJKSaG5LbS3vYyw0Y8RZv9K4zqVUxRh/CsyhU/cy1v5MKMJ2VfX5Ikx5GB/wdE+kUy/faBjDiaT2phF2tix/BB3o24udUzPrWZ4blWdPVPkTliAYXiBO6ad0lRu/KaJpR4cwHbTCf7+ieKKx29K5IknYdk4P9BCSEJXHBbCqMOnSClqItV8RfwUc71uLg2MqZPI1m5ejxqHmfasCvYZPDDU/03fNUqXtWkMNSeQ7XizKRXtrJi035H74okSecZGfh/QlpkGkNviuCCg8UklXSxImk8nx69FqNzEyMH1DAm1x3vyse5KX0gn/hMwJuFGFWtPCn6cqm6DLVi4ZYfS3nm/R8dvSuSJJ1HZOD/SYMTB9H3Gh/G7ysjvrSL5b0m8Fn2LAyGVoZnVDI+zw+fyid4LMGHpyPvwsO6AL2qkJvtvblNY8PHXMs/j8OMJz6VfX1Jkv4SMvBPw6i0kSRc4cTkPRXElnfxY9oklh2aiV5vYvCgUiblheJT9QzvhrRze+rLOJuewKjexgTCeUzrR1hXAVva3Bj3qOzrS5J05snAP02TMyYSdqGdqbuqiK7o4se0iSw7MBOdrpOMIUVMy43Eq/olVnrkMmPo+yhNH+Cq/pwUlQfPaRJJUo5QaXdm8itbWLlF9vUlSTpzZOB3g0szL8J3fBvTd9QSUWNhefoEPt93FRqNmX7D8rk4LxavmrfI1m3lwrFLaKnfh6fmRQLUGp5T+jNYfwyh2Lj5+1KeXbLc0bsjSdI5SgZ+N5k1/io6e2Vz8fZ6wmotLO83gS/2XoFKZaPv8FwuKYjDve4Tqm3fMnXSi5Q0avDRPIirqpNFXYOY4lKBt7meV48pXP3kp7TLvr4kSd1MBn43GhqXgXZwMZdsbSS43spPAyby5e7LQdjpM/wYlxXG49LwPV0dH3LhpHs42J6Br7gLJ1U1czrSuc6oENxZxKZWNy549AsKS6scvUuSJJ1DZOB3s3mX34ySnsNlW5oIaLSyImMCX+26DLsiSB12hCuK4jE2bkDX/BpXjructerr8bLejZMqm8m2WBY4BRJmP06l3ZlJ/9jMKtnXlySpm8jAPwPunjUPc0o2l29uxq/ZxspB4/lq16VY7BpShh/mitJ4nJr34lH3HLMzh/KF91M4tz6OUbWSdHx5Skkhxukowm5jzvel/H3JT47eJUmSzgEy8M8AIQT33ng7HTEHuXxzC94tdlZnjOPrnZdgtupIHnaIyyvi0TflEFD1JPcMiOaNuLdR17yHu+YdwlUGnunIoI93Ht7mel45ZmfWU5/S0Sn7+pIk/Xky8M8QtUrN/bfeQXvofq7c2Iqnyc7aQWP5ZufFdJgNJA05wGW1sTg1lRFa/gjPJrmxeOCHdBVvw0f3OJ5qOw83DiPLt5qgzlI2trhxwaIvKJJ9fUmS/iQZ+GeQRq3hvvnzaQ/czxUbTbh1wrrBY/hmx8WYulxIGnyASxtjMDY2ElrxAEsi7Nw39j0a8lrw09yNUdXMLc2DuMRTIdBygopTff3VWw44etckSeqBui3whRBqIcR+IcQPp557CSFWCyHyTj16dte2ehK9Vs89d95Ch89+rtzQhkuXYNPgkXy740JaOtxJyDjARW1RuNWZCal4gOX+Ldw+/SXKC8LxVd2OXlXExe3JzNP5EqA7AXY7N31fwnMfrHD0rkmS1MN05xH+fODYz54vBNYqihILrD31/LzkrHfmjrtuoMv9EFduaMfJqmLL4Ey+2zmdpnZPEgbuZ3pXGF61aoIqH2CHeyVzLnmQgorxeNsW4KTewTAlmEVd8YT6nsDb3MA/jtq4Rvb1JUn6A7ol8IUQIcBE4O2fvTwVWHLq+yXAtO7YVk/lYfTg1rtmYHE+xFUbOtDZNWzPGMb3O6dR1+pLXP/9TLYHEVDlgl/1InKd8rjuous53DYfl7ZncVEvI0HtxuN1A4iPKCaoo4wNLW6Me/RLisqqHb17kiT1AEJRlNNfiRBfAH8DXIG7FUWZJIRoUhTF42fLNCqK8l9tHSHEbGA2gK+vb/qyZctOux5HaWtrw8XF5VeXaehsIHtDBSpLCh9mGkDpov+ObYzr9wMB7lXk7U9jbVcDBSF1mDznYdD05R/r9tOn+X5UAYNpsMynzSZ4zWc9uyoDqNMF4qR0cXOiQnKU/xmv/2wm63ecnlw79Pz6s7Ky9iqK0u+3ljvtwBdCTAImKIpyixAikz8Y+D8XEROvFOXnnFY9jrRhwwYyMzN/c7n86nw+e34NNlsCH2UaUNs66LNzF+PSfyTIo5z8g6lsbGvhaHgVNrfrMTuP4PmtRQzKn4dLdCg15kew2A0sddvPtyY9Josn7WoDc1N03Dlz3Bmv/2wl63ecnlw79Pz6hRC/K/C7o6UzBJgihCgClgIjhRAfAdVCiMBTxQQCNb+1okqTnRlv72RbQR3d8ZfH2SrGP4Zptw1DI/K5amMHVq2BAwP6s3zvJEoawohJPUCmuyupJ0IRLe/g0vI9tw6LYEXvt2g8UkeAbj56dQ0zWvswW2XA07sWb3MDLx+xce3fltLZJfv6kiT9t9MOfEVR7lMUJURRlAjgcmCdoigzgO+AWacWmwV8+1vr8tQLjle1cuVbO7notW2sPVZ9zgZ/ckgyY+akorcXceXGTsw6I9n9+rJy/wQK66OI6XWQYT46BuZGYDF9TmD9xyxI9+PL4a9RdcANf9V89OqjXGCLYkFjKIExNQR1lLO+2ZVxi76guFz29SVJ+k9n8jz8xcAYIUQeMObU81/lrhdsuTeLx6cmU93SxfVL9jD+pc18d7ACm/3cC/7+0f0ZfEMUzpYSrtjUSbvBlWN9U1l9YBy5NXFEJR8mI1jF0GPRmDpXElXzTx7v5cn7E5+l7FAaXrYHMahX0Uflxb1FfYhPqca/s4IymwsTX9rE2u0HHb2LkiSdRbo18BVF2aAoyqRT39crijJKUZTYU48Nv2cdTlo1MwdFsOGeTJ67JBWLzc68T/cz6rkNfLa7BLPV3p0lO9zwxOGkXR2AW2c5V27qoNXFndw+vVh7aCxHqxOJTMimf4SNrOx4miw7Sax4mldjXXjxovspyZ+G0fRP3NVvE6bSc0fuMDJiyvFQGsBuZ/Y3Rbz4kTxfX5Kkk87aK221ahUXpYew+o4RvD6jLy5OGu798jAjnl3Pu1sKaTdbHV1itxmbNoa4K13wbK/iik0dNLl5cCI1kY3ZozlU2YvwuKP0ievkgoMJ1FmPklz2EEvDtDx60Q2U1s1F1P+Ej+5xPNR25paMZKJPDc5uHXiZG3kx28Z1i5fS2WV29G5KkuRgZ23g/4tKJRiXEsj3c4ey5LoBhHoZeeyHowx9ej2vrs+nucPi6BK7xdQBUwi6UMGnrZbLNndQ7+FNUa84thwbyb7yNMKij9Er2cTE/cnU2EtJKV3Aan+Fu6dOpFh5nM6iIwTo7sRJ3cK1DYOYYTXjFW0lqKOCdU2ujF/0OSWyry9J57WzPvD/RQjBiDhflt00iM/nDKJ3iDvPrsxh6OJ1PLPiOHVtPf/MlCtGXI7HxDYCmhu4dIuJGi9fypKi2JYzkp2l/QmJzCGpTxNT9qZQY28gqewO9rl3Mnd8OkU+L9B4pIkg3Vx06lKmdSVwXakbEX06COispORUX3/99kOO3k1JkhykxwT+z/WP8OL9awfww21DGR7ny2sbCxj69DoWfXeEiqYOR5d3Wq674Bp0o6sJbmzmki0mKnwCqEoIZVf+cLYWDyI4LI+EfvVM39OLBmsXMZV3UmRoZPaoWE4kvErFbiP+qtvRa/YwRARyw5FYUvvV4GZrQbEr3PBNIS99vMrRuylJkgP0yMD/l5Rgd169qi9r7hzBpN5BfLSjmBHPrmfBFwc5Udvm6PL+tFunzsE+rIiwhhYu3tZGqX8QtbFB7CsazIbCYQSFFhCXUcO0vcm0WRRCqu6hRZRzw7BAcjJepGR3FJ62J3DWfEmc2oXZBzIYlpiHzknBy9LEC4ctXC/7+pJ03unRgf8v0b4u/P2SVDbck8mVA8L49kAFo5/fyNxP9nG0osXR5f0pd1w6j/YBOUTUtnHRjlaKg0JpiAzgUGkGq09kERB8grghFUzdm4i5S49X3YNozLncNMiXA1lPUnxgMMb2D/DUvoCvRs3NOaMY652DS4iR4M4K1p7q65dW/Ob1cJIknSPOicD/lxBPI49OTWHLvSOZPTyaDTm1THh5M9e9v5u9xY2OLu8PEUKwcObdNKcdJrqyg+k7WykICac5zIejFQNYnj8Wv8BiYoeXM2V/DKoOF/SNT+HXuod5/b3ZMvZOSvMvgvrN+OsW4qw2M7syk3GtFYSkOf9fX//FjWzYcdjRuytJ0l/gnAr8f/F11bNwfAJb7x3JXWPi2F/SyEWvbePyN7ezOa+2x1y9K4Tg/uvuoT5pP7HlHUzb3UJueBSmIA9yq/rxbe5EfP1LiMksY+LBSAxtnnS1/oOYhjXcl+bJd2NnUNlwMx0lhQTpbkWraeCK9j6My7HTZ6gNV2srdrvC9V+f4OVPZF9fks5152Tg/4u7Uctto2LZcu9IHpyYSGGdiZnv7GLaq1tZeaQKew+4elej1vDAzfdQG7uXhJIupu5p5mhkDJ0BzhTWpvLl8an4+JYRk1XMuMOheDX70dj+IWk1X/BMsjvvjR5LAw/TeKyDYO0taLUFjFHCmbQziMwh5ag1arwtTTx/yMIba/NlX1+SzmHndOD/i7Neww3Doti0IIu/XdiLxnYLN324lwte3MTX+8uw2s7uq3d1ah33z72L6sg9JBWZmbq/mezoOMx+Bkoae7H06EV4+1QQO7qQ0ccCCGoIobzrOwZVvMk70S48O7IPLb7PUrHHiL+4CyfdJlLVHly2PZ2slANofH0I7qxguyWQCY9+TlllraN3WZKkM+C8CPx/0WvUXDEgjHV3jeCly9MQAu747CBZz23g453FdFpsji7xFxm0Bu6dN5/q0D0kF1iZfKiFg3FJ2Lw0VLSm8FH2ZXh4VhE7+gQj8ryIqoki37qZEWVP83WIEw8Mj6I9/kWKt4fgYfk7broPCFHruP7gKAa57sA/MZCAjkqKrS5MeGE9G3dlO3qXJUnqZudV4P+LRq1ialowK+YP562r++HlrOeBr7MZ/sx63tp0AlPX2Tltg6uTK3fMn0NN4B5651iYnN3E/oReCDeFmo4E3jt8Fe4eNcSOyWdQoSsJlfEctWeTVfYAG301zB8SQEe/5yjclYxT+5d4657EXa0wp2QMveuOEpZwqq+vwHVfFvDKp7KvL0nnkvMy8P9FpRKMSfLnm1sG8/ENA4nxc+HJ5ccY8vQ6XlqTR1P72dfP9nb25pb511ATsI/UI1YmHW1mT1IqGqONpq443jp0Na6udcSNzaV/mZ7eZSlkK6Vklt/BYXe4cbAXzYMepfDQMGjcS4B2PjpNB9e0ZDC8VM3EsTYUtRpvazN/P2hh9tOf0SX7+pJ0TjivA/9fhBAMifHhkxsz+OqWwfQL9+KFNbkMWbyOvy0/Rk1Lp6NL/A+B7oFcN/cyan0PknbYwsScZnam9EGn76LNGsMbB6/FxaWR+DE5pFWp6Vfch2ylgcEVt1LpZObGwZ5UD5lHadGFdJZXE6KdjUZXyzR7PH3WuTJteDE2d29CuipY1ehyqq9f5+jdliTpNMnA///0DfPk7Vn9+Gn+MEYl+vPW5hMMfWY9D35zmNKGdkeX928RPhFcecsk6r2z6XPAyoT8FranpuOkMdFONK8cuAGjczPxY46RVG9jSEF/jintpFXNxUwzN2Z4UjDoMqrbbqbxuJVg1RyE5hAZKl/GrEtmeMJ23KIiCeyspMjmwoQX1sm+viT1cDLwf0FioBsvX9GHdXdlclHfYD7bXUrm3zdw57ID5Ne0Oro8AOID45k6O5NGj2Ok77UyrqiVrWn9caYZqyaSF/fNwcnQSsLoY0S3msnKyyDfbiGy7g5cLVXcMsCLff2G02J8gMq9LgSKBzAaVhCjMTJz1yhi9WvoNSASV0srtlN9/VeXrnb0bkuS9CfJwP8NET7O/O3C3mxakMWsQREsP1zJmBc2MefDvRwua3Z0eaSFpzHmxn40uefQb6eFC0rb2Nx3IK6WetCH8/y+m9E5mUgYlU1YRwdjcgZTZhf4NN5PuCmXO9M9WZWWhCX0KUq2B+Bu/icehtfx0ai4OW88XqVbmH6BN3Yh8LY28+wBMzc98xlm87kxLbUknU9k4P9Oge4GHp6cxNZ7RzI3K4atBXVMfmULV7+7i50n6h1aW0Z0BoNnJdDqWkD/bWbGVrSxod8g3LuqURvDeXbPreh0nSSMPEygtY2xxzKotWpQt/6NtMZdLOrtwcepIaiSn+HE1ij0pp/w0z2IQWNjdn0Wxi0F3DDZToezN6FdlaxscGHComWUV8m+viT1JDLw/yBvFz13jY1n68KRLBgXz5HyZi57cweXvL6NQ7VWh03bkJWURe8rgjE5FzNgs5kxVSbW9R+Cp6kcvVsoi/fchkbbRWLWIfxo44IjA2i16GnpeJ0RNSt4Od6Vl1J90PVbTMHuNGg6RpD2ZtR6E5eZe+P7g8LlWQVogqMJ6qyk0ObChOfXsWn3EYfsryRJf5wM/D/JzUnLLZkxbLl3JIsmJ1He2MHze7uY+PIWfjxU6ZCbrk/sO4Hoi91oN5YxcGMXo2pNrBk4DO/mElw8gnly9+2o1BaSMvfjrW1lfHY/rF3OFHd9xsTKj/k40siiPm7o+j9M4bEsOisbCVVdh9q5ipEimPQfo+kXs5bY1ATcLC1YFcF1n+fzz8/W/OX7KknSHycD/zQZdGquGRLJhnuyuD5FR6fFxq2f7GPMCxv5fE8plr942oaLB11E4FQVnYZKBq3rIqvexKpBI/BuLMLTJ5gndt8BKoWkEfvwMLRwwZFUNJ0eHLSu5pKyl/kpSM/d/Zyh382UVV1KU54gyHYTetcDpGjcuHjLSFz5hskjo7EBXrZmntnfxRzZ15eks54M/G6i06gYFqJl9Z0jeOXKPug1au754hCZz25gybaiv3TahpkjZuA6zkSXUw2D13aR2WBixeBMvGvy8PYP5rGdd2JHkDR8L+6uLYw5koBLux+b7Xu4suxRdvhomDvAna7e06mz30z1AVd8zA/i6vYtQRoNc45OpiX3B+ZP96VZ70W4uZIVp/r6FbKvL0lnrdMOfCGEkxBilxDioBDiiBDi0VOvewkhVgsh8k49ep5+uWc/tUowqXcQy+cN5b1r+hPg7sQj3x1h6NPreG1DAa2df81R8E0XzEaTVYNFV8+QNV2MaDLx09CR+FYexz8okEd33oXVriZxyB7cPNsYfSQKn9Yw1tgLuKrsbvJd4aYMD+rjBmHyWkjpNl9cO97G2+UlXNUwt2ocNav2cM9FNsy+kQR3VXLC5sr459exZc/Rv2QfJUn6Y7rjCL8LGKkoSiqQBowTQmQAC4G1iqLEAmtPPT9vCCHISvDjizmDWDo7g8RAN55ecZwhi9fx3KocGkxnfrqC+VPnYR1chFXbyJBVnQxtNvHjsFH4lR0hOCyQR3beg9muI2nwblx8TWQeDSSkKYYVSg2XVN5Go97G7EFeFEbGYY95nBObQtC1rSfQcA8anZVZHQOwft7IJSNyCIhLxN3SjFURXLMsj9eWyb6+JJ1tTjvwlZP+dQNZ7akvBZgKLDn1+hJg2uluqycSQpAR5c2H1w/ku7lDGBztwz/W5TNk8Toe/+EoVc1ndtqGey69m7b+x7FrWhi2qoPBLW18P2IMPsWHCYsI4JEd99BpdSJp0C5cgzsZetyTyPoUVtpbGV99C0IxMSfDmwPhfuj6LCZ/Swz2xgJC1NejNpqYRDSRXwUTG76SMUMSsQHe9hae3tfFnGeXyb6+JJ1FuqWHL4RQCyEOADXAakVRdgL+iqJUApx69OuObfVkvUM8eH1mOqvvGM74lADe31bE8GfWc99XhymuN52RbQohuP+qhTT0PoiiMjF8VScDW0+Gvm/RAaJi/Hl4xwJMZmcSB+zEJayLQblGEqrTWWfvYkjDPHy7apnXz4uNYS44ZzxJ/v6+dFabCLHPRONRTn+NFxPWDqO1/SPmT42lQeNBhKWKFfXOTFy0jIpqx16nIEnSSaI7zxsXQngAXwO3AVsURfH42c8aFUX5rz6+EGI2MBvA19c3fdmyZd1Wz1+tra0NFxeX3718bbudnwotbCq3YrPDwEA1E6N0hLp2/2fpVruVn/avJqxwKApO/DBCcMTLi8mb11AcEE9FRQeLBjyLm1MLx/cOoO2EluwIJw4FbaWfSkW568Nku0Vx79FWLixVaNz/OqH+G/CIaqdCcy+KaRi1VhuvBn1NRGg6nx3yxr2tkgqdH672dm7tpSEuzKfb9+vn/uj4n216cv09uXbo+fVnZWXtVRSl328t162BDyCEeAQwATcCmYqiVAohAoENiqLE/9rvxsfHKzk5Od1az19pw4YNZGZm/uHfq2np5O0thXy0o5h2s43Rif7cmhVNn7Du/Zy7y9rFoleeICAvA0WjZ9UoPQecnZm8aTWN8QM4cqyWRwc+i7tTEzn7B9CaqyEvzJPdYWvprdFi1s9ls3df5uS3c02embacr/BTfYVfaiutLtfQWHcR7TZ4330jfab34/01LlQXF9KlNtAldNzdz5k5l47u1n36uT87/meLnlx/T64den79QojfFfjdcZaO76kje4QQBmA0cBz4Dph1arFZwLenu61zlZ+bE/dPSGTbwpHcPjqW3UUNTP/nNq56ewfb8uu67epdvUbPg7cspDJyByqLmbHru+jVbuL7YaPxzNlNeKiBh3csoLHDi/g+O3FLtBBXUs+QonEctlhRzP9kUs0aXo8x8lyKAWP8NOoNN1O23RvnliX4eTyLk1ZhdssIyj8uZsqQHDL6JmNXTvb1F+/r4pa/L8NiOTtvMCNJ57ru6B0EAuuFEIeA3Zzs4f8ALAbGCCHygDGnnku/wsOo4/bRcWxdOJL7JySQW93GlW/vZPo/t7H6aHW33HTdWefMwlvuoiJsB6ouC+M3mknqaOeHYSMJqThGau8AHt5xD/XtvsSn7sIt2UJEaSUjTkwi32ylxraMqyqW8nmonof6uKIOG0hHyAKKNgWgadlKsOF21EYrF9uScF/mjnfA98wZn0yN2p1IazXL65yZsOgzKmtkX1+S/mrdcZbOIUVR+iiK0ltRlBRFUR479Xq9oiijFEWJPfXYcPrlnh9c9BpmD49m84IsnpiWQl1bFzd+sIcJL2/m2wPlp33TdQ+jB3fdehuVwTtRm6xM3mIhrqODnwZn4Xp0O+lpATyy8x6q2/yJ67Ubj94WgsuLGZk3lTKzjWyxnjllr7A2QMud/Tzo9IuCpEWc2BiKvamEEK5G7dnOcFUAw1ZkcKz2dZ67Mo5a52AirNUUWF0Z//e1bN0rz9eXpL+SvNL2LOakVTMjI5z1d2fy/KWpWO0K85ceYNTzG1m6q4Qu65+/etfX1Zdbb7mBqsBdqFqsTN9hJaStlZ8GZ2I8sp2Bff15dOfdVLYGEZuyB4++Fvwq8xiTO4WGLsEG1WHuKH2cA54qbs3wpsHLC33/p8jfEU9nTRdBHVegCygjSePMdXum8+PWV3lxpj+KfyTetmYsimDWZ7m8+fnabhwxSZJ+jQz8HkCrVnFh3xBW3T6c12ek4+akZeFXhxnxzAbe2VJIu/nP9cRDPEO4bs6VVPvtRjRYuWi3hTBzF6sGDUefvY3B/QN4dNddlDaFEpu4F6/+ZjyrchiTO5HOLi3fqkq4s+Jeyox2bhrkQ4mXHrfBT5CfnU5LsR6/pjm4Bu8kQKPmzsKL+e7zZdwzzUZCYhx2BXyUVp7a23myr2+VfX1JOtNk4PcgKpVgXEoA380dwgfXDSDc28jjPxxl6NPreWVdHs0df/wip2i/aK6YPY1an704Naq48gAEWcysHjAM9eFtDB8QwOO776CoMZzo+P34DOzErfoYY3PHo+p0YSmNzKu6nXZNF7MHeHHUW+A16D4KK0ZSm+2Ce/3j+AR+hotGMK9xHPs/yCYr/RhXZSZRrXIn2nayrz/xkc+orJVdP0k6k2Tg90BCCIbH+fLZTYP4Ys4gUkPc+fuqkzddf3rFcWpbu/7Q+pKDk5l6/RhqvQ+glHVx9REVfjYL6/sNgsPbyBocyFN7bqegIYrI2IP4DmrHufoIY/JH4dzhzQdKO7PrbsNob+HmdA92+gt8+t5Mpe0SKnZ6Yqj/mEDvxWj0cLW5L/ZPtVhdvmTxxb0o1QcSba8h3+rK+GdXs33fsTM0apIkycDv4fpFePHetQP4cd5QRsT78vrGAoY+vY5F3x2hvKnjd68nPSKd4P7uNHgdxlbYyfV5GrwUOxv7DsR6cCujhwSxeO9t5NbFEBF9GL9hbeirsxlZMAxvUyjvW81c1nQHYZ2V3J7qysoQgUfcNBq8ZlOy2Q914w5CDHMR7nbGiTD6/tCHdUXP8861ibR5hOJra8GiqJi5NIe3vlh3BkdMks5fMvDPEclB7rx6ZV/W3DmCKalBfLSjmBHPrOeezw9yorbtt1cAxHnEkXFlMk0eR7HmdDCnUIcbClvT+tN5aAsXDA3mmX1zOVYbT1jEEQKHt6GtySazqD/BrXF8aLEyxnQ/fVvyeCDZhWVRelyCB9IZcxeFm4KxN1UQYrkKdaCJvhoPZu6Yxltr/8Y/rg4iKCoSmwI+ShtP7ulg7nOfy76+JHUzGfjnmGhfF569JJWNC7KYkRHOdwcrGPX8Rm79ZB9HKn77putjU8bS65Iwmt1y6Mpu47YKA0Yh2N4rnbZDW5kwIpjn9t/M4ZokQsKPEZTVgrr2GMNKexPVmMrSLhupXYsZU7eLv8cZeD3BGY1HJKL3g+RvjaCzzox/w+UYIsuI0Oq4M/dK3lr6FvMm2Zg8MI5qlRsx9hp+qDUy6ZHPqJJ9fUnqNjLwz1HBHgYWTUlm68KR3Dwimk05tUx8eQvXvreLvcW/HqJT06cSNdWTFtd8OvY2c2etM3q1ij1JqTQd3MrkrFBeOjCb/VW9CA7JISSrAVF7nEEV0STUDuSbLgv+yjtcVrWS9yOdWJzqgc3ZDePAR8nfl0BbqR7vyjm4x+zES6Pi7urprHp3C72TDvPQ5BQKNf7EU0Pev/r6+4//RaMmSec2GfjnOB8XPQvGJbBl4UjuHhvHwbJmLnptO5e9sZ1NubW/OG3D5UMuJ2C8ljbnItp2NHBvsxtCp2V/Qgq1B7cybWQYrx66nt2VfQgMziN0VB005JFeE0BqZSYrO9qxar7j5vKP+T5YxwPpvpj0KjyHPEZ+QTp1x1xwLXscv4hP0esEszuG0vhhF2X2T3jvmr5UGYIJpBUzaq7+9Dhvf7n+Lx45STr3yMA/T7gbtMwdGcuWe7N4eFISxfXtXP3uLqa8spUV2ZX/c9qGa0dei8uoTtqNpTRvrOWhdg8UvROHYxKoOLiN6aPDeePwLLaX9ycgsICwUTXQWEivelcGlI5nU3sLxU67WFj2Mlt81dzR3496Jwu+GQspbs6kcrc7+spPCPZ7CpWLimlKLPHfJfPp0Sf5eHYKzgEhKIqCN208sbud22RfX5JOiwz884xRp+G6oZFsXJDJ4gt70dJpYc5H+xj74ia2llv+66brt064BfWwejoMFTSsrmKR1QuLwcixyBjKDm7jkjGRvHt0BptLM/D3LyRiVAU0FxPbBEMLp7PH1MguQxGPlT3KMXe4ub83lW5W/NPmUKmfTulmb9R1OwnR3YzKXzBE48NFmyezeMXDvDAjmGG9I6kWbsRRy/e1RiYv+ozq2kbHDJ4k9XAy8M9Teo2ayweEsfbOEbx8RR80KsFbh81k/X0DH+4o/o+brt8x7Q7MGSV0OlVT92M5TwhvOlxcyQuL4sShbVwyNpIPjl/B+uKh+PqVEDG6HForCG9tJzP/co6YGlhuaOHJigXUOlm5sY8HBd52/OIupjHoWoo2BmBvrCSo7TJ0UZ0kaA3Mz57BMx8/x6yRFuaNTiAXX5JUdeRaXBn37Cp2HpB9fUn6o2Tgn+c0ahVTUoNYPm8Y8/vq8XHR89A32Qx7Zj1vbiqgrcuKEIL7Lr2P1r45dOlrqfmmjL9pvWlxc6cwKIzcQ9u5/IJoPsm9mNVFmfj4lBI9uhhMNQS11TAmZybF7c0s1dt5rOZ27CoTN6a6cihAhUfwEDrj51G4OQxzgxXfiktxTa4gQKtmQfkVLHt3OYEh+3l9Rl+KNP6Eqtowo2bGJ8d5R/b1JekPkYEvASenbejjp+HrWwbzyQ0DifN34anlJ2+6/sLqXJo7LDx41YPUpxymS9dA1eclPGPwodHDizK/II4e2s5V42P5PG8aP50YhadXBTFjChDmBnzai7jg2Cxq29t5R6PlgYb5eFnquTnFyLZQPQbPGOizgBO7o2gr1+FeMBvvXrtw0Qrmto6h6KN69jW/z+dzBmD3DEAtwAcTj+9uZ97zsq8vSb+XDHzpPwghGBzjw8c3ZPD1LYMZEOnFS2vzGLx4Hc+szOOWi++lNmEfFl0zVZ8W8ayLD7XePlR5+3Lw4A5mTozj64LJfF8wFg+PamJG5SOsLbi3H2fc0Wto77DzqnBiXutdxJmKuTvRwIpIZ9ROHjgNfIgTx5NoyHXGmPcYgbGfoTGoucKaQuC3Uby6dxEf35hKSkwIVcKVBFHHdzUn+/o1dU2OHjpJOuvJwJd+UZ8wT966uh8rbh/GmCR/3t58gqzntqD4X0h55E4smlaqPirk724+VPgFUOfhxe5DO7lmUjw/Fk3gm7zxuLvXEDsqB0EHzqYDXJA9Ezr0vGRzYlbXwwxqzubRBAOfxntiU4NHxoMUVvWlaq8b2qIPCfF/EuGtZYwqiIkbxrHw+7t5/JIAZg2J5rjiQ29NPbkWN8Y9s5KCMnlTFUn6NTLwpd+UEODGS5f3Yf3dmVzUN4Sv9tXyTd1A8oO2YlGbqP6wkOfcfSgLCKLZxZXth3dz7aREVpSM5Yucybi51RE/8ghCZ0bfvpsLjl6FsdOTF7q0TLA/z5TaTbwSo+efyT6Y1Z349V9AiXU4ZZu9UFXtIER1E+oIFalaV249OJMHP3qCCf1NPHNRb44rvkRp2+hEzeJDKt79Svb1JemXyMCXfrdwb2f+dmEvNi3IYtagBNZ3DiTPfwtW2qn54AR/9/SiKCgMk5OBzdm7uHFKImvKRrL0+DRcXBtJyMxGZbCjbtvG2GOX4d0RzMvtCn21H3NN1bd8Gq7nid5+tOvbCUqdQ5XnZIrX+2JvqMC//jIMvS2Ea7UsLJrJO+8vQ7jt5NPZGTQ5+eKkAi97G4/tamf+C19glX19SfovMvClPyzA3YmHJiWxecEUAvpMo9B/K1alk8p3C5iPnvyQcMxqLesO7+GmqUlsLB/Bx8cuxOjSTGLmQdSuAlo2MjpvGiGmWF5tacffuI67yt9mZaCWe3r50OJmJjD2EhrCr6JoUzCWBjNeeRfjMaAaT62KOxumcPijItZUvsU3czMICQmgRu1OkqqOb6sNJ/v69U2OHipJOqvIwJf+NC9nHQ9NGszsG2+gPGA7AguqLyvoW9DC8Ygo7MDK7D3cPC2FrZVD+eDIxRgMbSSO2IfaU4OtcT2ZRROIaU3l3aYGrO55PFbyDHu9VNya5EGDN/gFZ2JKnE3Rrkjaq7S4HLoe/357cdKrudY8ANevA/nbtvt597peZATrOGrzpq++gRyLG+OeXsnuQ7mOHiZJOmvIwJdOW1JgDNffOJNK/51ohY0Rhyz47akhOyoOtdXGj0f2ccv0XuyqHsR7Ry7FyamdpOF7UfvqsdatYWjZSJIaM/ikoYIS7zb+XvYQhS52rk90piJIg4dnIvbU+RQdTaCxwIj+8CMEJX6Gyk3HJMLJXDua276ex4WJrTwwMYkDZi8SnNrpRM2VHx3l/W82OHqIJOmsIANf6hbxAfFcft2FVPvuRm9XuLFaT3xOA4diEtB1dvBN9j7mTOvF3poBvJ19BXp9J8nDdqEOMGCpWc2AmoGk1WTxbX0hu731vFB5D626Tq6Ld+JEqBEngw+69LsoKUuh5oAr2rwlhPg9iTrEiUFaT2bvncknx7+mb3wj7107gDLhgatOhQ/tLNph4o4Xv8B2Gjd9l6RzwWkHvhAiVAixXghxTAhxRAgx/9TrXkKI1UKIvFOPnqdfrnQ2Sw1NZcqsC6jx2Yu9w8qsMh2zbIL9CSm4tLXyZfY+hvYL5UBdOq8fnoFOZyFpyE7UwS5YKteS1pzEwPKJrK7PY6WHF8/X3I5OaeaGeB2Hw91RqTW49VtASXtfyrZ6oirbQqD1JnTJemJ1eh4oms2LS96lUbWVb+cOxdXLmzqVC8nqer6uMjD50aXUyr6+dB7rjiN8K3CXoiiJQAZwqxAiCVgIrFUUJRZYe+q5dI4bGDWQkVcOos5rP22NXSRvNfGYh5E9yal4tzZzrLkI4WHgcF1vXjl4NTqdjeTB29GEuWAtW09SZzjDii9ma0Men7kF83TTfILMVdyaqGNrpBc20YVfn7so0w2meJ03Sn0pvmWX4DoY/LQq7qu9jC0fHebrwjf4Ys5AhicEcsTqxUBDE8fNJ/v6e2RfXzpPnXbgK4pSqSjKvlPftwLHgGBgKrDk1GJLgGmnuy2pZ8hKzGLgpb1o9DxEc3Unrj/W87CHgZ29+hDYVIenfzsdznqO1CXx8oFr0GgVkgdtQxPlhqV4AzE2L0bmz+BA4wneNkTwYPs9pLbmcm+Cnh+ifLCo2wjpdTOVfhdQtM4fS2Mnbgem0x59AGetijntI1C+dOXhzQt58YpE5o6MZWeHB71dOulAzZUfHWHJtxsdPUyS9Jfr1h6+ECIC6APsBPwVRamEk/8pAH7duS3p7DYhdQJJ08Jp9DxMfZkJn1XN3O+uZ2tqP0JrqwiNtBEd58ex+nhe3HctQgXJA7eijXXFUrSZcJWOsceuI7+lnFc0Ycy1PMbIhj0sjnfi/SgfzIZ2QmMvpy7qIoq2hdFRqyGh/CECh+xH66zlYiWe/quHMfubm7lyiCuvXNmH42YXvA0avEUHj2xvk3196bwjfumOR394RUK4ABuBJxVF+UoI0aQoisfPft6oKMp/9fGFELOB2QC+vr7py5Yt65Z6HKGtrQ0XFxdHl/GnnYn6t1ZshYNuuDcnYfRTOJxu51M3P7L2bKcwNJx0DKw+1EKk+wnmpb2NGjtH9w7BerwFJaQPzRZPViW/g1GjZ4G+npX2a/ncbyzTS83cWlCNR5cH1XX70B35hNCEEjzC26nwnUhL/Wxc2tUcNrfzavTbXBg9HrU5nJf3ddHWaSHCXkeOKoBYpYZ5wwNxdXbq1v3+M3ry+6cn1w49v/6srKy9iqL0+63luiXwhRBa4AdgpaIoz596LQfIVBSlUggRCGxQFCX+19YTHx+v5OTknHY9jrJhwwYyMzMdXcafdqbqf3PNm9StsuLekkB4Ly+ODNTyQquF0bu3kh8Vwy2JvXh66V4CjYXc3udNtCorR/cNwXK0GXVkf1rafFmd+g5anZq7jR3sbR3J60FXkFlt4alqM5oKK82mIqwH3yHEvwDflDbs4VlUK49gy22j1GLlheCPmDl1GoP9x3Lzx/vYdaKeIa4tbGt1xUsx8cbMdPr1iu32ff8jevL7pyfXDj2/fiHE7wr87jhLRwDvAMf+FfanfAfMOvX9LODb092W1DPNHj0blxFWWl3zKD7cQOp+G7c6a1jTfwjx+bm8evwo91/Zj6rOSJ7bNwcbOpLTt6BNccd6YjdGYwXj989B6VLxVJuaVI+d3F/6Chv9NNwaoqczzICrMRSnfndQ2ZZGxQ4PRNFG/E2zMWa4EqrV8GDF1fz0yXY+ynmdD67rz4xB4Wxtc2eAh5kOtFz5YTYffLvJ0UMlSWdUd/TwhwAzgZFCiAOnviYAi4ExQog8YMyp59J5at7EeYhBTbS5nKBgbx2DjgtuMqpYOXAYycezeen4MR64oh91XeE8vWcOCgaS+2xGn+oGxQex6fIYv/8mDGZ3nmo04edXyeLSJ8n2EFwXqqYlwg2dxoAxZR4V6nRK1ntCbSGeuRfhMVaDh1bF7a3jaPlKxX0bF3D/xGienJ7C3lYDQa5avEUHD29v5c6XvpR9femc1R1n6WxRFEUoitJbUZS0U1/LFUWpVxRllKIosaceG7qjYKnnWjB9AV3pFZici8nZVs3IQh3XOgmWD84kNfsAz+fk8NBV/Wkyh/LEzjkoOJOctgVDHxeMlcdo1xzmgn3X49bpz+KaCghW81L5A1QarMyKUFHsbwRhxSd5LhUegyhc5Y21sQ3nbdPwH9+C3knNDGsqiSv7c8O3cxiVoueTGzNoxIBJ606qtoGvKp2Y+uhS6hubHT1cktTt5JW20l9GCMGDlz9Ia68TtBvLyN5QwZRqZ2bo4fuhI0k/sJtnjufw8FUDaLWG8OiOW7DjQmKvLRjTnXGtyaNDt5dxB6/BszWc58qPURMcwCvVd2FWm7i5lwuF4b5YVSZCEm+iMiSLog2BdNar0K27iuARh1F7OjFeFcwl2y7nps/mYXAp57vbhhLq58ZBiyeZ7q0cMbtxwd9+Ym92nqOHTJK6lQx86S+lEioenvkw9QlH6TBUsm9FCZe1uHO5VuGb4WMYuHcHi3NyeGRGf9rtgTy87VZsuBGfvAXn/kaMNYWYjbuYcOxq/BoTebVsL+uco3ipbj7utnpuTNCyP9gXs76ViJirqImeQOneSForjGg230tI4ueoItzor3PljuybeOSTZzjUuIEv5gxmcmowG5pdGeZjxYSWKz7I5sPvZF9fOnfIwJf+clq1lkeve5Sa2IN0OlWz65tCZpm9uUhj58vMC8jYuZUnc/N4bMZAzPjz0Na52BQP4pO24jpQj6aqCJtxGxMLLyeivh8rWvbzd1MUzzbfQVR7MXck6Vnn743VrYPwkPE0xF5IZWESdcddUe1/kyCXJ9GnexGj0/FAyQ18vnQV7x55nRcv682CcfFsbtAR7eWEt+jkoW2t3C37+tI5Qga+5BB6jZ6Hr3+Yyqh9dOnr2LYsnznCnylqG5+PHM/QrRt5LDePx2cMxCZ8eXDrPCx2D2ITtuGeoUVVXUqXdhNjyqbSuyGTcn0O9zbG8UjHvQxoPsxDSU4s9XaDQDshvgMxxV5BbVsfKne7Q8FqfOpm4zrGCz+tirsbplP2tYkFG+7l2qHBvDOrH8UdGqwGd1J1jXzx775+i6OHTZJOiwx8yWFc9C48eOMDlIftpkvXwKaPjnO7UxATVFaWjp7IiE1reSQ3nydmDgS1Nw9snY/Z7k1s/A48BqkQdZWYWEdGxUgyaidics7j8c7ezLM8zsS6zbwY58RLrgaI1OHjGo0teib12v6UbvBEqcrBbf+F+E434qJTcUNnBmErU7jhu9kkh8LXtw7B1cVIttmT0Z5tp/r6y9l3JN/RwyZJf5oMfMmhPAweLLxpARUhOzFrm1n/3lEWuoUxRlj4eOxksjau5qHcEzx1dQZqtRf3b5lHl82H6LhdeA5WUDXW0mZZQ1rdQIaWTqea47zh1J+reJWZVd/xYaSOh51VWKKccdF7oY2+llrvgRSt8MDW2IR+zSQCp3WgddUxjSgmbb2YGz+/lS51MV/fOoShcb6saXRhTICCCS2XLznMR99vdvSwSdKfIgNfcjhfF1/umH075UE7sGhaWfNWNot8I8nCzEdjpzBq3Qruy83nqVkZ6LSeLNw8n06bL9Exe/AeYkW0t9FsWklCSzKjCmdwtPEIr5LMJONX3FH2LsuDtdzjKcj1dEKj0uAZdT3VIYMoWu2DuVGgWX4pwZlHUAc6M0zrxa2H5rDwk8fZWb2ed6/pz00jolhVrSbZ34iX6OTBLc0sePUb7PbumZZEkv4qMvCls0KwRzA333ATFQE7MKtaWfXaIf4WHMtwxcyHYyczas0KHj6YzSXTe2PQu3PvptupNIUTGbUPv0EtWNw9aFZ20sczhvGFN1DQdoKHq3zw1hzm0ZLn2OGj4ukEAwf1TlgVC0Ex11EWPoSiLUGYqoyo195JUNyXaJK9SNbpeaBwDl8t2cRbB95i4bgEXrwsjcONoHbzJE3XyLJSLZc/9wNtXfJm6VLPIQNfOmtE+0Yz67qZVPrvwoyJH1/ey7VlzqR1Wvlg3FRobubA99/iFWIn0NeFh7bewbaK/oSHHyJ94Aq8whsobDiIj7qBefrbsOrMvGSuY7X9BPeV3k++q517h7vxUbw3TWqIiZ5JZdx4io7H0ZDjjGrXq/jaHsI5K4BgnZpbW8dg/diN+75bRFaSC5/PGYRdpeW43YtR6iJ218HU51ZR39bl6KGTpN9FBr50VkkOTObyqy+m1m8fnZo68reXMe7HNibsbcNq9GPtwKEcDgpEpWsmLbqVfPtFvJM9iwbFmcSkzfTt+x06t0JKcosZf2IYV5umU4XCW5RzQels+jeu5s0YPRdluvNGtA6/sLEY+t3HsdaRVO33QnViLW4HL8LnQj88jBqmaYOZsjGTv738Mp2qXL6dO4TkEE/WW0O5oOsARU0Wpr+0jqrmTkcPnST9Jhn40lmnT3gf7rv7brwndnIs7jvqvXaSWJnPjHU1zF7RTN8TVnIjEtgek8QeN38a3XqzT/sSX+VdSZvKieTem+md+j1ubqWYagSjC0dxSckYittcyen8iOtyrmFw41reitEzebgz78a749HnJprDFnPicBI0laD7IZOAS2w4DQ4iWKfh5roR1LxUzrId7/LBdenMHRXPBqfeZDZuo7K5kwv/sYHShnZHD50k/SoZ+NJZydnozDVjruGFec8TfXkQq3t/wvq4dyFgI8OKD3PzygKmb2vCrcPA/vAIPtOr2ec2kt2aV9nfegtWnYrk9K3E9VqFp8sJbFYXBtWPZHLuaErqQilsXcLdR2YxtHkD70TrmTzMmfdTI2hPfoaSxruwdVpRL5uCs3kZwQv7Yw51IRVvJqxJ4esXP2BaHzVfz8+iMWk0mXUbqW82Me0fm8ivaXP00EnSL5KBL53V9Go9s1Jm8eUVX5I1NIuffL/n04TXsI7LZ8qATq6vzmPuj/sYfrieOjd3PnVS809dJp83vcih1jswOptJ6beDoJRtODkVIoQTMdZ0RhVOprSyF21Vn/Jo9tUMbt/M+1E6poxw5vV+Y8nXvUuXJQxD9mLanptG1KxYnGcl06hRyGxMwPxcMRXbNrD0liz6XTKLsXVraW9rY8o/NnGkXE68Jp2dZOBLPYKrzpX5fefzw/QfmBIzhc+qP+PR2kVYs1qY++iFPDkxgvsPbmL2iq3ElNezNcCZp92G8kTra3xfdS8+bq30H7AF15RsTJomVFjw1MbQq3Y8JWVDcT3yPQ8V3UiGbTsfROqYOCKYx2NepcmeiYd6Oy0Lh2Aq2Eefx4ZRluRKl1VD/G5Pjj21mmkBHsy59y6m1K2CjjamvrKFHw5WOHrIJOm/yMCXehR/Z38WDV7E11O+ZmDAQF498CrTfpjGPts+wgb14rGnb+XTCxJ5ZsMPXLlqLVbRydLAAdwuXuMfDYvRe6q5YPD3eCUVclw4YbVa0Oq98bcMpTY3k/AtG7i3dD7pYjvvRzmTNvQRVrrNxNO3GN/t93HwmkfBVYvbzER2K13o2o1YPitD91MZd865i0sbfsJgaWXup/tY+OUhLDa7o4dMkv5N4+gCJOnPiPKI4qWRL3Gg5gDP732ex3c8jq/GF1uRjTHhY5jx+ENc2NHJyk+WUfbdZnYl9WVjn37s0C0ioLORkX4/cFnAZ7SWJrApbwzO5ho8dV04qZJpzosj5cg+Bod+z/a4sVyTdh0XVYfzQs7TJEV9Rt3mYPJXeGML9mGvHjzKrcRWgeXrGq5Ln43n/k9532MQS3eVsK2gjn9elU5KsLujh0yS5BG+1LOl+aWxZNwS/jHyH6iFmrs23sVVy69id9VujAYnpl9/Nde+8xoXJITz8PuvctdHb+Fe1cgnmpncJt7iy9DhJGf9yNCYoxjbe9PSEooVKxhCaKkbSO81xdy9+SlaNB1MS3uRFl0Hnh4PExhUxIBWG6EVNqqsVlabOsnuNGMu6+BCn4t5srOZQKWNkoYOpryyhcU/HafTImfclBxLHuFLPZ4QgszQTGyBNlpCWnjlwCtct/I6hocMZ37f+cR5xjFj2gSaJ47l4627mfbxx0R+8jarM4axakgmO/WD8ImqYVjkBlKLTYjDw2nr1NPiegCrkyutSgrR+9pJtq7m78GTmdm5lYTOJ/g4cSJh5kRGViZQ3+rCAWsHpU1GAt0a6OWSyFt2K3fSSp6i5fWNBaw8UsXfLuxFRpS3o4dMOk/JwJfOGWqhZnrsdMZHjueT45/w9qG3ufi7i5kcPZm5aXMJdAnklsxB1A/uzzuHc/D4dClvL7qbI1Fx/DhhAt8EXcI3EQopEYdIqsrBd28SHrVudOrLMRlP0OEciaHBxretg9grQpiQu5ZvonW8P1zhwrbd9K+Jxl4VzfHqWFbQRqShlleJ4l5M7MNOi8nM5W/u4MqBYSwcn4Cbk9bRQyadZ2TgS+ccJ40T16Vcx0WxF/H24bf55NgnrChcwZWJV3JDrxvw1ruzID2ZmpSHeTX/Wiq++Za577yJU1cnyyeNZ8XAIRwOTMN1YgsJLXnoDnkTVd8X/1YbFqdCulxVFKr8eL09Dp/9ZYwv2sr9g+6kT9RepkQvIbjTjra6F3Vlffmpbg13KTG8ZvRnS4eFkS5GPttZwrpjNTwxLYXRSf6OHi7pPCIDXzpnuevduavfXVyZcCWvHniVJUeW8GXul1zf63quSrwKP70TjyZHUR59Ky9ddBG5a9dxyZrlzPrwI3YN78uP08aww603tmFairoqUAoUPGpciLdbCDI5o+gE1X6J1FjNzFjxBXVGTxb3fYgkl3Imh3yET/gmUqx6WqsiuOh4ICrrNNa1tXOpoqO13c4NH+xhYq9AHp2ajI+L3tHDJZ0HZOBL57xAl0CeGPoEVydfzUv7XuLFfS/yyfFPmJs2l8nRkwl20vFMQhjF4Vfw3MgsXttzgEvX/8Rj9/yd+kgnVlx/ASvdBlKRFEFdQhCNVTV0NB7BS9lDcnMMIWZPurwDcAUu3LUaTDZWB15Be4wv41UbCPXZgXZ0DjOsm1Htmseytkgutep4z2bk2UPVjDhaze2Do4h2kqdwSmeWDHzpvBHnGcero15ld9VuXtz7Ig9ve5glR5Ywv+98MkMzCTfoeTkxnLwwf57rm8rreYVcunkNlz25mks9Pmf/tb1YGZDB9sChdAVdAKYBbGs8hGL/Cs9OhbGVvbC7eGN3UxHWWYx28z5KbP4c8L+dQOdmevseYFbKUjTlg1lWOgKTxsxriiuHRDOLN+ejs6loKssmY3AIgdHuCJVw9JBJ55huCXwhxLvAJKBGUZSUU695AZ8BEUARcKmiKI3dsT1JOh39A/rz0YSPWFOyhpf3vcy89fPo69eXO9LvIM0vjVhnJ15PjuBouD/PRkcwZewUpu/ewox3V9JX9xY1F3/Kpoh01hsvIN95BGr7EOwtx/jEfQXupm1cVuuFWelLV2Ao2Gx4NW9HKbORnZ+CxW0UKcZ6FNcCfmyNpsU5n0fbQ1kqtBz1OcrH1UdZ8Wo8GU4exKX7E9vPH78IV4SQ4S+dvu46wn8feAX44GevLQTWKoqyWAix8NTze7tpe5J0WoQQjAkfQ2ZoJl/nfc1rB19j5k8zGRk6kvl95xPlEUWSi4H3ekVyINyfZwJ8mDwoiwuOH2LOylVc1LGRsdM2cCIiiA1MYIvLYFo97qHD0sCr3hsJbFpPVrUb9rZk8PDF4qVC3VaCtukASn0gafredLra2Iw/t2iredy9gaT63jyhMlOTuIq9LoUcKepF9pY0nF39iEn3J6afHz4hLjL8pT+tWwJfUZRNQoiI/+/lqUDmqe+XABuQgS+dZbQqLZfGX8qkqEl8ePRD3jvyHtO/m870mOncknYLfkY/0tyMfJIaza6mNhZ7uTEtKY0hVaXcvnU1vX9aQ+TEN7k87B32mYezxTyMfW5TKXCbSonfEQKaNyFOhJDeqCJQZ6QzxJ0uixlt4zp619vQOw9jnbs/d1V5cZmtln4GDcGFUxmrbaEh5msa0z7B1hpMcXE/srelYTRGEdPPn9h0f7yCnB09fFIPIxSle+7LeSrwf/hZS6dJURSPn/28UVEUz//xe7OB2QC+vr7py5Yt65Z6HKGtrQ0XFxdHl/Gnyfqh1dbKquZVbG7djEqoyHLNYrT7aAwqw7+XyVY0fIYTeWhIaKxh/oblhDWupW1sB9ZghfrWYA53DOBrjzHU6f0Rtlacm/djO+FC76pOUlXlqJy1KIqCtrme/C53NroNIshq5cI2FwI1apKMCj5qLW3qFiqivkGJWAcCrCZvmsv60laeht0ajXuoBrcw0Ls69qhfvnccKysra6+iKP1+azmHB/7PxcfHKzk5Od1SjyNs2LCBzMxMR5fxp8n6/09paymv7H+F5YXLcde7M7vXbC5PuBydWgeAoiisbWjl6ROVHG7rIEVl54Gje/A69j5NGdVYA0DToKOhMYw3vSZywHMQitCi7ShDFNvxKe1kuDUXV6MVNVDcaWSjkoCXtZkZrfV4q+MwqJ1IMmhwUwvq7V3kuB5CRC/H068IoVKwW4y0lqfSWp6GQTuQmL5hxKT74eZt+PWdOwPke8exhBC/K/DP5Fk61UKIQEVRKoUQgUDNGdyWJHWrUNdQnh7+NNckX8OL+17k2T3P8vGxj5nbZy4ToyaiEipGe7sxysuVn+qaebqwiisSBpDcZwgPlefgu/F1GpILcIvO56HKl3Db9zJ/D7iAtb7jMCeEURlr5qtKb9TlFvo05RHj1MxIWwHrRQxL3DVMqVmKrzoAVyUVoy6IAK2Owab+VO1NJ9vaTJt7AW7B+/AL34N7xHbs9rcoqY7n2Dt9MGiHEJ2aSEy6H87u8vx+6f+cySP8Z4H6n31o66UoyoJfW4c8wncsWf8v216xnRf2vsCxhmPEe8Zze/rtDAka8u8PUG2Kwnc1TTxbWMWJji5SXQ082F6D574XqA09hN1TwVhsx6+0gycDM1gfkEmncSCKyoC2pR2l0kxYRTkhKhP7OgMwCAvjLXvwbihGY9bSyzWNMOfe6NVaNECJWSGn00anyozduRwn/6P4Rm3Dxf3kcVV7Yyht5X0wqocTldKfmL5+GFx1Z2RsQL53HO33HuF3S+ALIT7l5Ae0PkA18AjwDbAMCANKgEsURWn4tfXIwHcsWf+vsyt2Vhat5OV9L1PWVsaAgAHckX4HKT4p/17Galf4vLqB54qqKOu0MNDdmYV6M+67FlPtuR27m4IxT8FYaOPhsESO+kVhNmbSZYgDux11TSdORU3QbEMnrIzTHsddMaGvrUDT0kS4IZE0z8E4qV2xA+XOWvKbLbSZrICCWl+D8CnEI2w33iGHEUKhq92D1speOIssIhOyiO4TiJNz987jI987jvWXtnQURbniF340qjvWL0lnA5VQMT5yPKPDRrMsdxlvHHyDK368ggsiLmBen3mEuYWhUQmuCPTmIn9PPqls4MWiaqY3WxjW5xHu9dZg3Pkk1UHraI/VcG92Hk2b6/hH7DGKPHUI7QjaPIdjCghCXdeO9UAT31pSGGQoITZAQ2dAGLnNDRTVvI+v8GWA9xjCTN4E6gT2foE0GHQUHzFSU+pHXXkGjao2VO4lOAUfwTdyOzrDZspMizn2VTxa2zBi46cQ0ycSnZO8/vJ80W0tne4gj/AdS9b/x7SZ21hydAlLjizBYrNwUdxFzEmdg4/B59/LdNjsfFhRx0vFNdRbrIz2duOeICd0a26i0u0oihMY9gnKc/15P8FOsXcHent/1NpMytxj0e6pB6uC3gP6WEqItdSgUeyorBZ0VSX4dOoY6DseT70/Vo3Ae0IkqiQfyo43ULznBCU5bXRZtIANvbECjX8BHuHbcfE9gaIIWhvDUDr6ExU1nZQB/dHq1H9qLOR7x7H+0pZOd5GB71iy/j+nrqOO1w++zhe5X6BT67gm+RpmJc/CWft/58mbrDbeLa/j1ZIamqw2Jvq6c69TKfbVs6iIFCgagWGXiqqj3nya0MnxUBvBtX74aEaytbkPdjNYenmi6bQQXVdOn84S3KwdCATq+iqCTWr6u2fhoffDphd4T4rBuZ8/igI1BY0Ub91H8ZFGaltPzsWv1bSi9SnBGLwHz9C9qHUdmFq96WpLwdd7LP2HTMHZ1fi7x0C+dxxLBr4D9PQ3jaz/9BQ1F/GP/f9gVfEqvJy8uKn3TVwSdwla9f/1y1usNt4oreGN0lpMNjtXuys8sOsWKjQlVITpUYTAuE1F415nvkixsSdWIT3HkwO6m2myuGFO9cbmb0A0m/GurKNvSyFRHbUnV27qILxVYaDzENy0Xli0FnynJOLSP/Df2zdV11KybgPFh2spbQzFrDgjhA29ex16/0N4hm9D716G2WyktSkGg34wyalTCA+PRqP55daPo8f+dPX0+mXgO0BPf9PI+rvH4drDPL/3efZU7yHUNZR5feYxNmIsKvF/dxRtsFh5raSGt8vqENYOlpW8TK/K5ZwINFIZokdRBM6bVXRt1fJtL8HmOGes5TfRrPgS66mlJNaXRg8N2Ozoq9qIqy9jQEshWpsFq0Uhok3FEF0/XNVutNKE7/gEvIfH/Me0DPbaAqrW/0TxwXKKm2Oot0YCoHXqROedh1voVlwCs1GEjabGYBRbKiFh44iJScXf3/8/1nW2jP2f1dPrl4HvAD39TSPr7z6KorClfAsv7HuBvMY8kr2TuSP9DgYGDvyP5WrNFv5RXMOS8lpuLVrCgqJ3aVHrKYvUUhmgB6vAeaMK1UYN3yUZ+CLgBrqsoYxxKsTXnMrWWBt5Ye7YVGqEyYJfVS2DGwvwN9VjswuCOvUMEsl44Ua1uQT3gcHEXDQclVbz82KhZAdtO7+i5EAFJaZESix9sNgNCKHgHJyDd/JHGNyrURRBc7MfbsYRDB52E87OEcDZNfZ/Rk+vXwa+A/T0N42sv/vZ7DZ+LPyRV/a/QqWpkiFBQ7gj/Q7iveL/Y7mKTjP3bt+Pc+1WXjz+N6yKEzZrB6WJeqr8nRBmcFmrQtlkYMGg6ygxxDLYcx/pZckoigtb+xZSFd2fQrsK7ArGqhb61BeT3FCICgV3i5EsWxIedgMn2g6jjXYldcYU3Px9/7NgSwfkLMe27zOqjpeT2zGUYx1j0OoUQnsrWJ1XYHfeiLNLEwBeXsNITFzMju3Hz7qx/yPOxvfOHyED3wF6+ptG1n/mdNm6WHp8KW8dfouWrhYyQzPp5dOLKPcoojyiCHUNZcumLUQOHMTbO37itm134GLvpKNEiybKxLaAAHThndClRrtKwytN17DLP5kk1RqSDAYiykdR7ZaP86hSCrymsLxFi00lUDd1EVFVwZCaY3RZ9UxqjyVR60mnzcTRxm3YfVSMmHMjnhHB/110azUc/pz6nWvZWjiMUnMarsY2ynRFuAoNfTObMes/RqXSY7ddRVbWnX/9wHaTs/m983vIwHeAnv6mkfWfeS3mFt49/C7LC5dTaar89+talRYftQ+9g3sT5RZNRYcfV259grj2Ykqb4vEyF7BKF4Wmdwdu4W3Y64y8v/oytnr3YVLxKgJdd2E1XIaLNZqqxGUMTTnBXuP1fNKaSptKjdZsJrGimOMVbqR1wJ1mFTEqJ9osTWQ3bcU7NIz+F12Ec5IvQq36r7qVykMUf/s5Kw/1w2p3xclNMG1+f5w8ajhy9A5aW7MJCryU2NgH0Wh63iyePeG982tk4DtAT3/TyPr/Wu2WdgqbCyloLqCgqYBdBbto0bRQ2lqKgoJWE8fTdU2MadzDZnUvwnaWscEzEoufjehxtSgqDUu+u4SNhnSm529k1rHv2ROn5VBMJJV+nqRGHKOfVx1HVIP5tH0WZQZvUEBV3YGqqoPRFjVzGu2EoKWhq5K8ln1EePciavBAjGm+6MLd/uOD2X0/HGX96oNYdW14WuKwmu2kDAsifWIIO/csROEnDIZwUpJfwM2ttwNH9o/rae+d/9/ZMHmaJEm/wqg1kuyTTLJPMgAbWk+GTqe1k+KWYgqaCtjdUEzVfj9mli9ny/BYElZWccgaTO6yQFImFjJj8hJ0O7r4OmYEtc56bt/3JUOO5VLqA2vSVDye7Iq7qorhttehZQDZvX3I9k/FGmhkZYuZtWUmJprszGn0YaB+IpXtJ9i6/BPitqTj6uuLMc0PY5ovOccb2f5DFZ6+5ZSo25hz2xUcWFVO9qZycndX4xl/IZkXXcmxnLvZs/cSoiJvJzx8NkL8uQu5pDNDBr4knWWcNE7Ee8UT7xXPhCiw9J3DDyueZ9yuJ9EMsxFnbie/0Zni7wKIGF7BZYOW4nGwic+ZgLa3lrv3LsW7Rcu1a8zMWG9mb2wja9LbiLW1oz3xJQnGD1F7pFPgfAHtSb58Y7axvKydKTUWbiOMAEMkx1v20iHySF2fRuv6UoRNoX+4CzU+KqrLu3C3FDHiihRShgez9Ys8Svc3sqJCkDHtU0xuz1Jw4u/UN2wiOek5nJyCHD2k0in/3ayTJOmsolUJJk24i22Xr+CFyBso9A2ln38lF4cdJbnARFiBjnG9V3Bh6I+sD02nsH8QLqpOAFR2hX65Ju77tAVhtzMsLxS3OnfqmtZjrLoT9+onMVoOYY505ouBnoztI1jibSPOLZ38knW8W/seH6hL6DKqCWruor7ECaNixLRqB/YOK97BLkyel0bYsJOtnxVv5FO0/npCfP5Oa+sRdu6aQHX1D44cPuln5BG+JPUQwxMGEBqWyg1HCilvrGFx6w6S935CcEkZviZQ4lazsjKT7xJG8nr0c3TWCMpOeNJRqkFns+NsMuHRauGJ93Kp8w7mQHwSR8KcqHFdR5nnlzR4DqTTO5NX/Fz5qdWGKmwCQ/duxi3/U16PSqEiYyIBzX4kltloPBpP4xM7cErwwpjmh1uAYMJlAzi8oYzdPxZRdtyDhMHv4hTxPNlH5lNXv4H4uEfQaFwdPYznNRn4ktSDRBr1/NA3jgfyjMypdGP4tElce2QTdRuWktDYyVDf/aytHsihKE9StSaiIuw8qe1kyFE77p2ttLq5IhQF37oyLqgrY9geLY3uvhRHDKTMz4Mq103UBgVQ65dGTlIG+fG9CC3byaht+4la9gw7eiXxfL8J/NRhJtPmztDiVqI/ridCo6K5MZ+END/iHh3I7h+LOLKpHO3eOUQNGk+l/Qmam/aSnPwc7u59HT2M5y0Z+JLUwxjUKp5PCGOAuzMLc8vIiRvKM8npHPzwdbxKyrEGaFmmGYrod4zEI/U80twAwbBS1LDbnoJP7zYajrqAFfRmK8E1FQTUfUtqZDidcbGUNusYWtDOu65t/BCloiQ0k/cuzySo6hiD9+4mZdlT7E6C98PjeSsuHvf4GIY0uDGqpJ6+e6txctWRmupL/OwUdm4s5/i6YNx838Cn14fs6bycyMjbiAi/GZVKxs9fTY64JPVQlwd609vVyPXZhVzXruW+eY+QteF71hytZWN+JqPCtrFnUDhfFLgxiQB86l2xNmrQJUFCUiWWVhVNRS405jqhWMFQUIhzQSFaoyuMGcjVId4M92ji81157PcMpDw8gS8mJePRVEd69g7Ct2/jYNwm2p1sbFJ5siY0DlVkLMntUYzJtzJ8h50+ekHUQF8O5rdyYt1VeIaNoav5FRoaNpGc9DwGQ6ijh/G8IgNfknqwJBcDK/vFc8fxEp4oayCm90gydGV8f8zMlu/7MnTiMcYEtrNXHc8lQx+H99+nbugTeLZsQ5f7E36uzfj1asZsUtNUZKDlhBHn9jaUzibsK4+hb9vEooX345zWl9d/2MrewiJ2BiaydugkdObRJOYexKe0gAavCszuxzBrdpIH5Hg68ZJ/DN5KLP1bopiiDyXUqCKnzI+m0kdpjtlCU/2lJPe+l8CAaY4exvOGDHxJ6uHcNGreTo74983UPw/wxpBTyWH38Xh+3EjE2DL6+u9i35pZ+PmHseV4CblOrjgFX4OHMOPTXoRn3T58jdX4JbdhaneltrketcEL930nqLv5Frw2bGBKn3BaP97ChUnpvHy0gWJXNYcT+2FPGUhwWRG2UjMurVZCdCV0eBVS7V5IizjCWieFNeFq1FHhRLcnMvxEfxryhtJc3J/GE18TP3gTiYmL0GrdHD2U5zwZ+JJ0DhBCMMHXgwt83Pmupon7c1vYXS0wZ91A0LFP6NdVQmRUGe6U/cfvNZ36IgbsFl+EFYx2K1jeQLF4Uddfi19NO6tfGU6Nqw9BQf507trDrcJAV1kT+UeD2BPUn5yo3piCXTG0myhtCiOqMZqJR8sxNmkpDoRjgQU0OheRq19FXuKPeLYHMqzoMmwHrqAyr4rstNsYN/4WvHwG/vfOSd1GBr4k9WBL7plLXUnRf73e360Xm7yHMuCnzzDYLTTnBHJY54dKb8PFO5AuHwMWJws6OtB11KO1N6PR21Dr7DTp7BgNrjThQoMxmNwwA+0YMNkNmCxG2p0jTj56GjFZjJgaDTjtqMJubcbVp5P2FF9yg2PJTYnF0GEisKaC0Cor/Y554tKewvFILfvi7HzutYmEii0MLL4A9dZZPHd0N5Y+bzI8fSoZQYP/41aRUveQgS9JPVhM/4z/Cny7Vo9FbwCgNSyCYo0LFqFHb/VAI1wxCzVtjYIaujCjpkvRYEZDl6LGjAbbr1yPKVDQCwt6YcFJbcag6sRbZSJUX4fKYGdvbSLe20u4Lf15Cmv7kGvvRbF/ECfCYwHQW21EttrI6BCEmKxc3GChzbmKwy4n8GuKQrUunh+ObmVRyJOE+4WQEZhBRmAG6f7pGLW//5aL0v8mA1+SerAhl85gyKUzyKk7zk0fX0VYow9R9MNkCQQbfCaywPazX1BOfmmx4YQTBmFDJcy404FOZcNJmPHWdeBlaMXFuQaV3Y6lU49GAbVNoLJqsFp1gAqVyoqbWy0Bgbm0tvhSVNobm7qSfZ0hvLD1JtxEB/HqJgYV1aDRCuqdXah2dqdKbSTfpkG02/hQZ0VjcSUrJJQ+Hnbcbc30KhxKcn06xfE7Wdq8lA+OfoBGpSHNN+3kfwBBGSR7J6ORp3X+YWd8xIQQ44CXADXwtqIoi8/0NiXpfFPdWUO9h5l69wqUukLiW7U42buIK8gn1CWcAL8UjNihKQfopEttpUFvwZx8CD+/Qpydm+nocMXJqZWfTZD5P7WanaloC6DSFMChtgBWFE+msi2AZrP7fyzXohjYbT35lwZWoAOoA+hCQ9e/l7MBK13hxwAdiiGA8CSFmDpXgmov4MkBPrgGDWBn9Q52VOzglQOv8MqBV3DVutI/oD8ZQRkMChxEuFv4f8zsKf1vZzTwxcmp8l4FxgBlwG4hxHeKohw9k9uVpPOFzdbJnr0XY+8o5vmQduxAnr6OH7JHMFBXgjbaSLO9gipdI4pKBQH/97tBQceIjjzw7+dWq5bOThcMhrZ/v2ZXBNsr+lPYEk5FWwAVbQG0Wv7E9AgGgdIFwv6/p2NXH2zgX/Nq1qoENZ46OtK8+NCUSezeNmbXCG7sM4vEzFSy67PZUbGDHZU7WFe6DoAA5wAyAk+G/4DAAbL//wvO9BH+ACBfUZQTAEKIpcBUQAa+JHUDIdR4eQ2hod6OzZSLCoVQ52rSNWUogE0tMHd0oDU1YzO6YTO6/Pt3q6piiI7ZA4DNokdjccHgWgJAh1VPTbsvOY0xfJZz4S9uXwEUgxrFqDn55az59/e923MwVSmUVXiis3VhtusAsLtqsca7Y3fXglpAlx3RaUV02BCdNkSHDXWZCd2eOszpPijaDgIbY4j4yo1q9W6ajIXYorXc4ncLAakBFFoL+fT4p3yT/w3f5H8DwEtZLzEybOQZGfOe7IzeAEUIcTEwTlGUG049nwkMVBRl7s+WmQ3MBvD19U1ftmzZGavnTGtra8PFxeW3FzxLyfod63TrV5QOoBCFAirKm6iqdsbU6Iq9pYwyQzEqu4LOrkcv3HFT+eGi9ceit+IdkENQ8HFM6Hj/yBVU/I/2zP9iN6ixe+pBK04GN+LkDdFPfU7AqWxRV3cguuy/uT6bpw7FS4+iVyPaLGhKTNgNaswZvqA7efwf0GEn1GQnY9ePZGxZisYOGhuoTz3+OFhP5egMxrmPw0Pj8bvHrqe/d7Kyss6KG6D8r6baf/wPoyjKm8CbcPKOVz35rjM9/a45sn7H6u76FUWhvb2Q5pZ9NDbtI692N4cbSskx1XHAXEQTVizCileXJynHBuDZFkVLVzAqNLhgxoIaMyqU//nPGFQdNlQd7d1Wr7rRDI3m/9qGbncd5kF+oBJUGVRUGVQcyRpJoy0Hld2CsJtR2S0YhBV6J/D81Bcwav7YjVd6+nvn9zrTgV8G/HyyjBCg4gxvU5IkTl6M5ewchbNzFEGBF5OcCJOsrbS0HKK5eT9NzfsobjhA2fYayqtNlITswTW4E09bCP1MKYSaA3CxanGx6cCup8Ouocjgjlk40Wx24rPOjr9kPxRn7X8dOrYbXPlo8oL/ufxrmw+z0G0DE6ImEOcZ9xdU2HOc6cDfDcQKISKBcuBy4MozvE1Jkn6BRuOKl9cQvLyGAJCmKNSZPqLpm6+w7anBx6sOkyqbPMNxclz1HHbVUmxU06YRmNRQGvggFqcosNjR7TCjarehaMT/9fCNGuxGDehVKGoBGhWoBYpGICx29FtqfledtoE+aIwa9IBVCCw/OwMnrcHKAa//jK5wnR2VSk2LxUKIZSvvH3mfd7Lf4ZK4S3h40MPdNn493RkNfEVRrEKIucBKTp6W+a6iKEfO5DYlSfr9hBD4Tp6J7+SZ/37N09ZBsGJnlNoAdgV7ezt2k4na0mx+LC8lp3otte3teHkEkOCWQIvaTJO6gzpdIxVOTXQ5tdGmseBUZyeyQk21p5ZaDy0mvZbOCHesNg02d1DcAJ0WRaVFQQNCi3LqC6HGItT8628IlWImuHEVoraUIkVNZIcTod7OuBp0qIQgRuuEVq1C6AQCQb3TZL7J/waDxuCQcT1bnfHz8BVFWQ4sP9PbkSSpe6jVPwtJNahdXVG7uhIUEMCN/f/vR4qiUGWq4kj9EZbvW4fi3EZjw1F0dh0x+hBCE0IJ6RdCiGsIoa6hBBsD8Tb6IoSg09ZJq7mVNnMbrZZWWrtaaTW3YbI2/Pu15q5WWswdNFvb6TC30KFrotmnFXtHM20WK0crFdQq0KgFBwQoKCiKgh07KOBr8OXalGv/+gE8i8lL1SRJ+lOEEAS6BBLoEoimUPOHPvQ0aAwYNAb8jH5/atvtZiuf7irlrU0nqGrppHeIO7dkxjA2yR+VSl6A9UvkTcwlSepxjDoN1w+NZOOCTBZf2IvmDgtzPtrLBS9u4uv9ZVhtv30a6PlIBr4kST2WXqPm8gFhrL1zBC9dnoZKCO747CBZz23gox3FdFpsv72S84gMfEmSejyNWsXUtGB+mj+Mt6/uh7eznge/yWb4M+t5a9MJTF1WR5d4VpCBL0nSOUOlEoxO8ufrWwbzyQ0DifV34cnlxxjy9DpeXJNLU7v5t1dyDpMf2kqSdM4RQjA4xofBMT7sK2nkn+sLeHFNHm9tOsGMjHCuHxqJn5uTo8v8y8nAlyTpnNY3zJO3Z/XjeFULr20o4K3NJ3hvWxGX9gvhpuHRhHqdPzdWkYEvSdJ5ISHAjZcu78Mdo+N4Y1MBn+0u5dNdpUxNC6Kf8fw4q0cGviRJ55UIH2f+dmFv5o2K5e3NhXyys4SvLTY2Ne7l1qwYeoX89kyhPZX80FaSpPNSoLuBhyYlseXeLCZFa9laUMfkV7Yw852d7DxRz5mcOt5RZOBLknRe83bRc1Gsjm0LR7JgXDxHK1q47M0dXPL6dtYfrzmngl8GviRJEuDqpOWWzBi23DuSR6ckU9HUwbXv72biy1v48VAltl+4PWNPIgNfkiTpZww6NbMGR7Dhniyevbg3nRYbt36yjzHPb2TZnlIsPXjaBhn4kiRJ/4NOo+KSfqGsvnMEr17ZFyetmgVfHCLz2Q0s2VbUI6dtkIEvSZL0K9QqwcTegfw4byjvXdufQHcnHvnuCEOfXsc/N+TT2mlxdIm/mzwtU5Ik6XcQQpAV70dWvB87T9Tz6oYCnlmRw2sbCrhmcATXDonEy1nn6DJ/lQx8SZKkP2hglDcDo7w5XNbMq+vz+ce6fN7eXMiVA8O4cVgUAe5n57QNMvAlSZL+pF4h7rw+M5286lZe21jA+9uK+GB7ERenn5y2IcLH2dEl/gfZw5ckSTpNsf6uPH9pGhvuzuSy/qF8ua+ckc9tYP7S/RyvanF0ef8mA1+SJKmbhHoZeWJaL7YsyOLGYVGsOVrNuBc3c8OSPewvaXR0eTLwJUmSupufmxP3TUhk68KR3D46lt1FDUz/5zauensH2/LrHHb1rgx8SZKkM8TDqOP20XFsXTiSByYkklvdxpVv72T6P7ex+mg19r/46t3TCnwhxCVCiCNCCLsQot//97P7hBD5QogcIcQFp1emJElSz+Wi13Dj8Cg2L8jiiWkp1LV1ceMHexj/0ma+PVD+l910/XSP8LOBC4FNP39RCJEEXA4kA+OAfwoh1Ke5LUmSpB7NSatmRkY4G+7O5IXLUrErCvOXHmDU8xv5dFcJXdYze/XuaQW+oijHFEXJ+R8/mgosVRSlS1GUQiAfGHA625IkSTpXaNQqpvcJYeXtw3l9RjruBi33fXWY4c+sp66t68xt9wytNxjY8bPnZade+y9CiNnA7FNPu4QQ2Weopr+CD1Dn6CJOg6zfsXpy/T25djhL6i8GfB/4U78a/3sW+s3AF0KsAQL+x48eUBTl21/6tf/x2v/8dEJRlDeBN09ta4+iKP3+13I9gazfsWT9jtOTa4dzo/7fs9xvBr6iKKP/xPbLgNCfPQ8BKv7EeiRJkqRucqZOy/wOuFwIoRdCRAKxwK4ztC1JkiTpdzjd0zKnCyHKgEHAj0KIlQCKohwBlgFHgRXArYqi/J6Pn988nXrOArJ+x5L1O05Prh3Ok/rFuXS/RkmSJOmXySttJUmSzhMy8CVJks4TZ13gCyEWCSHKhRAHTn1NcHRNf4YQ4m4hhCKE8HF0LX+EEOJxIcShU2O/SggR5Oia/gghxLNCiOOn9uFrIYSHo2v6vX5tqpKzmRBi3KkpVPKFEAsdXc8fIYR4VwhR01Ov/xFChAoh1gshjp1678z/teXPusA/5QVFUdJOfS13dDF/lBAiFBgDlDi6lj/hWUVReiuKkgb8ADzs4Hr+qNVAiqIovYFc4D4H1/NH/M+pSs5mp6ZMeRUYDyQBV5yaWqWneJ+T07/0VFbgLkVREoEM4NZfG/+zNfB7uheABfzCxWZnM0VRfn63Bmd62D4oirJKURTrqac7OHkNSI/wK1OVnM0GAPmKopxQFMUMLOXk1Co9gqIom4AGR9fxZymKUqkoyr5T37cCx/iFWQ3g7A38uaf+JH9XCOHp6GL+CCHEFKBcUZSDjq7lzxJCPCmEKAWuoucd4f/cdcBPji7iHBcMlP7s+S9OoyKdWUKICKAPsPOXlnHIPW1/bboG4DXgcU4eWT4OPMfJf7hnjd+o/35g7F9b0R/zW9NlKIryAPCAEOI+YC7wyF9a4G/4PdN9CCEe4OSfux//lbX9lj85VcnZ7HdPoyKdOUIIF+BL4Pb/76/0/+CQwP+90zUIId7iZB/5rPJL9QshegGRwEEhBJxsJ+wTQgxQFKXqLyzxV/2B6TI+AX7kLAv836pfCDELmASMUs6yC03+5FQlZzM5jYqDCSG0nAz7jxVF+erXlj3rWjpCiMCfPZ3OyQ+yegRFUQ4riuKnKEqEoigRnPzH0PdsCvvfIoSI/dnTKcBxR9XyZwghxgH3AlMURWl3dD3ngd1ArBAiUgih4+R9ML5zcE3nDXHyyPId4JiiKM//5vJn2QEQQogPgTRO/llYBNykKEqlI2v6s4QQRUA/RVEcPu3q7yWE+JKTU63aOTlb6xxFUcodW9XvJ4TIB/RA/amXdiiKMseBJf1uQojpwD8AX6AJOKAoyll/t7hTp06/CKiBdxVFedKxFf1+QohPgUxOTo9cDTyiKMo7Di3qDxBCDAU2A4c5+W8W4P5fOrvxrAt8SZIk6cw461o6kiRJ0pkhA1+SJOk8IQNfkiTpPCEDX5Ik6TwhA1+SJOk8IQNfkiTpPCEDX5Ik6Tzx/wBdp5B5iYc/TAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig1, ax1 = plt.subplots(num=1)\n",
    "\n",
    "plt.grid()\n",
    "plt.axis([-5, 2, -10, 60])\n",
    "\n",
    "# P_VT.plot(ax10, fill=False, edgecolor='r', linewidth=2)\n",
    "\n",
    "constr_x1 = [-5, 1.71]\n",
    "constr_x2 = [39.338, -10]\n",
    "\n",
    "\n",
    "plt.plot(constr_x1, constr_x2)\n",
    "\n",
    "for i in range(N_loop):\n",
    "    x_traj_array = np.hstack(x_diff_traj_list[i])\n",
    "    ax1.plot(x_traj_array[0,:],x_traj_array[1,:])\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "SMPC.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}