Cost on Est Error.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 674,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.linalg import inv as inv\n",
    "from scipy.linalg import solve_continuous_are\n",
    "from scipy.linalg import solve_discrete_are\n",
    "import matplotlib.pyplot as plt\n",
    "import cvxpy as cp\n",
    "import time\n",
    "import scipy.optimize as opt\n",
    "\n",
    "DIM = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 675,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computer K and P_t-1 for standard LQR controller\n",
    "def lqr(A, B, Q, R, P):\n",
    "    K = -1 * inv(B.T @ P @ B + R) @ B.T @ P @ A\n",
    "    Ptm1 = A.T @ P @ A + (B.T @ P @ A).T @ K + Q\n",
    "    return K, Ptm1\n",
    "\n",
    "# Compute terms when we are trying to be robust to estimation error\n",
    "def dist_lqr(A, B, Q, R, S, P):\n",
    "    I = np.eye(DIM)\n",
    "    G = inv(Q + A.T@P@A - S)@A@P\n",
    "    M = inv(Q + A.T@P@A - S)@(Q + A.T@P@A)\n",
    "    Qp = G.T@Q@G + (I-A@G).T@P@(I-A@G)\n",
    "    Ap = I - M\n",
    "    Qpp = G.T@Q - (I-A@G).T@P@A\n",
    "    \n",
    "    K = inv(B.T@Qp@B + R) @ B.T@Qpp@Ap \n",
    "    Ptm1 = (Ap-G@B@K).T@Q@(Ap-G@M@K) + K.T@R@K + ((A-A@M)-(I-A@G)@B@K).T@P@((A-A@M)-(I-A@G)@B@K)\n",
    "    \n",
    "    return K, Ptm1\n",
    "    \n",
    "# def get_opt_dist_uPOV(xhat, K, A, B, Q, R, S, P):\n",
    "#     dist = inv(Q + A.T@P@A - S) @ (A.T@P@B@K + (Q + A.T@Q@A))@xhat\n",
    "\n",
    "def get_opt_dist_distPOV(x, K, A, B, Q, R, S, P, d_max):\n",
    "    if d_max == 0: return 0\n",
    "    # helping...?\n",
    "    dist = - inv(K.T@R@K + K.T@B.T@P@B@K - S) @ (K.T@R@K + K.T@B.T@P@(A+B@K)) @ x\n",
    "    return dist\n",
    "    \n",
    "def get_opt_dist_SDP(x, K, A, B, Q, R, S, P, d_max):\n",
    "    if d_max == 0: return 0\n",
    "    \n",
    "    M = K.T@R@K + K.T@B.T@P@K@B - S\n",
    "    q = 2*(K.T@R@K + K.T@B.T@P@(A + B@K))@x\n",
    "    r = x.T@(Q + K.T@R@K + (A + B@K).T@P@(A + B@K))@x\n",
    "    # print(r)\n",
    "    loss = lambda d: -1* (d.T@M@d + d.T@q + r)\n",
    "    fun = lambda d: - d.T@M@d - d.T@q - r[0][0]\n",
    "    d0 = np.zeros(x.shape)+0.01\n",
    "    \n",
    "    res = opt.minimize(fun, x0=d0)\n",
    "#     print(res)\n",
    "#     print(res.x)\n",
    "#     print(get_opt_dist_distPOV(x, K, A, B, Q, R, S, P, d_max))\n",
    "    \n",
    "    return res.x@x\n",
    "#     Q0 = np.block([[M, q],[q.T, r]])\n",
    "#     t = cp.Variable()\n",
    "#     T = np.block([[np.zeros(S.shape),np.zeros(q.shape)],[np.zeros(q.T.shape),1]])\n",
    "    \n",
    "#     opt = cp.Minimize(t)\n",
    "#     c1 = Q0 << t*T\n",
    "#     prob = cp.Problem(opt, [c1])\n",
    "#     prob.solve()\n",
    "#     t = t.value\n",
    "    \n",
    "\n",
    "    t = loss(np.zeros(x.shape), M, q, r)\n",
    "    print(t)\n",
    "    d1 = 2*inv(M)*(-1*q + np.sqrt(q**2 - 4*M*(r-t)))\n",
    "    d2 = 2*inv(M)*(-1*q + np.sqrt(q**2 - 4*M*(r-t)))\n",
    "    print(d1)\n",
    "    print(d2)\n",
    "#     d = - inv(K.T@R@K + K.T@B.T@P@B@K - S) @ (K.T@R@K + K.T@B.T@P@(A+B@K)) @ x\n",
    "#     print()\n",
    "    return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 676,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loss\n",
    "J_lqr = lambda x,u,A,B,Q,R,P: (x.T@Q@x + u.T@R@u + (A@x + B@u).T@P@(A@x + B@u))[0][0]\n",
    "\n",
    "# Standard LQR with disturbance\n",
    "def run_lqr(x, A, B, Q, R, S, n_steps=5, d_max=1, noise=None, method='standard',seed=None):\n",
    "    if seed is not None:\n",
    "        np.random.seed(seed)\n",
    "    '''\n",
    "    x: starting position\n",
    "    S: paramterization of ellipse, S is PD\n",
    "    d_max: disturbance within a ball of radius d_max \n",
    "    '''\n",
    "    # Compute P_t for each t assuming P_T = Q\n",
    "    Ps = [Q]\n",
    "    Ps_o = [Q]\n",
    "    for i in range(n_steps-1):\n",
    "        K, Ptm1 = lqr(A, B, Q, R, Ps[-1])\n",
    "        K_o, Ptm1_o = dist_lqr(A, B, Q, R, S, Ps_o[-1])\n",
    "        Ps.append(Ptm1)\n",
    "        Ps_o.append(Ptm1_o)\n",
    "            \n",
    "    # flip around so P[t] = P_t\n",
    "    Ps = Ps[::-1] \n",
    "    Ps_o = Ps_o[::-1] \n",
    "    cost = []\n",
    "    d = 0\n",
    "    for i in range(n_steps):\n",
    "        P = Ps[i]\n",
    "        P_o = Ps_o[i]\n",
    "        K, _ = lqr(A, B, Q, R, P)\n",
    "        K_o, _ = dist_lqr(A, B, Q, R, S, P_o)\n",
    "\n",
    "        if noise is None:\n",
    "            d = 0\n",
    "        elif noise == 'random':\n",
    "            d = 2*d_max*(np.random.random(x.shape) - 0.5)\n",
    "        elif noise == 'standard':\n",
    "            d = get_opt_dist_SDP(x, K, A, B, Q, R, S, P, d_max)\n",
    "            print(d)\n",
    "        elif noise == 'robust_to_est_error':\n",
    "            d = get_opt_dist_SDP(x, K_o, A, B, Q, R, S, P_o, d_max)\n",
    "        \n",
    "        \n",
    "        xhat = x + d \n",
    "        if method == 'standard':\n",
    "            u = K@xhat\n",
    "            cost.append(J_lqr(x,u,A,B,Q,R,P))\n",
    "        elif method == 'robust_to_est_error':\n",
    "            u = K_o@xhat\n",
    "            cost.append(J_lqr(x,u,A,B,Q,R,P_o))\n",
    "\n",
    "        x = A@x + B@u\n",
    "\n",
    "    return cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 677,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = np.eye(DIM)\n",
    "R = np.eye(DIM) \n",
    "# A = np.random.normal(0,1,(DIM,DIM))\n",
    "# B = np.random.normal(0,1,(DIM,DIM))\n",
    "\n",
    "# A = np.array([[ 0.09114661,  0.2488551 ], [-2.04227927, -0.72706596]])\n",
    "# B = np.array([[-0.03184335,  0.39405137], [-0.57855091, -0.25887833]])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 709,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.03400208]\n",
      "[-0.26582252]\n",
      "[0.06833634]\n",
      "[-0.01756501]\n",
      "[0.00451035]\n",
      "[-0.00115023]\n",
      "[0.00027956]\n",
      "[-4.63302771e-05]\n",
      "[-2.7752773e-14]\n",
      "[-1.11007447e-14]\n",
      "[1.03400208]\n",
      "[0.15360217]\n",
      "[0.02009114]\n",
      "[0.0029966]\n",
      "[0.0003899]\n",
      "[5.82429633e-05]\n",
      "[7.41242748e-06]\n",
      "[9.1645026e-07]\n",
      "[-2.5045583e-15]\n",
      "[-1.45939911e-15]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2+0lEQVR4nO3dd3xUZfb48c+BBAKhRKoiJZTQ0yDUoIBYkSLo6iqIiIKoiP5ExLKKrrqrX1wbYsEFywrq0sReEBARkCUYughIxICFIpCEmuT8/riTIWUymYRMJmHO+/WaV+7c+9x7z4xyz9zynEdUFWOMMcGrUqADMMYYE1iWCIwxJshZIjDGmCBnicAYY4KcJQJjjAlyIYEOoLjq1aunkZGRgQ7DGGMqlKSkpH2qWt/TsgqXCCIjI1mzZk2gwzDGmApFRH4ubJldGjLGmCBnicAYY4KcJQJjjAlylgiMMSbIWSIwxpggZ4nAGGOCnCUCY4wJchWuH4HxXVa2knEik4zjziv9eJbrr/P+RGY2CqiCoq6/LqqnluWaBlzTmtPMve6pZfnm5WlbyP4CzcqxmwogIbIO57f22CfstFgiKEdUlSMnch+ss9wH7YwTpw7g6cezOOKel/fgnvuAf/RkVqA/UoUiEugIjPFubO+WlgjOBIs2/85/1/xC2rG8B/eM41lknMj0+Ydp9SqVCa8aQo2qIYRXrUx4lRDOrhVGeNUQ1/zcy3PNqxLinl8lpBIiIIjrL+B6D7jniYjr76lliOflziI5NS353uffnh19jQk4SwRl5GRWNlM+38r0ZT9xbkQ1GkWEUSe8Ck3qVKdGlbwH7zwH8Cr55zkH80qV7ABqjCkdlgjKwO+HjzFu9lr+l/In13dvxt8GtKNqSOVAh2WMMYAlAr9bsX0f49/9niMnsnj+r3EMjjs30CEZY0welgj8JDtbefnrHfzri620qF+Dd0Z3IqphzUCHZYwxBVgi8IM/M05w93+TWbJ1L4NiG/HPodGEV7Wv2hhTPtnRqZQl/3KQ22etZW/acR67oiPDuzW1J2OMMeWa3xKBiDQB3gIa4vQbmq6qz+drUxt4G2jqiuVpVX3dXzH5k6ry9qqf+ftHm2lQM4w5Y3sQ2yQi0GGZUnby5ElSU1M5duxYoEMxpoCwsDAaN25MaGhosdbz5xlBJjBBVdeKSE0gSUS+VNXNudrcDmxW1YEiUh/YKiKzVPWEH+MqdRnHM7l//gY+WLeHC9o24JmrY4moXiXQYRk/SE1NpWbNmkRGRtqZnilXVJX9+/eTmppK8+bNi7Wu3xKBqv4K/OqaThORLcC5QO5EoEBNcf5F1QAO4CSQCmPb72ncOmstP+1NZ+Ilbbi1d0t7xv8MduzYMUsCplwSEerWrcvevXuLvW6Z3CMQkUggHvgu36IXgQ+APUBN4BpVzS6LmErDwuTd3DdvA+FVK/P2zd3o2bJeoEMyZcCSgCmvSvr/pt8TgYjUAOYBd6nq4XyLLwGSgQuAlsCXIvJN/nYiMgYYA9C0aVN/h1yk45lZPPbRZt5etYuukXWYel08DWuFBTosY4wpEb+WoRaRUJwkMEtV53tociMwXx3bgZ1A2/yNVHW6qiaoakL9+qVfcKk4fjlwhL+8spK3V+3ilvNbMHt0N0sCpkzVqFHD4/zp06fTtm1b2rZtS0JCAkuXLnUv69OnD23atCE2NpYuXbqQnJxc5H5OnDjBXXfdRatWrYiKimLw4MGkpqYWud5zzz3HkSNH3O/79+/PwYMHi1yvKEuXLmXAgAFFtuvTpw9r1qzxad/5Y/XVI488wtNPP13s9corvyUC13X/GcAWVX2mkGa7gH6u9g2BNsBP/orpdH215XcGTF3Ozn0ZvHp9Z+7v346Qyjakgwm8jz76iFdffZXly5fzww8/MH36dIYPH87u3bvdbWbNmsW6deu47bbbmDhxYpHbfOCBB0hLS2Pr1q1s27aNK664gqFDh7rLihcm/8H1k08+ISIiosSf7XQUte+SJILMzAp1G9Mn/jyKJQLXAxeISLLr1V9ExorIWFebx4CeIrIB+AqYpKr7/BhTiWRmZfN/n/3ATW+uofFZ1fjojl5c0uHsQIdljNtTTz3FlClTqFfPuU/VqVMnbrzxRqZNm1agbY8ePfIkCE+OHDnC66+/zrPPPkvlyk5drBtvvJGqVauyePFiUlJSaNu2LcOGDaNdu3ZcddVVHDlyhBdeeIE9e/bQt29f+vbtC0BkZCT79u1zrzNy5Ehat27NsGHDWLRoEYmJiURFRbF69WoAVq9eTY8ePYiPj6dnz55s3brVa6xHjx7lr3/9K+3atWPIkCEcPXrUvSxn3xkZGVx++eXExsbSsWNH3nvvPY+x5j7bmjt3LiNHjgRg5MiRjB07lm7dunHvvfcCsG7dOnr06EFUVBSvvfYaAOnp6fTr149OnToRHR3NwoULAUhJSaFdu3aMHj2aDh06cPHFF7vj3L59OxdeeCGxsbF06tSJHTt2ADBlyhS6dOlCTEwMkydP9vodnC5/PjW0HPB650JV9wAX+yuG0vBH2jHGv/M9q346wLVdmzB5YAfCQq1gnIFHP9zE5j35b3udnvaNajF5YIdir7dp0yY6d+6cZ15CQgKvv16wW85nn33GFVdc4X4fFxdX4FLR9u3badq0KbVq1SqwzU2bNtGyZUu2bt3KjBkzSExMZNSoUbz00kvcc889PPPMMyxZssSdlPJvd86cOcycOZMuXbowe/Zsli9fzgcffMA//vEP3n//fdq2bcs333xDSEgIixYt4oEHHmDevHmFfvaXX36Z6tWrs2XLFtavX0+nTp08fuZGjRrx8ccfA3Do0CFq167tNdb8UlNTWbFiBZUrV+aRRx5h/fr1rFq1ioyMDOLj47n88stp0KABCxYsoFatWuzbt4/u3bszaNAgALZt28Y777zDa6+9xtVXX828efMYPnw4w4YN47777mPIkCEcO3aM7OxsvvjiC7Zt28bq1atRVQYNGsSyZcs4//zzi4yzJKxnsRff/bSfce98T9qxkzz9l1iu6tw40CEZU2LDhg3jxIkTpKen5znw+3K/wJMmTZqQmJgIwPDhw3nhhRe45557vK7TvHlzoqOjAejQoQP9+vVDRIiOjiYlJQVwDtI33HAD27ZtQ0Q4efKk120uW7aM8ePHAxATE0NMTEyBNtHR0UyYMIFJkyYxYMAAzjvvvOJ+XP7yl7+4z44ABg8eTLVq1ahWrRp9+/Zl9erVXH755TzwwAMsW7aMSpUqsXv3bn7//Xf3Z4+LiwOgc+fOpKSkkJaWxu7duxkyZAjgdAgD+OKLL/jiiy+Ij48HnDONbdu2WSIoS6rKq8t+YsrnW2lWpzr/uakrbc+uVfSKJqiU5Je7v7Rv356kpCQuuOAC97ykpCQSEhLc72fNmkXnzp2ZOHEid9xxB/Pne3p+w9GyZUt27dpFWloaNWueKpaYlJTkvmGb/1FFXx5drFq1qnu6UqVK7veVKlVyX3t/6KGH6Nu3LwsWLCAlJYU+ffoUud2itG7dmrVr1/LJJ5/wt7/9jX79+vHwww8XaJf7M+TvPR4eHl5o25z3s2bNYu/evSQlJREaGkpkZKR7O7k/e+XKlfNcwspPVbn//vu55ZZbfP+Qp8HudOZz6OhJRr+VxJOf/sClHc5m4bhESwKm3Lv33nuZNGkS+/fvB5xf+QsWLChwIBERHnvsMVatWsUPP/xQ6PbCw8O54YYbuPvuu8nKcoY8feuttzhy5Ig72ezatYuVK1cCMHv2bHr16gVAzZo1SUtLK/FnOXToEOee65Rrf+ONN4psf/755zN79mwANm7cyPr16wu02bNnD9WrV2f48OFMnDiRtWvXeoy1YcOGbNmyhezsbBYsWOB1vwsXLuTYsWPs37+fpUuX0qVLFw4dOkSDBg0IDQ1lyZIl/Pzzz163UbNmTRo3bsz7778PwPHjxzly5AiXXHIJM2fOJD09HYDdu3fzxx9/FPldlJQlglw27j7EgKnfsHTrH0we2J4Xr4unZljxanYY429HjhyhcePG7tczzzzDoEGDuOmmm0hMTKRVq1b06tWL999/H0+PW1erVo0JEyYwZcoUAPflivz++c9/EhYWRuvWrYmKimLOnDksWLDA/Uu4TZs2TJs2jXbt2vHnn39y6623AjBmzBguvfRS9w3Y4rr33nu5//77iY+P9+kJnVtvvZX09HTatWvHww8/XOBeCcCGDRvo2rUrcXFxPProo/ztb3/zGOuTTz7JgAED6NmzJ+ecc47X/cbExNC3b1+6d+/OQw89RKNGjRg2bBhr1qwhOjqat956i7ZtCzwNX8B//vMfXnjhBWJiYujZsye//fYbF198Mddddx09evQgOjqaq6666rSSa1GkqEfBypuEhATNeUa4tKgq76z+hUc+3ETd8Cq8eF0nOjc7q1T3Yc4MW7ZsoV27doEOw6vMzExuvPFGsrOzefvtt/3SEzolJYUBAwawcePGUt+2OT2F/T8qIkmqmuBhFbtHcOREJn9bsJH53+/mvKh6PP/XeOqEW8E4U3GFhITwn//8J9BhmAokqBPBjr3p3Pb2Wn78I43/d2Frxl3QispWMM6YIkVGRtrZwBkkaBPBR+v3MGnueqqGVuatUV05LyqwpSuMMSZQgi4RnMjM5h+fbOGNFSl0ahrBtGGdOKd2tUCHZYwxARNUiWD3waPcPmstyb8c5KZezbnvsraEWq0gY0yQC5pEsHLHfm6blcTJLOXlYZ24LNr7o2HGGBMsgubncJ3wKkTWC+eDcYmWBEyFlpqayuDBg4mKiqJly5bceeednDjh39FdU1JSEBGmTp3qnjdu3DifOnwVJqfA2549e7jqqqsKbXfw4EFeeumlEu0jd0lqU7igSQRtzq7J/Ft70qK+51ruxlQEqsrQoUO54oor2LZtGz/++CPp6ek8+OCDxdpOTm/h4mjQoAHPP/98qSedRo0aMXfu3EKXlzQRlOQzBqugSQRgQwyaim/x4sWEhYVx4403Ak7NmmeffZaZM2dy5MgR3njjDcaNG+duP2DAAPcANTVq1GDChAnExsaycuVK7rvvPtq3b09MTEyRxeIA6tevT79+/XjzzTcLLEtOTqZ79+7ExMQwZMgQ/vzzzwJtdu7c6e4pm9OzF5yzjY4dOwJOFdWcHsAxMTFs27aN++67jx07dhAXF8fEiRMLDFCT+8wkMjKSSZMm0alTJ+bMmQM4PXfj4uLo2LFjkaWu33jjDYYOHcqll15KVFSUu+Q0OBVMO3XqRGxsLP369QMgIyODUaNG0bVrV+Lj491lpyuaoLlHYEyp+/Q++G1D6W7z7Gi47MlCF3sqN12rVi2aNm3K9u3bvW46IyODbt268a9//Yv9+/dz00038cMPPyAi7lG8PvjgA9asWcPf//53j9uYNGkSl112GaNGjcozf8SIEUydOpXevXvz8MMP8+ijj/Lcc8/laXPnnXdy6623MmLECI/jJAC88sor3Hnnne5KqVlZWTz55JNs3LjRXSU198hrntStW9ddS+iVV17hyJEjJCcns2zZMkaNGsXGjRu9lrpOTk7m+++/p2rVqrRp04Y77riDsLAwRo8ezbJly2jevDkHDhwA4IknnuCCCy5g5syZHDx4kK5du3LhhRcWKFBX3lkiMCZIVK5cmSuvvBKA2rVrExYWxk033cSAAQPcv7AHDRrkrp/vSYsWLejWrZu7yBs4ReIOHjxI7969Abjhhhv4y1/+UmDdb7/91n2wvf7665k0aVKBNj169OCJJ54gNTWVoUOHEhUVVezPec011+R5f+211wJOcbrDhw9z8OBB0tLSCi113a9fP2rXrg04VV1//vln/vzzT84//3yaN28OQJ06dQCnXPQHH3zgHrby2LFj7Nq1q9yXIcnPEoExJeXll7u/tG/fvsD19MOHD7Nr1y5atWrF+vXryc7Odi/LXUo5LCzMXU8/JCSE1atX89VXXzF37lxefPFFFi9e7FMMDzzwAFdddZX7wF8cRV2eve666+jWrRsff/wx/fv359VXX6VFixZ52oSEhBT6GcG3ctHeSl3nLxftrfCdqjJv3jzatGnj9XOVd0F1j8CYiq5fv34cOXKEt956C3BuiE6YMIGRI0dSvXp1IiMjSU5OJjs7m19++cV9TTy/9PR0Dh06RP/+/Xn22WdZt26dzzG0bduW9u3b8+GHHwLO2cVZZ53FN998AzjX5D0licTERN59913AGRvBk59++okWLVowfvx4Bg8ezPr16wuUim7WrBmbN2/m+PHjHDx4kK+++sprvO+99x4Ay5cvp3bt2tSuXbvYpa67d+/OsmXL2LlzJ4D70tAll1zC1KlT3eM4f//990VuqzyyRGBMBSIiLFiwgDlz5hAVFUXr1q0JCwvjH//4B+AcbJs3b0779u0ZP368x2EbAdLS0hgwYAAxMTH06tWLZ555BnDuEXgasCW/Bx98kNTUVPf7N998k4kTJxITE0NycrLHbTz//PNMmzaN6OjoQsdM/u9//0vHjh2Ji4tj48aNjBgxgrp165KYmEjHjh2ZOHEiTZo04eqrr6Zjx45cffXV7lG8ChMWFkZ8fDxjx45lxowZQPFLXdevX5/p06czdOhQYmNj3ZefHnroIU6ePElMTAwdOnTgoYceKnJb5ZGVoTamGCpCGWoT3EpShtrOCIwxJsj5LRGISBMRWSIim0Vkk4jcWUi7PiKS7Grztb/iMcYY45k/nxrKBCao6loRqQkkiciXqro5p4GIRAAvAZeq6i4RaeDHeIwxxnjgtzMCVf1VVde6ptOALcC5+ZpdB8xX1V2udv4bndkYY4xHZXKPQEQigXjgu3yLWgNnichSEUkSkRGFrD9GRNaIyJq9e/f6OVpjjAkufk8EIlIDmAfcpaqH8y0OAToDlwOXAA+JSOv821DV6aqaoKoJ9evbSGLGGFOa/JoIRCQUJwnMUtX5HpqkAp+raoaq7gOWAbH+jMmYii5QZahzCsPlpqo8/vjj7j4NvXv3Zv369e7lkZGRREdHExMTQ+/evfn55589bv/QoUOMGDGCVq1a0bJlS0aMGMGhQ4eKjCun/0RRchfje+WVV9wd8jxZunQpK1as8Gm7uRX2HVUE/nxqSIAZwBZVfaaQZguBXiISIiLVgW449xKMMR4Esgy1J9OmTWPFihWsW7eOH3/8kQcffJBBgwaRkZHhbrNkyRLWr19Pnz59ePzxxz1u56abbqJFixZs376dHTt20Lx5c26++eYi9+9rIsht7NixjBjh8So0ULJE4EuntHJNVf3yAnoBCqwHkl2v/sBYYGyudhOBzcBGnMtHXrfbuXNnNSZQNm/eHND9L1q0SM8777w88w4dOqR16tTRjIwMff311/X22293L7v88st1yZIlqqoaHh6ud999t8bExOg333yjkyZN0nbt2ml0dLROmDDB63537typHTp0KDC/cePGumPHjjzzhg8frq+++qqqqjZr1kz37t2rqqqffvqpXnbZZQW2sW3bNo2MjNTMzEz3vMzMTI2MjNTt27frkiVL9LzzztP+/ftr69at9ZZbbtGsrCydNGmSVqpUSWNjY/W6664rsN2ZM2dqVFSUdunSRW+++Wb39zJ58mSdMmWKqqo+//zz7u/gmmuu0Z07d2rDhg21UaNGGhsbq8uWLdMbbrhB58yZ495ueHi4qqouWbJEe/XqpQMHDtSoqCjduXOntmnTRq+77jpt27atXnnllZqRkaGqqo8++qgmJCRohw4ddPTo0Zqdna2qqr1799Z7771Xu3TpolFRUbps2TL3558wYYJ26NBBo6Oj9YUXXlBV1TVr1uj555+vnTp10osvvlj37NlT4HMX9v8osEYLOa767fFRVV0OFDkAgKpOAab4Kw5j/OWp1U/xw4EfSnWbbeu0ZVLXglU5cwS6DHVuhw8fJiMjo0BRuISEBDZv3lyg/WeffcYVV1xRYP7mzZuJi4tzF8QDp9hbXFwcmzZtolatWqxevZrNmzfTrFkzLr30UubPn8+TTz7Jiy++6C5Pnduvv/7K5MmTSUpKonbt2vTt29djKYonn3ySnTt3UrVqVQ4ePEhERARjx46lRo0a7jEacspSeLJ27Vo2btxI8+bNSUlJYevWrcyYMYPExERGjRrFSy+9xD333MO4cePcZTeuv/56PvroIwYOHAg4ZxOrV6/mk08+4dFHH2XRokVMnz6dlJQUkpOTCQkJ4cCBA5w8eZI77riDhQsXUr9+fd577z0efPBBZs6cWWh8vrKexcYEicLKUM+fP5/q1asDThlqX5JAcfTt25dzzz2XTz/91F0Suri6du1KixYtqFy5Mtdeey3Lly/32v67776jT58+1K9fnypVqhQoTZ0jJiaGYcOG8fbbbxMSUvzfxV27dnWXpgZo0qQJiYmJAAwfPtwd55IlS+jWrRvR0dEsXryYTZs2udcZOnQoAJ07dyYlJQWARYsWccstt7hjqlOnDlu3bmXjxo1cdNFFxMXF8fjjj+ep93Q6ivzkItJcVXcWNc+YYOPtl7u/lIcy1Dlq1apFeHi4u2JojqSkJC6++GL3+yVLlhAREcGwYcOYPHmyu8Bd7s+UUzG1UiXnt2l2djbJycm0b9+e1NRUj6WkS8PHH3/MsmXL+PDDD3niiSfYsKHgQEO5y15nZ2fnuTHvS8nrY8eOcdttt7FmzRqaNGnCI488kue/S07Za19KXnfo0IGVK1cW/4MWwZczgnke5hU+wKgxxm/KQxnq3CZOnMj48eM5evQo4PyS3bRpU4HB6ENCQnjuued466233CWcc7Rq1Yr4+Pg8N5Iff/xxOnXqRKtWrQBnaMmdO3eSnZ3Ne++9R69evQAIDQ3NM6hMjm7duvH111+zf/9+Tp486R62Mrec76hv37489dRTHDp0iPT09AJlryMjI0lKSgKcS2ee9pdj165d7gP17Nmz6dWrl/ugX69ePdLT072Oz5zjoosu4tVXX3UnhgMHDtCmTRv27t3r3v7JkyfznFmcjkITgYi0FZErgdoiMjTXayQQVip7N8YUSyDLUG/dupXGjRu7X3PmzOGOO+6ga9euxMTEEBkZyYgRI/jyyy8JCyt4iDjnnHO49tprPQ5TOWPGDH788UdatmxJy5Yt+fHHH/Ncm+/SpQvjxo2jXbt2NG/enCFDhgAwZswY9+Wd/Pt65JFH6NGjB4mJiR6rcWZlZTF8+HCio6OJj49n/PjxREREMHDgQBYsWEBcXBzffPMNo0eP5uuvv3aP9extGMo2bdowbdo02rVrx59//smtt95KREQEo0ePpmPHjlxyySV06dKl0PVz3HzzzTRt2pSYmBhiY2OZPXs2VapUYe7cuUyaNInY2Fji4uJK9JirJ4WWoRaRwcAVwCDgg1yL0oB3VbV0IigmK0NtAsnKUBcuPT2dIUOG0KVLlxI91lmYpUuX8vTTT/PRRx+V2jbPZCUpQ13oPQJVXQgsFJEeqlr6F6WMMWeUGjVq8OWXXwY6DFMCvtwjGCIitUQkVES+EpG9IjLc75EZYwzQp08fOxvwM18SwcXq1AgaAKQArXA6gRkTlAq7nGpMoJX0/01fEkGo6+/lwBxVLboAiDFnqLCwMPbv32/JwJQ7qsr+/fs93qgvii89KD4UkR+Ao8CtIlIfOFbEOsackRo3bkxqaipWDt2UR2FhYTRu3LjY6xWZCFT1PhH5P+CQqmaJSAYwuAQxGlPhhYaG5ulJasyZwJeexaHAcOB8V6+5r4FX/ByXMcaYMuLLpaGXce4TvOR6f71rXtE1Yo0xxpR7viSCLqqae7CYxSJSsv7oxhhjyh1fnhrKEpGWOW9EpAVQOqNaGGOMCThfzggmAktE5Cec8QWaATf6NSpjjDFlxpenhr4SkSigjWvWVlU97t+wjDHGlJVCE4GrjISo6n9cB/71rvnXi0iWqs4uqyCNMcb4j7d7BHcACzzMnw9M8E84xhhjypq3RBCqqun5Z6pqBqfKThhjjKngvCWCaiJSYAQGEakJVClqwyLSRESWiMhmEdkkInd6adtFRDJF5KrC2hhjjPEPb4lgBjBXRJrlzBCRSOBd17KiZAITVLU90B24XUTa528kIpWBp4AvihG3McaYUuJtYJqnRSQdWCYiNVyz04EnVfXlojasqr8Cv7qm00RkC3AusDlf0ztwxkUuevw2Y4wxpc7r46Oq+grwiutyEKqa5q19YVxnEvHAd/nmnwsMAfriJRGIyBhgDEDTpk1LEoIxxphC+NKzGFVNO40kUAPnF/9drgFucnsOmKSq2UXsf7qqJqhqQv369UsShjHGmEL40rO4xFyVS+cBs1R1vocmCcC7rqqm9YD+IpKpqu/7My5jjDGn+C0RiHN0nwFsUdVnPLVR1ea52r8BfGRJwBhjypZPl4ZyE5GLRORLH5om4pSsvkBEkl2v/iIyVkTGFjtSY4wxfuGtxMQFOAPQNALex3nE83WcwnNPFLVhVV3uausTVR3pa1tjjDGlx9sZwb9wntSpC8wFVgJvqGrnQq73G2OMqYC83SNQVV3qmn5fRHar6otlEJMxxpgy5C0RRIjI0Nxtc7+3swJjjDkzeEsEXwMDc71fluu94lQhNcYYU8F5KzFho5AZY0wQ8Pr4qIjEi8jbIrLW9ZouIq1cy/zaGc0YY0zZKDQRiMiVwBxgMTDS9VqFU5G0B/B5GcRnjDHGz7z9qp8MXKiqKbnmrReRxcAPgMfewsYYYyoWb5eGQvIlAQBc835W1Qf8FZQxxpiy4+2M4KSINFXVXblnugaqOe7fsIwxxaLqemW7Xlm5prMhOyvXck/Lsj0sV5wHBDk1rZqzw3zLPc0rah31bR0KTnqcofkbqG/Lir1uUYrZvjjbr9sKGhYY3+u0FXVpaJGI/ANIcs1LAO4DJpV6JMaUB9nZkHUcMo9D1olTf3NPZx53tTmR7+9xD+3ytz/hfRtZJ5wYPB3Q88z3cEA3Z77Eu+CiR0t9s94eH31fRHYCE3BGEQPYBFytqutKPRJjSpMqHE+DI/sgYz8c2e+a3uea3p9reh8cOQAnj0B2ZunFIJUhpCpUruL6WxVCqhT8W7XmqXaVq0ClEBABqeS8KlU+NS2VvSzLt7zYyyrl3bYIIK6/5Jp2vfe23F1lrDjreNtPru3l+Y7zlzPLv9zbsqLW9bLfohTYdpEr+NYsvF4xt+ubokYoWweMyD9fRJ5W1Xv8EpExnmRlwtEDuQ7groN4xv5c064Des77rBOet1W5qvMPqnpd51WnOVSrA1Wq5ztI5z6IezuYe2pXxTnYGlMBlLQvwNWAJQJzelRh/3b48+d8B3PXAT33Af/oQQq99hpW23VQrwcRTaBRrDPtPti7/oa7pquEl+AXmzFnrpImAvtXZIovOxv+2Aw/fwspy+HnFc6BPrdKIacO3uF14ezoXAfyXL/iw+u55teByqGB+TzGnCG8jUdQp7BFWCIwvsjKhN/WOwf+n1c4r2MHnWW1m0LURdCsJ9Rve+oAH1bbfq0bU8a8nREk4ZyLe/pXWcjFVxPUMk/Anu9dB/5vYdd3cCLNWVanJbQfBM0SnYN/RNPAxmqMcfP21FDzwpYZA8DJo5C6xvVrfzn88j/IPOosq98OYq6GyETn4F/z7MDGaowplBWOM747ng6pqyHFdaln9xrXkzkCZ3eEziOdA3/THn57zM0YU/osEZjCHT0Iv3x36sbur8nOc/ZSGRrFQbexzq/9pt2hWkRgYzXGlJjfEoGINAHeAhri3GuYrqrP52szDKeXsgBpwK3WWS2AMvbDLtdN3ZTl8NsGQKFSKDROgMQ7nQN/k65OJyhjzBnBp0QgIr2AKFV9XUTqAzVUdWcRq2UCE1R1rYjUBJJE5EtV3ZyrzU6gt6r+KSKXAdOBbiX4HKYk0n53ru3/vMK53LN3izM/JAwad4E+9zk3dht3gdBqgY3VGOM3RSYCEZmMU2OoDfA6EAq8DSR6W09VfwV+dU2nicgW4Fxgc642K3KtsgpoXMz4TUmk/QYLb4fti5z3VWpAk24QfRVE9oJG8U7vWGNMUPDljGAIEA+sBVDVPa5f+D4TkUjXNr7z0uwm4NNC1h8DjAFo2tQeOzwt276EBWPhRAb0fRBa9YOzY6Gy3S4yJlj58q//hKqqiCiAiIQXZwciUgOYB9ylqocLadMXJxH08rRcVafjXDYiISGhuDVhDTjVLhc9CqumQcOOcNVMqN8m0FEZY8oBXxLBf0XkVSBCREYDo4DXfNm4iITiJIFZqjq/kDYxwL+By1R1v29hm2LZtx3mjYJf10HXMXDRYxAaFuiojDHlRJGJQFWfFpGLgMM49wkeVtUvi1pPRASYAWxRVY/DWopIU2A+cL2q/lisyE3RVGHdO/DxPU6lzL/OhraXBzoqY0w548vN4ruB93w5+OeTCFwPbBCRZNe8B4CmAKr6CvAwUBd4yckbZKpqQjH3Yzw5dhg+vhs2zIFmvWDodKh9bqCjMsaUQ75cGqoJfCEiB4D3gDmq+ntRK6nqcoooTqeqNwM3+xKoKYbUJOdS0MFfoO/f4Ly7rTa+MaZQ3gavB0BVH1XVDsDtwDnA1yKyyO+RmeLLzoblz8HMi51xaG/8BHpPtCRgjPGqOM8M/gH8BuwHGvgnHFNiab/DglvgpyXQbhAMegGqnRXoqIwxFYAv9whuwxmRrD4wBxidr3ewCbTcfQMGPg+dbrCa/sYYn/lyRtAEpw9Asp9jMcWVeQK+ehRWvggNOjh9Axq0DXRUxpgKxtsIZbVcHcCmuN7nGbFMVQ/4OTbjzf4dMHeUUxG0y2i4+DGrB2SMKRFvZwSzgQF4HqlMgRZ+jMt4k/wOfDzB+gYYY0qFtxHKBrj+2khl5cWxw/DJPbD+Pacc9NDXrG+AMea0Ffn4qIh85cs842e7k+DV850OYn0fhBs+tCRgjCkV3u4RhAHVgXoichanLg3VwiknbcpCdjasnApf/R1qnA0jP4FmPQIdlTHmDOLtHsEtwF1AI5z7BDmJ4DDwon/DMoDTN+D9sbBjsfUNMMb4jbd7BM8Dz4vIHao6tQxjMuAMGrNgLBxPgwHPOQPDW98AY4wf+FJ9dKqIdATaA2G55r/lz8CCVuYJWPx3WDEVGrR37gU0aBfoqIwxZzBfh6rsg5MIPgEuA5bjDExvSlOevgE3w8WPW98AY4zf+dKz+CogFvheVW8UkYY4Yxab0rTuPadsdKUQuOZtaDcw0BEZY4KEL4ngqKpmi0imiNTCKT7XxM9xBY/jac7AMevfhaY94crXoHbjQEdljAkiviSCNSISgTM8ZRKQDqz0Z1BBY/damHcT/JkCfR6A8++xktHGmDLny83i21yTr4jIZ0AtVV3v37DOcNnZziDyix6FGg1h5MfQrGegozLGBClvHco6eVumqmv9E9IZLmM/zB8NO76CtgNg0FSoXqfo9Ywxxk+8nRH8y8syBS4o5ViCw+cPQMo3cPkzkDDK+gYYYwLOW4eyvmUZSFA4tBs2zoWuY6DLTYGOxhhjAN/6EYzwNL+oDmUi0gSnr0FDnDOI6a7eyrnbCPA80B84Aow8oy85ffcyqEL3WwMdiTHGuPny1FCXXNNhQD9gLUV3KMsEJqjqWhGpCSSJyJf5hrm8DIhyvboBL7v+nnmOHYI1b0CHIRDRNNDRGGOMmy9PDd2R+73rUdJ3fVjvV+BX13SaiGzBqVqaOxEMBt5SVQVWiUiEiJzjWvfMsvYtOJEGPccFOhJjjMmjyPEIPMgAijVYjYhEAvHAd/kWnQv8kut9Kh5KXIvIGBFZIyJr9u7dW7xoy4Osk7DqZYg8DxrFBzoaY4zJw5d7BB/iXOMHJ3G0B/7r6w5EpAYwD7jLNQZysanqdGA6QEJCghbRvPzZtAAO73aqiBpjTDnjyz2Cp3NNZwI/q2qqLxsXkVCcJDBLVed7aLKbvOUqGrvmnTlU4dsXoH5baHVhoKMxxpgCfLlH8DWAq85QiGu6jqoe8Lae64mgGcAWVX2mkGYfAONE5F2cm8SHzrj7Az8thd83wKAXoVJJrsQZY4x/+XJpaAzwd+AYkI0zUpkCLYpYNRG4HtggIsmueQ8ATQFU9RWcstb9ge04j4/eWOxPUN6tmOqUkYi5OtCRGGOMR75cGpoIdFTVfcXZsKou59TwloW1UeD24my3Qvlto1NKot/DEFI10NEYY4xHvlyr2IHza90U18ppEBoOnc+8Ex1jzJnDlzOC+4EVIvIdcDxnpqqO91tUZ4LDe2DDHKeUhBWVM8aUY74kgleBxcAGnHsExhffvQKaZeUkjDHlni+JIFRV7/Z7JGeSY4dhzevQfjCcFRnoaIwxxitf7hF86urZe46I1Ml5+T2yiuz7/8Dxw9DTrp4ZY8o/X84IrnX9vT/XPF8eHw1OOeUkmvWCcwsd28cYY8oNXzqUFauuUNDbvBAO/QL9ny66rTHGlAN+G48gKKnCt89DvdYQdXGgozHGGJ/4czyC4LNzGfy2Hga+YOUkjDEVht/GIwhKK6ZCeAOIuSbQkRhjjM/KZDyCoPD7Ztj+JXQbA6FhgY7GGGN85vfxCILGymkQWh0SbFB6Y0zF4tfxCIJG2m+w/j1IuNHKSRhjKpxCE4GItAIa5oxHkGt+oohUVdUdfo+uorByEsaYCszbPYLnAE9DSx52LTMAx9NgzUxoNxDqWB87Y0zF4y0RNFTVDflnuuZF+i2iiub7t+HYISsnYYypsLwlgggvy6qVchwVU1YmrHwJmvaExgmBjsYYY0rEWyJYIyKj888UkZuBJP+FVIFsWQiHdkHPO4pua4wx5ZS3p4buAhaIyDBOHfgTgCrAED/HVf6pwrcvQN1W0PrSQEdjjDElVmgiUNXfgZ4i0hfo6Jr9saouLpPIyruU5fBrMgx4zspJGGMqNF9KTCwBlhR3wyIyExgA/KGqHT0srw28DTR1xfG0qr5e3P0EzIqpUL0exP410JEYY8xp8edP2TcAb9dMbgc2q2os0Af4l4hU8WM8peePH2Db59B1DITafXNjTMXmt0SgqsuAA96aADVFRIAarraZ/oqnVK18EUKqQZebAx2JMcactkBe3H4RaAfsATYAd6pqtqeGrqEy14jImr1795ZljAXllJOIHwbhdQMbizHGlAJfag35yyVAMnAB0BL4UkS+UdUCvZlVdTowHSAhIUHzLy9Tq6c7w1F2v61Yq53IOkHK4RR2HNzBbxm/ke3Keeqq56eq7vd5plF3yT/3e1f73G3zb8s97Vq/0GXGmAqj29nd6N2kd6lvN5CJ4EbgSXWOZNtFZCfQFlgdwJi8O54O/5sB7QZA3ZYem5zMPsmuw7vYfnA72w9uZ8fBHWw/uJ1dh3eRpVmlGo4giAiCuN+7Jp1lOfNdbZyrcKeW5ywzxlQM4aHhZ1wi2IUz2tk3ItIQaAP8FMB4ipY8C44dhJ7jyczO5Je0X9hxcAfbDm5jx8Ed7Di4g5TDKWRmO7c6KkklmtRsQsvaLbmw6YW0imhFq7NacW6Nc6kklQoeqE8dxfMc3PMc7O3AbYwpZX5LBCLyDs7TQPVEJBWYDIQCqOorwGPAGyKyAed37CRV3eeveEoqKzuL3em72X5gK9u/f5HtzVqz4/un2Ll0JyezT7rbNa7RmFYRrejduDctI1rSKqIVzWs3JyzEBqkxxpRvfksEqnptEcv3AOVmhPdszWZP+h73pZycyzo/HfqJ41nHnUbV4JwqYbSq3pDERol5DvjVQ6sH9gMYY0wJBfLSUECoKr9l/OY+0Odc1vnp0E8czTzqbtegegNaRbTimrOvoVXtlrT8+hlaHs0gfNxSqFQ5cB/AGGNKWdAkglW/rmLq2qnsOLSDjJMZ7vn1qtWjZURLroy60v0Lv0VEC2pVqXVq5ZRvYfcGuPwZSwLGmDNO0CSCqpWrUjWkKgNbDCTqrCj3Qb921dpFr7xiKlSvC7Fer3YZY0yFFDSJIL5BPDMvmVn8Fff+CD9+Cr3vgyp2H8AYc+axsplFWfkihIRB1wJDMxhjzBnBEoE36X/Aunch7joIrxfoaIwxxi8sEXizejpknYDutwc6EmOM8RtLBIU5kQH/+ze0vRzqtQp0NMYY4zeWCAqTPBuO/mnjERtjzniWCDzJznJuEjfuAk26BToaY4zxK0sEnvzwEfyZAj3HgxV5M8ac4SwR5KcK374AZzV37g8YY8wZzhJBfrtWwe410ON2KydhjAkKlgjyWzEVqtWBuGGBjsQYY8qEJYLc9m2DrZ84g9JbOQljTJCwRJDbymlQuYqVkzDGBBVLBDnS9zp9B+KuhRoNAh2NMcaUGUsEOf73GmQdhx7jAh2JMcaUKUsEACeOwOrXoE1/qBcV6GiMMaZMWSIAWDcbjh6wchLGmKDkt0QgIjNF5A8R2eilTR8RSRaRTSLytb9i8So7y7lJfG5naNojICEYY0wg+fOM4A3g0sIWikgE8BIwSFU7AH/xYyyF++FjOPCTczZg5SSMMUHIb4lAVZcBB7w0uQ6Yr6q7XO3/8FcsXq2YChHNoN2ggOzeGGMCLZD3CFoDZ4nIUhFJEpERZR7Bru8gdbXzpJCVkzDGBKlADl4fAnQG+gHVgJUiskpVf8zfUETGAGMAmjZtWnoRrHgBwiIg3spJGGOCVyDPCFKBz1U1Q1X3AcuAWE8NVXW6qiaoakL9+vVLZ+/7dzj3B7rcDFXCS2ebxhhTAQUyESwEeolIiIhUB7oBW8ps7yunQeVQ6DqmzHZpjDHlkd8uDYnIO0AfoJ6IpAKTgVAAVX1FVbeIyGfAeiAb+LeqFvqoaanK2AfJsyDmGqjZsEx2aYwx5ZXfEoGqXutDmynAFH/FUKj//Rsyj1kHMmOMIRh7Fp88CqunQ+tLoX6bQEdjjDEBF3yJYN07cGS/nQ0YY4xLcCWC7GxY8SI0iodmiYGOxhhjyoXgSgRbP4EDO6ychDHG5BJciWDFVIhoCu0GBzoSY4wpN4InEfyyGn5ZBd1vh8qB7FBtjDHlS/AkAoCW/SB+eKCjMMaYciV4fho36QrXzw90FMYYU+4E1xmBMcaYAiwRGGNMkLNEYIwxQc4SgTHGBDlLBMYYE+QsERhjTJCzRGCMMUHOEoExxgQ5UdVAx1AsIrIX+DnQcZymesC+QAdRjtj3kZd9H6fYd5HX6XwfzVTV46DvFS4RnAlEZI2qJgQ6jvLCvo+87Ps4xb6LvPz1fdilIWOMCXKWCIwxJshZIgiM6YEOoJyx7yMv+z5Ose8iL798H3aPwBhjgpydERhjTJCzRGCMMUHOEkEZEpEmIrJERDaLyCYRuTPQMQWaiFQWke9F5KNAxxJoIhIhInNF5AcR2SIiPQIdUyCJyP9z/TvZKCLviEhYoGMqSyIyU0T+EJGNuebVEZEvRWSb6+9ZpbEvSwRlKxOYoKrtge7A7SLSPsAxBdqdwJZAB1FOPA98pqptgViC+HsRkXOB8UCCqnYEKgN/DWxUZe4N4NJ88+4DvlLVKOAr1/vTZomgDKnqr6q61jWdhvMP/dzARhU4ItIYuBz4d6BjCTQRqQ2cD8wAUNUTqnowoEEFXghQTURCgOrAngDHU6ZUdRlwIN/swcCbruk3gStKY1+WCAJERCKBeOC7AIcSSM8B9wLZAY6jPGgO7AVed10q+7eIhAc6qEBR1d3A08Au4FfgkKp+EdioyoWGqvqra/o3oGFpbNQSQQCISA1gHnCXqh4OdDyBICIDgD9UNSnQsZQTIUAn4GVVjQcyKKXT/orIde17ME6CbASEi8jwwEZVvqjz7H+pPP9viaCMiUgoThKYparzAx1PACUCg0QkBXgXuEBE3g5sSAGVCqSqas4Z4lycxBCsLgR2qupeVT0JzAd6Bjim8uB3ETkHwPX3j9LYqCWCMiQignMNeIuqPhPoeAJJVe9X1caqGolzE3CxqgbtLz5V/Q34RUTauGb1AzYHMKRA2wV0F5Hqrn83/Qjim+e5fADc4Jq+AVhYGhu1RFC2EoHrcX79Jrte/QMdlCk37gBmich6IA74R2DDCRzXmdFcYC2wAedYFVTlJkTkHWAl0EZEUkXkJuBJ4CIR2YZz1vRkqezLSkwYY0xwszMCY4wJcpYIjDEmyFkiMMaYIGeJwBhjgpwlAmOMCXKWCEzQEJG6uR7b/U1Edrum00XkJT/sr42ILHXtY4uITHfNj7PHhk15EhLoAIwpK6q6H+f5fETkESBdVZ/24y5fAJ5V1YWufUa75scBCcAnfty3MT6zMwIT9ESkT854CCLyiIi8KSLfiMjPIjJURP5PRDaIyGeuEiGISGcR+VpEkkTk85xu//mcg1M6AgBV3SAiVYC/A9e4zhSuEZFwV+351a6Cc4Nd+xgpIgtdZxXbRGSy/78NE4wsERhTUEvgAmAQ8DawRFWjgaPA5a5kMBW4SlU7AzOBJzxs51lgsYh86hpkJUJVTwAPA++papyqvgc8iFNioyvQF5iSq/JoV+BKIAb4i4gk+OtDm+Bll4aMKehTVT0pIhtwBkT5zDV/AxAJtAE6Al86ZXCojFMqOQ9VfV1EPscZXGQwcIuIxHrY38U4Bfjucb0PA5q6pr90XdJCROYDvYA1p/0JjcnFEoExBR0HUNVsETmpp+qwZOP8mxFgk6oWOZSkqu7BOWOY6RpysKOHZgJcqapb88wU6UbBMsNWE8aUOrs0ZEzxbQXq54wpLCKhItIhfyMRuTTXPYWzgbrAbiANqJmr6efAHa4qm4hIfK5lF7nGqa2GMxrVt374PCbIWSIwpphc1/mvAp4SkXVAMp5r5V8MbHS1+RyY6Co3vQRon3OzGHgMCAXWi8gm1/scq3HGr1gPzFNVuyxkSp1VHzWmnBKRkTiDt48LdCzmzGZnBMYYE+TsjMAYY4KcnREYY0yQs0RgjDFBzhKBMcYEOUsExhgT5CwRGGNMkPv/oVqA9jqE4lgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "S = np.eye(DIM)*.95\n",
    "S_o = S\n",
    "\n",
    "# Starting position, goal at 0\n",
    "A = np.array([[1]])\n",
    "B = np.array([[1]])\n",
    "x = [[1]] # ,[1]] # ,[1],[1],[1]]\n",
    "# LQR disturbance, max norm -> find sigma for \n",
    "n_steps = 10\n",
    "d_max = 0.5 # Maximum disturbance\n",
    "d_max_opt = 0.5# .1245\n",
    "# Standard LQR, d_max = 0\n",
    "no_dist_y = run_lqr(np.array(x), A, B, Q, R, S, n_steps=n_steps, d_max=0, noise=None, method='standard')\n",
    "# Standard LQR with optimal estimation error, d_max = d_max\n",
    "opt_dist_y = run_lqr(np.array(x), A, B, Q, R, S, n_steps=n_steps, d_max=d_max, noise='standard', method='standard')\n",
    "# Standard LQR with random estimation error, d_max = d_max\n",
    "rnd_dist_y = run_lqr(np.array(x), A, B, Q, R, S, n_steps=n_steps, d_max=d_max, noise='random', method='standard')\n",
    "\n",
    "# Our thing \n",
    "no_dist_y_opt = run_lqr(np.array(x), A, B, Q, R, S_o, n_steps=n_steps, d_max=0, noise=None, method='robust_to_est_error')\n",
    "# Our thing\n",
    "opt_dist_y_opt = run_lqr(np.array(x), A, B, Q, R, S_o, n_steps=n_steps, d_max=d_max_opt, noise='robust_to_est_error', method='robust_to_est_error')\n",
    "# Our thing\n",
    "rnd_dist_y_opt = run_lqr(np.array(x), A, B, Q, R, S_o, n_steps=n_steps, d_max=d_max_opt, noise='random', method='robust_to_est_error')\n",
    "\n",
    "lqr_dist_y_opt = run_lqr(np.array(x), A, B, Q, R, S_o, n_steps=n_steps, d_max=d_max_opt, noise='standard', method='robust_to_est_error')\n",
    "\n",
    "# opt_dist_y_lqr = run_lqr(np.array(x), A, B, Q, R, S, n_steps=n_steps, d_max=d_max_opt, noise='robust_to_est_error', method='standard')\n",
    "\n",
    "\n",
    "x = [i for i in range(1,n_steps+1)]\n",
    "\n",
    "# plt.figure()\n",
    "# plt.plot(x, np.cumsum(no_dist_y), label=\"LQR: No disturbance\")\n",
    "plt.plot(x, np.cumsum(opt_dist_y), label=\"LQR: Optimal disturbance\")\n",
    "# plt.plot(x, np.cumsum(rnd_dist_y), label=\"LQR: Random disturbance\")\n",
    "plt.plot(x, np.cumsum(no_dist_y_opt), label=\"Ours: No disturbance\")\n",
    "# plt.plot(x, np.cumsum(opt_dist_y_opt), label=\"Ours: Optimal disturbance\")\n",
    "# plt.plot(x, np.cumsum(rnd_dist_y_opt), label=\"Ours: Random disturbance\")\n",
    "plt.ylabel(\"Cumulative LQR Cost\")\n",
    "plt.xlabel(\"Time Step\")\n",
    "plt.plot(x, np.cumsum(lqr_dist_y_opt), label=\"Ours: LQR Opt disturbance\")\n",
    "\n",
    "\n",
    "# plt.plot(x, np.cumsum(opt_dist_y_lqr), label=\"LQR: Our Opt disturbance\") # not super informative\n",
    "\n",
    "plt.legend(loc='upper center', bbox_to_anchor=(0.73,0.94))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 680,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optimization not loss \n",
    "\n",
    "# get a game, game is not trivial\n",
    "# control uses dist, dist uses control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}