{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 1 System Modelling\n", "\n", "This notebook contains the theory on using the vehicle Kinematics Equations to derive the linearized state space model" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "ExecuteTime": { "end_time": "2024-10-22T09:43:55.338138Z", "start_time": "2024-10-22T09:43:55.319761Z" } }, "outputs": [], "source": [ "import numpy as np\n", "from scipy.integrate import odeint\n", "from scipy.interpolate import interp1d\n", "import cvxpy as cp\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "plt.style.use(\"ggplot\")\n", "\n", "import time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## kinematics model equations\n", "\n", "The variables of the model are:\n", "\n", "* $x$ coordinate of the robot\n", "* $y$ coordinate of the robot\n", "* $v$ velocity of the robot\n", "* $\\theta$ heading of the robot\n", "\n", "The inputs of the model are:\n", "\n", "* $a$ acceleration of the robot\n", "* $\\delta$ steering of the robot\n", "\n", "These are the differential equations f(x,u) of the model:\n", "\n", "$\\dot{x} = f(x,u)$\n", "\n", "* $\\dot{x} = v\\cos{\\theta}$ \n", "* $\\dot{y} = v\\sin{\\theta}$\n", "* $\\dot{v} = a$\n", "* $\\dot{\\theta} = \\frac{v\\tan{\\delta}}{L}$\n", "\n", "Discretize with forward Euler Integration for time step dt:\n", "\n", "${x_{t+1}} = x_{t} + f(x,u)dt$\n", "\n", "* ${x_{t+1}} = x_{t} + v_t\\cos{\\theta}dt$\n", "* ${y_{t+1}} = y_{t} + v_t\\sin{\\theta}dt$\n", "* ${v_{t+1}} = v_{t} + a_tdt$\n", "* ${\\theta_{t+1}} = \\theta_{t} + \\frac{v\\tan{\\delta}}{L} dt$\n", "\n", "----------------------\n", "\n", "The Model is **non-linear** and **time variant**, but the Numerical Optimizer requires a Linear sets of equations. To approximate the equivalent **LTI** State space model, the **Taylor's series expansion** is used around $\\bar{x}$ and $\\bar{u}$ (at each time step):\n", "\n", "$ f(x,u) \\approx f(\\bar{x},\\bar{u}) + \\frac{\\partial f(x,u)}{\\partial x}|_{x=\\bar{x},u=\\bar{u}}(x-\\bar{x}) + \\frac{\\partial f(x,u)}{\\partial u}|_{x=\\bar{x},u=\\bar{u}}(u-\\bar{u})$\n", "\n", "This can be rewritten usibg the State Space model form Ax+Bu :\n", "\n", "$ f(\\bar{x},\\bar{u}) + A|_{x=\\bar{x},u=\\bar{u}}(x-\\bar{x}) + B|_{x=\\bar{x},u=\\bar{u}}(u-\\bar{u})$\n", "\n", "Where:\n", "\n", "$\n", "A =\n", "\\quad\n", "\\begin{bmatrix}\n", "\\frac{\\partial f(x,u)}{\\partial x} & \\frac{\\partial f(x,u)}{\\partial y} & \\frac{\\partial f(x,u)}{\\partial v} & \\frac{\\partial f(x,u)}{\\partial \\theta} \\\\\n", "\\end{bmatrix}\n", "\\quad\n", "=\n", "\\displaystyle \\left[\\begin{matrix}0 & 0 & \\cos{\\left(\\theta \\right)} & - v \\sin{\\left(\\theta \\right)}\\\\0 & 0 & \\sin{\\left(\\theta \\right)} & v \\cos{\\left(\\theta \\right)}\\\\0 & 0 & 0 & 0\\\\0 & 0 & \\frac{\\tan{\\left(\\delta \\right)}}{L} & 0\\end{matrix}\\right]\n", "$\n", "\n", "and\n", "\n", "$\n", "B = \n", "\\quad\n", "\\begin{bmatrix}\n", "\\frac{\\partial f(x,u)}{\\partial a} & \\frac{\\partial f(x,u)}{\\partial \\delta} \\\\\n", "\\end{bmatrix}\n", "\\quad\n", "= \n", "\\displaystyle \\left[\\begin{matrix}0 & 0\\\\0 & 0\\\\1 & 0\\\\0 & \\frac{v \\left(\\tan^{2}{\\left(\\delta \\right)} + 1\\right)}{L}\\end{matrix}\\right]\n", "$\n", "\n", "are the *Jacobians*.\n", "\n", "\n", "\n", "So the discretized model is given by:\n", "\n", "$ x_{t+1} = x_t + (f(\\bar{x},\\bar{u}) + A|_{x=\\bar{x}}(x_t-\\bar{x}) + B|_{u=\\bar{u}}(u_t-\\bar{u}) )dt $\n", "\n", "$ x_{t+1} = (I+dtA)x_t + dtBu_t +dt(f(\\bar{x},\\bar{u}) - A\\bar{x} - B\\bar{u}))$\n", "\n", "The LTI-equivalent kinematics model is:\n", "\n", "$ x_{t+1} = A'x_t + B' u_t + C' $\n", "\n", "with:\n", "\n", "$ A' = I+dtA|_{x=\\bar{x},u=\\bar{u}} $\n", "\n", "$ B' = dtB|_{x=\\bar{x},u=\\bar{u}} $\n", "\n", "$ C' = dt(f(\\bar{x},\\bar{u}) - A|_{x=\\bar{x},u=\\bar{u}}\\bar{x} - B|_{x=\\bar{x},u=\\bar{u}}\\bar{u}) $" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "-----------------\n", "[About Taylor Series Expansion](https://courses.engr.illinois.edu/ece486/fa2017/documents/lecture_notes/state_space_p2.pdf):\n", "\n", "In order to linearize general nonlinear systems, we will use the Taylor Series expansion of functions.\n", "\n", "Typically it is possible to assume that the system is operating about some nominal\n", "state solution $\\bar{x}$ (possibly requires a nominal input $\\bar{u}$) called **equilibrium point**.\n", "\n", "Recall that the Taylor Series expansion of f(x) around the\n", "point $\\bar{x}$ is given by:\n", "\n", "$f(x)=f(\\bar{x}) + \\frac{df(x)}{dx}|_{x=\\bar{x}}(x-\\bar{x})$ + higher order terms...\n", "\n", "For x sufficiently close to $\\bar{x}$, these higher order terms will be very close to zero, and so we can drop them.\n", "\n", "The extension to functions of multiple states and inputs is very similar to the above procedure.Suppose the evolution of state x\n", "is given by:\n", "\n", "$\\dot{x} = f(x1, x2, . . . , xn, u1, u2, . . . , um) = Ax+Bu$\n", "\n", "Where:\n", "\n", "$ A =\n", "\\quad\n", "\\begin{bmatrix}\n", "\\frac{\\partial f(x,u)}{\\partial x1} & ... & \\frac{\\partial f(x,u)}{\\partial xn} \\\\\n", "\\end{bmatrix}\n", "\\quad\n", "$ and $ B = \\quad\n", "\\begin{bmatrix}\n", "\\frac{\\partial f(x,u)}{\\partial u1} & ... & \\frac{\\partial f(x,u)}{\\partial um} \\\\\n", "\\end{bmatrix}\n", "\\quad $\n", "\n", "Then:\n", "\n", "$f(x,u)=f(\\bar{x},\\bar{u}) + \\frac{df(x,u)}{dx}|_{x=\\bar{x}}(x-\\bar{x}) + \\frac{df(x,u)}{du}|_{u=\\bar{u}}(u-\\bar{u}) = f(\\bar{x},\\bar{u}) + A_{x=\\bar{x}}(x-\\bar{x}) + B_{u=\\bar{u}}(u-\\bar{u})$\n", "\n", "-----------------" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## ODE Model\n", "Motion Prediction: using scipy intergration" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "ExecuteTime": { "end_time": "2024-10-22T09:43:55.339159Z", "start_time": "2024-10-22T09:43:55.327823Z" } }, "outputs": [], "source": [ "# Define process model\n", "# This uses the continuous model\n", "def kinematics_model(x, t, u):\n", " \"\"\"\n", " Returns the set of ODE of the vehicle model.\n", " \"\"\"\n", "\n", " L = 0.3 # vehicle wheelbase\n", " dxdt = x[2] * np.cos(x[3])\n", " dydt = x[2] * np.sin(x[3])\n", " dvdt = u[0]\n", " dthetadt = x[2] * np.tan(u[1]) / L\n", "\n", " dqdt = [dxdt, dydt, dvdt, dthetadt]\n", "\n", " return dqdt\n", "\n", "\n", "def predict(x0, u):\n", " \"\"\" \"\"\"\n", "\n", " x_ = np.zeros((N, T + 1))\n", "\n", " x_[:, 0] = x0\n", "\n", " # solve ODE\n", " for t in range(1, T + 1):\n", "\n", " tspan = [0, DT]\n", " x_next = odeint(kinematics_model, x0, tspan, args=(u[:, t - 1],))\n", "\n", " x0 = x_next[1]\n", " x_[:, t] = x_next[1]\n", "\n", " return x_" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "ExecuteTime": { "end_time": "2024-10-22T09:43:55.355431Z", "start_time": "2024-10-22T09:43:55.333525Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.37 ms, sys: 7 μs, total: 1.38 ms\n", "Wall time: 1.42 ms\n" ] } ], "source": [ "%%time\n", "\n", "u_bar = np.zeros((M, T))\n", "u_bar[0, :] = 0.2 # m/ss\n", "u_bar[1, :] = np.radians(-np.pi / 4) # rad\n", "\n", "x0 = np.zeros(N)\n", "x0[0] = 0\n", "x0[1] = 1\n", "x0[2] = 0\n", "x0[3] = np.radians(0)\n", "\n", "x_bar = predict(x0, u_bar)" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "ExecuteTime": { "end_time": "2024-10-22T09:43:55.414551Z", "start_time": "2024-10-22T09:43:55.354098Z" } }, "outputs": [ { "data": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAEDCAYAAABTZPIVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAi0lEQVR4nO3deXwU9f3H8ddMDghHDo6QYDiCIVzKIVY8sCCWQ0yLsRar9UA51FCrtlTF+KtgOUqtihXwAgStWASJKKBcXiCtVCvEGuSQUyAkMYSEhJBjvr8/VqIhCSQhyWQ37+fjwcPsHLufT8advHf2OzOWMcYgIiIiIl7NdrsAERERETl3CnUiIiIiPkChTkRERMQHKNSJiIiI+ACFOhEREREfoFAnIiIi4gMU6kRERER8gEKdiIiIiA9QqBMRERHxAf5uF+CWo0ePUlRUdNblWrduTXp6eh1UVPfUm3fy5d6g8v35+/sTFhZWBxV5F+3b1Js38+X+6mLf1mBDXVFREYWFhWdcxrKskmV97W5q6s07+XJv4Pv91QXt29Sbt/Ll/uqqN339KiIiIuIDFOpEREREfIBCnYiIiIgPUKgTERER8QEKdSIiIiK1yOTn1cnJHwp1IiIiIrXEfLuX4sn3kZP0Wq2/lkKdiIiISC0wn3+C85cHIT2V3HffxBScrNXXa7DXqRMRERGpDcYpxixfhFm1BACrWy/C//QkabknavVrWIU6ERERkRpi8o7jzH0KvvwMAGvIddi/HIVfcCjknqjV11aoExEREakB5vABnFlTIe0QBARi3fZb7EsHltxRorYp1ImIiIicI7Pl3zjznob8E9CiFXZCIlaH8+u0BoU6ERERkWoyjoNZsRjzzuueCbEXYN/1IFZwaJ3XolAnIiIiUg3mRB7O/Kdhy6cAWIPisH51J5a/O/FKoU5ERESkikzqQZw50+DwAfD3x7olAfuKn7lak1dfpy4pKYmRI0eyYMECt0sRERGRBsJ8+RnOtAmeQBfaAvuP010PdODFR+p27drFunXr6NChg9uliIiISANgjMG8uxTz1j/AGDi/K/Y9E7FCwtwuDfDSI3X5+fk8++yz3HXXXTRt2tTtckRERMTHmfwTOC/MwCS9CsZg/XQY9oSp9SbQgZceqZs7dy59+vShZ8+eLFu27IzLFhYWUlhYWPLYsiyCgoJKfj6TU/Pr6voydUm9eSdf7g18vz8R8U4mPRVn9lQ4uA/8/LFuGoc9YJjbZZXhdaHuk08+Yc+ePUyfPr1SyyclJbF06dKSx9HR0cyYMYPWrVtX+jUjIiKqXKe3UG/eyZd7A9/vT0S8h0n5AueFJyDvOISEYd/9EFZMd7fLKpdXhbqMjAwWLFhAYmIigYGBlVonPj6euLi4ksenjgCkp6dTVFR0xnUtyyIiIoLU1NRavVebG9Sbd/Ll3qBq/fn7+1fpw5mISFUYYzBr3sK8uRCMA9GxnvFzYS3dLq1CXhXqdu/ezbFjx3j44YdLpjmOw7Zt23jvvfdYtGgRtl16mGBAQAABAQHlPl9l/ygaY3zyDyioN2/ly72B7/cnIvWbOXkS88oszOaPALCuuBrrN/dgBVTugJJbvCrUXXjhhfztb38rNe25556jbdu2jBgxokygExEREakK812a5/pz+3eDnx/WjWOwBg73irG+XhXqgoKCaN++falpjRo1onnz5mWmi4iIiFSF2f4lzvMz4Hg2NAvGvvthrC4XuF1WpXlVqBMRERGpacYYzPsrMG/MA8eB9udjJzyC1dK7xu16faibNGmS2yWIiIiIlzKFBZhX52D+9T4A1qUDsW4djxXYyOXKqs7rQ52ISH1WXFzMkiVL2LBhA1lZWYSFhTFw4ECuv/76knHAxhiWLFnC+vXrOX78OJ07d2b06NG0a9fO5epFfJvJzMB5bjrs3QmWjXXDKKzBI7xi/Fx5FOpERGrR8uXLWbt2LePHjycqKordu3czZ84cmjRpwvDhw0uWWblyJQkJCURGRrJs2TKmTJnCzJkzSy6WLiI1y+xM8QS6nGPQtDn2uD9ide/tdlnnRKeLiojUoh07dnDxxRdz0UUXER4ezqWXXkrPnj355ptvAM9RulWrVhEfH0+/fv1o374948eP5+TJk2zcuNHl6kV8jzEG58N3cZ5M9AS6qI7YiU96faADHakTEalVXbt2Ze3atRw6dIi2bduyd+9etm/fzu233w5AWloaWVlZ9OrVq2SdgIAAunfvzvbt2xk8eHC5z6tbIJZPvXmvuujPFBZiXn8B8/Fqz2td3B/7jvuwGjWutdeEutt2CnUiIrVoxIgR5OXl8cADD2DbNo7j8Otf/5r+/fsDkJWVBUBISEip9UJCQsjIyKjweXULxDNTb96rtvorzswg42+JFH+dDJZFyO3jaX7D7XUakmt72ynUiYjUok2bNrFhwwZ+97vf0a5dO/bu3cuCBQtKTpg45fQ/LGe7o4ZugVg+9ea9arM/883XFM+ZDscyoUlT7HF/JPeCvuSmptbo61Skrm6BqFAnIlKL/vGPfzBixAiuuOIKANq3b096ejpvvfUWAwcOJDQ0FKDkzNhTsrOzyxy9+zHdAvHM1Jv3qun+nA1rMIueh6IiaNsee/wjWOFtXfkd1va204kSIiK16OTJk2VuYWjbdsmOPTw8nNDQUJKTk0vmFxUVkZKSQpcuXeq0VhFfYooKcV57HvPKLE+g63Mp9sS/YoW3dbu0WqMjdSIitahv374sW7aMVq1aERUVxd69e1mxYgVXXXUV4PlaZvjw4SQlJREZGUlERARJSUk0atSoZNydiFSNyc7Cef4vsDMFAGvEb7CG/wrLx+8Rr1AnIlKL7rzzThYvXszcuXM5duwYLVq0YPDgwdxwww0ly4wYMYKCggLmzp1Lbm4uMTExJCYm6hp1ItVg9u3CmT0NjmZA4yDsMX/A6nWJ22XVCYU6EZFaFBQUxKhRoxg1alSFy1iWxciRIxk5cmTdFSbig5x/fYB5dTYUFkCb87DHJ2JFRrldVp1RqBMRERGvZoqLMUtfxqx72zOh50+wR/8eq0lTdwurYwp1IiIi4rVMTjbOi3+Frz0nG1lxN2L9/CafHz9XHoU6ERER8Upm/26cOdPguzRo1Bj7zvuxLrrc7bJco1AnIiIiXsfZ/DFm4d+hoABaR3jGz53Xwe2yXKVQJyIiIl7DOMWYZa9iVi/zTOjRB3vsH7GaNnO3sHpAoU5EBLj99turtd7kyZPp2LFjzRYjIuUyuTk4L/0NvvoCAGvYL7Hib8Gy/VyurH5QqBMRAfLz8+nTpw/BwcGVWt5xHDZs2IDjOLVcmYgAmIP7cGZPhfRUCAzEGnUf9k+udLusekWhTkTkezfccAMxMTGVWra4uJgNGzbUckUiAmD+uwln/kw4mQ8twz3j59pFu11WvaNQJyICXHPNNYSGhlZ6edu2q7yOiFSNcRzM24swK9/wTOjaE3vcg1jNK3dEvaFRqBMRgTPe8aE8lmVVeR0RqTyTl4sz7ylI/g8A1uARWL8cheWn8XMVUagTERGResUc/hZnzlRIPQgBgVi3jce+9Cq3y6r3FOpERE6TkpJS4TzbtmnSpAnnnXcefjpiIFLjnC2bceb+DfJPQItW2AmPYHWo3FjXhk6hTkTkNJMnTz7rMo0bNyYuLo5f/epXdVCRiO8zjsOx1+fi/ON5z4TYHth3PYQVHOpqXd5EoU5E5DQPPfQQ8+fPp23btlxxxRWEhISQlZXFJ598wqFDh7jxxhvZtm0bb775Js2aNeOaa65xu2QRr2by8zDznyH7i38BYF11LdbI0Vj+iilVod+WiMhptmzZQrdu3Rg/fnyp6QMHDmTWrFl8/fXX3HXXXQCsX79eoU7kHJi0QzizpsLhA+AfgH3LPVhX/MztsryS7XYBIiL1zaZNm7jiiivKnde/f38+/fRTAPr27cvhw4frsjQRn2L+9znO1D94Al1oC8JnvITdf7DbZXktHakTETnNyZMnyc7OLnfesWPHKCgoADzj6nSyhEjVGWMw772JSXoVjIHzu+J3z0Qade0B+qBUbQp1IiKn6dKlC4sXLyYmJoa2bduWTD948CBvvPEGXbt2BSAtLY2WLVu6VaaIVzIn8zEL/o75bCMA1k+HYv16HFZgoMuVeT+FOhGR04waNYrHHnuM3//+97Rr147Q0FCysrI4cOAATZs25fbbbwcgMzOTAQMGuFytiPcw6ak4c6bBt3vBzx/rpnHYA4a5XZbPUKgTETlNVFQUTz75JCtWrGDbtm0cOXKE5s2bExcXx7XXXktYWBjguVesiFSO2bYV54W/Qm4OBIdi3/0wVufubpflUxTqRETKERoayi233OJ2GSJezxiDWfc2ZsnLYBzo2Bn7nolYLVq5XZrPUagTEalAXl4eO3bsICcnhz59+tCsWTO3SxLxKqbgJObV2Zh/fwiAdfnVWLfcgxWg8XO1watCXVJSEps3b+bgwYMEBgYSGxvLLbfcUmogs4hITVi6dCnLly8vOdN1+vTpNGvWjMcff5yePXty3XXXuVugSD1nvkv3jJ/b/w3YNtaNYzwXFbYst0vzWV51nbqUlBSGDh3K1KlTefTRR3EchylTppCfn+92aSLiQ1avXs3SpUu56qqrePjhh0vNu+iii/jvf//rUmUi3sFs/x/O1N97Al2zYOzf/xl7UJwCXS3zqiN1iYmJpR4nJCQwZswYdu/eTffuGmwpIjXjvffeIy4ujltuuQXHcUrNi4yM1AWHRSpgjMF8sBLzxjwoLob2nbATHsFqGe52aQ2CV4W60+Xl5QGccZxLYWEhhYWFJY8tyyIoKKjk5zM5Nd8XP1moN+/ky71B/ekvLS2NXr16lTsvKCioZN8jIj8whQWY157HfLIOAKvfAKxbf4vVqJHLlTUcXhvqjDEsXLiQrl270r59+wqXS0pKYunSpSWPo6OjmTFjBq1bt670a0VERJxTrfWZevNOvtwbuN9fkyZNOHbsWLnz0tLSCA4OruOKROo3c/Q7nOemw54dYNlYN4zCGjzC9Q9oDY3Xhrp58+axf/9+Hn/88TMuFx8fT1xcXMnjU/+DpaenU1RUdMZ1LcsiIiKC1NRUjDHnXnQ9ot68ky/3BlXrz9/fv0ofzqriggsuYPny5Vx88cUEfn+Ve8uyKC4uZu3atRUexRNpiMyuFJzn/gLZWdC0Ofa4P2J17+12WQ2SV4a6+fPn8/nnnzN58uSz3qInICCAgICAcudV9o+iMcYn/4CCevNWvtwbuN/fjTfeyMSJE/n973/PJZdcAnjG2e3du5eMjAweeOCBKj1fZmYm//jHP9iyZQsFBQVERkZyzz330KlTJ8DT75IlS1i/fj3Hjx+nc+fOjB49mnbt2tV4byI1yfn4PcyiF6G4CM7rgD0+Eau1b3+TUJ95VagzxjB//nw2b97MpEmTCA/XwEsRqXkRERH8+c9/ZuHChaxevRqAjz/+mB49enDvvffSqlXlL5p6/Phx/u///o8ePXrwyCOPEBwczJEjR2jSpEnJMsuXL2flypUkJCQQGRnJsmXLmDJlCjNnziwZAyxSn5iiQszrL2E+fg8Aq+8VWHfch9WoscuVNWxeFermzZvHxo0befDBBwkKCiIrKwvwjH8J1I2ARaQGRUVFkZiYSGFhITk5OTRr1qxa+5nly5fTsmVLEhISSqb9+AOpMYZVq1YRHx9Pv379ABg/fjxjx45l48aNDB48+NybEalB5thRnOf/Aru2gWVhxd+KNeyXGj9XD3hVqFuzZg0AkyZNKjU9ISGBgQMH1n1BIuLzAgICaNGiRbXX/+yzz+jVqxdPPfUUKSkptGjRgiFDhvCzn/0M8Jx4kZWVVWqcXkBAAN27d2f79u0Vhjqd2V8+9Va7zO4dOHOmQlYmBDXFHjcB+8KLa+S560N/taWuevOqUPfGG2+4XYKI+KgfnyVfGTfccEOllktLS2Pt2rVce+21xMfHs2vXLl5++WUCAgIYMGBAyTcOISEhpdYLCQkhIyOjwufVmf1npt5q3vG1b3N01nQoKsS/XTSt/u9JAs6r+OoT1aVtV31eFepERGrLkiVLqrR8ZUOd4zicf/753HzzzYAnfB04cIA1a9YwYMCAkuVO/wR/thNFdGZ/+dRbzTNFRThvzMO8v8JTR59LMaMfIMMOgBq8ELe2nce5nNmvUCciAixevLjk58OHDzNt2jQGDRpE//79CQ0NJSsriw0bNvDBBx/wyCOPVPp5w8LCiIqKKjUtKiqKTz/9FIDQ0FAAsrKyCAsLK1kmOzu7zNG7H9OZ/Wem3mrotbKzcF74K+z4HwDWL27GunYk2Hat1aBtV31ede9XEZG6sGDBAgYMGEB8fDytW7cmICCA1q1bc/311/PTn/6Ul19+udLP1aVLFw4dOlRq2qFDh0o+iYeHhxMaGkpycnLJ/KKiIlJSUujSpUvNNCRSDWbfLs/9W3f8DxoHYY9PxP75r7FsRYf6SltGROQ027ZtqzBQdenSha+//rrSz3Xttdeyc+dOli1bRmpqKhs3bmT9+vUMHToU8HwtM3z4cJKSkti8eTP79+9n9uzZNGrUiP79+9dIPyJV5fz7Q5wZD0NmBrQ5D/uRv2H17ud2WXIW+vpVROQ0AQEBfPPNN1x44YVl5u3evRt//8rvOmNiYpgwYQKLFi3izTffJDw8nNtvv50rr7yyZJkRI0ZQUFDA3Llzyc3NJSYmhsTERF2jTuqcKS7GvLkAs3a5Z8KFF2OP+QNWk6buFiaVolAnInKan/zkJyxdupTGjRvTv39/mjVrxvHjx9m4cSNLly6t8hG0vn370rdv3wrnW5bFyJEjGTly5LmWLlJt5ng2zotPwLatAFjDR2KNuFlft3oRhToRkdPcfvvtHDlyhJdffpmXX34ZPz8/iouLAejWrRu33367yxWK1CxzYA/O7KnwXRo0aox9x31Yfa9wuyypIoU6EZHTBAUF8dhjj7Flyxa++uorcnJyaN68OT169KBXr14+eXFUabic/2zELHgGCk5C6wjshEewojq6XZZUg0KdiEgFevfuTe/evd0uQ6RWGKcYk/QPzHtveiZ074097o9YTZu7W5hUm0KdiIhIA2Nyj+O89AR89QUA1tB4rPjbsPz8XK5MzoVGP4qIABMmTGD//v2VXt5xHCZMmMC3335bi1WJ1DxzcL/n+nNffQGBgVhjJ2DfcIcCnQ/QkToREeDAgQMUFBTU+joibjL//RfO/Jlw8gS0DPeMn2vfye2ypIYo1ImIfO+JJ56o8NZbIt7MOA7mndcxK76/HV6XC7HvegirebC7hUmNUqirgHGKoagIJz8fU3ASgwWWBbYFlq2z30R8zIABA6q1XnCw/ihK/WbycnHmPQXJ/wHA+tkvsPR1q09SqKvI7h0Uz3iIgxXNt2zw8wM/f/Czv//v948DAiEwEAIbff9zI6wf/VwyL6gJBDX1XKn7+58JagpNvp8eEFiXHYs0aAkJCW6XIFLjTOq3nuvPpR4E/wCsW8djXz7I7bKklijUVZdxoMiBosLKLV6def4BP4S9Jk2haTOs5qEQ/MM/69TPzUOgeYg+eYmICABm639w5j0JJ/IgrBV2wkSsjp3dLktqkUJdRTrF4jd7CREREaQePoRxjCfIGQOOA8XF3/8rAudHPxcVQWGh5yKOhScxBQXf//z9f089LsiHE3mYvFzPG+5E7vf/8jz/wBMYc455/n3v9ABY6rFlQdPmZQPfqcctWkPLcGjRSkcBRUR8lHEczKolmLcXef5mxXTHvuchrOAwt0uTWqZQVwHL9sNq5I/dOAircRPPG6M6z1ONdYzjQP6JH4Jenif0meM5kJ1V8s/k/PAzx7M9NR7P9vw7tL/iAGhZOKEtOdI2iuLmoZ6g1zIcq1WbH0KfvwaLi4h4G5Ofh/PyM/DffwFgDRyOdeNo7dMbCIW6esiybc/XrU2alp5+hnWMU+wJc6cC34/CH9lZmGNZkJnuua9fYQEczaDgaEbp5/ihAAhr4Ql6LdtAq+9DX8twaNUGWrTW17wiIvWMSTuEM3saHNoP/v5YN9+NfeUQt8uSOqRQ5yMs2w+Cwzz/qDgAGmMgJwvru3RCiws4ums7JuMI5rs0yEj7IfRlZkBmBmZnyg/rnvrB3x/anIfVtj20bff9f9tD60iFPRERF5j//ddzh4i8XAhpgX3Pw1jnd3W7LKljCnUNjGVZEByGFdKCJpGRHOt8oSfofe9U6CMj7UdBr5zQd3Af5uA+zzqnVlbYExGpU8YYzHvLMEmvesZ9d+riCXShLd0uTVygUCelnAp9BIdhdepSZr5xHE+wO3wAc2i/Z+zeoQNw+IDnBBCFPfEhhw4dIjMzk4KCAoKDg2nbti1NmjRxuywRAMzJfMyCv2M+2wiA1X8w1s13Y+kC2g2WQp1UiWXb0DoCWkdg9fxJyfRzDnsdzseKjoUOMViNGtd9YyLf27FjB2vXrmXLli1kZ2eXmmfbNh07duTKK69k4MCBCnjiGpOeijNnGny7F/z8sH49FmvANbowfgNXrVC3e/duOnXSveLkB+cc9v6zwRP2bBvadsDqFAvRsVjRXSAyyvP8IrVo7969LFiwgG3btnHeeefRr18/OnXqRHBwMIGBgRw/fpwjR46wc+dOXn/9dRYvXkx8fDxxcXH4++vzsdQdZ9tWnOdnQG4ONA/BvvthrNgebpcl9UC19kQTJ04kJiaGYcOGcdlll2mHJhWqVNg7uA+zdyfs3gFZ38G3ezDf7oGPV3uCXuOg7wOe5x+dYnW9JalxjzzyCP379+e2224764fW/Px8Nm3axPLlyykuLuaXv/xlHVUpDZkxhpyk13DmPeMZP9chxnNB4Rat3S5N6olqpbGEhARWr17NrFmzeOWVV7j66qsZPHgwLVtqYKZUToVh7+h3sGc7Zvd2zJ4dsHeX55p927Zitm394avbluGegBcd6zmq1/58z63YRKrpySefJDIyslLLNm7cmEGDBjFw4EAyMjLOvoLIOTIFJzGvzibr3x8CYF12FdYtCdrvSSnVCnUDBgxgwIAB7Nq1i/fee4933nmH5cuX07dvX6655hp69NBhYKkeK6wlhF2OddHlAJjiYs9Xtnu2w+4dnqB3+AB89/3ZuZ9t9AQ9Pz+IivYczevUhcJ+V2CsAM9dNkQqobKB7sds2yY8PLwWqhH5gfku3TN+bv83YPthj7wTBsVp/JyUYRlTzVsl/Eh2djbr1q1j3bp1fPfdd0RFRTFs2DAGDBhAYGD9vB1Veno6hYVnvm+rZVlERkZy+PBhauDXVK94c2/mRB7s3YnZ833I273dc5Hl0zVpCh1jsc7vitW9l+dnLx8q4M3brTKq0l9AQACtW9fO106//e1vmTBhAh07diwzb//+/fz1r39l1qxZtfLa50r7Nt/qzez4n2f8XM4xaNac1o88wdHw83yit9P52rb7sbrat9XIXzh/f38aNWpUMrbu5MmTzJ07l2XLlvHAAw8QGxtbEy8jAoAV1AS69cLq1gv4/tp6memY3Ts8X93u2YG1f7fnvropX2BSvsC887pnbF6XC7F69MHq1hvatNUnXSlXeno6RUVF5c4rLCwkPT29jiuShsYYg/lwFWbxXM+9xdtF4zc+kcYX9obDh90uT+qpcwp1+/btY/Xq1WzcuJGioiIuvfRSfve73xETE8O+fft48cUXeemll3jiiSdqql6RMizL+uE2Zj/pj2VZRLRuzeHP/oWzezvs+AqzbQscz4GtmzFbN3u+sm0ZjtW9N3TrjdWtJ1azYHcbEa9w5MgRgoKC3C5DfJgpLMC89jzmk3UAWJf8FOu2e7Ea63JPcmbVCnWbNm1i9erVfP311wQHBxMXF8eQIUMIDQ0tWaZDhw7cdNNNTJ06taZqFak0y98fq0MMdvvzYeBwz9m2B3ZjUrZgvvoCvtnmGZe3YQ1sWIOxLM/JFj36eILe+V11A+wG5sMPP+Sjjz4qeTx37twy4a2goIB9+/bRvXv3ui5PGghz9Duc56bDnh1g2Vi/vB1ryHX6VkEqpVqh7plnnqFjx47cc8899O/fv8JLmrRu3Zorr7zynAoUqQmWbXsubNwhBq65AXMyH3Z+hflqi+co3sF9sG8XZt8uzKolENjI81Vt915Y3fp47oShnapPKygoKHWx4dzc3DJj0wICArj88ssZOXJkXZcnDYDZtQ3n+b/AsaPQpBn2uD9i9ejjdlniRaoV6iZPnkzXrme/UXCbNm1ISEiozkuI1CqrUWO4oC/WBX0BMFnfYVK2wrYtmJQtnhMvvvwM8+Vnnq9qQ1t4xuH16IPVrRdWcKhrtUvtGDJkCEOGDAFg/Pjx/OEPfyj3RAmR2uB8vBqz6AUoLoLzOmAnPIIVXvUzsqVhq1aoq0ygq02rV6/m7bffJisri6ioKEaNGkW3bt1crUm8mxXaEuvyQXD5IM+JFwf3fv9V7RbY+RVkZWL+9T78631PyIuKxurR2/NVbUx3XSvKx8yePbvWnjspKYnXX3+d4cOHM2rUKMAzKH7JkiWsX7+e48eP07lzZ0aPHk27du1qrQ6pH0xRIeb1lzAfv+eZcNHl2Hfch9VY4zal6rzu+g6bNm1iwYIFjBkzhi5durBu3TqmTZvG008/TatWrdwuT3yAZVme0BYVDUPiMYUFsGsb5qsvPF/V7t9dctcLszoJAgKhc3es7r2x+g3ECm3hdgtSDRkZGdXah2RmZtKiReW2+a5du1i3bh0dOnQoNX358uWsXLmShIQEIiMjWbZsGVOmTGHmzJk6KcOHmWNHPV+37toGloU14jdYw3+loR5SbV4X6lasWMGgQYO4+uqrARg1ahRbt25lzZo13HzzzTX2OsbAiRMWubmQl2fhY5fMwbLUW+U1gg69Pf+Gg8nOwmz/ErYnY7YlQ1YmJH8NyV9jR/fFCqy9O6v48naDH/pzo7f77ruPn/3sZ1xzzTVERESccdmioiL+85//sGzZMvr168cNN9xw1ufPz8/n2Wef5a677mLZsmUl040xrFq1ivj4ePr16wd4vv4dO3YsGzduZPDgwefWmNRLZs8OnDnTPbdGDGqKPfYPWBde7HZZ4uW8KtQVFRWxe/durrvuulLTe/bsyfbt28tdp7CwsNRgZ8uySj75nunT0IkTFjExp3bsZ97Bezf1VnWRQDegnMHya2rpJcvw5e0G33xjExTk1OlrPvrooyxcuJD33nuPmJgYevToQXR0NCEhIQQEBHD8+HGOHDnCjh072Lp1K/n5+QwfPpy4uLhKPf/cuXPp06cPPXv2LBXq0tLSyMrKolevXiXTAgIC6N69O9u3b68w1FV33/bj+b54RMgbenM+WYfz6hwoKoTIdviNT8SKOO+s63lDb+fCl/urq968KtRlZ2fjOA4hISGlpoeEhJCVlVXuOklJSSxdurTkcXR0NDNmzDjr1Zpzc8+5XBGppjZt2tC0ad2+Zrdu3fjLX/7CF198wdq1a3n33XcpKCgos1x4eDhDhw5l8ODBhIWFVeq5P/nkE/bs2cP06dPLzDu17ypvv3am+8pWd9/2Y2c7IunN6mNvpqiIrHkzOf72PwFo3O+ntJzwOHaTZlV6nvrYW03y5f5quzevCnWnlJd0K0q/8fHxpT5Jn1ruTFeMB8/XP998Y9OmTRuOHDnik7csUW/ex5d7gx/6y85OJTv7zP35+/vXym3C+vTpQ58+fSgqKmLv3r0cPXqUgoICmjdvTlRUVKXHz52SkZHBggULSExMPONtE0/fh51t+1Z333Zq2YiICFJTU33u/6P62pvJOUbxc3+BHf8DwP7FzRTG3ciRYzlwLKdSz1Ffe6spvtxfVXo7l32bV4W64OBgbNsuc1Tu2LFjZT7lnhIQEEBAQPkXkT3bLzYoyKFpU89/ffF/MPXmfXy5N/ihv+xs43p//v7+xMTEnPPz7N69m2PHjvHwww+XTHMch23btvHee+8xc+ZMwHPE7sdH/rKzsyvcr8G57dt+vJzbv+faUp96M/u+wZkzDTLToVEQ9pgHsHpf6plXjRrrU2+1wZf7q+3evCrU+fv706lTJ5KTk7nkkktKpicnJ/OTn/zExcpERMp34YUX8re//a3UtOeee462bdsyYsQI2rRpQ2hoKMnJyURHRwOe8cMpKSn85je/caNkqUHOvz/EvDILCgsgvC32+Eew2rZ3uyzxUV4V6gDi4uJ49tln6dSpE7Gxsaxbt46MjAydISYiNerw4cOsXbuWgwcPlhlbZ1kWf/rTnyr1PEFBQbRvX/qPeKNGjWjevHnJ9OHDh5OUlERkZCQREREkJSXRqFEj+vfvXzPNSJ0zxcWYZQsxa97yTLjwYuwxv8eq4vg5karwulB3+eWXk5OTw5tvvsnRo0dp164dEydOrJWxNSLSMO3fv5/ExERatGhBamoqHTp0ICcnh8zMTFq2bEmbNm1q9PVGjBhBQUEBc+fOJTc3l5iYGBITE3WNOi9ljmfjvPgEbNsKgDV8JNaIm7BsP5crE1/ndaEOYOjQoQwdOtTtMkTER73++uv06tWLBx54gJtvvpm7776bTp068d///pfnnnuOX//61+f0/JMmTSr12LIsRo4cqXvK+gDz7R6c2dMg4wgENvLcHeJiHXGVumG7XYCISH2zZ88eBg4cWHJG6amBzRdddBE///nPWbRokZvlST1lPtuIM/1BT6Br1QZ74hMKdFKnFOpERE6Tm5tLs2bNsG0bPz8/cn904cpOnTqxZ88eF6uT+sY4xTjLXsF54a9QcBK698Z+9CmsqI5ulyYNjFd+/SoiUptatGhBdnY24LlYaEpKCj179gQ84+0aN27sZnlSj5i84zgvPQn/+xwAa0g81vW3Yflp/JzUPYU6EZHTdOnShR07dnDJJZfQv39/lixZQlZWFv7+/nz44YdceeWVbpco9YA5tN8zfi7tEAQGYt12L3a/AW6XJQ2YQp2IyGmuv/56jh49CsB1111HVlYWGzduxLIsLrvsMm699VaXKxS3mS/+jTPvaTh5Alq09lx/rv35bpclDZxCnYjIaSIiIkru0WjbNnfeeSd33nmny1VJfWAcB/POPzErPPdvpcuF2Hc9iNW84rt/iNQVnSghInKaOXPmkJaWVu689PR05syZU8cVSX1gTuThzJlWEuisq3+Off9kBTqpNxTqRERO89FHH5WcKHG6nJwcPvroozquSNxmUg/iTJsAWzeDfwDWqPuwfz0Wy19feEn9of8bRUSq4Pjx4wQEBLhdhtQhk/wfnLlPwok8CG2JnfAIVnRnt8sSKUOhTkQESElJISUlpeTx+vXr2bJlS6llCgoK+M9//kNUVFQdVyduMMZgVi3BLH8NjIGY7tj3PIQVHOZ2aSLlUqgTEQG++uorli5dWvL4/fffL3e5Vq1aMXr06LoqS1xi8k/gvPwM/HcTANbAa7BuHIPlr6O0Un8p1ImIACNGjGDYsGEYYxg7diyJiYlER0eXWiYgIEAXHm4ATNphnDnT4OA+8PPHuvku7J/qfuNS/ynUiYgAgYGBBAYGAjBr1izCwsLw1yD4Bsd89QXOi09A3nEICcO++2GsmG5ulyVSKdpjiYicpnXr1gAcPHiQlJQUcnJyGDRoEKGhoWRmZtKsWbOSACi+wRiDWfMW5s2FYByIjsVOmIgV2tLt0kQqTaFOROQ0juPwwgsv8OGHH5ZM6927N6Ghobz44otER0dz4403uleg1Chz8iTmlWcxmz8GwOo/GOvmu7F0lrN4GV2nTkTkNMuWLWPjxo3ceuutPPnkk6Xm9enTp8xZseK9TMYRnBkPegKdn58nzN32WwU68Uo6UicicpoPP/yQX/7yl8TFxeE4Tql54eHhFd5tQryL+ToZ54UZcDwHmodg3/0QVuwFbpclUm0KdSIip8nMzCQ2NrbceQEBAeTn59dxRVKTjDGY9e9glswHx4EOMZ7xcy1au12ayDlRqBMROU1ISEiFR+MOHTpEixYt6rgiqSmm4CTmH3Mw//oAAOvSq7BuTcAKbORyZSLnTqFOROQ0ffr0YdmyZSUnRwBYlkVeXh7vvvsuffv2dbdAqRaTmY4zZzrs2wW2jfWrO7Gu/jmWZbldmkiNUKgTETnNyJEj+eKLL3jggQfo0aMHAK+//joHDhzAz8+PG264weUKparMjq9wnv8L5ByDZs2xxz2I1a2X22WJ1CiFOhGR04SGhjJ9+nTeeOMNvvjiC2zbZt++fVx00UXceOONNGvWzO0SpZKMMZiP3sX88yUoLoaoaM/4udYRbpcmUuMU6kREyhEaGsq4cePcLkPOgSksxLz+AmbDGgCsn1yJdfu9WI10qzfxTQp1IiLic0zWd57xc7u3g2VjXX8r1tDrNX5OfJpCnYhIOb7++ms2btxIeno6BQUFpeZZlsWf/vQnlyqTszm5LZniP0+AY5nQpCn22D9iXXCR22WJ1DqFOhGR03zwwQc8//zzNGvWjMjISAJOu7uAMcalyuRsnI9Xk7boeSgqgrbtscc/ghXe1u2yROqEQp2IyGnefvttLrvsMsaPH18m0En9ZIoKMYvnYT5cBYB10WVYd9yH1biJy5WJ1B2FOhGR06Snp3PHHXco0HkJk30U5/kZsDMFLIvgW+4i98prQOPnpIFRqBMROc15553HsWPHauS5kpKS2Lx5MwcPHiQwMJDY2FhuueUW2rb94StBYwxLlixh/fr1HD9+nM6dOzN69GjatWtXIzX4MrNnJ85z0+FoBgQ1wR7zB0KGjSDv8GF9TS4Nju12ASIi9c1NN93EW2+9RWZm5jk/V0pKCkOHDmXq1Kk8+uijOI7DlClTSt0/dvny5axcuZI777yT6dOnExoaypQpUzhx4sQ5v74vcza9j/PXhz2BLuI87Ef+ht3rErfLEnGNjtSJiAAzZswo9TgvL4/77ruPjh07lrnYsGVZPPjgg5V63sTExFKPExISGDNmDLt376Z79+4YY1i1ahXx8fH069cPgPHjxzN27Fg2btzI4MGDz6Er32SKijBLX8asf8czodcl2Hc+gNWkqbuFibhMoU5EBNi/f3+px7ZtExwcTGZmZo0csTslLy8PoCQopqWlkZWVRa9eP9yyKiAggO7du7N9+/YKQ11hYSGFhYUljy3LIigoqOTnMzk13xuv2WZyjmFe+Cvm62QArJ//GvvnN2HZni+evLm3s/Hl3sC3+6ur3hTqRESA2bNn1/prGGNYuHAhXbt2pX379gBkZWUBEBISUmrZkJAQMjIyKnyupKQkli5dWvI4OjqaGTNm0Lp160rXExHhXbfKKvhmOxnT/4hJO4wV1IQWv59Mk8uvKndZb+utKny5N/Dt/mq7N68JdWlpabz55pv873//IysrixYtWnDllVdy/fXX4+/vNW2IiBdISUmhU6dONG5c9nZS+fn5JV+dVtW8efPYv38/jz/+eJl5p3+CP9sg//j4eOLi4sqsn56eTlFR0RnXtSyLiIgIUlNTveZkAufTj3AW/h0KCiA8Ent8IsfO68Cxw4dLLeeNvVWWL/cGvt1fVXrz9/ev0oezUutWay0XHDp0CGMM48aNIyIiggMHDvDCCy+Qn5/Pbbfd5nZ5IuJDJk+ezNSpU4mJiSkz79ChQ0yePJnFixdX6Tnnz5/P559/zuTJk2nZsmXJ9NDQUMBzxC4sLKxkenZ2dpmjdz8WEBBQ4SVXKvsH0RhT7/94GqcYs+wVzOokz4QLLsIeMwGaNjtj7d7QW3X5cm/g2/3Vdm9eE+p69+5N7969Sx63adOGQ4cOsWbNGoU6EakzRUVF2HblLxxgjGH+/Pls3ryZSZMmER4eXmp+eHg4oaGhJCcnEx0dXfIaKSkp/OY3v6nR2r2Nyc3BefFvkPIFANY1v8S67hYs28/lykTqJ68JdeXJy8src1ba6RrqYOKzUW/eyZd7A3f7y8vLKzmJATxHzk4f01ZQUMBHH31UcnStMubNm8fGjRt58MEHCQoKKhlD16RJEwIDA7Esi+HDh5OUlERkZCQREREkJSXRqFEj+vfvXxOteSXz7V6cOdMgPRUCG2GNug/7Jw339yFSGV4b6lJTU3n33XfPepSuIQ4mrgr15p18uTdwp7+VK1eW2lc88cQTFS4bHx9f6edds2YNAJMmTSo1PSEhgYEDBwIwYsQICgoKmDt3Lrm5ucTExJCYmFjyAbShMZ9/gvPyM3AyH1qGY49PxGoX7XZZIvWeZVz+4vqNN94otSMtz/Tp0zn//PNLHmdmZjJp0iS6d+/O3XfffcZ1KzpS56uDiStLvXknX+4N6m4wcXl27NjB9u3bMcbw2muvMWzYMFq1alVqmYCAANq3b1+tkyTqSnp6eql9XnksyyIyMpLD9eyuC8YpxixfhFm1xDOhWy/scX/EahZc6eeor73VBF/uDXy7v6r0FhAQ4L0nSgwbNowrrrjijMv8uLnMzEwmT55MbGws48aNO+vzN5TBxNWl3ryTL/cG7vQXGxtLbGwsACdPnuTqq6+mRYsWdVpDQ2byjuPMfQq+/AwAa8h1WNffjuWn8XMileV6qAsODiY4uHKfwk4FuujoaBISEqo0WFlEpLJ+9atfuV1Cg2IOH8CZPQ2OHISAQKzbxmNfWv7150SkYq6Huso69ZVrq1atuO2228jOzi6ZV5VByyIiUn+YLZ/izHsK8k9Ai1bYCYlYHc4/+4oiUobXhLrk5GRSU1NJTU0tM47ujTfecKkqERGpDuM4mBWLMe+87pkQewH2XQ9iBYe6WpeIN/OaUDdw4MCSM8VERMR7mRN5OPOfhi2fAmANisP61Z1YujuQyDnRO0hEROqMST3ouf7c4QPg7491SwL2FT9zuywRn6BQJyIidcJ8+RnOS0/CiVwIbYF9z0SsTl3cLkvEZyjUiYhIrTLGYN5dinnrH2AMnN/VE+hCws6+sohUmkKdiIjUGpN/AmfBM/D5JgCsnw7Dumksln/51w8VkepTqBMRkVph0lNxZk+Fg/vAzx/rpnHYA4a5XZaIz1KoExGRGmdSvsB54QnIOw7Bodj3PIwVU39vrybiCxTqRESkxhhjMGvfwixdCMaB6FjP+Lmwlm6XJuLzFOpERKRGmJMnMa/Mwmz+CADriquxfnMPVkCgy5WJNAwKdSIics7Md2me68/t3w1+flg3jsEaOBzLstwuTaTBUKgTEZFzYrZ/ifP8DDieDc2Cse9+GKvLBW6XJdLgKNSJiEi1GGMw76/AvDEPHAfad8JOSMRq2drt0kQaJIU6ERGpMlNYgHl1DuZf7wNgXToQ69bxWIGNXK5MpOFSqBMRkSoxmRk4z02HvTvBsrF+dQfWz36h8XMiLlOoExGRSjM7UzyBLucYNG2OPe6PWN17u12WiKBQJyIileR8+C7mny9CcTFEdcROeASrdYTbZYnI9xTqRETkjExhIeb1FzAb1gBgXdwfa9TvsBo1drkyEfkxhToREamQycrEef4v8M3XYFlY8bdhDbte4+dE6iGFOhERKZf55muc5/4CxzKhSVPssROwLujrdlkiUgGFOhERKcPZuBbz2nNQVASR7bDHJ2K1aet2WSJyBgp1IiJSwhQVYd6Yi/lglWdC70uxR9+P1biJu4WJyFkp1ImICAAmO8szfm5nCgDWiJuxho/Esm2XKxORylCoExERzL5dOHOmQWYGNA7CHvMHrF6XuF2WiFSBQp2ISAPn/OsDzKuzobAA2pznGT8XGeV2WSJSRQp1IiL1xOrVq3n77bfJysoiKiqKUaNG0a1bt1p7PVNcjLNkPmbd254JF17sOULXpGmtvaaI1B4NlBARqQc2bdrEggULuP7665kxYwbdunVj2rRpZGRk1MrrFR/Lwnn6TyWBzrp2JPZvH1WgE/FiCnUiIvXAihUrGDRoEFdffXXJUbpWrVqxZs2aGn8ts383R+6/FfN1MjRqjH3Pw9jX3aITIkS8nL5+FRFxWVFREbt37+a6664rNb1nz55s37693HUKCwspLCwseWxZFkFBQSU/V8Ts+IrimY9BwUkIj8RvfCLWeR3OvYl64lTvvnjHC1/uDXy7v7rqTaFORMRl2dnZOI5DSEhIqekhISFkZWWVu05SUhJLly4teRwdHc2MGTNo3br1GV/LCQ0hLTIKv5bhtHxwKnbz4HOuvz6KiIhwu4Ra48u9gW/3V9u9KdSJiNQT5X2Kr+iTfXx8PHFxcWWWS09Pp6io6MwvdN8kWsXEciQ9HXM8t/oF10OWZREREUFqairGGLfLqVG+3Bv4dn9V6c3f3/+sH84qXLdaa4mISI0JDg7Gtu0yR+WOHTtW5ujdKQEBAQQEBJQ772x/NKzgUCw/P4wxPvfH8xT15r18ub/a7k2jYkVEXObv70+nTp1ITk4uNT05OZkuXbq4VJWIeBsdqRMRqQfi4uJ49tln6dSpE7Gxsaxbt46MjAwGDx7sdmki4iUU6kRE6oHLL7+cnJwc3nzzTY4ePUq7du2YOHFitcfWiEjDo1AnIlJPDB06lKFDh7pdhoh4KY2pExEREfEBDfZInb9/5VuvyrLeRr15J1/uDSrXn6//DqpL+zYP9ea9fLm/2t63WcZXzxsWERERaUD09esZnDhxgoceeogTJ064XUqNU2/eyZd7A9/vr77w5d+zevNevtxfXfWmUHcGxhj27NnjkxdBVG/eyZd7A9/vr77w5d+zevNevtxfXfWmUCciIiLiAxTqRERERHyAQt0ZBAQEcMMNN1R4f0Vvpt68ky/3Br7fX33hy79n9ea9fLm/uupNZ7+KiIiI+AAdqRMRERHxAQp1IiIiIj5AoU5ERETEByjUiYiIiPgA373BWiWtXr2at99+m6ysLKKiohg1ahTdunWrcPmUlBQWLlzIt99+S1hYGL/4xS8YMmRIHVZ8dklJSWzevJmDBw8SGBhIbGwst9xyC23btq1wna+++orJkyeXmf70009z3nnn1Wa5VfLGG2+wdOnSUtNCQkJ46aWXKlzHG7YZwPjx40lPTy8zfciQIYwZM6bM9Pq+zVJSUnj77bfZs2cPR48eZcKECVxyySUl840xLFmyhPXr13P8+HE6d+7M6NGjadeu3Rmf99///jeLFy/myJEjtGnThptuuqnU84qH9m0e9f19coov79vAt/Zv9Xnf1qBD3aZNm1iwYAFjxoyhS5curFu3jmnTpvH000/TqlWrMsunpaUxffp0rr76au699162b9/O3LlzCQ4O5tJLL3Whg/KlpKQwdOhQzj//fIqLi/nnP//JlClTeOqpp2jcuPEZ1505cyZNmjQpeRwcHFzb5VZZu3bt+L//+7+Sx7Zd8QFnb9lmANOnT8dxnJLH+/fvZ8qUKVx22WVnXK++brOTJ0/SsWNHrrrqKp588sky85cvX87KlStJSEggMjKSZcuWMWXKFGbOnElQUFC5z7ljxw5mzpzJjTfeyCWXXMLmzZt5+umnefzxx+ncuXNtt+Q1tG8rq76+T37MV/dt4Fv7t/q8b2vQoW7FihUMGjSIq6++GoBRo0axdetW1qxZw80331xm+TVr1tCqVStGjRoFQFRUFN988w3vvPNOvXoTJSYmlnqckJDAmDFj2L17N927dz/juiEhITRt2rQ2yztntm0TGhpaqWW9ZZtB2Z3VW2+9RZs2bbx2m/Xp04c+ffqUO88Yw6pVq4iPj6dfv36A55P82LFj2bhxI4MHDy53vZUrV9KzZ0/i4+MBiI+PJyUlhZUrV3L//ffXSh/eSPu2surr++THfHXfBr61f6vP+7YGG+qKiorYvXs31113XanpPXv2ZPv27eWus3PnTnr27FlqWu/evfnggw8oKirC379+/jrz8vIAaNas2VmXffDBByksLCQqKorrr7+eCy64oLbLq7LU1FTuuusu/P396dy5MzfddBNt2rQpd1lv3WZFRUVs2LCBa6+9FsuyzrisN2yz06WlpZGVlUWvXr1KpgUEBNC9e3e2b99e4Y5vx44dXHvttaWm9erVi1WrVtVqvd5E+7byecP7pCHs28C3929u79vq71avZdnZ2TiOQ0hISKnpISEhZGVllbtOVlZWucsXFxeTk5NDWFhYbZVbbcYYFi5cSNeuXWnfvn2Fy4WFhTFu3Dg6depEUVERH3/8MX/+85957LHHzvpJqi517tyZ8ePH07ZtW7Kysli2bBmPPvooTz31FM2bNy+zvDduM4DNmzeTm5vLwIEDK1zGW7ZZeU69x8rbNhkZGWdc7/QjGaGhoRW+Zxsi7dtK85b3SUPZt4Fv79/c3rc12FB3SnmfEs70yeH0eaduyHG2TxtumTdvHvv37+fxxx8/43Jt27YtNdg4NjaWjIwM3nnnnXr1BvrxIe/27dsTGxvLvffey0cffURcXFy563jbNgP44IMP6N27Ny1atKhwGW/ZZmdS0bapCmNMvd6WbtG+zcNb3icNZd8GDWP/5ta+rcFe0iQ4OBjbtsuk4GPHjpVJ2KeUl5qzs7Px8/Or1OH/ujZ//nw+//xzHnvsMVq2bFnl9WNjY0lNTa2FympO48aNad++PYcPHy53vrdtM4D09HSSk5NLxkNVhTdsM6DkE2l526ai99+p9arynm2ItG87O294n/jivg18f//m9r6twYY6f39/OnXqRHJycqnpycnJdOnSpdx1OnfuXGb5rVu30qlTp3o1fsEYw7x58/j000/505/+RHh4eLWeZ8+ePZUetOuWwsJCDh48WOFXDd6yzX7sgw8+ICQkhIsuuqjK63rDNgMIDw8nNDS01LYpKioiJSWlwvcfeHbqX375ZalpycnJxMbG1lqt3kb7trPzhveJL+7bwPf3b27v2xpsqAOIi4tj/fr1vP/++3z77bcsWLCAjIyMkoGMixYtYtasWSXLDxkyhIyMjJLrAr3//vu8//77/PznP3erhXLNmzePDRs2cN999xEUFERWVhZZWVkUFBSULHN6bytXrmTz5s0cPnyYAwcOsGjRIj799FOGDRvmRgsVeuWVV0hJSSEtLY2dO3fy5JNPcuLECQYMGAB47zY7xXEcPvzwQwYMGICfn1+ped62zfLz89m7dy979+4FPAOI9+7dS0ZGBpZlMXz48JLrju3fv5/Zs2fTqFEj+vfvX/Ics2bNYtGiRSWPhw8fztatW3nrrbc4ePAgb731Fl9++WWZAcYNnfZt3vM+OcXX923gO/u3+rxvq79xvg5cfvnl5OTk8Oabb3L06FHatWvHxIkTad26NQBHjx4tNbAxPDyciRMnsnDhQlavXk1YWBh33HFHvTt9fM2aNQBMmjSp1PSEhISSgamn91ZUVMSrr75KZmYmgYGBtGvXjocffrhan6ZqU2ZmJs888wzZ2dkEBwfTuXNnpk6d6vXb7JQvv/ySjIwMrrrqqjLzvG2bffPNN6UuHvrKK68AMGDAAMaPH8+IESMoKChg7ty55ObmEhMTQ2JiYqnrOJ3aSZ7SpUsX7r//fv75z3+yePFiIiIiuP/++3WNutNo3+Y975NTfH3fBr6zf6vP+zbLVGf0noiIiIjUKw3661cRERERX6FQJyIiIuIDFOpEREREfIBCnYiIiIgPUKgTERER8QEKdSIiIiI+QKFORERExAco1ImIiIj4AIU6ERERER+gUCciIiLiAxTqRERERHyAQp34hIKCAh588EHuvfde8vLySqZnZWUxduxYJk2ahOM4LlYoIiJSuxTqxCcEBgbywAMPkJ2dzZw5cwBwHIe///3vANx3333Ytv53FxER36W/cuIzIiMjueuuu9i8eTOrVq1i6dKlfPXVV9x7772EhYW5XZ6IiEit8ne7AJGadPnll5OSksKrr76K4zjEx8fTs2dPt8sSERGpdTpSJz7nqquuori4GD8/P4YPH+52OSIiInVCoU58Sn5+PrNmzSIyMpLAwECef/55t0sSERGpEwp14lNeeuklMjIymDBhAnfffTefffYZK1ascLssERGRWqdQJz5j/fr1bNiwgdGjR9OuXTsuvfRShg0bxmuvvcauXbvcLk9ERKRWKdSJT9i/fz8vv/wyAwYMYODAgSXTb731Vjp06MDTTz9Nbm6uewWKiIjUMssYY9wuQkRERETOjY7UiYiIiPgAhToRERERH6BQJyIiIuIDFOpEREREfIBCnYiIiIgPUKgTERER8QEKdSIiIiI+QKFORERExAco1ImIiIj4AIU6ERERER+gUCciIiLiA/4fPu/8mGRyMxYAAAAASUVORK5CYII=" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot trajectory\n", "plt.subplot(2, 2, 1)\n", "plt.plot(x_bar[0, :], x_bar[1, :])\n", "plt.plot(np.linspace(0, 10, T + 1), np.zeros(T + 1), \"b-\")\n", "plt.axis(\"equal\")\n", "plt.ylabel(\"y\")\n", "plt.xlabel(\"x\")\n", "\n", "plt.subplot(2, 2, 2)\n", "plt.plot(np.degrees(x_bar[2, :]))\n", "plt.ylabel(\"theta(t) [deg]\")\n", "# plt.xlabel('time')\n", "\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## State Space Linearized Model" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "ExecuteTime": { "end_time": "2024-10-22T09:43:55.415384Z", "start_time": "2024-10-22T09:43:55.412871Z" } }, "outputs": [ { "data": { "text/plain": "'\\n控制问题描述\\n'" }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Control problem statement.\n", "\"\"\"\n", "\n", "N = 4 # number of state variables\n", "M = 2 # number of control variables\n", "T = 20 # Prediction Horizon\n", "DT = 0.2 # discretization step" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "ExecuteTime": { "end_time": "2024-10-22T09:43:55.419986Z", "start_time": "2024-10-22T09:43:55.419039Z" } }, "outputs": [], "source": [ "def get_linear_model(x_bar, u_bar):\n", " \"\"\"\n", " Computes the LTI approximated state space model x' = Ax + Bu + C\n", " \"\"\"\n", "\n", " L = 0.3 # vehicle wheelbase\n", "\n", " x = x_bar[0]\n", " y = x_bar[1]\n", " v = x_bar[2]\n", " theta = x_bar[3]\n", "\n", " a = u_bar[0]\n", " delta = u_bar[1]\n", "\n", " A = np.zeros((N, N))\n", " A[0, 2] = np.cos(theta)\n", " A[0, 3] = -v * np.sin(theta)\n", " A[1, 2] = np.sin(theta)\n", " A[1, 3] = v * np.cos(theta)\n", " A[3, 2] = v * np.tan(delta) / L\n", " A_lin = np.eye(N) + DT * A\n", "\n", " B = np.zeros((N, M))\n", " B[2, 0] = 1\n", " B[3, 1] = v / (L * np.cos(delta) ** 2)\n", " B_lin = DT * B\n", "\n", " f_xu = np.array(\n", " [v * np.cos(theta), v * np.sin(theta), a, v * np.tan(delta) / L]\n", " ).reshape(N, 1)\n", " C_lin = DT * (\n", " f_xu - np.dot(A, x_bar.reshape(N, 1)) - np.dot(B, u_bar.reshape(M, 1))\n", " )\n", "\n", " return np.round(A_lin, 4), np.round(B_lin, 4), np.round(C_lin, 4)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "ExecuteTime": { "end_time": "2024-10-22T09:43:55.424811Z", "start_time": "2024-10-22T09:43:55.422385Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.81 ms, sys: 361 μs, total: 2.17 ms\n", "Wall time: 576 μs\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:17: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n", " A[0, 2] = np.cos(theta)\n", "/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:18: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n", " A[0, 3] = -v * np.sin(theta)\n", "/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:19: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n", " A[1, 2] = np.sin(theta)\n", "/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:20: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n", " A[1, 3] = v * np.cos(theta)\n", "/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:21: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n", " A[3, 2] = v * np.tan(delta) / L\n", "/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:26: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n", " B[3, 1] = v / (L * np.cos(delta) ** 2)\n" ] } ], "source": [ "%%time\n", "\n", "u_bar = np.zeros((M, T))\n", "u_bar[0, :] = 0.2 # m/s\n", "u_bar[1, :] = np.radians(-np.pi / 4) # rad\n", "\n", "x0 = np.zeros(N)\n", "x0[0] = 0\n", "x0[1] = 1\n", "x0[2] = 0\n", "x0[3] = np.radians(0)\n", "\n", "x_bar = np.zeros((N, T + 1))\n", "x_bar[:, 0] = x0\n", "\n", "for t in range(1, T + 1):\n", " xt = x_bar[:, t - 1].reshape(N, 1)\n", " ut = u_bar[:, t - 1].reshape(M, 1)\n", "\n", " A, B, C = get_linear_model(xt, ut)\n", "\n", " xt_plus_one = np.dot(A, xt) + np.dot(B, ut) + C\n", "\n", " x_bar[:, t] = np.squeeze(xt_plus_one)" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "ExecuteTime": { "end_time": "2024-10-22T09:43:55.586059Z", "start_time": "2024-10-22T09:43:55.425290Z" } }, "outputs": [ { "data": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAEDCAYAAABTZPIVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8QElEQVR4nO3deXQV9f3/8edMNsKSDQgJJEAQgqKyaCsuWFDKUqTFWIt1R1HUUKv+qlZNF/CLUOpGK+AGCi64gASsoGyKlWKlWhVrFERWgZDEcAkkhORmPr8/LlyNSSCBJJM7eT3O4SR37sy973eGTF535jMzljHGICIiIiIhzXa7ABERERE5cQp1IiIiIh6gUCciIiLiAQp1IiIiIh6gUCciIiLiAQp1IiIiIh6gUCciIiLiAQp1IiIiIh6gUCciIiLiAeFuF+CWvXv34vf7jzlf+/btyc/Pb4SKGp96C01e7g1q3194eDjx8fGNUFFo0bZNvYUyL/fXGNu2Zhvq/H4/5eXlR53HsqzgvF67m5p6C01e7g28319j0LZNvYUqL/fXWL3p8KuIiIiIByjUiYiIiHiAQp2IiIiIByjUiYiIiHiAQp2IiIhIAzKlJY1y8odCnYiIiEgDMd9spWLibezPfrHB30uhTkRERKQBmI/+hfOXuyE/l+I3X8OUHWrQ92u216kTERERaQjGqcAsnodZOh8A65Q+JP7pYfKKDzboYViFOhEREZF6YkoO4Mx6BD77EABr6MXYvxxDWEwcFB9s0PdWqBMRERGpB2b3DpwZk2HPToiIxLpmPPbZFwTvKNHQFOpERERETpD55AOc2Y9A6UFIaIedmYXV5aRGrUGhTkREROQ4GcfBvPEK5h8vBSakn4Z9091YMXGNXotCnYiIiMhxMAdLcJ55FD75AADrwpFYv7oeK9ydeKVQJyIiIlJHJncnzszJsHsHhIdjXZWJfd5PXa0ppK9Tl52dzejRo5kzZ47bpYiIiEgzYT77EGfynYFAF5eAfdcU1wMdhPCeuk2bNrFy5Uq6dOnidikiIiLSDBhjMG8uwCx6AYyBk07GvuVerNh4t0sDQnRPXWlpKY899hg33XQTrVq1crscERER8ThTehDnyamY7OfBGKyfDMe+84EmE+ggRPfUzZo1i379+tG7d28WLlx41HnLy8spLy8PPrYsi+jo6OD3R3Pk+ca6vkxjUm+hycu9gff7E5HQZPJzcWY8ADu3QVg41uXjsAcOd7usKkIu1P3rX/9iy5YtTJkypVbzZ2dns2DBguDjtLQ0pk6dSvv27Wv9nklJSXWuM1Sot9Dk5d7A+/2JSOgwOR/jPPkglByAmDjsW+7B6t7L7bKqFVKhrqCggDlz5pCVlUVkZGStlsnIyGDkyJHBx0f2AOTn5+P3+4+6rGVZJCUlkZub26D3anODegtNXu4N6tZfeHh4nT6ciYjUhTEGs3wR5rW5YBxISw+Mn4tv63ZpNQqpULd582b27dvHPffcE5zmOA5ffPEFb731FvPmzcO2Kw8TjIiIICIiotrXq+0fRWOMJ/+AgnoLVV7uDbzfn4g0bebQIcxz0zHr3gXAOm8w1pW3YEXUboeSW0Iq1J1++uk89NBDlaY9/vjjdOzYkVGjRlUJdCIiIiJ1Yb7NC1x/bvtmCAvDuuwGrEEjQmKsb0iFuujoaDp37lxpWlRUFG3atKkyXURERKQuzIbPcJ6YCgeKoHUM9s33YPU8ze2yai2kQp2IiIhIfTPGYN5egnl1FjgOdO6GnZmF1Ta0xu2GfKibMGGC2yWIiIhIiDLlZZjnZ2LefxsAq/9ArKt/gxUV5XJldRfyoU5EpCmrqKhg/vz5vPfee/h8PuLj4xk0aBCXXHJJcBywMYb58+ezatUqDhw4QI8ePRg7diypqakuVy/ibaawAOfxKbD1K7BsrEvHYA0ZFRLj56qjUCci0oAWL17MihUrGD9+PCkpKWzevJmZM2fSsmVLRowYEZxnyZIlZGZmkpyczMKFC5k0aRLTpk0LXixdROqX+SonEOj274NWbbDH3YXVq6/bZZ0QnS4qItKANm7cyI9+9CPOOOMMEhMTOfvss+nduzdff/01ENhLt3TpUjIyMujfvz+dO3dm/PjxHDp0iDVr1rhcvYg3OavfxHk4KxDoUrpiZz0c8oEOtKdORKRBnXzyyaxYsYJdu3bRsWNHtm7dyoYNG7j22msByMvLw+fz0adPn+AyERER9OrViw0bNjBkyJBqX1e3QKyeegtdjdGfKS/HeelJzD+XBd7rRwOwr7sNK6pFg70nNN66U6gTEWlAo0aNoqSkhDvuuAPbtnEch1//+tcMGDAAAJ/PB0BsbGyl5WJjYykoKKjxdXULxKNTb6GrofqrKCyg4KEsKr5cD5ZF7LXjaXPptY0akht63SnUiYg0oLVr1/Lee+/x29/+ltTUVLZu3cqcOXOCJ0wc8cM/LMe6o4ZugVg99Ra6GrI/8/WXVMycAvsKoWUr7HF3UXzamRTn5tbr+9SksW6BqFAnItKAXnjhBUaNGsV5550HQOfOncnPz2fRokUMGjSIuLg4gOCZsUcUFRVV2Xv3fboF4tGpt9BV3/05a1ZgXnwc/H5ITsUen4XVoaMrP8OGXnc6UUJEpAEdOnSoyi0MbdsObtgTExOJi4tj/fr1wef9fj85OTn07NmzUWsV8RLj9+PMewIz97FAoOt7NvZ9D2J16Oh2aQ1Ge+pERBrQmWeeycKFC2nXrh0pKSls3bqVN954gwsuuAAIHJYZMWIE2dnZJCcnk5SURHZ2NlFRUcFxdyJSN6bIh/PEX+CrHACsUVdgjRiN5fF7xCvUiYg0oOuvv55XXnmFWbNmsW/fPhISEhgyZAiXXnppcJ5Ro0ZRVlbGrFmzKC4upnv37mRlZekadSLHwWzbhDNjMuwtgBbR2GP/H1bf/m6X1SgU6kREGlB0dDRjxoxhzJgxNc5jWRajR49m9OjRjVeYiAc577+DeX4GlJdBh06B8XPJKW6X1WgU6kRERCSkmYoKzII5mJWLAxNO/xH2Db/DatnK3cIamUKdiIiIhCyzvwjnqb/Cl4GTjayLRmP94grPj5+rjkKdiIiIhCSzfTPOzMnwbR5EtcC+/nasM851uyzXKNSJiIhIyHHW/RMz9+9QVgbtkwLj5zp1cbssVynUiYiISMgwTgVm4fOYZQsDE07th33jXVitWrtbWBOgUCciIiIhwRTvx3n6Ifj8YwCs4b/EyrgKyw5zubKmQaFOREREmjyzcxvOjAcgPxciI7Gu/S32WT9xu6wmRaFOREREmjTz37U4z0yDQ6XQNjEwfi41ze2ymhyFOhEREWmSjONgXp+HWfJqYMLJvbHH3Y3VJsbdwpoohToRERFpckxJMc7sR2D9fwCwfjoK69IxWGEaP1cThToRERFpUszub3BmPgC5OyEiEuua8dhnX+B2WU2eQp2IiIg0Gc4n63BmPQSlByGhHXbmfVhdurtdVkhQqBMRERHXGcdh30uzcF54IjAh/VTsm36PFRPnal2hRKFOREREXGVKSzDP/I2ij98HwLrgIqzRY7HCFVPqQj8tERERcY3J24Uz/QHYvQPCI7CvugXrvJ+6XVZIUqgTERERV5j/fRS4Q0RJMcQlkPjHRyiMbYsxxu3SQpJCnYiIiDQqYwzmrdcw2c+DMXDSyYTdci9RJ58Ku3e7XV7IUqgTERGRRmMOlWLm/B3z4RoArJ8Mw/r1OKzISJcrC30KdSIiItIoTH4uzszJ8M1WCAvHunwc9sDhbpflGQp1IiIi0uDMF5/iPPlXKN4PMXHYN9+D1aOX22V5ikKdiIiINBhjDGbFYsyCOWAcSEsPBLqEdm6X5jkKdSIiItIgTNkhzHPTMR+8C4B17mCsq27BitD4uYYQUqEuOzubdevWsXPnTiIjI0lPT+eqq66iY8eObpcmIiIi32O+zQ+Mn9v+Ndg21mU3BC4qbFlul+ZZIRXqcnJyGDZsGCeddBIVFRW8/PLLTJo0iUceeYQWLVq4XZ6IiIgAZsP/cJ6cCvv3QesY7Jt/j9XzdLfL8ryQCnVZWVmVHmdmZnLDDTewefNmevXSYEsRERE3GWMw7yzBvDobKiqgczfszPuw2ia6XVqzEFKh7odKSkoAaN26dY3zlJeXU15eHnxsWRbR0dHB74/myPNe3FWs3kKTl3sD7/cn4mWmvAzz4uOYf60CwOo/EOvq32BFRblcWfMRsqHOGMPcuXM5+eST6dy5c43zZWdns2DBguDjtLQ0pk6dSvv27Wv9XklJSSdUa1Om3kKTl3sD7/cn4jVm77c4j0+BLRvBsrEuHYM1ZJQ+oDWykA11s2fPZvv27dx///1HnS8jI4ORI0cGHx/5D5afn4/f7z/qspZlkZSURG5urufuQ6feQpOXe4O69RceHl6nD2ci0jDMphycx/8CRT5o1QZ73F1Yvfq6XVazFJKh7plnnuGjjz5i4sSJtG3b9qjzRkREEBERUe1ztf2jaIzx5B9QUG+hysu9gff6Kyws5IUXXuCTTz6hrKyM5ORkbrnlFrp16wYE+p0/fz6rVq3iwIED9OjRg7Fjx5Kamupy5SJH5/zzLcy8p6DCD526YI/PwmqvPe1uCalQZ4zhmWeeYd26dUyYMIHERA28FJGm7cCBA/zxj3/k1FNP5b777iMmJoY9e/bQsmXL4DyLFy9myZIlZGZmkpyczMKFC5k0aRLTpk0LjgEWaUqMvxzz0tOYf74FgHXmeVjX3YYVpStRuCmkQt3s2bNZs2YNd999N9HR0fh8PgBatmxJpG4ELCJN0OLFi2nbti2ZmZnBad//QGqMYenSpWRkZNC/f38Axo8fz4033siaNWsYMmRIo9cscjTGV4jzxF/g6y/BsrAyrsYa/kuNn2sCQirULV++HIAJEyZUmp6ZmcmgQYMavyAR8bxdu3ZRWFhIWVkZMTExdOzYsdJetmP58MMP6dOnD4888gg5OTkkJCQwdOhQfvrTnwKQl5eHz+ejT58+wWUiIiLo1asXGzZsqDHU6cz+6qm3hmU2b8SZ+QD4CiG6Ffa4O7FP/1G9vHZT6K+hNFZvIRXqXn31VbdLEJFmYOPGjaxYsYJPPvmEoqKiSs/Ztk3Xrl05//zzGTRo0DEDXl5eHitWrOCiiy4iIyODTZs28eyzzxIREcHAgQODRxxiY2MrLRcbG0tBQUGNr6sz+49OvdW/AyteZ+/0KeAvJzw1jXZ/fJiITjVffeJ4ad0dv5AKdSIiDWnr1q3MmTOHL774gk6dOtG/f3+6detGTEwMkZGRHDhwgD179vDVV1/x0ksv8corrwTPsA8Pr35z6jgOJ510EldccQUQCF87duxg+fLlDBw4MDjfDz/BH+tEEZ3ZXz31Vv+M34/z6mzM228E6uh3NmbsHRTYEbB7d729j9ZdwImc2a9QJyJy2H333ceAAQO45pprgmem1qS0tJS1a9eyePFiKioq+OUvf1ntfPHx8aSkpFSalpKSwgcffABAXFwcAD6fj/j4+OA8RUVFVfbefZ/O7D869VZP71XkC9zua+PnAFi/uALrotFg2w1Wg9bd8VOoExE57OGHHyY5OblW87Zo0YILL7yQQYMGHfUwac+ePdm1a1elabt27Qp+Ek9MTCQuLo7169eTlpYGgN/vJycnhyuvvPI4OxE5cWbbJpyZk6GwAFpEY4/9f1h9+7tdlhyF7XYBIiJNRW0D3ffZtn3UyytddNFFfPXVVyxcuJDc3FzWrFnDqlWrGDZsGBA4LDNixAiys7NZt24d27dvZ8aMGURFRTFgwIDj7kXkRDj/Xo0z9Z5AoEvsiH3fQwp0IUB76kREqvGb3/yGO++8k65du1Z5bvv27fz1r39l+vTpx3yd7t27c+eddzJv3jxee+01EhMTufbaazn//POD84waNYqysjJmzZpFcXEx3bt3JysrS9eok0ZnKiowr83BrFgcmHD6j7Bv+H9YLWu+x7o0HQp1IiLVONoJB+Xl5eTn59f6tc4880zOPPPMGp+3LIvRo0czevToOtcpUl/MgSKcpx6ELz4FwBoxGmvU5Vh2mMuVSW0p1ImI1NGePXu0F008xezYgjPjAfg2D6JaYF93G9aZ57ldltSRQp2IyGGrV6/m3XffDT6eNWtWlfBWVlbGtm3b6NWrV2OXJ9IgnP+swcz5G5QdgvZJ2Jn3YaV0dbssOQ4KdSIih5WVlVW62HBxcXGluzZA4FIi5557rg6VSsgzTgVm0QuYN18LTOjVF3vcXVit2rhbmBw3hToRkcOGDh3K0KFDgcD9V3/3u99Ve6KESKgzxQdwZj0E//svANbQDKxLrsEK0/i5UKZQJyJSjRkzZrhdgkiDMDu348yYBPm5EBmJdc2t2P0HHntBafIU6kREDisoKKBdu3Z1Xq6wsJCEhIQGqEikfpn/vo/zzDQ4dBDaJmJn3ovV+SS3y5J6oosPi4gcdtttt/Hss8+Sm5t7zHn9fj/vv/8+d911F2+//XYjVCdy/Izj4Cyeh/P4lECg63k6dtYjCnQeoz11NTBOBfj9OKWlmEOlGANYgGUF/mGBbYFlV7kRt4iEpj/84Q/MnTuXt956i+7du3PqqaeSlpZGbGwsERERHDhwgD179rBx40Y+/fRTSktLGTFiBCNHjnS7dJEamZJinGcehU/XAWAN/jnWpddhhSsCeI3WaE02b6Ri6u/ZWZt5bRvsMAgLhzA78NUO++778AiIiITIyMNfo7AiIr83Leq751q0gOhWWNEtIboV/OCrVcMNvEXkxJ1yyin85S9/4eOPP2bFihW8+eablJWVVZkvMTGRYcOGMWTIEOLj412oVKR2TO43gevP5e6E8Aisq8djn3uh22VJA1Goqw+OE/jnLz/2vIeZ430+PKJy0GsZ+Gq1bA1tYiEmDmLisNrEQkw8xMRCqzZYto60i9RWv3796NevH36/n61bt7J3717Kyspo06YNKSkpGj8nIcF8+h+c2Q/DwRKIbxcYP9e1h9tlSQNSqKtJt3TCZswnKSmJ3NxcjDHww3+YQJirqAj8cyqgwg8VTuCrc3h6eVngX1kZpvwQlH33OPhceVngwo+HSjEHS+BgceAXsaQ48H3pwUBd/nLYvy/w73t+GAIrPbZsaBMTCHxtYrEOB7+iTqk4hAWCX5tYiEuAmHgFQJHDwsPD6d69u9tliNSJcRzM0vmY1+cF/lZ174V9y++xYrRX2esU6mpg2WFYUeHYLaKxolocDnH18LrHuZxxKgLB7kjgKymBgyWYI+GveH8g6BX5MIe/st8HB/aDcQKPi3yB1zr8mvuqe6PwcEhoD20TsdomQttEaJeIlRD4SlyC7gMoItJEmdISnGf/Bv99HwBr0Aisy8ZihWvoTnOgUBciLDsMWrYO/Pv+9GMsZ/x+OFB0OOTtwxwJe0X7iPaXcXDPLkzR4RBYtBf8fsjbDXm7K+3tC34fFgbx7aqGviPfx7fTxSvFM3bv3s2KFSvYuXNnlbF1lmXxpz/9yaXKRKoyebtwZkyGXdshPBzripuxzx/qdlnSiBTqPM4KDw8cVo0LjAE6EgIty6JtcjK7d+8OHFoGTEUF+L6Fb/MwBXmBGzt/m4c5/JXCgsBh5YI9ULCn+tBn2RDfNhD0OnSCjqlYyZ0hORXi2+pMYQkZ27dvJysri4SEBHJzc+nSpQv79++nsLCQtm3b0qFDB7dLFAky//svztMPBobsxCZg33IP1kknu12WNDKFOgmywsICe9vaJmKlV33eOBXg21s56B35viAPCvMCe/oK86EwH7Px88ByR14guiUkp2IlpwbCXsfOkNwZEtop7EmT89JLL9GnTx/uuOMOrrjiCm6++Wa6devGf//7Xx5//HF+/etfu12iCMYYzLKFmIXPB4badOsZCHRxbd0uTVygUCe1ZtlhkNAuEMJ69KryvHEOj937Ng9TsAdyv8Hs2hE4FJC3KzD2b/MGzOYNgfmPLBgVfXiPXur3vnaGhPY6aUNcs2XLFm644YbgB44je7TPOOMMfv7znzNv3jwmTpzoZonSzJlDpZi5j2H+8x4A1vlDsS6/SZe+asYU6qTeWLYdPNT7w93+xl8Oe3Zjdm2H3dth1w7M7h2wZ1fg6uZbNmK2bAzMe2ShyKjv7dnrjNUxNXAYt10HhT1pcMXFxbRu3RrbtgkLC6O4uDj4XLdu3ViwYIGL1UlzZ/JzcWZOgW+2QFgY1q/HYQ0crqMezdxxhbrNmzfTrVu3+q5FPMwKj4BOnbE6da403fj9kL/7cMg7HPZ2bYc9OwOXeNm2CbNtU2DeIwu1iIYu3bG6pWN1TYdu6TrUIPUuISGBoqIiAJKSksjJyaF3795AYLxdixYt3CxPmjHni09xnpgauOpBm1jsm+/BSj/V7bKkCTiuUHfvvffSvXt3hg8fzjnnnEO4bjUix8kKDw/sfUtOxeLc4HRTUfG9sLfju7CX+03g0i4bPsNs+Oy7oBffDtLSsdJ6YHfriRMX60o/4h09e/Zk48aNnHXWWQwYMID58+fj8/kIDw9n9erVnH/++W6XKM2MMYb9i+bhzJ4WuEZql+6BCwontHe7NGkijiuNZWZmsmzZMqZPn85zzz3H4MGDGTJkCG3bam+J1A8rLAySUiApBYtzgtONUxEIeFs2fnfIdud22FsAewsw/11LBbDTtgOHbNPSD4e99MB4PV1jT2rpkksuYe/evQBcfPHF+Hw+1qxZg2VZnHPOOVx99dUuVyjNiSk7hHlhJr733wHAOucCrKsysSKjXK5MmhLLmOO/qu6mTZt46623eP/993EchzPPPJOf/exnnHpq098NnJ+fT3n50W/rZVkWyT+47IdXeKk3U3oQtn2N2bIBs+Ur2LIxEPJ+KKpF4LBtWjpWt3Tomo6V0K7xCz4BXlpv1alLfxEREbRvrz0UP6Rtm/d6M9/m4zw+BbZtAjsMe/T1cOFIz42f8+K6O6Kxtm0ndNy0e/fu/OY3v+Gaa65h5cqVrFy5kvvvv5+UlBSGDx/OwIEDiYyMPJG3EDkmq0U09DwNq+dpgceWRWJkOLn/fi8Q9DZvhK2bAidkbPwfZuP/vjtsG5cQ3JNnpaVD1+5YLVq61os0HTNnzuTSSy8lMTGxynP5+fnMnz+fzMxMFyqT5sRs/F9g/Nz+fdC6De3ve5C9iZ08F3qkftTLYLjw8HCioqKCY+sOHTrErFmzWLhwIXfccQfp6dVc9EykAYW1bY99xjmYfmcDhw/b7v7mu8O2mzfCzm3gK4SP/435+N+BoGdZgcO2XXtgXXmLLg3QjL377rsMHTq02lC3f/9+3n33XYU6aTDGGMzqpZhXZgXuIZ6aRtj4LFqc3hd273a7PGmiTijUbdu2jWXLlrFmzRr8fj9nn302v/3tb+nevTvbtm3jqaee4umnn+bBBx+sr3pFjotlh0GnLlidusCAIUDgGk+Bw7YbMVs2BA7bFhbAzm2Y4gPYCnRSgwMHDhCh/x/SQEx5OebFxzH/WgmAddZPsK65FUtnXMsxHFeoW7t2LcuWLePLL78kJiaGkSNHMnToUOLi4oLzdOnShcsvv5wHHnigvmoVqVdWVAtIP7XSpQCMrxC2bgyM05NmJycnh5ycnODjVatW8cknn1Sap6ysjP/85z+kpKQ0cnXSHJi93wbGz23ZCJaNdem1WEMu9tz4OWkYxxXq/va3v9G1a1duueUWBgwYUOMlTdq3b6/T/iWkWHEJ0PdstPlsnj7//PNKFxV+++23q52vXbt2jB07trHKkmbCbPoC54m/wL690LI19ri7sE7t53ZZEkKOK9RNnDiRk08+9o2CO3TooDEnIhIyRo0axfDhwzHGcOONN5KVlUVaWlqleSIiInThYal3zj+XYeY9CRV+6NQFO/M+rMRkt8uSEHNcoa42ga4hLVu2jNdffx2fz0dKSgpjxozhlFNOcbUmEQl9kZGRwTP2p0+fTnx8fL1fXD07O5uXXnqJESNGMGbMGCAwKH7+/PmsWrWKAwcO0KNHD8aOHUtqamq9vrc0PcZfjnn5acy7bwUmnHEu9nW3Bc7qF6mjkLuB5tq1a5kzZw6XXHIJU6dO5ZRTTmHy5MkUFFRzXTIRkePUvn17wsPD2blzJytWrGDhwoX4fD4ACgsLKSsrq/Nrbtq0iZUrV9KlS5dK0xcvXsySJUu4/vrrmTJlCnFxcUyaNImDBzW208vMvr04D/8hEOgsC+viq7Bv/r0CnRy3kLu/1xtvvMGFF17I4MGDARgzZgyffvopy5cv54orrqi39zEGDh60KC6GkhILr10SyLLUWyjycm/wXX9NoTfHcXjyySdZvXp1cFrfvn2Ji4vjqaeeIi0tjcsuu6zWr1daWspjjz3GTTfdxMKFC4PTjTEsXbqUjIwM+vfvD8D48eO58cYbWbNmDUOGDKm3nqTpMFs24sycAr5vIboV9o2/wzr9R26XJSEupEKd3+9n8+bNXHzxxZWm9+7dmw0bNlS7THl5eaWrq1uWRXR0dPD7mhw8aNG9e9LhR0k1zhf61Fto8nJv8PXXNtHRjqs1LFy4kDVr1nD11VfTt29ffve73wWf69evH6tXr65TqJs1axb9+vWjd+/elUJdXl4ePp+PPn36BKdFRETQq1cvNmzYUGOoO95t2/ef9+IZlaHQm/OvlTjPzwR/OSSnEjY+Cyup0zGXC4XeToSX+2us3kIq1BUVFeE4DrGxlW/WHhsbGzws8kPZ2dmVzmZLS0tj6tSpx7wFR3HxCZcrIsepQ4cOtGrlbg2rV6/ml7/8JSNHjsRxKgfMxMRE8vLyav1a//rXv9iyZQtTpkyp8tyRbVd127WjDSs53m3b9yUleffDQVPszfj9+GZP48DrLwPQov9PaHvn/dgtW9fpdZpib/XJy/01dG8hFeqOqC7p1pR+MzIyGDlyZJX58vPz8fv9Nb6HMYG9BR06dGDPnj2euyWLZVnqLQR5uTf4rr+iolyKio7eX3h4eIPe+7WwsLDGu+FERERQWlpaq9cpKChgzpw5ZGVlHfW2iT/chh1r/R7vtu3IvElJSeTm5nru/1FT7c3s34fzxFTMhs8AsH9xBeUjL2PPvv2wb3+tXqOp9lZfvNxfXXo7kW1bSIW6mJgYbNuusldu3759VT7lHhEREVHjld+P9YONjnZo1Srw1Yv/wdRb6PFyb/Bdf0VFxvX+YmNja9wbt2vXLhISEmr1Ops3b2bfvn3cc889wWmO4/DFF1/w1ltvMW3aNCCwxy4+Pj44T1FRUY3bNTixbdv353P759xQmlJvZtvXODMnQ2E+REVj33AHVt/DtzA8jhqbUm8Nwcv9NXRvIRXqwsPD6datG+vXr+ess84KTl+/fj0//vGPXaxMRLymX79+LFy4MHhyBARCZ0lJCW+++SZnnnlmrV7n9NNP56GHHqo07fHHH6djx46MGjWKDh06EBcXx/r164PXxPP7/eTk5HDllVfWa0/S+Jx/r8Y8Nx3KyyCxI/b4+7A6dna7LPGokAp1ACNHjuSxxx6jW7dupKens3LlSgoKCnSGmIjUq9GjR/Pxxx9zxx13cOqpgVvJvfTSS+zYsYOwsDAuvfTSWr1OdHQ0nTtX/iMeFRVFmzZtgtNHjBhBdnY2ycnJJCUlkZ2dTVRUFAMGDKjfpqTRmIoKzMK5mOWLAhNOOzNwhmsdx8+J1EXIhbpzzz2X/fv389prr7F3715SU1O59957G3RsjYg0P3FxcUyZMoVXX32Vjz/+GNu22bZtG2eccQaXXXYZrVvX3x/nUaNGUVZWxqxZsyguLqZ79+5kZWUFz2aV0GIOFOE89SB88SkA1ohfYY26AssOc7ky8bqQC3UAw4YNY9iwYW6XISIeFxcXx7hx4+r9dSdMmFDpsWVZjB49mtGjR9f7e0njMt9swZkxGQr2QGRU4O4QP9IeV2kcIRnqREREmhrz4RqcZ/8GZYegXQfs8VlYKV3dLkuaEYU6EZEafPnll6xZs4b8/PwqtwWzLIs//elPLlUmTYlxKjCLXsS8efi6gb36Yo+7C6tVG3cLk2ZHoU5EpBrvvPMOTzzxBK1btyY5ObnK5UO8eskFqRtTcgDn6Yfhfx8BYA29GOuSa7HCNH5OGp9CnYhINV5//XXOOeccxo8fX+P14KR5M7u2B8bP5e2CyEisa27F7j/Q7bKkGVOoExGpRn5+Ptddd50CnVTLfPxvnNmPwqGDkNA+cP25zie5XZY0cwp1IiLV6NSpE/v27XO7DGlijONg3ngZ84/A/VvpeTr2TXdjtan57h8ijcV2uwARkabo8ssvZ9GiRRQWFrpdijQR5mAJzszJwUBnDf459u0TFeikydCeOhGRw6ZOnVrpcUlJCbfddhtdu3atcrFhy7K4++67G7M8cZHJ3Ykz4wHI/QbCI7CuzsQ+d7DbZYlUolAnInLY9u3bKz22bZuYmBgKCwur7LGzLKsxSxMXmfX/wZn1MBwsgbi22Jn3YaX1cLsskSoU6kREDpsxY0bw+5ycHNLS0qq9VVdpaSmbN29uzNLEBcYYzNL5mMUvgjHQvRf2Lb/Hiol3uzSRamlMnYhINSZOnMjOnTurfW7Xrl1MnDixkSuSxmRKD+I8ORWz6AUwBmvQz7B/938KdNKkaU+diEgd+f1+bFufib3K5O3GmTkZdm6DsHCsK27C/onuNy5Nn0KdiMhhJSUllJSUBB/7fD4KCgoqzVNWVsa7775LXFxcI1cnjcF8/jHOUw9CyQGITcC+5R6sk052uyyRWlGoExE5bMmSJSxYsCD4+MEHH6xx3oyMjMYoSRqJMQazfBHmtblgHEhLx868FyuurdulidSaQp2IyGF9+vShRYsWGGN48cUXGT58OO3atas0T0REBJ07d6ZXr14uVSn1zRw6hHnuMcy6fwJgnfdTrCtvwdLdRCTEKNSJiByWnp5Oeno6AIcOHWLw4MEkJCS4XJU0JPNtXuD6czu2QFgY1mU3Yg36mS5ZIyFJoU5EpBq/+tWv3C5BGpj5cj3Ok3+FA0XQJhb75nuw0k91uyyR46ZQJyIizYoxBrPqH5j5z4DjQJfugfFzCe3dLk3khCjUiYhIs2HKDmFemIl5/x0ArLMvwLo6EysyyuXKRE6cQp2IiDQLpjAfZ+YU2LYJbBvrV9djDf65xs+JZyjUiYiI55mvcnAenwL790HrNtjj7sY6pY/bZYnUK4U6ERHxLGMM5t03MS8/DRUVkJKGPf4+rHYd3C5NpN4p1ImIiCeZ8nLMS09i3lsOgPXj87Gu/S1WlMbPiTcp1ImIiOcY37eB8XObN4BlY/3yGqyhGRo/J56mUCciIp5y6MvPqPi/38G+QmjZCvvGu7BOO8PtskQanEKdiIh4hvPPZeTNexL85dCxc2D8XGJHt8sSaRQKdSIiEvKMvxzzymzM6qUAWGecg3XdbVgtWrpcmUjjUagTEZGQZor24jwxFb7KAcsi9qqbOXD+cND4OWlmFOpERBpQdnY269atY+fOnURGRpKens5VV11Fx47fHRI0xjB//nxWrVrFgQMH6NGjB2PHjiU1NdXFykOD2fpV4ISIvQUQ3RL7xjuJGfYLinfvxhjjdnkijcp2uwARES/Lyclh2LBhPPDAA/zhD3/AcRwmTZpEaWlpcJ7FixezZMkSrr/+eqZMmUJcXByTJk3i4MGDLlbe9Dlr38aZek8g0CV1wr7vIezeP3a7LBHXKNSJiDSgrKwsBg0aRGpqKl27diUzM5OCggI2b94MBPbSLV26lIyMDPr370/nzp0ZP348hw4dYs2aNS5X3zQZvx/n5acxz04LnBDR5yzsex/CSkpxuzQRV+nwq4hIIyopKQGgdevWAOTl5eHz+ejT57tbVkVERNCrVy82bNjAkCFDqn2d8vJyysvLg48tyyI6Ojr4/dEceT4Ur9lm9u/DPPlXzJfrAbB+/mvsn1+OZQf2UYRyb8fi5d7A2/01Vm8KdSIijcQYw9y5czn55JPp3LkzAD6fD4DY2NhK88bGxlJQUFDja2VnZ7NgwYLg47S0NKZOnUr79u1rXU9SUlIdqndf2dcbKJhyFyZvN1Z0SxL+30RanntBtfOGWm914eXewNv9NXRvIRPq8vLyeO211/jf//6Hz+cjISGB888/n0suuYTw8JBpQ0SasdmzZ7N9+3buv//+Ks/98BP8sQb5Z2RkMHLkyCrL5+fn4/f7j7qsZVkkJSWRm5sbMicTOB+8izP371BWBonJ2OOz2NepC/t27640Xyj2Vlte7g283V9degsPD6/Th7NKyx7XUi7YtWsXxhjGjRtHUlISO3bs4Mknn6S0tJRrrrnG7fJERI7qmWee4aOPPmLixIm0bds2OD0uLg4I7LGLj48PTi8qKqqy9+77IiIiiIiIqPa52v5BNMY0+T+exqnALHwOsyw7MOG0M7BvuBNatT5q7aHQ2/Hycm/g7f4aureQCXV9+/alb9++wccdOnRg165dLF++XKFORJosYwzPPPMM69atY8KECSQmJlZ6PjExkbi4ONavX09aWhoAfr+fnJwcrrzySjdKbjJM8X6cpx6CnI8BsH72S6yLr8Kyw1yuTKRpCplQV52SkpLgYOOaNNfBxMei3kKTl3sDb/Y3e/Zs1qxZw9133010dHRwDF3Lli2JjIzEsixGjBhBdnY2ycnJJCUlkZ2dTVRUFAMGDHC3eBeZb7bizJwM+bkQGYU15jbsHzffn4dIbYRsqMvNzeXNN9885l665jiYuC7UW2jycm/grf6WL18OwIQJEypNz8zMZNCgQQCMGjWKsrIyZs2aRXFxMd27dycrKyv4AbS5MR+txXl2GhwqhbaJ2OOzsFLT3C5LpMmzjMsHrl999dVKoas6U6ZM4aSTTgo+LiwsZMKECfTq1Yubb775qMvWtKfOq4OJa0u9hSYv9waNN5jYy/Lz8ytt86pjWRbJycnsbmJ3XTBOBWbxS5ilrwYmnNIHe9xdWK1jav0aTbW3+uDl3sDb/dWlt4iIiNA9UWL48OGcd955R53n+80VFhYyceJE0tPTGTdu3DFfv7kMJj5e6i00ebk38H5/UpUpOYAz6xH47EMArKEXY11yLVaYxs+J1JbroS4mJoaYmNp9CjsS6NLS0sjMzMS2dUMMEZFQZ3bvwJkxGfbshIhIrGvGY59d/fXnRKRmroe62jpyyLVdu3Zcc801FBUVBZ87ckkAEREJLeaTD3BmPwKlByGhHXZmFlaXk469oIhUETKhbv369eTm5pKbm1tlHN2rr77qUlUiInI8jONg3ngF84+XAhPST8O+6W6smDhX6xIJZSET6gYNGhQ8U0xEREKXOViC88yj8MkHAFgXjsT61fVYujuQyAnRb5CIiDQak7szcP253TsgPBzrqkzs837qdlkinqBQJyIijcJ89iHO0w/DwWKIS8C+5V6sbj3dLkvEMxTqRESkQRljMG+9hsl+HoyBk04OBLrY+GMvLCK1plAnIiINxpQexMz5O+ajfwFg/WQ41uU3YoVXf/1QETl+CnUiItIgTH4uzowHYOc2CAvHunwc9sDhbpcl4lkKdSIiUu9Mzsc4Tz4IJQcgJg77lnuwuvdyuywRT1OoExGRemOMwaxYhFkwF4wDaemB8XPxbd0uTcTzFOpERKRemEOHMM9Nx6x7FwDrvMFYV96CFRHpcmUizYNCnYiInDDzbV7g+nPbN0NYGNZlN2ANGoFlWW6XJtJsKNSJiMgJMRs+w3liKhwogtYx2Dffg9XzNLfLEml2FOpEROS4GGMwby/BvDoLHAc6d8POzMJq297t0kSaJYU6ERGpM1NehnnhcczaVQBYZw/Cuno8VmSUy5WJNF8KdSIiUiemsADn8Smw9SuwbKxfXYf1019o/JyIyxTqRESk1sxXOYFAt38ftGqDPe4urF593S5LRFCoExGRWnJWv4l5+SmoqICUrtiZ92G1T3K7LBE5TKFORESOypSXY15+CvPPZQBYPxqANea3WFEtXK5MRL5PoU5ERGpkfIU4T/wFvv4SLAsr4xqs4Zdo/JxIE6RQJyIi1TJff4nz+F9gXyG0bIV9451Yp53pdlkiUgOFOhERqcJZswLz4uPg90NyKvb4LKwOHd0uS0SOQqFORESCjN+PeXUW5p2lgQn9zsa+/nasFi3dLUxEjkmhTkREADBFPpwnp8LGzwGwRl2JNeJXWLbtcmUiUhsKdSIigtm2CWfmZCgsgBbR2Df8DqvPWW6XJSJ1oFAnItLMOe+/g3l+BpSXQYdOgfFzySlulyUidaRQJyLSRCxbtozXX38dn89HSkoKY8aM4ZRTTmmw9zMVFTjzn8WsXByYcPqPAnvoWrZqsPcUkYajgRIiIk3A2rVrmTNnDpdccglTp07llFNOYfLkyRQUFDTI+1Xs8+E8+qdgoLMuGo39mz8o0ImEMIU6EZEm4I033uDCCy9k8ODBwb107dq1Y/ny5fX+XmbHFvbccQ3my/UQ1QL7lnuwL75KJ0SIhDgdfhURcZnf72fz5s1cfPHFlab37t2bDRs2VLtMeXk55eXlwceWZREdHR38viZm4+dUTPszlB2CxGTCxmdhdepy4k00EUd69+IdL7zcG3i7v8bqTaFORMRlRUVFOI5DbGxspemxsbH4fL5ql8nOzmbBggXBx2lpaUydOpX27dsf9b2cuFjyklMIa5tI27sfwG4Tc8L1N0VJSUlul9BgvNwbeLu/hu5NoU5EpImo7lN8TZ/sMzIyGDlyZJX58vPz8fv9R3+j2ybQrns6e/LzMQeKj7/gJsiyLJKSksjNzcUY43Y59crLvYG3+6tLb+Hh4cf8cFbjsse1lIiI1JuYmBhs266yV27fvn1V9t4dERERQURERLXPHeuPhhUThxUWhjHGc388j1BvocvL/TV0bxoVKyLisvDwcLp168b69esrTV+/fj09e/Z0qSoRCTXaUyci0gSMHDmSxx57jG7dupGens7KlSspKChgyJAhbpcmIiFCoU5EpAk499xz2b9/P6+99hp79+4lNTWVe++997jH1ohI86NQJyLSRAwbNoxhw4a5XYaIhCiNqRMRERHxgGa7py48vPat12XeUKPeQpOXe4Pa9ef1n8Hx0rYtQL2FLi/319DbNst49bxhERERkWZEh1+P4uDBg/z+97/n4MGDbpdS79RbaPJyb+D9/poKL/+c1Vvo8nJ/jdWbQt1RGGPYsmWLJy+CqN5Ck5d7A+/311R4+ees3kKXl/trrN4U6kREREQ8QKFORERExAMU6o4iIiKCSy+9tMb7K4Yy9RaavNwbeL+/psLLP2f1Frq83F9j9aazX0VEREQ8QHvqRERERDxAoU5ERETEAxTqRERERDxAoU5ERETEA7x7g7VaWrZsGa+//jo+n4+UlBTGjBnDKaecUuP8OTk5zJ07l2+++Yb4+Hh+8YtfMHTo0Eas+Niys7NZt24dO3fuJDIykvT0dK666io6duxY4zKff/45EydOrDL90UcfpVOnTg1Zbp28+uqrLFiwoNK02NhYnn766RqXCYV1BjB+/Hjy8/OrTB86dCg33HBDlelNfZ3l5OTw+uuvs2XLFvbu3cudd97JWWedFXzeGMP8+fNZtWoVBw4coEePHowdO5bU1NSjvu6///1vXnnlFfbs2UOHDh24/PLLK72uBGjbFtDUf0+O8PK2Dby1fWvK27ZmHerWrl3LnDlzuOGGG+jZsycrV65k8uTJPProo7Rr167K/Hl5eUyZMoXBgwdz6623smHDBmbNmkVMTAxnn322Cx1ULycnh2HDhnHSSSdRUVHByy+/zKRJk3jkkUdo0aLFUZedNm0aLVu2DD6OiYlp6HLrLDU1lT/+8Y/Bx7Zd8w7nUFlnAFOmTMFxnODj7du3M2nSJM4555yjLtdU19mhQ4fo2rUrF1xwAQ8//HCV5xcvXsySJUvIzMwkOTmZhQsXMmnSJKZNm0Z0dHS1r7lx40amTZvGZZddxllnncW6det49NFHuf/+++nRo0dDtxQytG2rqqn+nnyfV7dt4K3tW1PetjXrUPfGG29w4YUXMnjwYADGjBnDp59+yvLly7niiiuqzL98+XLatWvHmDFjAEhJSeHrr7/mH//4R5P6JcrKyqr0ODMzkxtuuIHNmzfTq1evoy4bGxtLq1atGrK8E2bbNnFxcbWaN1TWGVTdWC1atIgOHTqE7Drr168f/fr1q/Y5YwxLly4lIyOD/v37A4FP8jfeeCNr1qxhyJAh1S63ZMkSevfuTUZGBgAZGRnk5OSwZMkSbr/99gbpIxRp21ZVU/09+T6vbtvAW9u3prxta7ahzu/3s3nzZi6++OJK03v37s2GDRuqXearr76id+/elab17duXd955B7/fT3h40/xxlpSUANC6detjznv33XdTXl5OSkoKl1xyCaeddlpDl1dnubm53HTTTYSHh9OjRw8uv/xyOnToUO28obrO/H4/7733HhdddBGWZR113lBYZz+Ul5eHz+ejT58+wWkRERH06tWLDRs21Ljh27hxIxdddFGlaX369GHp0qUNWm8o0bateqHwe9Ictm3g7e2b29u2prvWG1hRURGO4xAbG1tpemxsLD6fr9plfD5ftfNXVFSwf/9+4uPjG6rc42aMYe7cuZx88sl07ty5xvni4+MZN24c3bp1w+/3889//pP/+7//489//vMxP0k1ph49ejB+/Hg6duyIz+dj4cKF/OEPf+CRRx6hTZs2VeYPxXUGsG7dOoqLixk0aFCN84TKOqvOkd+x6tZNQUHBUZf74Z6MuLi4Gn9nmyNt2yoLld+T5rJtA29v39zetjXbUHdEdZ8SjvbJ4YfPHbkhx7E+bbhl9uzZbN++nfvvv/+o83Xs2LHSYOP09HQKCgr4xz/+0aR+gb6/y7tz586kp6dz66238u677zJy5Mhqlwm1dQbwzjvv0LdvXxISEmqcJ1TW2dHUtG7qwhjTpNelW7RtCwiV35Pmsm2D5rF9c2vb1mwvaRITE4Nt21VS8L59+6ok7COqS81FRUWEhYXVavd/Y3vmmWf46KOP+POf/0zbtm3rvHx6ejq5ubkNUFn9adGiBZ07d2b37t3VPh9q6wwgPz+f9evXB8dD1UUorDMg+Im0unVT0+/fkeXq8jvbHGnbdmyh8HvixW0beH/75va2rdmGuvDwcLp168b69esrTV+/fj09e/asdpkePXpUmf/TTz+lW7duTWr8gjGG2bNn88EHH/CnP/2JxMTE43qdLVu21HrQrlvKy8vZuXNnjYcaQmWdfd8777xDbGwsZ5xxRp2XDYV1BpCYmEhcXFyldeP3+8nJyanx9w8CG/XPPvus0rT169eTnp7eYLWGGm3bji0Ufk+8uG0D72/f3N62NdtQBzBy5EhWrVrF22+/zTfffMOcOXMoKCgIDmScN28e06dPD84/dOhQCgoKgtcFevvtt3n77bf5+c9/7lYL1Zo9ezbvvfcet912G9HR0fh8Pnw+H2VlZcF5ftjbkiVLWLduHbt372bHjh3MmzePDz74gOHDh7vRQo2ee+45cnJyyMvL46uvvuLhhx/m4MGDDBw4EAjddXaE4zisXr2agQMHEhYWVum5UFtnpaWlbN26la1btwKBAcRbt26loKAAy7IYMWJE8Lpj27dvZ8aMGURFRTFgwIDga0yfPp158+YFH48YMYJPP/2URYsWsXPnThYtWsRnn31WZYBxc6dtW+j8nhzh9W0beGf71pS3bU03zjeCc889l/379/Paa6+xd+9eUlNTuffee2nfvj0Ae/furTSwMTExkXvvvZe5c+eybNky4uPjue6665rc6ePLly8HYMKECZWmZ2ZmBgem/rA3v9/P888/T2FhIZGRkaSmpnLPPfcc16ephlRYWMjf/vY3ioqKiImJoUePHjzwwAMhv86O+OyzzygoKOCCCy6o8lyorbOvv/660sVDn3vuOQAGDhzI+PHjGTVqFGVlZcyaNYvi4mK6d+9OVlZWpes4HdlIHtGzZ09uv/12Xn75ZV555RWSkpK4/fbbdY26H9C2LXR+T47w+rYNvLN9a8rbNsscz+g9EREREWlSmvXhVxERERGvUKgTERER8QCFOhEREREPUKgTERER8QCFOhEREREPUKgTERER8QCFOhEREREPUKgTERER8QCFOhEREREPUKgTERER8QCFOhEREREPUKgTTygrK+Puu+/m1ltvpaSkJDjd5/Nx4403MmHCBBzHcbFCERGRhqVQJ54QGRnJHXfcQVFRETNnzgTAcRz+/ve/A3Dbbbdh2/rvLiIi3qW/cuIZycnJ3HTTTaxbt46lS5eyYMECPv/8c2699Vbi4+PdLk9ERKRBhbtdgEh9Ovfcc8nJyeH555/HcRwyMjLo3bu322WJiIg0OO2pE8+54IILqKioICwsjBEjRrhdjoiISKNQqBNPKS0tZfr06SQnJxMZGckTTzzhdkkiIiKNQqFOPOXpp5+moKCAO++8k5tvvpkPP/yQN954w+2yREREGpxCnXjGqlWreO+99xg7diypqamcffbZDB8+nBdffJFNmza5XZ6IiEiDUqgTT9i+fTvPPvssAwcOZNCgQcHpV199NV26dOHRRx+luLjYvQJFREQamGWMMW4XISIiIiInRnvqRERERDxAoU5ERETEAxTqRERERDxAoU5ERETEAxTqRERERDxAoU5ERETEAxTqRERERDxAoU5ERETEAxTqRERERDxAoU5ERETEAxTqRERERDzg/wPCSlzg3E/bqwAAAABJRU5ErkJggg==" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot trajectory\n", "plt.subplot(2, 2, 1)\n", "plt.plot(x_bar[0, :], x_bar[1, :])\n", "plt.plot(np.linspace(0, 10, T + 1), np.zeros(T + 1), \"b-\")\n", "plt.axis(\"equal\")\n", "plt.ylabel(\"y\")\n", "plt.xlabel(\"x\")\n", "\n", "plt.subplot(2, 2, 2)\n", "plt.plot(np.degrees(x_bar[2, :]))\n", "plt.ylabel(\"theta(t)\")\n", "# plt.xlabel('time')\n", "\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results are the same as expected, so the linearized model is equivalent as expected." ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "ExecuteTime": { "end_time": "2024-10-22T09:43:55.586267Z", "start_time": "2024-10-22T09:43:55.558478Z" } }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "name": "python3", "language": "python", "display_name": "Python 3 (ipykernel)" }, "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.5" } }, "nbformat": 4, "nbformat_minor": 4 }