1198 lines
127 KiB
Plaintext
1198 lines
127 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {},
|
|
"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",
|
|
"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",
|
|
"* $v$ linear velocity of the robot\n",
|
|
"* $\\delta$ angular velocity 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 \\theta} \\\\\n",
|
|
"\\end{bmatrix}\n",
|
|
"\\quad\n",
|
|
"=\n",
|
|
"\\quad\n",
|
|
"\\begin{bmatrix}\n",
|
|
"0 & 0 & -v\\sin{\\theta} \\\\\n",
|
|
"0 & 0 & v\\cos{\\theta} \\\\\n",
|
|
"0 & 0 & 0\\\\\n",
|
|
"\\end{bmatrix}\n",
|
|
"\\quad\n",
|
|
"$\n",
|
|
"\n",
|
|
"and\n",
|
|
"\n",
|
|
"$\n",
|
|
"B = \n",
|
|
"\\quad\n",
|
|
"\\begin{bmatrix}\n",
|
|
"\\frac{\\partial f(x,u)}{\\partial v} & \\frac{\\partial f(x,u)}{\\partial \\delta} \\\\\n",
|
|
"\\end{bmatrix}\n",
|
|
"\\quad\n",
|
|
"= \n",
|
|
"\\quad\n",
|
|
"\\begin{bmatrix}\n",
|
|
"\\cos{\\theta} & 0 \\\\\n",
|
|
"\\sin{\\theta} & 0 \\\\\n",
|
|
"\\frac{\\tan{\\delta}}{L} & \\frac{v}{L{\\cos{\\delta}}^2}\n",
|
|
"\\end{bmatrix}\n",
|
|
"\\quad\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": [
|
|
"### Kinematics Model"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Control problem statement.\n",
|
|
"\n",
|
|
"N = 4 #number of state variables\n",
|
|
"M = 2 #number of control variables\n",
|
|
"T = 20 #Prediction Horizon\n",
|
|
"dt = 0.25 #discretization step"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def get_linear_model(x_bar,u_bar):\n",
|
|
" \"\"\"\n",
|
|
" \"\"\"\n",
|
|
" L=0.3\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([v*np.cos(theta), v*np.sin(theta), a,v*np.tan(delta)/L]).reshape(N,1)\n",
|
|
" C_lin = dt*(f_xu - np.dot(A,x_bar.reshape(N,1)) - np.dot(B,u_bar.reshape(M,1)))\n",
|
|
" \n",
|
|
" return np.round(A_lin,4), np.round(B_lin,4), np.round(C_lin,4)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Motion Prediction: using scipy intergration"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Define process model\n",
|
|
"# This uses the continuous model \n",
|
|
"def kinematics_model(x,t,u):\n",
|
|
" \"\"\"\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" #[x,y,v,th]\n",
|
|
" #[a,d]\n",
|
|
" \n",
|
|
" L=0.3\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,\n",
|
|
" dydt,\n",
|
|
" dvdt,\n",
|
|
" dthetadt]\n",
|
|
"\n",
|
|
" return dqdt\n",
|
|
"\n",
|
|
"def predict(x0,u):\n",
|
|
" \"\"\"\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" x_bar = np.zeros((N,T+1))\n",
|
|
" \n",
|
|
" x_bar[:,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,\n",
|
|
" x0,\n",
|
|
" tspan,\n",
|
|
" args=(u[:,t-1],))\n",
|
|
"\n",
|
|
" x0 = x_next[1]\n",
|
|
" x_bar[:,t]=x_next[1]\n",
|
|
" \n",
|
|
" return x_bar"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Validate the model, here the status w.r.t a straight line with constant heading 0"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"CPU times: user 4.28 ms, sys: 163 µs, total: 4.45 ms\n",
|
|
"Wall time: 3.86 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": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Check the model prediction"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAChCAYAAACSyiu8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de1xVdb7/8dfaIAIisAXUxDQxrfRYqZBNZZSgv1Iz6jiO9yEtJ0kpLUvPaUzHLKeJofDaFKkxNaNlUjZO46CJTuYjFG9BaTpolheEzUVuclnf3x975EiAXGTvtTZ8no8HDzd7rb32e6/Hcn9Y37W+36+mlFIIIYQQJmMxOoAQQghRFylQQgghTEkKlBBCCFOSAiWEEMKUpEAJIYQwJSlQQgghTMnd6ABCiOYrLi5mzZo1nD59Gk3TmDlzJt26dSM+Pp4LFy4QFBTEnDlz8PHxMTqqEE2mST8oIVzXihUruOWWW4iIiKCyspJLly6xefNmfHx8iIqKIjk5maKiIiZPnmx0VCGaTJr4hHBRJSUlfPvttwwbNgwAd3d3OnToQFpaGuHh4QCEh4eTlpZmZEwhmk2a+IRwUdnZ2fj6+rJq1SpOnTpFSEgI0dHRFBQUYLVaAbBarRQWFhqcVIjmMdUZlK7rPP/88yxbtszoKEKYXlVVFVlZWYwYMYLXXnuN9u3bk5yc3OjXp6SkMH/+fObPn+/AlEI0n6nOoLZu3UpwcDClpaWNWv/MmTN1Ph8YGEhOTk5LRrtmZsskeRp2tUzdunVzcpraAgICCAgIoE+fPgDceeedJCcn4+fnR15eHlarlby8PHx9fet8fWRkJJGRkdW/u8r/J7PlAfNlcrU89f1/Ms0ZVG5uLunp6URERBgdRQiX4O/vT0BAQHVhOXLkCN27dyc0NJTU1FQAUlNTCQsLMzKmEM1mmjOodevWMXny5KuePaWkpJCSkgLAsmXLCAwMrHM9d3f3epcZxWyZJE/DzJjp56ZNm0ZCQgKVlZV07tyZmJgYlFLEx8ezY8cOAgMDmTt3rtExRRunlxQ363WmKFD79+/Hz8+PkJAQMjIy6l3v500S9Z0ymu30FsyXSfI0zOxNfAA33HBDnddsFy5caEAaIWpSFRWoLX8hd892+O0baH7WJr3eFAXq6NGj7Nu3jwMHDlBeXk5paSkJCQnExsYaHU0IIUQzqNNZ6O/Gw48naT9sFOXtPJq8DVMUqIkTJzJx4kQAMjIy2LJlixQnIYRwQaqqCvX5JtSWv0IHHyyzXsQvYmSzWkhMUaCEEEK4PnXuR/R334CsY2ih96BNehLNp+67SBvDdAWqf//+9O/f3+gYQgghGknpOmrHZ6iP3wOP9mgz5mEJG3rN2zVdgRJCCOE6VM559HUJcPQIDAjFMnUWmn+nFtm2FCghhBBNppRC/eufqA2JoIH269lod0eiaVqLvYcUKCGEEE2i8nPR31sJR/bBTQOwPPY0WkDnFn8fKVBCCCEaRSmFStuNen8NVJajjZ+Bdv9INItjBiWSAiWEEKJB6mIh6v3VqP1fQshNWB57Bq1rsEPfUwqUEE6k63qj1tM0rUXb8oW4FurQ1+jvrYDiIrRHp6KNeATNzc3h7ysFSggnmjBhQqPW8/DwICkpycFphLg6VVKM2vAOas926N4Ly5zFaN17Oe39pUAJ4UQeHh788Y9/vOo6Simef/55JyUSom7q20Po696EPBvayHFoD/0Kzb2dUzNIgRLCiR566CGCgoIaXG/06NFOSCNEbepSGWrTOtQXW6FrMJb5v0cLucmQLFKghHCicePGNWq9sWPHOjiJELWp49+ir30Dss+iRY5Be2QKmkd7w/JIgRLCIOfPn6/z+Xbt2uHv74/FQbfuCvFzqqIc9ckHqG3J0CkQy3NL0W4aYHQsKVBCGOVqI/ZbLBYGDx7M448/jr+/vxNTibZGnTphnxbjzA9oQ0egjZuG5ultdCxACpQQhvnNb35DZmYmY8eOrZ4c8aOPPuKmm26iX79+vP/++yQmJvLss88aHVW0QqqyEvX3j1B/2wA+flhiF6INCDU6Vg3ShiCEQTZu3MiMGTPo2rUr7u7udO3alSeeeIJNmzYRHBxMTEwMmZmZRscUrZA68wP6sudRn36ANvgeLIuXm644gZxBCWEYpRQXLlwgOPj/euPn5ORUd+b19PSkqqrKqHiiFVJ6FSplC2pzEnh6YnnyBbTBdxsdq15SoIQwyMiRI/nd737HfffdR0BAADabjS+++IKRI0cCkJ6eTt++fQ1OKVoLdeGc/Q697zPhtjuwTH0KzddqdKyrkgIlhEEefvhhevbsyVdffUVWVhb+/v7MnDmT22+/HYA77riDO+64w+CUwtUppVC7/oH68F2wWNAeexrtF8NcYigtKVBCGOj222+vLkhCtDRly0FfvxwyD8Att2GJjkXr1HBHcbOQAiWEQSoqKvjoo4/48ssvuXjxIuvXr+fQoUOcPXuWBx54wOh4woUppdD3foH64E9QVYk26Um08Add4qzpSnIXnxAGWb9+PadPnyY2Nrb6i+P6669n27ZtBicTrkwV5lPw2v+iEuMhuAeWl97Ect9IlytOYJIzqJycHFauXEl+fj6aphEZGVl9oViI1urrr78mISEBT0/P6i+PTp06YbPZDE4mXJVK/wr9z6u4VFqCNvYxtOFj0CyOnxbDUUxRoNzc3JgyZQohISGUlpYyf/58br31Vrp37250NCEcxt3dvdb8UIWFhXTs2NGgRMJVqeIi1F//hNq7E3r0JuDlxeR7+xod65qZokBZrVasVvvtjl5eXgQHB2Oz2aRAiVbtzjvvZMWKFURHRwOQl5fHunXruOuuuxq9DV3XmT9/Pp06dWL+/PkUFRURHx/PhQsXCAoKYs6cOfj4+DjoEwgzUN+k22+EKMxDe2g82shxuHftCjk5Rke7Zqa7BpWdnU1WVhY33nij0VGEcKiJEyfSuXNnnn32WUpKSoiNjcVqtfLLX/6y0dvYunVrjY6+ycnJDBgwgISEBAYMGEBycrIjogsTUGWl6Emr0N9cBF7eWBb8AcuYiWjupjjvaBGm+iRlZWXExcURHR2Nt3ftwQpTUlJISUkBYNmyZQQGBta5HXd393qXGcVsmSRPwxydyd3dnejoaKKjo6ub9ppyITs3N5f09HQeffRRPvvsMwDS0tJYtGgRAOHh4SxatIjJkyc7Ir4wkDqWYe90m5ttn349ahJaOw+jY7U40xSoyspK4uLiGDp0KEOGDKlzncjISCIjI6t/z6nnFPbywJtmYrZMkqdhV8vUrVu3Zm2zvik2AEpLS6sfd+nSpcFtrVu3jsmTJ9d4XUFBQXVzudVqpbCwsFk5hTmpinJU8p9R//wEArtgee4VtL79jY7lMKYoUEop1qxZQ3BwsMwkKlq1q02xcaUNGzZcdfn+/fvx8/MjJCSEjIyMZmVx1RYJs+UB52SqOP4tBW8uQf/xJF4PPILP1KeweNU9LYbZ9lFz85iiQB09epRdu3bRo0cP5s2bB8CECRMYNGiQwcmEaFlXFp4vvviCI0eO8Mtf/pKgoCAuXLjARx99xIABDU8Ud/ToUfbt28eBAwcoLy+ntLSUhIQE/Pz8yMvLw2q1kpeXh69v/XdyuWqLhNnygGMzqcoK1N8+RG3dCL5WLM8sprz/QGzFJVBc4vQ8zdFQnvpaJExRoG6++WY2btxodAwhnGrDhg0kJCTg4WG/dnDdddcxY8YMnn76ae67776rvnbixIlMnDgRgIyMDLZs2UJsbCxJSUmkpqYSFRVFamoqYWFhjv4YwoHUT6fskwn+8G+0X9yPNv4JNO+2c1emKQqUEG2RUors7Owa3SkuXLhQq29UU0RFRREfH8+OHTsIDAxk7ty5LRFVOJnSq1DbklGfvA9eHbDE/A/awDuNjuV0UqCEMMioUaOqp9u43ASSmprKqFGjmrSd/v3707+//UJ5x44dWbhwoSPiCidR58/Y79A78R0M+gWWyTFoHf2MjmUIKVBCGGTMmDH06NGDr776ipMnT9aabkO0LUrXUal/R320Dtzd0abPRRsS7pJj6LUUKVBCGEim2xAAKvcC+voE+PYQ/NcgLFNno1kDjI5lONONJCFEa7Z9+/ZGrbdjxw4HJxFmoJRC/3I7+uLZ8O+jaFNisMS+JMXpP6RACeFE7733nv1LSdev+pOUlGR0VOFgqiAPfeVS1Lo34fpeWF5KwHLvA226Se/npIlPCCcqKytj/PjxDa7Xrl07J6QRRlH7/oX+/mooK0MbNx0t4iE0i5wv/JwUKCGcaMWKFY1aT/6Kbp1U8UXU+2tQabvhhj5Ypj2Ddt31RscyrVZXoPSvviA/Yz9V5eWgaWiaBSyXf9zA3d3+084D3D2gXTvw8ACP9tDeC629J3h62X+8vMGrg/3Hw0O+NMQ1CwoKMjqCMIg6nIb+3gooKkR7eBLag2PR3Fx3MkFnaHUFiuJCqs79BJWVoBRK10Ep0KugqhKqqqCiAior7P+qmp0iVX3bdW8HHXzA2wc6+oKPL5qPL3T0g47+4OuP5ucP/p3Ar5O90Akh2jxVWoLamIj61z8huKf9JogeIUbHcgmtrkBZIh8mYPz0Ro9DpSoroaIcLpVBeRmUlUFZKZSVoEpLoLQESouhuAhKilDFF6HoIpz7CXUx0/74P0WuRnHz8gb/ALAGolkDqBw/DTzbzhAlQghQ3x1GX5cAthz7GdNDE9Dk+mKjtboC1VTa5Sa/OkYFbkyDntKr7MWrMB8KbKj8PCiwQb4NlZcDebmoM6dQlyZJgRKijVCXLqE2v4favgU6d8PywjK03jcbHcvltPkCda00i9t/mvn8ILhnvUWtXWBgq5iCWQhxderEd+hr34TzP9nvzntkKlr79kbHcklSoIRwssLCQnbt2kV6ejqnTp2ipKQEb29vevbsye23385999131WkyhDmpigrUZ39F/X0TWAOwzF2CdsttRsdyaY0uUOvXryc8PJwbbrjBgXGEaN0++OADdu/ezcCBAxk2bBjBwcF4eXlRWlrKTz/9RGZmJi+88AL33HMPkyZNMjquaKSKrO/R//gS/HgS7e5Ie98m7w5Gx3J5jS5QVVVVLF26FF9fX4YOHcrQoUMJCJDhOIRoCqvVSkJCQp0dcXv16sU999xDeXm5DHXkIlRVFerzTdg++yt4+2CZ9SLabXcYHavVaHSBmjZtGtHR0Rw4cIDdu3fz8ccf06dPH+69916GDBmCp6fcVi1EQx588MHqx/n5+fj7+9dap6SkhAceeMCZsUQzqHM/or/7BmQdo/3dEVSMfcze9US0mCaNrWGxWBg8eDDPPPMMS5cupbCwkFWrVvHEE0+wZs0abDabo3IK0eo8/fTTdT4/Z84cJycRTaF0HT3lU/TfPQPZZ9FmzMP/uSVSnBygSTdJlJSUsHfvXnbv3s2pU6cYMmQI06dPJzAwkM8++4xXXnmF119/3VFZhWhVlKrdLbykpASLjMlmWirnvL1f09EjMCAUy9RZaP6djI7VajW6QMXFxXHo0CFuueUWhg8fTlhYWI129KlTpxIdHe2IjEK0KjNnzgSgvLy8+vFlRUVF3H333UbEElehlEL965+oDYmggfbr2fabIWT4M4dqdIHq06cP06dPr7PNHOzNf2+//Xazgxw8eJC1a9ei6zoRERFERUU1e1tCmNns2bNRSvHqq68ye/bsGsv8/f3p1q2bQclEXVS+zT6G3pF9cNMALI89jRbQ2ehYbUKjC9SYMWMaXKd9Mzuj6bpOYmIiL774IgEBASxYsIDQ0FC6d+/erO0JYWb9+vUDIDExsdn/Z4TjKaVQabtR76+BynK08TPQ7h8p02I4kSn29PHjx+natStdunTB3d2du+66i7S0NKNjCdHitm7dSkVFBVD/H3QVFRVs3brVmbHEz6iLhai3XkO9/Tp0Dcby2zexRIyW4uRkphhJwmaz1ehTFRAQwPfff9+sbS1c6Mv337tTUWGuPlrt2pkrk+Rp2ODBbixY0LLbzM/PJzY2loEDB9KvXz+6deuGp6cnZWVlnDlzhszMTA4cOEB4eHjLvrFoNHXoa3uTXnER2qNT0UY8ItNiGMQUBaquu5nquviYkpJCSkoKAMuWLSMwMLDWOl5ebmiaZroZSc2WSfI0zGLR6jzGrsXEiRMZPXo0O3fuZMeOHfzwww8UFxfj4+NDjx49GDhwIBMmTKBjx44t+r6iYaqkGLXhHdSe7dC9F5Y5i9G69zI6VptmigIVEBBAbm5u9e+5ublYrdZa60VGRhIZGVn9e11TaixYAIGBgY2ebsNZzJZJ8jTsapmu5UYGX19fxowZ06jrusI51LeH0Ne9Cfk2tFHj0Eb/Cs3dXH8wtUWmaFDt3bs3Z8+eJTs7m8rKSvbs2UNoaKjRsYQQrZy6VIb+wRr0P/4WPNpjmf8alqjJUpxMwhRnUG5ubkybNo2lS5ei6zr3338/119/vdGxhHCokpISPvzwQzIzM7l48WKNpu7Vq1c3+PqcnBxWrlxJfn4+mqYRGRnJyJEjKSoqIj4+ngsXLhAUFMScOXPw8ZG5yH5OHf8Wfe0b9tEgIsegPTIFzUPuqjQTUxQogEGDBjFo0CCjYwjhNO+88w42m42xY8eyfPlyZs+ezaeffsqQIUMa9Xo3NzemTJlCSEgIpaWlzJ8/n1tvvZWdO3cyYMAAoqKiSE5OJjk5mcmTJzv407gOVVGO+uQD1LZk6BSI5bmlaDcNMDqWqIMpmviEaIsOHz7Ms88+S1hYGBaLhbCwMObMmcPu3bsb9Xqr1UpISAgAXl5eBAcHY7PZSEtLq74LMDw8XLpsXEGdOoH+8lzUPz5GuycSy6IEKU4mZpozKCHaGqUU3t7eAHh6elJcXIy/vz/nzp1r8rays7PJysrixhtvpKCgoPomI6vVSmFhYYvmdkWqshL1949Qf9sAPn5YYl9CGzDY6FiiAVKghDBIz549yczMZMCAAdx8880kJibi6enJdddd16TtlJWVERcXR3R0dHXBa4zGdNsAcHd3b/Hb7a9FU/NUns6iIGEJlce/w3PocDo+8SyWji078rir7yNHa24eKVBCGOQ3v/lN9Y0R06ZN44MPPqC4uJhZs2Y1ehuVlZXExcUxdOjQ6mtXfn5+5OXlYbVaycvLq3f6+MZ02wDzdQFobB6lV6FSPkVt/jN4emJ58gUqBt+N7VI5XGrZz+Oq+8hZGspTX7cNKVBCGKSwsJA+ffoA9r5RTz75JGAf+qsxlFKsWbOG4OBgRo8eXf18aGgoqampREVFkZqaSlhYWMuHNzmVfdber+n7TLjtDixTn0Lzrd23Upib3CQhhEFefvnlOp9funRpo15/9OhRdu3axTfffMO8efOYN28e6enpREVFcfjwYWJjYzl8+HCbmhlAKYW+8+/ov3safjyJ9tgzWJ76XylOLkrOoIRwMl3Xgf+Mlv2fn8vOnz+PWyPHfbv55pvZuHFjncsWLlx47UFdjLLloK9fDpkH4JbbsETHonUKMjqWuAZSoIRwsgkTJlQ/Hj9+fI1lFouFRx55xNmRXJpSCrV3J+ovf4KqSrSJT6Ld96BMJtgKSIESwslWrFiBUopFixaxePFilFJomoamafj6+uLh4WF0RJehCvPRk1bBwb1w4y32yQQ7y4SPrYUUKCGcLCjI3uy0atUqwN7kd2XfJdE4Kn2PvTiVlaCNjUYb/jCaRabFaE2kQAlhkOLiYt555x327t2Lu7s7SUlJ7Nu3j+PHj9dq+hP/Ry8qRE/8I2rvTujRG8u0OWjBPYyOJRxA7uITwiBvv/023t7erFq1Cnd3+9+Kffv2Zc+ePQYnMy/1TTq5z0xBfb0L7aHxWBb8QYpTKyZnUEIY5MiRI7z11lvVxQns/aEKCgoMTGVOqqwU9eFa1K7Pcbu+F5aZC9B63mh0LOFgUqCEMIi3tzcXL16sce0pJydHrkX9jDr2DfraNyE3G23EIwRMjyW38KLRsYQTSIESwiARERHExcUxfvx4lFIcO3aMv/zlLwwfPtzoaKagyi+hkv+MSvkUArtgee4VtL79/zNnkxSotkAKlBAGefjhh2nXrh2JiYlUVVWxevXq6kkH2zqV9b19MsGzp+19mv47Gs3Ty+hYwsmkQAlhEE3TGDVqFKNGjTI6immoygrU3zaitn4IvlYszyxG6z/Q6FjCIFKghDDQmTNnOHnyJGVlZTWeHzZsmEGJjKN+OoX+bjz88G+0O+9Hm/AEmrdMVd+WSYESwiAff/wxmzZtomfPnrRv377GsrZUoJRehdqWjPrkffDqgCXmf9AG3ml0LGECUqCEMMjWrVt55ZVX6Nmzp9FRDKPOn7FfazrxHQz6BZbJMWgd/YyOJUzC8AKVlJTE/v37cXd3p0uXLsTExNChQwejYwnhcB4eHgQHBxsdwxBK11E7t6I2rQP3dmiPP4t2x70ywKuowfCRJG699Vbi4uJ4/fXXue6669i8ebPRkYRwGF3Xq39+9atf8e6775KXl1fj+cvTcbRWKvcC+hsv2Ucf7/tfWBatwDIkXIqTqMXwM6jbbrut+nHfvn3Zu3evgWmEcKwrp9q4bPv27bWe27BhgzPiOJVSCrVnB2rD26ArtCkxaEP/nxQmUS/DC9SVduzYwV133WV0DCEcZsWKFUZHMIQqyENPWgmHvrafNUXHogV1NTqWMDmnFKglS5aQn59f6/nx48cTFhYG2O9ocnNzY+jQofVuJyUlhZSUFACWLVtGYGBgneu5u7vXu8woZsskeRrmiEyXp9oA+PTTTxkzZkytdT777DNGjx7dou9rJLXvX+jvr4ayMrRx09EiHkKzGH51QbgApxSo3/72t1ddvnPnTvbv38/ChQuverofGRlJZGRk9e85OTl1rhcYGFjvMqOYLZPkadjVMnXrdu2T4m3atKnOArVp06ZWUaBU8UXU+2tQabuh541Yps9Bu+56o2MJF2J4E9/Bgwf55JNPWLx4ca2+IEK0Rt988w1gv2Hi8uPLzp8/j5eX6w/po47sQ1+/HIoK0R6ehPbgWDQ3mUxQNI3hBSoxMZHKykqWLFkCQJ8+fZgxY4bBqYRwnNWrVwNQXl5e/RjsQx/5+/szbdo0o6JdM1VagvrwXdTubRDcE0vsQrQevY2OJVyU4QVq+fLlRkcQwqlWrlwJ2G+YmDVrlsFpWo767jD6ugSw5aA9+N9oD01Ea9fO6FjChRleoIRoq1pLcVKXLqE2v4favgU6d8PywjK03jcbHUu0AlKghBDNpk58Z59M8PxPaMNGoz06Fa29p9GxRCshBUoI0WSqsgK15a+ov28Caycsc5eg3XJbwy8UogmkQAkhmkSdzrJPi/HjSbS7I9DGPY7mLeNnipYnBUqIVujgwYOsXbsWXdeJiIggKirqmrepqqpQn29CbfkrdPDBMutFtNvuaIG0QtRNCpQQrYyu6yQmJvLiiy8SEBDAggULCA0NpXv37s3epjr3I/q7b0DWMbTQe9AmPYnm49uCqYWoTQqUEK3M8ePH6dq1K126dAHgrrvuIi0trVkFSuk6JVs2oCetBo/2aDPmYQmrfzgyIVqSFCghWhmbzUZAQED17wEBAXz//fdN3o7Sq9DfWMTFbw/BgFAsU2eh+XdqyahCXJUUKCFaGaVUrefqGuOyMYMvFw/+Be7DH8LjvgdNMy1GWxlY+Fq0ljxSoIRoZQICAsjNza3+PTc3F6vVWmu9Rg2+HD7SdAP5mi0PmC+Tq+Wpb/BlGfNeiFamd+/enD17luzsbCorK9mzZw+hoaFGxxKiyeQMSohWxs3NjWnTprF06VJ0Xef+++/n+utlmgvheqRACdEKDRo0iEGDBhkdQ4hroqm6rqgKIYQQBmuV16Dmz59vdIRazJZJ8jTMjJmMYLb9YLY8YL5MrSVPqyxQQgghXJ8UKCGEEKbktmjRokVGh3CEkJAQoyPUYrZMkqdhZsxkBLPtB7PlAfNlag155CYJIYQQpiRNfEIIIUzJpftBNTTnjVKKtWvXcuDAAdq3b09MTIzDTntzcnJYuXIl+fn5aJpGZGQkI0eOrLFORkYGr732Gp07dwZgyJAhjB071iF5Lnvqqafw9PTEYrHg5ubGsmXLaix35j46c+YM8fHx1b9nZ2czbtw4Ro0aVf2cM/bRqlWrSE9Px8/Pj7i4OACKioqIj4/nwoULBAUFMWfOHHx8fGq91hHzLJmVGT9rQ8ezo13LsePMTBs3bmT79u34+tqnRJkwYYLT+sXV913YrP2kXFRVVZWaNWuWOnfunKqoqFDPPfecOn36dI119u/fr5YuXap0XVdHjx5VCxYscFgem82mTpw4oZRSqqSkRMXGxtbK880336hXX33VYRnqEhMTowoKCupd7sx9dKWqqir1+OOPq+zs7BrPO2MfZWRkqBMnTqi5c+dWP5eUlKQ2b96slFJq8+bNKikpqc7MDR1zrYVZP2tDx7OjNffYcXamDRs2qE8++cSpOS6r77uwOfvJZZv4rpzzxt3dvXrOmyvt27ePe++9F03T6Nu3L8XFxeTl5Tkkj9VqrT7z8PLyIjg4GJvN5pD3aknO3EdXOnLkCF27diUoKMjh7/Vz/fr1q/WXW1paGuHh4QCEh4fXOpagccdca9GWPmtTNPfYcXYmI9X3Xdic/eSyTXyNmfPGZrPVGOI9ICAAm81W58jOLSk7O5usrCxuvPHGWsuOHTvGvHnzsFqtTJkyxSljpC1duhSA4cOH1xi9GozbR19++SV33313ncuM2EcFBQXVn9lqtVJYWFhrnZaaZ8kVmPmzXu14NkJjjh0j/OMf/2DXrl2EhIQwdepUQ4rYld+FzdlPLlugVCPmvGnMOi2trKyMuLg4oqOj8fb2rrGsV69erFq1Ck9PT9LT0/nDH/5AQkKCQ/MsWbKETp06UVBQwMsvv0y3bt3o169f9XIj9lFlZSX79+9n4sSJtZYZsY8ay4h9ZRSzftaGjmdhN2LEiOprtxs2bOC9994jJibGqRmu9l3YWC7bxNeYOW8CAgJqzEFS37w4LaWyspK4uDiGDh3KkCFDai339vbG09MTsA/mWVVV5fC/tjp1ss+A6jTcCygAAAMdSURBVOfnR1hYGMePH6+x3Nn7CODAgQP06tULf3//WsuM2Edg3z+Xmzbz8vKqLy5fqbHzLLUGZv2sDR3PRmjMseNs/v7+WCwWLBYLERERnDhxwqnvX9d3YXP2k8sWqMbMeRMaGsquXbtQSnHs2DG8vb0d9p9MKcWaNWsIDg5m9OjRda6Tn59f/Zfp8ePH0XWdjh07OiQP2P+CKS0trX58+PBhevToUWMdZ+6jy67WvOfsfXRZaGgoqampAKSmphIWFlZrnbY0z5IZP2tjjmcjNObYcbYrryN//fXXTp1upb7vwubsJ5fuqJuens769eur57x59NFH2bZtG2A/xVVKkZiYyKFDh/Dw8CAmJobevXs7JMt3333HwoUL6dGjR3VTyIQJE6rPTkaMGMHnn3/Otm3bcHNzw8PDg6lTp3LTTTc5JA/A+fPnef311wGoqqrinnvuMXQfAVy6dImZM2eyYsWK6tP+K/M4Yx+98cYbZGZmcvHiRfz8/Bg3bhxhYWHEx8eTk5NDYGAgc+fOxcfHB5vNxltvvcWCBQuAuo+51spsn7W+49mZmnLsGJkpIyODkydPomkaQUFBzJgxw2lnwPV9F/bp06fJ+8mlC5QQQojWy2Wb+IQQQrRuUqCEEEKYkhQoIYQQpiQFSgghhClJgRJCCGFKUqCEEEKYkhQoIYQQpiQFSgghhClJgWpDzp07x2OPPca///1vwD5i9fTp08nIyDA4mRBC1CYFqg3p2rUrkyZNYvny5Vy6dInVq1cTHh5O//79jY4mhBC1yFBHbdDvf/97srOz0TSNV199lXbt2hkdSQghapEzqDYoIiKC06dP88ADD0hxEkKYlhSoNqasrIz169czbNgwPvzwQ4qKioyOJIQQdZIC1casXbuWXr168eSTTzJo0CD+9Kc/GR1JCCHqJAWqDUlLS+PgwYPMmDEDgF//+tdkZWWxe/dug5MJIURtcpOEEEIIU5IzKCGEEKYkBUoIIYQpSYESQghhSlKghBBCmJIUKCGEEKYkBUoIIYQpSYESQghhSlKghBBCmJIUKCGEEKb0/wHaxLzQ/z4cswAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 2 Axes>"
|
|
]
|
|
},
|
|
"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": [
|
|
"Motion Prediction: using the state space model"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 7,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"CPU times: user 4.49 ms, sys: 267 µs, total: 4.75 ms\n",
|
|
"Wall time: 1.92 ms\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": 8,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAChCAYAAACSyiu8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3wU9b3/8dd3EkISQpIlCSJBkAAqULRgUuolppDAUQQbrT8UBIRgqSJQoaVCf4p4MIrWNBquXiI39QgWiWKptZESsBzacBMIyq2IKELIndzIZb7njy0pNAm5kN2ZTT7PxyMPk53dyXvHYd+Z2Z3vV2mtNUIIIYTNGFYHEEIIIeoiBSWEEMKWpKCEEELYkhSUEEIIW5KCEkIIYUtSUEIIIWzJ2+oAQojmKykpYfny5Zw8eRKlFI899hhdu3YlOTmZs2fPEhYWxsyZMwkICLA6qhBNpuQ6KCE81+LFi+nbty+xsbFUVVVx/vx5NmzYQEBAAPHx8aSlpVFcXMy4ceOsjipEk8kpPiE8VGlpKV9++SVDhw4FwNvbmw4dOpCZmUlMTAwAMTExZGZmWhlTiGaTU3xCeKjs7GwCAwNZunQpJ06cICIigokTJ1JYWIjD4QDA4XBQVFRkcVIhmsdWR1CmafKb3/yGhQsXWh1FCNurrq7m+PHjDB8+nJdeeon27duTlpbW6Menp6czZ84c5syZ48KUQjSfrY6gNm3aRHh4OGVlZY26/6lTp+q8PTQ0lJycnJaMdsXslknyNOxymbp27ermNLWFhIQQEhJCnz59APjxj39MWloaQUFB5Ofn43A4yM/PJzAwsM7Hx8XFERcXV/Ozp/x7slsesF8mT8tT378n2xxB5ebmsnv3bmJjY62OIoRHCA4OJiQkpKZY9u/fT7du3YiMjCQjIwOAjIwMoqKirIwpRLPZ5ghq5cqVjBs37rJHT+np6aSnpwOwcOFCQkND67yft7d3vcusYrdMkqdhdsz0nxISEkhJSaGqqorOnTszdepUtNYkJyezefNmQkNDmTVrltUxRRtnlpY063G2KKhdu3YRFBREREQEWVlZ9d7vP09J1HfIaLfDW7BfJsnTMLuf4gO49tpr63zPdt68eRakEeJSurISvfF/yN3+GTz9CirI0aTH26KgDh06xM6dO9mzZw8VFRWUlZWRkpLCjBkzrI4mhBCiGfTJ45hvJcO3X9N+6N1UtPNp8jpsUVBjx45l7NixAGRlZbFx40YpJyGE8EC6uhr9yXr0xvegQwDGtKcIih3RrDMktigoIYQQnk+f/hbzrVfg+GFU5O2ohx5FBdT9KdLGsF1B9e/fn/79+1sdQwghRCNp00Rv/hj9wWrwaY+aMhsjKvqK12u7ghJCCOE5dM4ZzJUpcGg/DIjEmDANFdypRdYtBSWEEKLJtNboz/+CXpsKCtTD01G3xaGUarHfIQUlhBCiSXRBLubqJbB/J1w/AGPSL1EhnVv890hBCSGEaBStNTpzG/qd5VBVgXpwCmrICJThmkGJpKCEEEI0SJ8rQr+zDL3rbxBxPcakJ1Bdwl36O6WghBBCXJbe+3fM1YuhtAR13wTU8HtRXl4u/71SUEIIIeqkS0vQa99Eb/8Mul2LMfO/Udf0dNvvl4ISQghRiz64F3NVCuTnoUaMRo16AOXdzq0ZpKCEEELU0OfL0etXov+6CbqEY8x5ERVxvSVZpKCEEEIAoI8edA5VdPY0Ku4e1L3jUT7tLcsjBSWEEG2crqxAf/gu+tMN0CkM49eJqOsHWB1LCkoIIdoyfeKYc1qMU9+gooejRiegfP2tjgVIQQkhRJukq6rQf/oD+o9rISAIY8YzqAE3Wx3rElJQQgjRxuhT3zjfazpxFDU4BjVmCqpDR6tj1SIFJYQQbYQ2q9HpG9Eb1oCvL8ajT6Juvs3qWPWSghJCiDZAnz2NueIVOHIQbvoRxoTHUYEOq2NdlhSUEEK0Ylpr9NY/o99/CwwDNekJ1C1DWnRaDFeRghJCiFZK5+VgrloEB/dA35swJs5AdQqzOlajSUEJIUQro7XG3PFX9LuvQ3UV6qFHUTF3ecRR08WkoIQQohXRRQUUpv4evWML9O7rnEywc1erYzWLLQoqJyeHJUuWUFBQgFKKuLg4RowYYXUsIYTwKHr3/2K+vZTzZaWo+yehht2DMlw/LYar2KKgvLy8GD9+PBEREZSVlTFnzhxuvPFGunXrZnU0IYSwPV1SjH7vdedRU/dehDz3LAX+gVbHumK2KCiHw4HD4fy4o5+fH+Hh4eTl5UlBCdEA0zSZM2cOnTp1Ys6cORQXF5OcnMzZs2cJCwtj5syZBAQEWB1TuJA+sNv5QYiifNSoB1EjRuPdpQvk5Fgd7Yq5ZiL5K5Cdnc3x48fp3bu31VGEsL1NmzYRHv7vabfT0tIYMGAAKSkpDBgwgLS0NAvTCVfS5WWYby/FfHU++Plj/PZljHvGorxtcdzRImz1TMrLy0lKSmLixIn4+9cerDA9PZ309HQAFi5cSGhoaJ3r8fb2rneZVeyWSfI0zI6ZLpabm8vu3bu57777+PjjjwHIzMxk/vz5AMTExDB//nzGjRtnYUrhCvpwFubKVyHnjHP69fiHUO18rI7V4mxTUFVVVSQlJREdHc3gwYPrvE9cXBxxcXE1P+fUcwgbGhpa7zKr2C2T5GnY5TJ17dr8T0UVFRWxdetWdu/ezYkTJygtLcXf358ePXrwwx/+kJ/85CcEBjb8/sHKlSsZN24cZWVlNbcVFhbWnC53OBwUFRU1O6ewH11ZgU57G/2XDyH0KozZL6D69LM6lsvYoqC01ixfvpzw8HBGjhxpdRwhXObdd99l27ZtDBw4kKFDhxIeHo6fnx9lZWV89913HDx4kCeffJLbb7+dhx56qN717Nq1i6CgICIiIsjKympWFk89I2G3POCeTJVHv6Tw1QWY336N3533EjDhcQy/uqfFsNs2am4eWxTUoUOH2Lp1K927d2f27NkAjBkzhkGDBlmcTIiW5XA4SElJoV27drWW9ezZk9tvv52Kigo2b9582fUcOnSInTt3smfPHioqKigrKyMlJYWgoCDy8/NxOBzk5+df9kjMU89I2C0PuDaTrqpC/3EdetM6CHRgPPEsFf0HkldSCiWlbs/THA3lqe+MhC0K6oYbbmDdunVWxxDC5e66666a7wsKCggODq51n9LSUu68887Lrmfs2LGMHTsWgKysLDZu3MiMGTNYs2YNGRkZxMfHk5GRQVRUVMs+AeFW+rsTzskEv/mnc/y8B3+O8m87n8q03af4hGgrfvnLX9Z5+8yZM5u9zvj4ePbt28eMGTPYt28f8fHxzV6XsI42qzE/WY/53EzIz8WY+luMhJltqpzAJkdQQrRFWutat5WWlmIYTfu7sX///vTv3x+Ajh07Mm/evBbJJ6yhs085JxM89hUMugVj3FRUxyCrY1lCCkoIN3vssccAqKioqPn+guLiYm67zb4TyAnX0aaJzvgT+g8rwdsb9civUD+6w+MGeG1JUlBCuNn06dPRWvPCCy8wffr0S5YFBwdf0UfYhWfSuWcxV6XAl1/ADwZhTJiOcoRYHctyUlBCuFm/fs7rVlJTU2nfvr3FaYSVtNbo7ZvRa98AU6PGP46KHt6mj5ouJh+SEMKNNm3aRGVlJUC95VRZWcmmTZvcGUtYQBfmYy5JRK98Fa6JwHjmVYw7/kvK6SJyBCWEGxUUFDBjxgwGDhxIv3796Nq1K76+vpSXl3Pq1CkOHjzInj17iImJsTqqcCG983PMd5bB+fOoByajho5CNfHDMW2BFJQQbjR27FhGjhzJli1b2Lx5M9988w0lJSUEBATQvXt3Bg4cyJgxY+jYsaPVUYUL6JJz6HeWozO3Qc/rMCY9gbpaZm2oT6srKPN//0pB1i6qKypAKUA5/zIxDDC8wMsLvL3Bu53zq50PtGsHPj7g0x7a+6HaO/+Lrz/4+YFfB/Dv0CoHYxTuFxgYyD333MM999xjdRThRnpfJubqxVBchIofh7rzZygvz51M0B1aXUFRUkT16e+gqgq0Bq3R2gTThOrqf31VOZdXVkBVZa1V1L465V+820GHjuDfAQI6QkAgKiAQOgZBYDB0DEIFBkNQJwhygJ+/nE8Woo3TZaXodanoz/8C4T0wZjyD6h5hdSyP0OoKyoj7KSEPTm70OFRaa2dJVVTA+XKoKHf+t7wcysvQ5aVQVgKlF76K0SXnoKQYsr9HH/sKioucBch/lJuPD+r+BIwhMn29qK20tJT333+fgwcPcu7cuUsu3F22bJmFyURL0V/tw1yZAnk5qLvuR40ag6pjHEZRt1ZXUE2llPrXaT4f6FB7GJHGHP9o04TSYjhXCIX56MJ8KMyDwnxUePeWDy1ahTfffJO8vDzuv/9+Fi1axPTp0/noo4/qnW5GeA59/jx6w2r0Zxuhc1eMJxeiet1gdSyP0+YLqiUow4CAQOfX1dc0qtSE2LdvH8nJyXTs2BHDMIiKiqJXr168+OKLMu2MB9PHvsJc8Sqc+Q4VOwp17wTn+9qiyaSghLCI1rpm5mhfX19KSkoIDg7m9OnTFicTzaErK9Efv4f+03pwhGDMWoDqe5PVsTxaowtq1apVxMTEcO2117owjhBtR48ePTh48CADBgzghhtuIDU1FV9fX66++mqro4kmqjx+BPP3z8C3X6Nui0ONnozy72B1LI/X6IKqrq4mMTGRwMBAoqOjiY6OJiRExooSorl+8Ytf1HwwIiEhgXfffZeSkhKmTZtmcTLRWLq6Gv3JevI+fg86dMSY9jTqJpmDq6U0uqASEhKYOHEie/bsYdu2bXzwwQf06dOHO+64g8GDB+Pr6+vKnEK0OkVFRfTp0wdwXhv16KOPAnD06FErY4lG0t9/i7niFTh+mPa3x1L5s0nOy05Ei2nS2BqGYXDzzTfzxBNPkJiYSFFREUuXLuXnP/85y5cvJy8vz1U5hWh1nnvuuTpvT0xMdHMS0RTaNDHTP8Jc8ARkf4+aMpvgXy2QcnKBJn1IorS0lB07drBt2zZOnDjB4MGDmTx5MqGhoXz88cc8//zzvPzyy67KKkSrYF64Zk7rmq8Lzpw5g5eMLmBbOueM87qmQ/thQCTGhGmo4E5Wx2q1Gl1QSUlJfPHFF/Tt25dhw4YRFRVFu4suOJswYQITJ050RUYhWpUxY8bUfP/ggw9esswwDO699153RxIN0FqjP/8Lem0qKFAPT3d+GEJGinGpRhdUnz59mDx5MsHBwXUuNwyDN954o9lB9u7dy4oVKzBNk9jYWOLj45u9LiHsbPHixWitmT9/Ps8++yxaa5RSKKUIDAzEx0fGfLQTXZCLuXoJ7N8J1w/AmPRLVEhnq2O1CY0uqMYMbNncyddM0yQ1NZWnnnqKkJAQ5s6dS2RkJN26ySi/ovUJCwsDYOnSpYBz/y8sLMThcFgZS/wHrTX6H1vR774GVRWoB6eghoyQaTHcyBYX6h49epQuXbpw1VVXAXDrrbeSmZkpBSVatZKSEt5880127NiBt7c3a9asYefOnRw9erTWqT/hXvpcoXO+pl3bIeJ657QYXcKtjtXm2KKg8vLyLrmmKiQkhCNHjjRrXfPmBXLkiDeVlfa6RqtdO3tlkjwNu/lmL+bOdd3633jjDTp06MDSpUuZNWsWANdddx2rV6+WgrKQ3vt357QYpSWo+yaght8r02JYxBYFdfGnmC6o683H9PR00tPTAVi4cCGhoaG17uPn54VS6pIPcNiB3TJJnoYZhqpzH2sp+/fv57XXXsPb+9//DAMDAyksLHTZ7xT106Ul6LVvord/Bt16Ysz6b1S3nlbHatNsUVAhISHk5ubW/Jybm1vn+fi4uDji4uJqfq5rSo25cyE0NLTR0224i90ySZ6GXS5T165dr3j9/v7+nDt37pJ9PScnR96LsoA+uBdzVQoU5KHuHo0a+QDK215/MLVFtni3r1evXnz//fdkZ2dTVVXF9u3biYyMtDqWEC4VGxtLUlISBw4cQGvN4cOHWbJkCcOGDbM6Wpuhz5djvrscM3ke+LTHmPMSRvw4KSebsMURlJeXFwkJCSQmJmKaJkOGDOGaa66xOpYQLvXTn/6Udu3akZqaSnV1NcuWLSMuLo4RIxo3wWVOTg5LliyhoKAApVTNY4uLi0lOTubs2bOEhYUxc+ZMAgJqz3XW1umjXzqHKsr+HhV3D+re8SgfmRbDTmxRUACDBg1i0KBBVscQwm2UUtx9993cfffdzXq8l5cX48ePJyIigrKyMubMmcONN97Ili1bGDBgAPHx8aSlpZGWlsa4ceNaOL3n0pUV6A/fRX+aBp1CMX6diLp+gNWxRB1sU1BCtEWnTp3i66+/pry8/JLbhw4d2uBjHQ5HzftVfn5+hIeHk5eXR2ZmJvPnzwcgJiaG+fPnS0H9iz5xDPOtZDj1DSp6OGp0AsrX3+pYoh5SUEJY5IMPPmD9+vX06NGj1kXujSmoi2VnZ3P8+HF69+59yUW/DoeDoqKiFsvsqXRVFfpPf0D/cS0EBGHMeAY14GarY4kGSEEJYZFNmzbx/PPP06NHjytaT3l5OUlJSUycOLFmht7GaMxlGwDe3t4u/bh9UzU1T9XJ4xSmLKDq6Ff43jGcjo/MwujYsiOPe/o2crXm5pGCEsIiPj4+hIdf2egEVVVVJCUlER0dzeDBgwEICgoiPz8fh8NBfn4+gYF1vxg35rINsN8lAI3No81qdPpH6A1vg68fxqNPUnnzbeSdr4DzLft8PHUbuUtDeeq7bMMWHzMXoq0wTbPm64EHHuCtt94iPz//ktsvTMfREK01y5cvJzw8nJEjR9bcHhkZSUZGBgAZGRlERbW9GV712dOYL/9/9Psr4AeDMJ5dhLr5NqtjiSaSIygh3OjiqTYu+Oyzz2rdtnbt2gbXdejQIbZu3Ur37t2ZPXt2zfrj4+NJTk5m8+bNhIaG1gyj1BZordFb/4x+/y0wDNSkJ1C3DJFpMTyUFJQQbrR48WLA+UK6Y8cObrnllkuWa635+9//3qh13XDDDaxbt67OZfPmzbuyoB5I5+Vgrl4EWXug700YE2egOoVZHUtcASkoIdzowlQbAOvXr69zGpsPPviAUaNGuTOWR9Nao3dsQf/P61BdhXroUVTMXXLU1ApIQQnhZgcOHACgurq65vsLzpw5g5+fnxWxPJIuKsB8eyns2QG9+zonE+x85eMkCnuQghLCzZYtWwZAZWVlzffgHFkiODiYhIQEq6J5FL17O+aapVBeirp/EmrYPShDpsVoTaSghHCzJUuWAM73o6ZNm2ZxGs9jFhdhpv4evWMLdO+FkTATFd7d6ljCBaSghLCIlFPT6QO7yX17CbogDzVqDGrE/0N5y8tYayX/Z4UQtqfLy9Dvr0Bv/QSva3piPDYX1aO31bGEi0lBCSFsTR8+gLniVcjNRv3XvYQkzCC36JzVsYQbSEEJIWxJV5xHp72NTv8IQq/CmP0Cqk+/f83ZJAXVFkhBCSFsRx8/4pxM8PuTqJ+MQP3sYZSvfPy+rZGCEkLYhq6qRP9xHXrT+xDowHjiWVT/gVbHEhaRghJC2IL+7oRzMsFv/on68RDUmJ+j/GWq+rZMCkoIYSltVqM/TUN/+A74dcCY+lvUwB9bHUvYgBSUEMIy+swp53tNx76CQbdgjJuK6hhkdSxhE5YX1Jo1a9i1axfe3t5cddVVTJ06lQ4dOlgdSwjhQto00Vs2odevBO92qEd+hfrRHTLAq7iE5QV14403MnbsWLy8vHj77bfZsGED48aNszqWEMJFdO5ZzFUp8OUXzskEJ0xHOUKsjiVsyPKCuummm2q+v+6669ixY4eFaYQQrqK1Rm/fjF77BpgaNf5xVPRwOWoS9bK8oC62efNmbr31VqtjCCFamC7Mx1yzBL74B1z3A+dkgmFdrI4lbM4tBbVgwQIKCgpq3f7ggw8SFRUFOCdp8/LyIjo6ut71pKenk56eDsDChQsJDQ2t837e3t71LrOK3TJJnobZMZMn0js/x3xnGZSXo0ZPRsWOQhmG1bGEB3BLQT399NOXXb5lyxZ27drFvHnzLnu4HxcXR1xcXM3POTk5dd4vNDS03mVWsVsmydOwy2Xq2lUmxWuILjmHfmc5OnMbXNsHI+EJ1NXXWB1LeBDLT/Ht3buXDz/8kGeffZb27dtbHUcI0QL0vkzM1YuhuAj104dQd92P8pLJBEXTWF5QqampVFVVsWDBAgD69OnDlClTLE4lhGgOXVaKXpeK/vwvEN4DY8YzqO4RVscSHsryglq0aJHVEYQQLUB/tQ9zZQrk5TiPmEaNQbVrZ3Us4cEsLyghhGfT58+jN6xGf7YROnfFeHIhqtcNVscSrYAUlBCi2fSxr5yTCZ75zvnpvHsnoOS9ZNFCpKCEEE2mKyvRH7+H/tN6cIRgzFqA6ntTww8UogmkoIQQTaJPHndOi/Ht16jb4lAPPILy87c6lmiFpKCEaIX27t3LihUrME2T2NhY4uPjr3iduroa/cl69Mb3IKAjxrSnUTdFtUBaIeomBSVEK2OaJqmpqTz11FOEhIQwd+5cIiMj6datW7PXqb//1nnU9PURVFQ0auwvUAGBLZhaiNqkoIRoZY4ePUqXLl246qqrALj11lvJzMxsVkFp06R041rMNcvApz1qymyMqPqHIxOiJUlBCdHK5OXlERLy7+krQkJCOHLkSJPXo81qzFfmc+7LL2BAJMaEaajgTi0ZVYjLkoISopXRWte6ra4xLhsz+HLJzbfgPWwUPj+5yzbTYthxEF+7ZWoteaSghGhlQkJCyM3Nrfk5NzcXh8NR636NGnw5ZoTtBvK1Wx6wXyZPy1Pf4Msy5r0QrUyvXr34/vvvyc7Opqqqiu3btxMZGWl1LCGaTI6ghGhlvLy8SEhIIDExEdM0GTJkCNdcI9NcCM8jBSVEKzRo0CAGDRpkdQwhrojSdb2jKoQQQlisVb4HNWfOHKsj1GK3TJKnYXbMZAW7bQe75QH7ZWoteVplQQkhhPB8UlBCCCFsyWv+/PnzrQ7hChER9ptm2m6ZJE/D7JjJCnbbDnbLA/bL1BryyIckhBBC2JKc4hNCCGFLHn0dVENz3mitWbFiBXv27KF9+/ZMnTrVZYe9OTk5LFmyhIKCApRSxMXFMWLEiEvuk5WVxUsvvUTnzp0BGDx4MPfff79L8lzw+OOP4+vri2EYeHl5sXDhwkuWu3MbnTp1iuTk5Jqfs7OzGT16NHfffXfNbe7YRkuXLmX37t0EBQWRlJQEQHFxMcnJyZw9e5awsDBmzpxJQEBArce6Yp4lu7Ljc21of3a1K9l33Jlp3bp1fPbZZwQGOqdEGTNmjNuui6vvtbBZ20l7qOrqaj1t2jR9+vRpXVlZqX/961/rkydPXnKfXbt26cTERG2apj506JCeO3euy/Lk5eXpY8eOaa21Li0t1TNmzKiV58CBA/qFF15wWYa6TJ06VRcWFta73J3b6GLV1dX6kUce0dnZ2Zfc7o5tlJWVpY8dO6ZnzZpVc9uaNWv0hg0btNZab9iwQa9Zs6bOzA3tc62FXZ9rQ/uzqzV333F3prVr1+oPP/zQrTkuqO+1sDnbyWNP8V085423t3fNnDcX27lzJ3fccQdKKa677jpKSkrIz893SR6Hw1Fz5OHn50d4eDh5eXku+V0tyZ3b6GL79++nS5cuhIWFufx3/ad+/frV+sstMzOTmJgYAGJiYmrtS9C4fa61aEvPtSmau++4O5OV6nstbM528thTfI2Z8yYvL++SId5DQkLIy8urc2TnlpSdnc3x48fp3bt3rWWHDx9m9uzZOBwOxo8f75Yx0hITEwEYNmzYJaNXg3Xb6G9/+xu33XZbncus2EaFhYU1z9nhcFBUVFTrPi01z5InsPNzvdz+bIXG7DtW+POf/8zWrVuJiIhgwoQJlpTYxa+FzdlOHltQuhFz3jTmPi2tvLycpKQkJk6ciL+//yXLevbsydKlS/H19WX37t387ne/IyUlxaV5FixYQKdOnSgsLOS5556ja9eu9OvXr2a5FduoqqqKXbt2MXbs2FrLrNhGjWXFtrKKXZ9rQ/uzcBo+fHjNe7dr165l9erVTJ061a0ZLvda2Fgee4qvMXPehISEXDIHSX3z4rSUqqoqkpKSiI6OZvDgwbWW+/v74+vrCzgH86yurnb5X1udOjlnQA0KCiIqKoqjR49estzd2whgz5499OzZk+Dg4FrLrNhG4Nw+F05t5ufn17y5fLHGzrPUGtj1uTa0P1uhMfuOuwUHB2MYBoZhEBsby7Fjx9z6++t6LWzOdvLYgmrMnDeRkZFs3boVrTWHDx/G39/fZf/ItNYsX76c8PBwRo4cWed9CgoKav4yPXr0KKZp0rFjR5fkAedfMGVlZTXf79u3j+7du19yH3duowsud3rP3dvogsjISDIyMgDIyMggKiqq1n3a0jxLdnyujdmfrdCYfcfdLn4f+R//+Idbp1up77WwOdvJoy/U3b17N6tWraqZ8+a+++7j008/BZyHuFprUlNT+eKLL/Dx8WHq1Kn06tXLJVm++uor5s2bR/fu3WtOhYwZM6bm6GT48OF88sknfPrpp3h5eeHj48OECRO4/vrrXZIH4MyZM7z88ssAVFdXc/vtt1u6jQDOnz/PY489xuLFi2sO+y/O445t9Morr3Dw4EHOnTtHUFAQo0ePJioqiuTkZHJycggNDWXWrFkEBASQl5fHa6+9xty5c4G697nWym7Ptb792Z2asu9YmSkrK4uvv/4apRRhYWFMmTLFbUfA9b0W9unTp8nbyaMLSgghROvlsaf4hBBCtG5SUEIIIWxJCkoIIYQtSUEJIYSwJSkoIYQQtiQFJYQQwpakoIQQQtiSFJQQQghbkoJqQ06fPs2kSZP45z//CThHrJ48eTJZWVkWJxNCiNqkoNqQLl268NBDD7Fo0SLOnz/PsmXLiImJoX///lZHE0KIWmSoozboxRdfJDs7G6UUL7zwAu3atbM6khBC1CJHUG1QbGwsJ0+e5M4775RyEkLYlhRUG1NeXs6qVasYOnQo77//PsXFxVZHEkKIOklBtf60KKMAAACESURBVDErVqygZ8+ePProowwaNIjXX3/d6khCCFEnKag2JDMzk7179zJlyhQAHn74YY4fP862bdssTiaEELXJhySEEELYkhxBCSGEsCUpKCGEELYkBSWEEMKWpKCEEELYkhSUEEIIW5KCEkIIYUtSUEIIIWxJCkoIIYQtSUEJIYSwpf8DGL248MtlBqYAAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 2 Axes>"
|
|
]
|
|
},
|
|
"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": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## PRELIMINARIES"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 9,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def compute_path_from_wp(start_xp, start_yp, step = 0.1):\n",
|
|
" final_xp=[]\n",
|
|
" final_yp=[]\n",
|
|
" delta = step #[m]\n",
|
|
"\n",
|
|
" for idx in range(len(start_xp)-1):\n",
|
|
" section_len = np.sum(np.sqrt(np.power(np.diff(start_xp[idx:idx+2]),2)+np.power(np.diff(start_yp[idx:idx+2]),2)))\n",
|
|
"\n",
|
|
" interp_range = np.linspace(0,1,np.floor(section_len/delta).astype(int))\n",
|
|
" \n",
|
|
" fx=interp1d(np.linspace(0,1,2),start_xp[idx:idx+2],kind=1)\n",
|
|
" fy=interp1d(np.linspace(0,1,2),start_yp[idx:idx+2],kind=1)\n",
|
|
" \n",
|
|
" final_xp=np.append(final_xp,fx(interp_range))\n",
|
|
" final_yp=np.append(final_yp,fy(interp_range))\n",
|
|
"\n",
|
|
" return np.vstack((final_xp,final_yp))\n",
|
|
"\n",
|
|
"def get_nn_idx(state,path):\n",
|
|
"\n",
|
|
" dx = state[0]-path[0,:]\n",
|
|
" dy = state[1]-path[1,:]\n",
|
|
" dist = np.sqrt(dx**2 + dy**2)\n",
|
|
" nn_idx = np.argmin(dist)\n",
|
|
"\n",
|
|
" try:\n",
|
|
" v = [path[0,nn_idx+1] - path[0,nn_idx],\n",
|
|
" path[1,nn_idx+1] - path[1,nn_idx]] \n",
|
|
" v /= np.linalg.norm(v)\n",
|
|
"\n",
|
|
" d = [path[0,nn_idx] - state[0],\n",
|
|
" path[1,nn_idx] - state[1]]\n",
|
|
"\n",
|
|
" if np.dot(d,v) > 0:\n",
|
|
" target_idx = nn_idx\n",
|
|
" else:\n",
|
|
" target_idx = nn_idx+1\n",
|
|
"\n",
|
|
" except IndexError as e:\n",
|
|
" target_idx = nn_idx\n",
|
|
"\n",
|
|
" return target_idx\n",
|
|
"\n",
|
|
"def road_curve(state,track):\n",
|
|
" \n",
|
|
" #given vehicle pos find lookahead waypoints\n",
|
|
" nn_idx=get_nn_idx(state,track)-1\n",
|
|
" LOOKAHED=6\n",
|
|
" lk_wp=track[:,nn_idx:nn_idx+LOOKAHED]\n",
|
|
"\n",
|
|
" #trasform lookahead waypoints to vehicle ref frame\n",
|
|
" dx = lk_wp[0,:] - state[0]\n",
|
|
" dy = lk_wp[1,:] - state[1]\n",
|
|
"\n",
|
|
" wp_vehicle_frame = np.vstack(( dx * np.cos(-state[3]) - dy * np.sin(-state[3]),\n",
|
|
" dy * np.cos(-state[3]) + dx * np.sin(-state[3]) ))\n",
|
|
"\n",
|
|
" #fit poly\n",
|
|
" return np.polyfit(wp_vehicle_frame[0,:], wp_vehicle_frame[1,:], 3, rcond=None, full=False, w=None, cov=False)\n",
|
|
"\n",
|
|
"def f(x,coeff):\n",
|
|
" return round(coeff[0]*x**3 + coeff[1]*x**2 + coeff[2]*x**1 + coeff[3]*x**0,6)\n",
|
|
"\n",
|
|
"# def f(x,coeff):\n",
|
|
"# return round(coeff[0]*x**5+coeff[1]*x**4+coeff[2]*x**3+coeff[3]*x**2+coeff[4]*x**1+coeff[5]*x**0,6)\n",
|
|
"\n",
|
|
"def df(x,coeff):\n",
|
|
" return round(3*coeff[0]*x**2 + 2*coeff[1]*x**1 + coeff[2]*x**0,6)\n",
|
|
"# def df(x,coeff):\n",
|
|
"# return round(5*coeff[0]*x**4 + 4*coeff[1]*x**3 +3*coeff[2]*x**2 + 2*coeff[3]*x**1 + coeff[4]*x**0,6)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### MPC Problem formulation\n",
|
|
"\n",
|
|
"**Model Predictive Control** refers to the control approach of **numerically** solving a optimization problem at each time step. \n",
|
|
"\n",
|
|
"The controller generates a control signal over a fixed lenght T (Horizon) at each time step."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"\n",
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Linear MPC Formulation\n",
|
|
"\n",
|
|
"Linear MPC makes use of the **LTI** (Linear time invariant) discrete state space model, wich represents a motion model used for Prediction.\n",
|
|
"\n",
|
|
"$x_{t+1} = Ax_t + Bu_t$\n",
|
|
"\n",
|
|
"The LTI formulation means that **future states** are linearly related to the current state and actuator signal. Hence, the MPC seeks to find a **control policy** U over a finite lenght horizon.\n",
|
|
"\n",
|
|
"$U={u_{t|t}, u_{t+1|t}, ...,u_{t+T|t}}$\n",
|
|
"\n",
|
|
"The objective function used minimize (drive the state to 0) is:\n",
|
|
"\n",
|
|
"$\n",
|
|
"\\begin{equation}\n",
|
|
"\\begin{aligned}\n",
|
|
"\\min_{} \\quad & \\sum^{t+T-1}_{j=t} x^T_{j|t}Qx_{j|t} + u^T_{j|t}Ru_{j|t}\\\\\n",
|
|
"\\textrm{s.t.} \\quad & x(0) = x0\\\\\n",
|
|
" & x_{j+1|t} = Ax_{j|t}+Bu_{j|t}) \\quad \\textrm{for} \\quad t<j<t+T-1 \\\\\n",
|
|
"\\end{aligned}\n",
|
|
"\\end{equation}\n",
|
|
"$\n",
|
|
"\n",
|
|
"Other linear constrains may be applied,for instance on the control variable:\n",
|
|
"\n",
|
|
"$ U_{MIN} < u_{j|t} < U_{MAX} \\quad \\textrm{for} \\quad t<j<t+T-1 $\n",
|
|
"\n",
|
|
"The objective fuction accounts for quadratic error on deviation from 0 of the state and the control inputs sequences. Q and R are the **weight matrices** and are used to tune the response.\n",
|
|
"\n",
|
|
"Because the goal is tracking a **reference signal** such as a trajectory, the objective function is rewritten as:\n",
|
|
"\n",
|
|
"$\n",
|
|
"\\begin{equation}\n",
|
|
"\\begin{aligned}\n",
|
|
"\\min_{} \\quad & \\sum^{t+T-1}_{j=t} \\delta x^T_{j|t}Q\\delta x_{j|t} + u^T_{j|t}Ru_{j|t}\n",
|
|
"\\end{aligned}\n",
|
|
"\\end{equation}\n",
|
|
"$\n",
|
|
"\n",
|
|
"where the error w.r.t desired state is accounted for:\n",
|
|
"\n",
|
|
"$ \\delta x = x_{j,t,ref} - x_{j,t} $"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Problem Formulation: Study case\n",
|
|
"\n",
|
|
"In this case, the objective function to minimize is given the sum of the **cross-track error** plus **heading error**\n",
|
|
"\n",
|
|
"$\n",
|
|
"\\begin{equation}\n",
|
|
"\\begin{aligned}\n",
|
|
"\\min_{} \\quad & \\sum^{t+T-1}_{j=t} cte_{j|t,ref}^TQcte_{j,t,ref} + \\psi_{j|t,ref}^TP\\psi_{j|t,ref} + u^T_{j|t}Ru_{j|t} \\\\\n",
|
|
"\\textrm{s.t.} \\quad & x(0) = x0\\\\\n",
|
|
" & x_{j+1|t} = Ax_{j|t}+Bu_{j|t}) \\quad \\textrm{for} \\quad t<j<t+T-1 \\\\\n",
|
|
"\\end{aligned}\n",
|
|
"\\end{equation}\n",
|
|
"$\n",
|
|
"\n",
|
|
"where R,P,Q are the cost matrices used to tune the response.\n",
|
|
"\n",
|
|
"## Error Formulation\n",
|
|
"\n",
|
|
"The track can be represented by fitting a curve trough its waypoints, using the vehicle position as reference\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"A fitted cubic poly has the form:\n",
|
|
"\n",
|
|
"$\n",
|
|
"f = K_0 * x^3 + K_1 * x^2 + K_2 * x + K_3\n",
|
|
"$\n",
|
|
"\n",
|
|
"The derivative of a fitted cubic poly has the form:\n",
|
|
"\n",
|
|
"$\n",
|
|
"f' = 3.0 * K_0 * x^2 + 2.0 * K_1 * x + K_2\n",
|
|
"$\n",
|
|
"\n",
|
|
"\n",
|
|
"* **crosstrack error** cte: desired y-position - y-position of vehicle -> this is the value of the fitted polynomial\n",
|
|
"\n",
|
|
"* **heading error** epsi: desired heading - heading of vehicle -> is the inclination of tangent to the fitted polynomial\n",
|
|
"\n",
|
|
"Then using the fitted polynomial representation in vehicle frame the errors can be easily computed as:\n",
|
|
"\n",
|
|
"$\n",
|
|
"cte = f(px) \\\\\n",
|
|
"\\psi = -atan(f`(px)) \\\\\n",
|
|
"$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"-----------------------------------------------------------------\n",
|
|
" OSQP v0.6.0 - Operator Splitting QP Solver\n",
|
|
" (c) Bartolomeo Stellato, Goran Banjac\n",
|
|
" University of Oxford - Stanford University 2019\n",
|
|
"-----------------------------------------------------------------\n",
|
|
"problem: variables n = 282, constraints m = 343\n",
|
|
" nnz(P) + nnz(A) = 948\n",
|
|
"settings: linear system solver = qdldl,\n",
|
|
" eps_abs = 1.0e-04, eps_rel = 1.0e-04,\n",
|
|
" eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,\n",
|
|
" rho = 1.00e-01 (adaptive),\n",
|
|
" sigma = 1.00e-06, alpha = 1.60, max_iter = 10000\n",
|
|
" check_termination: on (interval 25),\n",
|
|
" scaling: on, scaled_termination: off\n",
|
|
" warm start: on, polish: on, time_limit: off\n",
|
|
"\n",
|
|
"iter objective pri res dua res rho time\n",
|
|
" 1 0.0000e+00 1.06e+00 1.06e+02 1.00e-01 3.10e-04s\n",
|
|
" 200 5.3723e+00 4.66e-04 2.38e-05 1.11e+00 1.71e-03s\n",
|
|
" 350 5.3725e+00 2.13e-04 9.37e-06 1.11e+00 2.73e-03s\n",
|
|
"\n",
|
|
"status: solved\n",
|
|
"solution polish: unsuccessful\n",
|
|
"number of iterations: 350\n",
|
|
"optimal objective: 5.3725\n",
|
|
"run time: 2.95e-03s\n",
|
|
"optimal rho estimate: 4.64e+00\n",
|
|
"\n",
|
|
"CPU times: user 234 ms, sys: 568 µs, total: 235 ms\n",
|
|
"Wall time: 231 ms\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"%%time\n",
|
|
"\n",
|
|
"MAX_SPEED = 1.25\n",
|
|
"MIN_SPEED = 0.75\n",
|
|
"MAX_STEER = 1.57\n",
|
|
"MAX_ACC = 1.0\n",
|
|
"REF_VEL=1.0\n",
|
|
"\n",
|
|
"track = compute_path_from_wp([0,3,6],\n",
|
|
" [0,0,0],0.5)\n",
|
|
"\n",
|
|
"# Starting Condition\n",
|
|
"x0 = np.zeros(N)\n",
|
|
"x0[0] = 0\n",
|
|
"x0[1] = -0.25\n",
|
|
"x0[2] = 0.0\n",
|
|
"x0[3] = np.radians(-0)\n",
|
|
" \n",
|
|
"#starting guess\n",
|
|
"u_bar = np.zeros((M,T))\n",
|
|
"u_bar[0,:]=0.5\n",
|
|
"u_bar[1,:]=0.1\n",
|
|
"\n",
|
|
" \n",
|
|
"K=road_curve(x0,track)\n",
|
|
"\n",
|
|
"# dynamics starting state w.r.t vehicle frame\n",
|
|
"x_bar=np.zeros((N,T+1))\n",
|
|
"x_bar[2]=x0[2]\n",
|
|
"#prediction for linearization of costrains\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",
|
|
" A,B,C=get_linear_model(xt,ut)\n",
|
|
" xt_plus_one = np.squeeze(np.dot(A,xt)+np.dot(B,ut)+C)\n",
|
|
" x_bar[:,t]= xt_plus_one\n",
|
|
" \n",
|
|
"\n",
|
|
"x = cp.Variable((N, T+1))\n",
|
|
"u = cp.Variable((M, T))\n",
|
|
"\n",
|
|
"#CVXPY Linear MPC problem statement\n",
|
|
"cost = 0\n",
|
|
"constr = []\n",
|
|
"\n",
|
|
"for t in range(T):\n",
|
|
" \n",
|
|
" #cost += 1*cp.sum_squares(x[2,t]-1.0) # move car to \n",
|
|
" cost += 100*cp.sum_squares(x[3,t]-np.arctan(df(x_bar[0,t],K))) # psi\n",
|
|
" cost += 1*cp.sum_squares(f(x_bar[0,t],K)-x[1,t]) # cte\n",
|
|
" # Actuation effort\n",
|
|
" cost += cp.quad_form( u[:, t],1*np.eye(M))\n",
|
|
" \n",
|
|
" # Actuation rate of change\n",
|
|
" if t < (T - 1):\n",
|
|
" cost += cp.quad_form(u[:, t + 1] - u[:, t], 1*np.eye(M))\n",
|
|
" \n",
|
|
" # KINEMATICS constrains\n",
|
|
" A,B,C=get_linear_model(x_bar[:,t],u_bar[:,t])\n",
|
|
" constr += [x[:,t+1] == A@x[:,t] + B@u[:,t] + C.flatten()]\n",
|
|
" \n",
|
|
"# sums problem objectives and concatenates constraints.\n",
|
|
"constr += [x[:,0] == x_bar[:,0]] #<--watch out the start condition\n",
|
|
"constr += [x[2, :] <= MAX_SPEED]\n",
|
|
"# constr += [x[2, :] >= MIN_SPEED]\n",
|
|
"constr += [cp.abs(u[0, :]) <= MAX_ACC]\n",
|
|
"constr += [cp.abs(u[1, :]) <= MAX_STEER]\n",
|
|
"\n",
|
|
"prob = cp.Problem(cp.Minimize(cost), constr)\n",
|
|
"solution = prob.solve(solver=cp.OSQP, verbose=True)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 11,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXgUVb7/8ffpBLKQpJN0AggkSFgEXEYgyCIYMCHjgoper4oLg+gwyCag/kYYFBSjuDAwqMCMICLqHXFG0ZmrF4wooMgQDKCCSsKOAUL2jUCSOr8/mkRi9q2ruvN9PU+e3qqqP2lSfLtO1TlHaa01QgghhMXYzA4ghBBCVEcKlBBCCEuSAiWEEMKSpEAJIYSwJClQQgghLEkKlBBCCEvyNjtAU6WlpVV6HBYWRkZGhklpms6d87tzdqg+f6dOnUxKYy2/3s/KWfHf3GqZJE/datrP5AhKCCGEJUmBEkIIYUlu38QnhBA10QV5cPQg+thBOHYIbDYI6wBhHVBhHZ33g0NRNvmubkWWKlAZGRm8+uqr5OTkoJQiLi6OG264wexYQgg3oEtLYN9uCk6nUfbDd3DsIGRdcK4lJMx5m/MFaE3FGG/e3hDaHhXVC3XTXaj2ct7RKixVoLy8vLjvvvuIiorizJkzPP7441xxxRV06dLF7GhCCIvSp0+it2xAf5UI+bkUKgUdOqN69IXIKFRElPM2IMi5fEkJZKZDxil0xqnztyfRu7ajk75EjbwBdeMdFcsL81iqQIWEhBASEgKAn58fnTt3JisrSwqUEKISXVoK3+7A2LwB9u1yNt1dcRW2a+IJGxJDZkFhjeuqNm2gY2fo2Bl14TZzs9Efvo3+7N/obZ85i9TI0c7lhSksVaAulJ6ezqFDh+jRo0el5xMTE0lMTARg4cKFhIWFVXrd29u7ynPuxJ3zu3N2cP/8rYHOyUJv/gS99VPIzYKQMNTNd6OGjUKFOABQvn5QS4GqibKHoMZNRcfehPGPN9DvrUZ//jHqtt+hoq9GKVX3RkSzUlacbqO4uJh58+Zx2223MWjQoFqXlX5Q1uHO2UH6QdXG7H5QOjMd/X/vo7/8FMpK4bIB2GKug8sGoLy8WiST3rcL473VcPwwdO+NbdLjqODQBm/HavuF1fJAzfuZ5Y6gSktLWbRoEcOHD6+zOAkhPJtOT0N//A/09s8BhRp6Leq6/0K1v6jF31v17YftiSvQ2zah//4axpJ52B57FtUusMXfWzhZqkBprVmxYgWdO3dm9OjRZscRQphEpx1Ff/weesdW8PZGxVyP+u2tqNBwl+ZQNi/UsFFoR3uMpU9hLH0a26wFKB9fl+ZorSxVoH766Se2bNlCZGQkjz32GABjx46lf//+JicTQrQ0rTWk/oDx6XrYtR18fFHxt6BGjUHZQ0zNpvr8BtvvH8NY8TzGsmexTX1CLp5wAUsVqN69e7Nu3TqzYwjhtnbv3s3q1asxDIPY2FjGjBlT7XKpqan86U9/YubMmQwePNjFKSvTpaXob75Cf/ohHEmFdoGo0XeiYm+y1KXeqv8Q1O+mot9YirFqEbaJj6FsXnWvKBrNUgVKCNF4hmGwatUq5s6di8PhYPbs2URHR1fppmEYBm+//TZXXnmlSUmddGEBeusG9Kb/hewM52Xf9zyEGnItysfH1Gw1sV0dh1FUiF63Cv3Wcrhvilzd14KkQAnhIVJTU+nYsSMdOnQAYOjQoSQlJVUpUJ988gmDBg3iwIEDLs+oS0og5Xtnp9ivP4ezxdD7Cmz3PuS8Is8NhhyyjboFo6gA/e93wT8Adft4syN5LClQQniIrKwsHA5HxWOHw0FKSkqVZXbs2MG8efNYvnx5jduqq79hufr0HSvLOs3Zb77m3DfbOLcnCV18Btq2xffqOPxvupM23XrW91esF1f0Z9MTppNfVsqZT/6JX/sOtLvtPlPzNITV8tRGCpQQHqK6Lo2/bn564403uOeee7DVcaQSFxdHXFxcxeOa+s2U96nRpaWQnwt5OZCXg87LgfQ09PffwNGDzoVDw1CDYrBdPhB6X0GJjw+5zo036Pesi8v6Zo25D5WZQcHa5RT6+GMbFGNqnvqyWh5wo35QQojGcTgcZGZmVjzOzMysGDqs3IEDB/jLX/4CQF5eHrt27cJms3HVVVc1+P2MD9aSsWcHZTlZUJhfdQFlg+69UbeNQ10eDZ27etT5GmWzwf0PozNPof/+Gvqy/tJHqplJgRLCQ3Tv3p0TJ06Qnp5OaGgo27ZtY/r06ZWWefXVVyvdHzBgQKOKEwD+7fCO7IbR81IItENQMCooGMp/7KGWvdihuShvb2z3PISxYCZ6/Vuoex4yO5JHkQIlhIfw8vJiwoQJJCQkYBgGI0eOJCIigo0bNwIQHx/frO9n++1tBFuwucjVVEQ31LU3ojf9G311HOri5j2n1ppJgRLCg/Tv379Kx/aaCtOUKVNcEalVUDffjd75JcbbK7DNftEtrkZ0B/IpCiFEEyn/dqjb74fDKegvN5odx2NIgRJCiGagBsVAr8vQ769F5+eZHccjSIESQohmoJTCdvckOFOI/uBNs+N4BClQQgjRTFTnSFTczegvP0Uf/MnsOG5PCpQQQjQjddNdYA/BeHsF2igzO45bkwIlhBDNSPn6o+54AI4eQG/eYHYctyYFSgghmpmKHga9r0CvX4uRk2V2HLdlqQK1bNkyHnzwQR555BGzowghRKNVXDBx9iz5b//V7Dhuy1IFasSIEcyZM8fsGEII0WTqoi6oEddT/PnH6OzMulcQVViqQPXt25eAgACzYwghRLNQsTeB1ujP/212FLfkdkMd1TVPjTvNdVIdd87vztnB/fML61HhHfEZdA1nN29A33gnysfX7Ehuxe0KVF3z1FhxrpOGcOf87pwdqs9f0zw1QtSX/013cfbrL9Bfb0KNuMHsOG7FUk18Qgjhadr0vhwu7on+9CO0YZgdx61IgRJCiBaklEKNugXS0+C7b8yO41YsVaCWLFnC3LlzSUtLY9KkSWzatMnsSEII0WSq/1AICcNI/NDsKG7FUuegZsyYYXYEIYRodsrbGxU7Gv2PN9BHD6Iio8yO5BYsdQQlhBCeSg2PBx9fdOJHZkdxG1KgRIuRgTKF+IXyD0ANjUXv2IKW4Y/qRQqUaBF6/16MJ6fKlANCXEDF3QRGGfqLj82O4hakQIlmpc+dxXh3FcZLc8AoA63NjiSEZaj2neA3V6E3f4I+d9bsOJYnBUo0m3M/fY/x9Ax04oeomOuxPfkXVPfeZscSwlJscbdAQT56++dmR7E8S13FJ9yTLilBf/QO2Rs/gGAHtplPo/peaXYsIayp16UQ2R2d+C/08N+ilDI7kWVJgRJNovfvxXh7OaQdxS/uJs7efA/Kz9/sWG6ptLSUlJQUjhw5QmFhIe3ataNr16707NkTb2/ZVT2Fs+PuzehVi2FvMlw2wOxIliV/9aJR9OmTGP98A77ZBiFh2KbPI2jkb916LD6z5OXlsX79ejZv3kxAQACdO3fG19eX4uJiPvnkEwoKCoiJiWHMmDEEBQWZHVc0AxU9DP2PNRiffoiXFKgaSYESDaLPFKH/dx36s4/A5oW6+W5U/K0oHx+zo7mtefPmMXLkSF588UVCQ0OrvJ6VlcWXX37JvHnzWLx4sQkJRXNT3m1QI29Ar38LffJnVMfOZkeyJClQol60UYb+8lP0+rchPxc15FrUrfehQhxmR3N7L774Yq1NeKGhodx8883ccIOMhO1J1NWx6A/fQf/nC9Qt95gdx5LkKj5RK33uLMaW/8OYNw29dhl06IztT4uwTZghxamZXFicXn/99WqXeeONN+Q8lIdRwQ7ocwX6689llPMaSIES1dK52Rgfvo3xxwechamtD7ZJj2P7f8+hLu5pdjyPtXnz5mqf37Jli4uTCFdQQ0ZCZjqk/mB2FEtq0FeyNWvWEBMTw8UXX9xCcYTZ9PHD6MQP0f/ZDGVlcMVAbKPGQK9L5XLYFlQ+cn9ZWVmVUfzT09MJDAw0I5ZoYarfELTPcvT2z1G9LjU7juU0qECVlZWRkJBAUFAQw4cPZ/jw4Tgc0szj7nR+LnrnV+gdWyB1H7T1QQ2PR8XejOogM8q6wtatWwHnpebl98vZ7XamTJliRizRwpSPr7NI7fwKPXYiqk1bsyNZSoMK1IQJExg/fjy7du1i69atvP/++/Ts2ZNrrrmGQYMG4evr21I5RTPTRYXo3dudRemHPWAY0CkSddvvUNfEo9rJN3ZXmjdvHgB///vfueuuu0xOI1xJDRnpHFVizw6IHmZ2HEtp8Dkom83GgAEDmDFjBgkJCeTl5bFs2TJ+//vfs2LFCrKyGj9K7+7du3n44YeZNm0a69evb9Q2FizwatR6ixY17j/kxq5nxntOuB2M7Z9TtuxZjEfGoVf/BU6loX57G7Z5S/F66hVs1/9XtcXJ1VnN+Fwb+7fTVKWlpRX3aytOJSUlrogjXK335RAcirH9C7OTWI7SumGjeRYVFbF9+3a2bt3KkSNHGDRoEDExMYSFhfHvf/+b77//npdeeqnBQQzD4OGHH2bu3Lk4HA5mz57Nww8/TJcuXWpdLy0trdLjzp078fPPaTUsXTNXr1fTumFhYXV2dq3ve+qzZyHle/S+3eh9u+HnI84X7CGo6GGogcMh6pJ6nVuqz3tWl90qn2tj1+vUqeWbOGfMmMHIkSMZPnx4tf2gsrOz2bJlC1988UWd/aB2797N6tWrMQyD2NhYxowZU+n1rVu38uGHzlldfX19efDBB+t1TvnX+1m5+vy9uprVMtUnj/GPN9CJH2J78Q1UoN30PK5W037WoCa+RYsWsWfPHvr06cOoUaMYOHAgbdq0qXh93LhxjB8/vlEBU1NT6dixIx06dABg6NChJCUl1VmgxC90YQEcSUEf3I/+6Tvn+aTSUvD2hh59WfjjVOa8FQURUSibXMBpFU8//TTr16/nscceIyAggIsuugg/Pz/OnDnDiRMnKCoqIiYmhqeeeqrW7RiGwapVqyp9yYuOjq60D7Vv35758+cTEBDArl27+Nvf/sazzz7b0r+iqIMaMhK94X30jq2o2NFmx7GMBhWonj178sADDxAcHFzt6zabjddee61RQbKysipdcOFwOEhJSamyXGJiIomJiQAsXLiQsLAwFizw4plnfmme6dzZWY3nzi3jiSdqnjTP1evVZ11vb2/CwsLqXC8qIpRLg/bz/+74jmuj9lKS8gPGiWMVr3t37U7bG/+bJ/9nMK9vGUCx4Tw/uGyo8/Xhww0SE0upTUN/z/LsVvxcW+I9m0tQUBDjxo3j7rvvJiUlhaNHj1JYWEhAQACRkZH06NGjXn2g6vMl75JLLqm437NnTzIzM5v/FxINpjp3hYhu6K83gRSoCg1u4mspX3/9NXv27GHSpEmAs99HamoqEyZMqHU9T2/i02eL4eRx9IljkHYMfeI4B75Mo3vQceeFDQDBDri4J6pbT1S3XtC1O8o/oMWz/po08TWPnJwcUlNTyc/P58Ld89prr611ve3bt7N79+5K+1BKSgoPPPBAtct/9NFHpKWlVSx/oV9/ETx37ly12/D29q50Ds0KrJapvnkKP/o7BauX4nj5Hby7XGx6Hldq27b6qxct0zXd4XBU+jaXmZlJSEiIiYlcQ2sNeTmQcQqdcYrCMwUYRw+hM05B+glnJ75yXl7QvhM/5fekx11DUF27OwtTsFzq7yl27NjBK6+8QseOHTl27BgREREcO3aM3r1711mgqvuuWdP5xe+//57PP/+cp59+utrX4+LiiIuLq3hc0zkLK57PsFqm+ubRlw4AZSPrkw+w3Xqf6XlcqVnOQbWk7t27c+LECdLT0wkNDWXbtm1Mnz69wduZO7dxzTKzZuU363raMKCwAAryIC8HnZMJudmQmwU5WeicLPbckoMxJR1Kfvl2WgAQaIewDqioS2BYHOqiSOgUAeEXoby92b8okJtuaXjewYMbP4Nnc38+LbVeU9Zt7N9Oc3r33Xd56KGHGDJkCPfffz8vvPACn3/+OceOHatz3fp+yTty5Ah//etfmT17tnQAthBlD4FL+6G3f4G+5R45T4yFCpSXlxcTJkwgISEBwzAYOXIkERERDd7OE0+U0ZgvB488Uvk/NW0YcO4sFJ+Bs8XO2/M/+kwhnCmEwgJmRhZirClEF50vRvl5ztuCfNDVjK/Vpi0Eh4I9lNDfXAwhA5zFyNEBwjoQdkkfMgsKG5S1vv75z8afb2jse7p6vaas29i/neaUkZHBkCFDKj0XExPDxIkTGTduXK3r1udLXkZGBi+99BJTp041pflS1E4NHoFeuQhS9sIll5sdx3SWKVAA/fv3p3///o1e31i3ijylMM4UnX/CcP5oA8rK0GVlUFbq/CkthdISKCn55fbcWedPyblKRzW1ausD/u3Arx0EBsFFEaiAIOf9gCAItDsvGw0Odf74tav1sm7l6wd1FCjhuYKCgsjJySE4OJjw8HD2799PYGAgRj0GE63pS97GjRsBiI+P5x//+AcFBQWsXLmyYp2FCxe26O8k6k9dORjt64f+ehNKCpS1ClRT6T1JnD17BmdLvAKbDZRynruxeTkfe3uDl7fzuTZtwdcf2rRBebdxFpuKn7bOW18/8PVzFg4f53382kG7APDzd64nRDOJjY3lxx9/ZPDgwdx444089dRTKKUYPbp+V3ZV9yUvPj6+4v6kSZOqvShCWIPy8UENGIr+Zht67KRWP8+aRxUor4QVljwBKER9XdixNiYmhksvvZTi4mLpD9iKqMEj0V99ht7zH9RV15gdx1RyFk4ICwsLC5Pi1Nr0ugxCw9Bff252EtNJgRJCCAtRNhtq0AjYtwudm212HFNJgRJCCItRQ0aCYThnG2jFpEAJIYTFqIsioGsPdCsf4VwKlBBCWJAaHANHD6DTjpodxTRSoIQQwoLUVdeAzYb+z2azo5hGCpQQQliQCgqBvlc6hz6qR0dtTyQFSgghLEoNHglZp51zu7VCUqCEEMKi1JWDwMe31V4sIQVKCCEsSvn4ovoNQe/8Cl3f8UE9iBQoIYSwMDVkhHP2hG93mh3F5aRACSGElfW+AuyhGNtb39BHUqCEEMLClM0LNega+O4bdEGe2XFcyjIF6uuvv2bWrFnceeedHDhwwOw4QghhGWrQCCgrRe/80uwoLmWZAhUREcGjjz5Knz59zI4ihBDWEtENOkW2uqv5LFOgunTpIlNQCyFENZRSzj5RB35Enz5pdhyXcbsJCxMTE0lMTARg4cKFhIWFVXrd29u7ynPuxJ3zu3N2cP/8wrOpQdegP3gTvf0L1E13mR3HJVxaoBYsWEBOTk6V5++66y4GDhxYr23ExcURFxdX8fjXs+e6+4y67pzfnbND9fnlqF5YhQoNh16XOYc+Gn0nSimzI7U4lxaoJ554wpVvJ4QQHkUNHoFe8zIc2g9Rl5gdp8VZ5hyUEEKI2qn+Q6FN21ZzsYRlCtSOHTuYNGkS+/fvZ+HChSQkJJgdSQghLEX5t0P95ip00lZ0aanZcVqcZS6SuOqqq7jqqqvMjiGEEJamBo9w9ofauwt+U79z9+7KMkdQQggh6uHSfhAQhPHZR2itzU7ToqRACSGEG1HebVA3j4Uf9qC//NTsOC1KCpQQQrgZFXM9XHI5et0qdOZps+O0GClQQgjhZpTNhu1300BrjDVLPbapTwqUEEK4IRXeEXX7/c6mvi0bzI7TIqRACSGEm1Ix10Gf36DfW43OOGV2nGYnBUoIIdyUUsrZ1KfAWPMy2jDMjtSspEAJIYQbU472qP+eAD9+i978f2bHaVZSoIQQws2p4fHQtx/6H6s9ajoOKVBCeJDdu3fz8MMPM23aNNavX1/lda01r7/+OtOmTePRRx/l4MGDJqQUzc3Z1DcVvLww3viLxzT1SYESwkMYhsGqVauYM2cOixcv5quvvuL48eOVltm1axcnT55k6dKlTJw4kZUrV5qUVjQ3FRqOuuMB2L8X/fE69Kk0dEEe2iir1/paa/TZs+jiInTxGfTZYvS5s+iSc+iSEnRpCbqszKWXtFtmLD4hRNOkpqbSsWNHOnToAMDQoUNJSkqiS5cuFcvs3LmTa665BqUUvXr1orCwkOzsbEJCQsyKLZqRujoOnfw1+sN30B++c/5JBf4B0C4Q2gWQbQ+mrKgQiovhbDEUn4GzZ5z3G1J8lA1syrl9ZTt/++s5qn55rP77fmwx1zXo95ECJYSHyMrKwuFwVDx2OBykpKRUWebCWYMdDgdZWVlVClRdM1eXs+IsxFbL5Oo8+omXOPftNxj5ORj5eej8PIz8XIyCXHR+HrqwgDZt2qLC7Sg/P5SvPzY/f5SvH8rHF7y8zxcq7bw9/6Mr7htglN8azucNw/m4UpDKxc6nz+W0beDnIAVKCA9RXdPLr2ddrc8yUPfM1eWsOIuy1TKZkqdrzxpfMuvzOQtQw/vWNHO1nIMSwkM4HA4yMzMrHmdmZlY5MnI4HJX+c6puGSGsQgqUEB6ie/funDhxgvT0dEpLS9m2bRvR0dGVlomOjmbLli1ordm/fz/+/v5SoIRlSROfEB7Cy8uLCRMmkJCQgGEYjBw5koiICDZu3AhAfHw8/fr1Izk5menTp9O2bVsmT55scmohaqa0pw6DK4QQwq15XBPf448/bnaEJnHn/O6cHdw/vxms+JlZLZPkaTyPK1BCCCE8gxQoIYQQluQ1f/78+WaHaG5RUVFmR2gSd87vztnB/fObwYqfmdUySZ7GkYskhBBCWJI08QkhhLAkj+kHtXv3blavXo1hGMTGxjJmzBizI9VbRkYGr776Kjk5OSiliIuL44YbbjA7VoMZhsHjjz9OaGioW10pBFBYWMiKFSs4duwYSikeeughevXqZXYsS7PaPjdlyhR8fX2x2Wx4eXmxcOFCl2dYtmwZycnJ2O12Fi1aBEBBQQGLFy/m9OnThIeHM3PmTAICAkzLs27dOj777DOCgoIAGDt2LP3793dJngbTHqCsrExPnTpVnzx5UpeUlOhHH31UHzt2zOxY9ZaVlaUPHDigtda6qKhIT58+3a3yl/vXv/6llyxZop977jmzozTYyy+/rBMTE7XWWpeUlOiCggKTE1mbFfe5yZMn69zcXFMz7N27Vx84cEDPmjWr4rm1a9fqDz74QGut9QcffKDXrl1rap53331Xf/jhhy7L0BQe0cR34TQD3t7eFdMMuIuQkJCKk5Z+fn507tyZrKwsk1M1TGZmJsnJycTGxpodpcGKior44YcfuPbaawHn6NPt2rUzOZW1ufs+11L69u1b5egoKSmJmJgYAGJiYlz6OVWXx514RBNffaYZcBfp6ekcOnSIHj16mB2lQd544w3uvfdezpw5Y3aUBktPTycoKIhly5Zx5MgRoqKiGD9+PL6+vmZHsyyr7nMJCQkAjBo1qtJo7GbKzc2tGO8wJCSEvLw8kxPBhg0b2LJlC1FRUYwbN86yRcwjjqB0PacQsLri4mIWLVrE+PHj8ff3NztOvX3zzTfY7Xa3uXT118rKyjh06BDx8fG88MIL+Pj4VDtduviFFfe5BQsW8PzzzzNnzhw2bNjAvn37TM1jVfHx8bz88su88MILhISE8Oabb5odqUYeUaDqM82A1ZWWlrJo0SKGDx/OoEGDzI7TID/99BM7d+5kypQpLFmyhO+//56lS5eaHaveHA4HDoeDnj2dc+gMHjyYQ4cOmZzK2qy4z4WGhgJgt9sZOHAgqamppuYpZ7fbyc7OBiA7O7vi4gSzBAcHY7PZsNlsxMbGcuDAAVPz1MYjClR9phmwMq01K1asoHPnzowePdrsOA129913s2LFCl599VVmzJjBZZddxvTp082OVW/BwcE4HA7S0tIA+O677ypNky6qsto+V1xcXNG8XFxczLfffktkZKRpeS4UHR3N5s2bAdi8eTMDBw40NU95sQTYsWMHERERJqapncd01E1OTmbNmjUV0wzcdtttZkeqtx9//JEnn3ySyMjIimYSS1/6WYu9e/fyr3/9y+0uMz98+DArVqygtLSU9u3bM3nyZMu2y1uFlfa5U6dO8dJLLwHOJtthw4aZkmfJkiXs27eP/Px87HY7d9xxBwMHDmTx4sVkZGQQFhbGrFmzXPa3VV2evXv3cvjwYZRShIeHM3HiRNOPfmviMQVKCCGEZ/GIJj4hhBCeRwqUEEIIS5ICJYQQwpKkQAkhhLAkKVBCCCEsSQqUEEIIS5ICJYQQwpKkQAkhhLAkKVAe7uTJk9x///0cPHgQcI5C/cADD7B3716TkwkhRO2kQHm4jh07cs899/Dyyy9z9uxZli9fTkxMDJdeeqnZ0YQQolYy1FEr8fzzz5Oeno5Siueee442bdqYHUkIIWolR1CtRGxsLMeOHeO6666T4iSEcAtSoFqB4uJi1qxZw7XXXst7771HQUGB2ZGEEKJOUqBagdWrV9OtWzcmTZpE//79+dvf/mZ2JCGEqJMUKA+XlJTE7t27mThxIgC/+93vOHToEFu3bjU5mRBC1E4ukhBCCGFJcgQlhBDCkqRACSGEsCQpUEIIISxJCpQQQghLkgIlhBDCkqRACSGEsCQpUEIIISxJCpQQQghLkgIlhBDCkqRACSGEsCQpUEIIISxJCpQQQghLkgIlhBDCkqRACSGEsCRvswM0VVpaWrXPh4WFkZGR4eI01bNKFqvkAPfJ0qlTJxensSbZz9wvB7hPlpr2MzmCEkIIYUlufwQlrEOfKYL0NPTJnyH9BOTngAbQoHWl+3m+vhjFxaAUoEAByua8rbTRGh80m7Mjr4dOF7fItoX70KUlUJAPhfkVt/pcMZw76/w5e/aX++fOUfH3+Ks5X/P8/DAMwNcX2vo6b318UT6+4OsP9mCwh0K7QJRNjhFqIwVKNIouLUHv+g/sTUafSoP0NMjL+WUBpcCvHdjKC5D65XmlOGvzQhtl5wvXBTu61lWL1IVPqCovNlnZJZdKgfJQ2jAgPxdyMiEnG53rvCU3C52bDTlZUJDnLErFZ+reoM0GPr7g3cZ5v8Ivf5dn0ejiYjhXXKl4Vfl65eUFQX5AUrEAAB/oSURBVCFgd/6okDDo0AnVsTN07AKh4a2+gEmBEg2is06jt2xAb93oLEgBQXBRF9QVA6F9J1SHTtChE4R3RLX1qXE7Vmob9w8Lo8giWUTD6cJ8OPkzZ/YVYRxKhazT6Ix0yDz/U1pSdaVA+y+F4aIu0C7Q+bccEAjtglABgdAuAHz9oK3PLz9e3qg6viSV/21rrZ1HWmfPwNli521RIeTlOItjbhbk5qBzsyAzHZ2yF4oKfylkbdpC+4tQHbtApwhU155wcQ+UPaTZP0OrkgIl6qQNA37cg/H5J7BnB6Dh8mhsI66HS/uhbF5mRxQeTmsNp0/Az0edTcinjp+//dnZHAfklS8caAdHe1SXi+E3V4Ej3Hl0EhzqLEpBwSjvNi2eWSkFPj7On1+/Vs3yWmvn0d7Jn9Enj8Opn9Enf0YfOwjJX6O14Vww2OEsVF17oC7uARf3RAUEtejvYhYpUKJWeueXGB+85WzCCwhCXXcravhvUeEdzY4mPJguKYEjqegDP6BTf4ADPzr/8y4XFAwdO6P6DXHeduhCSK8+ZNu8ned63JBSyvl7BQWjel1a6TV9thiOHUQfTnV+LodT0Xt2OIsawEURqJ59oUdf562jfZ1Heu5ACpSols7PQ7+zAr3zS4johnpgFmrA1ag2Lf/NU7Q+2jDgcAp6TxJ6/3dwOPWXprn2F6EuGwA9eqMioqBDZ5R/uyrb8A4LQ3loU63y8XUWnx59K57TZ4qcxergT+jUH9BJX8KWDc4mwpAwVI8+FF15FbpjBHTp6pYtHVKgRBV693aMtcugsAA15l7Udf+F8nK/P25PtXv3blavXo1hGMTGxjJmzJhKr2utWb16Nbt27cLHx4fJkycTFRVVr3VdSZecgx+/Re/+D3pPkvOcjM3mbLK69kZU9z7OohTUes65NITy84feV6B6XwHgvOgo7Sg6ZR+k7EOn7CM/aatzYV8/iOqN6tEH1aMPdOuF8vUzMX39SIESFXRhAfrvr6G3fw5dumGb8RQqopvZscQFDMNg1apVzJ07F4fDwezZs4mOjqZLly4Vy+zatYuTJ0+ydOlSUlJSWLlyJc8++2y91m1JuuQc/HwEfeQAet8u2LvLefGAjx/qsv5w5SDU5QNQ7QJdksfTKJsXdOmG6tINRt6I1ppQXUrmjq/gfFOp/tf/OJsFbTZns2Dni51HV526QpeuzisHLdQ0KAVKAKC//wZjzcuQl4MafSfqxjtcciJZNExqaiodO3akQ4cOAAwdOpSkpKRKRWbnzp1cc801KKXo1asXhYWFZGdnc/r06TrXbQidn0uZUYLOyqrmRe28lPvIQTh6AH30AJw4BmVlzteDQ1GDR6CuHASXXCFNxy1AKYVX+EXYBo+AwSMA0EWFcPBHZ7E6ehCdug92bP7lykE/f+gU6byoxM8ffPzAz895BObrD75+qKBgCAmDEEetV+o2BylQrZzWGr3+bfTH6+CiCGxT/oS6uKfZsUQNsrKycDgcFY8dDgcpKSlVlgkLC6u0TFZWVr3WbQj999fI2LGl7gUD7dC1O+qKgajI7hAZBWEdLPVNvbVQ/u3gsgHOc3rn6aJCSDuCPn7EefvzEfSxQ85+YcVFzqPcC1TqzxUQCMFhEBqGCg1DjRqDan9Rs+WVAtWK6dJS9JuvoL/ehBo2CnX3H1Bt2podS9RC66qjafz6P/qalqnPuuUSExNJTEwEYOHChZUKXrlzN92JHjICwyirdhu2QDve3S/BFhLmkmLk7e1dbU5Xs0oOqG+WMIjsWuOruqwMfbYYXVSIPlOIkZ1JWeZpjMx0yjLSKctMx8hMp/TbJPwd4QTcO6kJWX61ToOWFh5DF5/B+Ovz8H0y6ua7nc168o3W8hwOB5mZmRWPMzMzCQkJqbLMhZ2gy5cpLS2tc91ycXFxxMXFVTyutlN1xwjCLutXe4drA7jgPVuSVTp/WyUHNHcWG/gFOn+qG3nlkXEUpZ+kuIb3k8FiRb3ovByMl/4Ee3ejxk3FdtNdUpzcRPfu3Tlx4gTp6emUlpaybds2oqOjKy0THR3Nli1b0Fqzf/9+/P39CQkJqde6QjRaoB19YV+1ZiBHUK2MTk/DWDIfcrOc55t+M9DsSKIBvLy8mDBhAgkJCRiGwciRI4mIiGDjxo0AxMfH069fP5KTk5k+fTpt27Zl8uTJta4rRLMItFfuTN0MpEC1IiUp+zAW/hG0ge2RBFTUJWZHEo3Qv39/+vfvX+m5+Pj4ivtKKR588MF6rytEc1CBdvSRA826TSlQrYT+6TuyXl4AAUHO/k0dO5sdSQjhSQKCoECOoEQD6QM/Yry8AK/2F6Efnt+qRkMWQrhIoN05GntpSbP1oXRJgVq2bBnJycnY7XYWLVpU5fXahmYRTaOPHsRY+hTYQwiZ/xeyDbMTCSE8UuD5EdUL8pwjrjcDl1zFN2LECObMmVPj6xcOzTJx4kRWrlzpilgeT584jrFkHvj6YZu1AK9Qa/TNEEJ4HhVod97Jz6t9wQZwSYHq27cvAQEBNb5e09AsovH06ZMYf34ClMI2cwHK0d7sSEIITxZQXqCa7zyUJc5B1TQ0S3WdCOvTwx3csTd38ynLOk32X+ajSs8R8swy2nTtbkqO2kgWITzM+SY+nZ9b7YSMjWGJAtWQIVjq1cMdT+7NXTudn4vx4hzIycb2yAJy29nh/Hu31s+kLo3p4S6E+JXAYOdtgZs18dWlpqFZRMPookLnOaeMU9imPYHq1svsSEKI1qJdAChbszbxWaJA1TQ0i6g/XVqKsWIh/HwE20OzUZdcZnYkIUQromw25+jmzXiRhEua+JYsWcK+ffvIz89n0qRJ3HHHHZSWlgK1D80i6kdrjX57OfywBzX+YdTlA+peSQghmltAELoZO+u6pEDNmDGj1tdrG5pF1E1//B76y09Ro+/EdnWs2XGEEK1VoB3yPKyJTzSe8Z/N6PVvOWcnvflus+MIIVqzwOYd7kgKlBvTKfvQb/wFel2KGjdNpswQQphKBdrd7xyUaH765M8YryZAWAdsk+eg2jTP2Fei8fLy8tiyZQvJyckcOXKEoqIi/P396dq1K1deeSUjRowgKCjI7JhCtJwAOxTmo8vKUF5eTd6cFCg3pPNznePr2WzYps9DtQs0O1Kr984777B161b69evHtddeS+fOnfHz8+PMmTP8/PPP7Nu3jz/+8Y8MGzaMe+65x+y4QrSMoPOjSRTmQVDTr8SWAuVmdMk555FTTha2R55BhXc0O5IAQkJCWLp0KW2qOZLt1q0bw4YN49y5c2zatMmEdEK4SMAF4/E1Q4GSc1BuRGuNfvMVOPAjtgdmorr3NjuSOO/666+vKE45OTnVLlNUVMR1113nylhCuJQqH9G8mTrrSoFyI3rD++jtX6DG3IsacLXZcUQNHn744WqfnzlzpouTCOFi50c01810oYQUKDeh9+xAv/8mauBw1A3/bXYcUYvqxpYsKirCZpPdTXi4ijmhmucISs5BuQH98xGM1xZBZHfU+OlyOblFPfTQQwCcO3eu4n65goICrr5ajnqFh2t3vkA1U2ddKVAWp/PzMF5e4Jx0cMqfUG19zI4kajBt2jS01jz33HNMmzat0mvBwcEyMrrweMrLC9oFyhFUa6BLS5wDwOZmY/t/z6FCmmcaZdEy+vbtC8CqVavw8ZEvEqKVCgxCy0USnk1rjX7nr7D/e2eznkydYWkff/wxJSUlADUWp5KSEj7++GNXxhLC9QKabzQJOYKyKL3pf9FbN6Kuvx3boBiz44g65OTkMH36dPr160ffvn3p1KkTvr6+FBcXk5aWxr59+9i1axcxMfJvKTxckB1OHG+WTUmBsiC9bzd63Uq4chBqzL1mxxH1cPfddzN69Gi++OILNm3axNGjRyksLCQgIIDIyEj69evH2LFjCQyUUT+EZ1MBdnTBvmbZVoMKVHFxMYWFhbRr1w5fX99mCSAq0+lpGH99ATp2cXbGlUuT3UZQUBA333wzN998c4u9R0FBAYsXL+b06dOEh4czc+ZMAgICqiy3e/duVq9ejWEYxMbGMmbMGADWrVvHZ599VjEm4NixY+nfv3+L5RWtUGAQFOSjDaPJ/3/VWaCOHj1KYmIiycnJnD59uuL59u3bc+WVVzJq1CgiIyObFEI46TNFGK8kgFLYps5F+fqbHUlYzPr167n88ssZM2YM69evZ/369dx7b+WjbMMwWLVqFXPnzsXhcDB79myio6Pp0qULADfeeGOLFlHRygXaQRtQWPBLv6hGqrVALVmyhOPHjzN06FCmTZtW7QCYS5cupUuXLnVOSihqpw0DY9Wf4dTP2GY8JWPsubGioiLee++9ilmkL+y4u3z58iZtOykpifnz5wMQExPD/PnzqxSo1NRUOnbsSIcOHQAYOnQoSUlJFQVKiBYVUD7cUU7LFqhhw4YRHR1d9f0DArjkkku45JJLuPXWW/nmm2+aFEKA/vBt2LMDdfcfUH1+Y3Yc0QQrV64kKyuL22+/nZdffplp06bx0UcfMWjQoCZvOzc3l5AQ5yCcISEh5OVVvVoqKysLh+OXLgkOh4OUlJSKxxs2bGDLli1ERUUxbty4apsIExMTSUxMBGDhwoWEhYVVm8fb27vG11zNKlmskgPMyXK2SyQ5gN2maHvBezcmS60F6sLilJKSQs+ePassk5qayoABAxr0pqIyY8cW9MfvoYbHo0bcYHYc0UTffvstixcvJjAwEJvNxsCBA+nevTvPP/88o0ePrnP9BQsWVDvg7F133VWv969uqKXy0Ufi4+O5/fbbAXj33Xd58803mTx5cpXl4+LiiIuLq3ickZFR7XuFhYXV+JqrWSWLVXKAOVm04bzN/fkoqmNEvbLU1Im93hdJPPPMM6xZs6bK8wkJCaxevbq+mxG/oo+kotcshR59nUdPMoyR29Na4+/vPH/o6+tLYWEhwcHBnDx5sl7rP/HEEzW+Zrfbyc7OJiQkhOzs7GonQHQ4HGRmZlY8zszMrDjqCg4Orng+NjaW559/vl6ZhKi38816Oj+Xpv5vVuclFoZhYBiGs+Oo1hWPDcPgxIkTeDXDrImtlc7Nxnj1WQiwY3vocZS3zIrrCbp27cq+fc7LbHv37s2qVatYuXIlF110UZO3HR0dzebNmwHYvHkzAwcOrLJM9+7dOXHiBOnp6ZSWlrJt27aK1pDs7OyK5Xbs2EFERESV9YVokopzUE3vrFvnEdTYsWMr7v+6icFms3Hrrbc2OURrpEtKMJY/B4X52P74PCoouO6VhFv4wx/+UNHMNmHCBN555x0KCwuZOnVqk7c9ZswYFi9ezKZNmwgLC2PWrFmA87zTX//6V2bPno2XlxcTJkwgISEBwzAYOXJkRSF66623OHz4MEopwsPDmThxYpMzCXEh5d0G/Ns1y5xQdRaoV155Ba018+fP56mnnvolhFIEBQXRtm3bJodobbTW6LeWOSce/MP/Q0VGmR1JNKO8vLyK87VBQUFMmjQJcJ6vbarAwECefPLJKs+HhoYye/bsisf9+/evtn/TrwexFaJFBNihoOlHUHU28YWHh9O+fXuWLVtGeHh4xU9YWJgUp0bSG9ejt32GumksKnqY2XFEM3vmmWeqfT4hIcHFSYQwSTMNGFtrgVqzZk2N01eXy8nJqfbiCVE9vScJ/c83UAOuRo2+0+w4ohnJ+Vohzgu0t3wTX6dOnZg9ezZdunShT58+dOrUqaKj7okTJ9i3bx9paWncdtttTQ7SGjgnHnwJIqJQ98+QYYw8jJyvFcJJBdrRB39q8nZqLVCjRo1i5MiR7Ny5k127dpGUlERRURHt2rUjMjKSUaNGMWDAAPlmWA86Pw/jlWd+mXhQ5gvyOL8+X6u1Rikl52tF6xMQBAV5TR6Pr86LJLy9vRk8eDA//PADDzzwAD169Gj0m7VWuqQEY8VzzokHH3sWFWqNXuaieYWHhwOwbNkywNnkd+HID0K0GoF2MAw4U+icYbeR6t1RVynFiy++iI+PD8OGDWPYsGENmsK6ptGVy+3du5cXXniB9u3bAzBo0KCKHu/uTGtN3t9egv17UQ8+IhMPtgKFhYWsXLmS7du34+3tzdq1a9m5cyepqan1Hg1CCLcWaHfe5uc2qUDV+9hr/PjxLF++nAcffJCMjAz+9Kc/8cc//pF///vfda5bPrrynDlzWLx4MV999RXHj1ed0KpPnz68+OKLvPjiix5RnAD0Zx9RnPgv1A13yMSDrcRrr72Gv78/y5Ytw9vb+R2wV69ebNu2zeRkQriGqihQTbvUvEGNgzabjSuuuILJkyezaNEiAgMDWbt2bZ3rXTi6sre3d8Xoyp5O70lCr1uNz6AY1C13mx1HuMh3333H/fffX6lpLygoiNzcpl/VJIRbKB/FvIlX8jV4wsIdO3bw1VdfsW/fPvr27cuUKVPqXK+u0ZXL7d+/n8cee4yQkBDuu+++aodhcZdRlksO/kT2ypfwjupF6KNPY1hgGCOzP5MLeXIWf39/8vPzKxWojIwMORclWo8A5xGULmjaeHz1LlB//vOf2bVrF1FRUVx99dVMmTKl2oEqq1Pb6MrlunXrxrJly/D19SU5OZkXX3yRpUuXVlnPHUZZ1lkZGM89Cv7tMCY9juHdxhKjG7f2UZZr0phRlmsTGxvLokWLuOuuu9Bas3//fv7nf/6HUaNGNTWqEO6hmZr46l2gyueOacw3zdpGVy5XPvozOIdpWbVqFXl5efUuglahi4swXl4AxWecY+wFh5odSbjYLbfcQps2bVi1ahVlZWUsX76cuLg4brhBplIRrYNq0wZ8/VzXxPfrq+4a4sLRlUNDQ9m2bRvTp0+vtExOTg52ux2lFKmpqRiGQWBg46/+MIMuK8P420uQdgTbtCdRXS42O5IwgVKKG2+8kRtvvNHsKEKYpxlGk2jQOajGqml05Y0bNwLOSdS2b9/Oxo0b8fLyom3btsyYMcOt5kbSWqP//hp8txN132TUZVUH6hStR1paGocPH6a4uLjS89dee61JiYRwsYCmj8fnkgIF1Y+uHB8fX3H/uuuu47rrrnNVnGanP/sI/cXHqPhbsV3jvr+HaLr333+ff/7zn3Tt2hWfX40YIgVKtBqBdshq2nlmlxUoT6Z3b0evex36D0H91+/MjiNM9vHHH/Pss8/StWtXs6MIYRoVaEcfPdCkbchopU2kD/7kHAC2aw9sE2bJALCCtm3b0rlzZ7NjCGGuQDvk51V7FXd9yf+mTaBP/ozx8tNgD8U2ba4MANuKXTi1xp133snrr79OdnZ2pecNwzA7phCuExgEZaVwpqjRm5AmvkbSudkYS+aBsmGbMR8VJJ0wW7MLp9oo99lnn1V57t1333VFHCHMd76zLgW5zingG0EKVCPoM0UYS5+C/Fxsjz6Lat/wzpzCs7zyyiuA82rO7du3M2TIkEqva635z3/+Y0Y0IUyhAu1ogLxcaOT/kdLE10C6tARj+XNw/DC2SY+juvU0O5KwgPDwcMLDw2nfvj3//Oc/Kx5f+Pz7779vdkwhXKd8PL6Cxl9qLkdQDaANA/3GUvhhD2r8w6jLB5gdSVjI999/D0BZWVnF/XKnTp3Cz8/PjFhCmOP8cEc6P6/R4/FJgWoA/f4a9H82o8bci+3qWLPjCItZvnw5ACUlJRX3wTmyRHBwMBMmTDArmhCuF3DBnFCNJAWqnozEj9AbPkCNuAF1w3+bHUdY0Kuvvgo4z0dNnTrV5DRCmEv5+ICPb5MGjJVzUPVgbP8c/e5K6DcYNfb3bjUEk3A9KU5CnBcQJOegWpLeswO9+i9wyeXYfv8oyuZldiTRihUUFLB48WJOnz5NeHg4M2fOJCAgoMpyy5YtIzk5GbvdzqJFixq8vhDNItDepPH45AiqFvqn7zBWPA+R3bFN/ROqTVuzI4lWbv369Vx++eUsXbqUyy+/nPXr11e73IgRI5gzZ06j1xeiWZwfTaKxpEDVQB9JxXjlGQjviG36PJSvf90rCdHCkpKSiImJASAmJoakpKRql+vbt2+1R0b1XV+I5qACgpp0kYQUqGroE8cxlsyHdoHYZjyFCnSvSROF58rNza2Y7DMkJIS8vIZ9O23q+kI0yPk5oRo7Hp+cg/oVnZmOsfhJUArbzKdRoQ2fQViIpliwYAE5OTlVnr/rrrtcliExMZHExEQAFi5cWONM2t7e3o2aZbslWCWLVXKA+VkKO15EQWkJjgD/RmWRAnUBnZeD8ecnndO1P5qA6iBDGAnXe+KJJ2p8zW63k52dTUhICNnZ2QQFNezovr7rx8XFERcXV/E4I6P6eX3CwsJqfM3VrJLFKjnA/CyGzVliMg8fIrzPZTVm6dSp+v9rpYnvPF1Y4Bz8NScD27QnUJFRZkcSooro6Gg2b94MwObNmxk4cKBL1xeiIVRg0zrrSoECdPH5wV/TjmF7aDaqZ1+zIwlRrTFjxvDtt98yffp0vv32W8aMGQNAVlYWzz33XMVyS5YsYe7cuaSlpTFp0iQ2bdpU6/pCtIiKAtW4c52tvolPnzuL8UoCHE7B9oc/oi6T8fWEdQUGBvLkk09WeT40NJTZs2dXPJ4xY0aD1heiRQQ4m5B1IzvrtuojKF1SgrF8Iez/HnX/DFT/IXWvJIQQon6kia9xdFkZxsqX4PtvUPdOxjZ4hNmRhBDCs/j4Qpu2UqAawjltxl8g+WvUnQ9gu+a3ZkcSQgiPo5RyzgslBap+tNbot1egt3/hnDYj7hazIwkhhOcKsKMbeZFEqypQWmv0e6+jt/wf6vrbsd14h9mRhBDCs8kRVP3o9W+jP/0Qde1o1K33mR1HCCE8ngq0Q4EcQdXK+N916I/XoYbHo+58UOZ0EkIIVzg/Hl9jtIoCZXz6IXr9W6jBI1D3PoSytYpfWwghzBdoh3Nn0WeLG7yqx/9PbXzxCXrdKhgwFDX+YZlwUAghXOl8Z10jN7vBq7psJIndu3ezevVqDMMgNja2yhArWmtWr17Nrl278PHxYfLkyURFNW08PGPbZ+i3l8MVA7E9+AjKS4qTEEK4kgq0owEjLweCwxu0rkuOoAzDYNWqVcyZM4fFixfz1Vdfcfz48UrL7Nq1i5MnT7J06VImTpzIypUrm/SexV8mot94GfpeiW3SH1HebZq0PSGEEI1wfjSJxhxBuaRApaam0rFjRzp06IC3tzdDhw6tMpPnzp07ueaaa1BK0atXLwoLC8nObvgvBKB3bSd38VPQsw+2yTJVuxBCmCawvImv6hxndXFJE19WVhYOh6PiscPhICUlpcoyF05m5XA4yMrKqpj9s1x9JlLLO5JCac8+BM9bjM2vXXP+Ko1i9qRhVssBkkWIViMw2HkeyjAavKpLClR10/3++jLv+iwD9ZtITd9yL2GBAWQWFELhmcZEblZmTxpmtRzgPllqmkhNCFE/ys8fr8Vv4RcWRmED93mXNPE5HA4yMzMrHmdmZlY5MnI4HJX+k6humfpSSqF8/RoXVgghhCW4pEB1796dEydOkJ6eTmlpKdu2bSM6OrrSMtHR0WzZsgWtNfv378ff37/RBUoIIYT7c0kTn5eXFxMmTCAhIQHDMBg5ciQRERFs3LgRgPj4ePr160dycjLTp0+nbdu2TJ482RXRhBBCWJTS1Z38EUIIIUzmsSNJPP7442ZHqGCVLFbJAZLFU1jps7NKFqvkAPfP4rEFSgghhHuTAiWEEMKSvObPnz/f7BAtpalj+TUnq2SxSg6QLJ7CSp+dVbJYJQe4dxa5SEIIIYQlSROfEEIIS5ICJYQQwpJcNh+UK9U195SrTJkyBV9fX2w2G15eXixcuNBl771s2TKSk5Ox2+0sWrQIgIKCAhYvXszp06cJDw9n5syZBAQEmJJl3bp1fPbZZwQFOUc6Hjt2LP3792/RHBkZGbz66qvk5OSglCIuLo4bbrjBtM/F3cl+JvtZdZp1P9MepqysTE+dOlWfPHlSl5SU6EcffVQfO3bMlCyTJ0/Wubm5prz33r179YEDB/SsWbMqnlu7dq3+4IMPtNZaf/DBB3rt2rWmZXn33Xf1hx9+6JL3L5eVlaUPHDigtda6qKhIT58+XR87dsy0z8WdyX7mJPtZVc25n3lcE1995p5qDfr27Vvl20lSUhIxMTEAxMTEuOxzqS6LGUJCQiquIvLz86Nz585kZWWZ9rm4M9nPnGQ/q6o59zOPa+Krz9xTrpSQkADAqFGjKk0TYobc3NyKAXhDQkLIy8szNc+GDRvYsmULUVFRjBs3zqU7V3p6OocOHaJHjx6W+1zcgexnNbPa35M772ceV6B0PeeVcoUFCxYQGhpKbm4uzzzzDJ06daJv376mZLGa+Ph4br/9dgDeffdd3nzzTZcNEFxcXMyiRYsYP348/v7+LnlPTyP7mXtw9/3M45r46jP3lKuEhoYCYLfbGThwIKmpqabkKGe328nOzgYgOzu74sSpGYKDg7HZbNhsNmJjYzlw4IBL3re0tJRFixYxfPhwBg0aBFjrc3EXsp/VzEp/T+6+n3lcgarP3FOuUFxczJkzZyruf/vtt0RGRro8x4Wio6PZvHkzAJs3b2bgwIGmZSn/QwXYsWMHERERLf6eWmtWrFhB586dGT16dMXzVvpc3IXsZzWz0t+Tu+9nHjmSRHJyMmvWrKmYe+q2225zeYZTp07x0ksvAVBWVsawYcNcmmPJkiXs27eP/Px87HY7d9xxBwMHDmTx4sVkZGQQFhbGrFmzXNIeXV2WvXv3cvjwYZRShIeHM3HixBb/Bv7jjz/y5JNPEhkZWdEcNXbsWHr27GnK5+LuZD+T/aw6zbmfeWSBEkII4f48rolPCCGEZ5ACJYQQwpKkQAkhhLAkKVBCCCEsSQqUEEIIS5ICJYQQwpKkQAkhhLCk/w/om/HOnlJ+igAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 4 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"x_mpc=np.array(x.value[0, :]).flatten()\n",
|
|
"y_mpc=np.array(x.value[1, :]).flatten()\n",
|
|
"v_mpc=np.array(x.value[2, :]).flatten()\n",
|
|
"theta_mpc=np.array(x.value[3, :]).flatten()\n",
|
|
"a_mpc=np.array(u.value[0, :]).flatten()\n",
|
|
"w_mpc=np.array(u.value[1, :]).flatten()\n",
|
|
"\n",
|
|
"#simulate robot state trajectory for optimized U\n",
|
|
"x_traj=predict(x0, np.vstack((a_mpc,w_mpc)))\n",
|
|
"\n",
|
|
"#plt.figure(figsize=(15,10))\n",
|
|
"#plot trajectory\n",
|
|
"plt.subplot(2, 2, 1)\n",
|
|
"plt.plot(track[0,:],track[1,:],\"b+\")\n",
|
|
"plt.plot(x_traj[0,:],x_traj[1,:])\n",
|
|
"plt.axis(\"equal\")\n",
|
|
"plt.ylabel('y')\n",
|
|
"plt.xlabel('x')\n",
|
|
"\n",
|
|
"#plot v(t)\n",
|
|
"plt.subplot(2, 2, 2)\n",
|
|
"plt.plot(a_mpc)\n",
|
|
"plt.ylabel('a(t)')\n",
|
|
"#plt.xlabel('time')\n",
|
|
"\n",
|
|
"\n",
|
|
"plt.subplot(2, 2, 4)\n",
|
|
"plt.plot(theta_mpc)\n",
|
|
"plt.ylabel('theta(t)')\n",
|
|
"\n",
|
|
"plt.subplot(2, 2, 3)\n",
|
|
"plt.plot(v_mpc)\n",
|
|
"plt.ylabel('v(t)')\n",
|
|
"\n",
|
|
"plt.tight_layout()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"full track demo"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 18,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"CVXPY Optimization Time: Avrg: 0.1990s Max: 0.3323s Min: 0.1631s\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"track = compute_path_from_wp([0,3,4,6,10,13],\n",
|
|
" [0,0,2,4,3,3],0.5)\n",
|
|
"\n",
|
|
"# track = compute_path_from_wp([0,5,7.5,10,12,13,13,10],\n",
|
|
"# [0,0,2.5,2.5,0,0,5,10],0.5)\n",
|
|
"\n",
|
|
"sim_duration = 80 #time steps\n",
|
|
"opt_time=[]\n",
|
|
"\n",
|
|
"x_sim = np.zeros((N,sim_duration))\n",
|
|
"u_sim = np.zeros((M,sim_duration-1))\n",
|
|
"\n",
|
|
"MAX_SPEED = 1.25\n",
|
|
"MIN_SPEED = 0.75\n",
|
|
"MAX_STEER = 1.57\n",
|
|
"MAX_ACC = 1.0\n",
|
|
"\n",
|
|
"# Starting Condition\n",
|
|
"x0 = np.zeros(N)\n",
|
|
"x0[0] = 0\n",
|
|
"x0[1] = -0.25\n",
|
|
"x0[2] = 0\n",
|
|
"x0[3] = np.radians(-0)\n",
|
|
"x_sim[:,0]=x0\n",
|
|
" \n",
|
|
"#starting guess\n",
|
|
"u_bar = np.zeros((M,T))\n",
|
|
"u_bar[0,:]=0.5\n",
|
|
"u_bar[1,:]=0.01\n",
|
|
"\n",
|
|
"for sim_time in range(sim_duration-1):\n",
|
|
" \n",
|
|
" iter_start=time.time()\n",
|
|
" K=road_curve(x_sim[:,sim_time],track)\n",
|
|
" \n",
|
|
" # dynamics starting state w.r.t vehicle frame\n",
|
|
" x_bar=np.zeros((N,T+1))\n",
|
|
" x_bar[2,0]=x_sim[2,sim_time]\n",
|
|
" \n",
|
|
" #prediction for linearization of costrains\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",
|
|
" A,B,C=get_linear_model(xt,ut)\n",
|
|
" xt_plus_one = np.squeeze(np.dot(A,xt)+np.dot(B,ut)+C)\n",
|
|
" x_bar[:,t]= xt_plus_one\n",
|
|
" \n",
|
|
" #CVXPY Linear MPC problem statement\n",
|
|
" cost = 0\n",
|
|
" constr = []\n",
|
|
" x = cp.Variable((N, T+1))\n",
|
|
" u = cp.Variable((M, T))\n",
|
|
" \n",
|
|
" #Prediction Horizon\n",
|
|
" for t in range(T):\n",
|
|
"\n",
|
|
" #cost += 30*cp.sum_squares(x[2,t]-np.arctan(df(x_bar[0,t],K))) # psi\n",
|
|
" cost += 50*cp.sum_squares(x[3,t]-np.arctan(df(x_bar[0,t],K))) # psi\n",
|
|
" cost += 20*cp.sum_squares(f(x_bar[0,t],K)-x[1,t]) # cte\n",
|
|
"\n",
|
|
" # Actuation rate of change\n",
|
|
" if t < (T - 1):\n",
|
|
" cost += cp.quad_form(u[:, t + 1] - u[:, t], 100*np.eye(M))\n",
|
|
" \n",
|
|
" # Actuation effort\n",
|
|
" cost += cp.quad_form( u[:, t],1*np.eye(M))\n",
|
|
" \n",
|
|
" # Kinrmatics Constrains (Linearized model)\n",
|
|
" A,B,C=get_linear_model(x_bar[:,t],u_bar[:,t])\n",
|
|
" constr += [x[:,t+1] == A@x[:,t] + B@u[:,t] + C.flatten()]\n",
|
|
"\n",
|
|
" # sums problem objectives and concatenates constraints.\n",
|
|
" constr += [x[:,0] == x_bar[:,0]] #<--watch out the start condition\n",
|
|
" constr += [x[2, :] <= MAX_SPEED]\n",
|
|
" constr += [x[2, :] >= 0.0]\n",
|
|
" constr += [cp.abs(u[0, :]) <= MAX_ACC]\n",
|
|
" constr += [cp.abs(u[1, :]) <= MAX_STEER]\n",
|
|
" \n",
|
|
" # Solve\n",
|
|
" prob = cp.Problem(cp.Minimize(cost), constr)\n",
|
|
" solution = prob.solve(solver=cp.OSQP, verbose=False)\n",
|
|
" \n",
|
|
" #retrieved optimized U and assign to u_bar to linearize in next step\n",
|
|
" u_bar=np.vstack((np.array(u.value[0, :]).flatten(),\n",
|
|
" (np.array(u.value[1, :]).flatten())))\n",
|
|
" \n",
|
|
" u_sim[:,sim_time] = u_bar[:,0]\n",
|
|
" \n",
|
|
" # Measure elpased time to get results from cvxpy\n",
|
|
" opt_time.append(time.time()-iter_start)\n",
|
|
" \n",
|
|
" # move simulation to t+1\n",
|
|
" tspan = [0,dt]\n",
|
|
" x_sim[:,sim_time+1] = odeint(kinematics_model,\n",
|
|
" x_sim[:,sim_time],\n",
|
|
" tspan,\n",
|
|
" args=(u_bar[:,0],))[1]\n",
|
|
" \n",
|
|
"print(\"CVXPY Optimization Time: Avrg: {:.4f}s Max: {:.4f}s Min: {:.4f}s\".format(np.mean(opt_time),\n",
|
|
" np.max(opt_time),\n",
|
|
" np.min(opt_time))) "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 19,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABC8AAALICAYAAABfINo9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzde3xcdYH///eZmcxMkuY2M72YtAgtpdCVWyhQC6WWhqp8USMoIqLW/vwpUC0gv/5sFZQVq/0JadkKXfRHt/jddRVX1/pdFtYaQcBWJND2C6UsNFAL9EKae3OZSZPz+f4xydCQSZpJZs6ZJK/n49HHZCZnznn3NA3Nm8/FMsYYAQAAAAAAZCmP2wEAAAAAAACGQnkBAAAAAACyGuUFAAAAAADIapQXAAAAAAAgq1FeAAAAAACArEZ5AQAAAAAAsprP7QCjcejQIcevGYlEVF9f7/h1Jyrut7O4387jnjuL++0s7rezuN+pKS0tdTuC41L9t3M2fU2RJTmyJEeW5MZKlsG+P4/p8gIAAAAYT3bv3q0tW7bItm0tWbJElZWV/T7/8ssv60c/+pGmTJkiSbr44ov1qU99yo2oAOAoygsAAAAgC9i2rc2bN+uOO+5QOBzWmjVrNG/ePE2fPr3fcWeddZZWr17tUkoAcAdrXgAAAABZoLa2VtOmTdPUqVPl8/m0YMEC1dTUuB0LALIC5QUAAACQBRobGxUOhxPPw+GwGhsbBxz32muvadWqVfrBD36gt956y8mIAOAapo0AAAAAWcAYM+A1y7L6PT/ttNO0adMmBYNB7dy5U/fcc482btyY9HzV1dWqrq6WJK1bt06RSCSlPD6fL+X3ZApZkiNLcmRJbqxnobwAAAAAskA4HFZDQ0PieUNDg0pKSvodk5eXl/i4vLxcmzdvVmtrqwoLCwecr6KiQhUVFYnnqe4yMFZ2JnAaWZIjS3JkSW4ku40wbQQAAADIArNmzdLhw4dVV1en7u5u7dixQ/Pmzet3THNzc2KERm1trWzbVkFBgRtxAcBRjLwAAAAAsoDX69Xy5cu1du1a2batxYsXa8aMGdq2bZskaenSpXr22We1bds2eb1e+f1+3XrrrQOmlgDAeER5AQAAAGSJ8vJylZeX93tt6dKliY8/8pGP6CMf+YjTsQDAdUwbAQAAAAAAWY3yAgAAAAAAZDXKCwAAAAAAkNUoLwAAAAAMytg96nlwnczr/+12FAATGAt2AgAAABhcW6v0wg6ZgiJZs850Ow2ACYqRFwAAAAAGF+2UJJk3XnM5CICJjPICAAAAwOCi0fjjwb/JdMXczQJgwqK8AAAAADC4WG950dMjvfm6u1kATFiUFwAAAAAG1zttRJLM/n0uBgEwkVFeAAAAABhcrLe8sDzSfta9AOAOygsAAAAAgzJ900beP0vmjVfdDQNgwqK8AAAAADC43gU7rbPOkRrqZFqbXQ4EYCKivAAAAAAwuGiHJMk689z4c6aOAHAB5QUAAACAwcWiktcrzTpL8nhk3qC8AOA8ygsAAAAAg4tFpUCurEBAKnu/zH7WvQDgPMoLAAAAAIOLdkrBoCTJmjlH+ts+Gdt2ORSAiYbyAgAAAMCgTKxTCuTGn5w2R+rskN456G4oABMO5QUAAACAwcWiUqBv5MUZksS6FwAcR3kBAAAAYHDRTinYO/JiapmUmy+x7gUAh1FeAAAAABhc9ISRFx6PdOrpMmyXCsBhlBcAAAAABhfrlNW35oUka87Z0lv7ZeoOuRgKwERDeQEAAABgcCdOG5FkXVIheTwyTz7mYigAEw3lBQAAAIDBxaKJrVIlySoOybrgUpnt1TLRDheDAZhIKC8AAAAAJGXsHqkrlljzoo+15Cqps0PmL0+6lAzAREN5AQAAACC5WCz+eMK0EUmyZs6RTjtD5olHZWzbhWAAJhrKCwAAAADJxTrjj4HcAZ+yLr9KOnJQ2rvb4VAAJiLKCwAAAADJRaPxx/dMG5Eka94lUlGJ7CcedTgUgImI8gIAAABAcr0jL6xgkpEXvhxZl31Eeul5mXfYNhVAZlFeAAAAAEhuiJEXkmQtXCpJMv/7r04lAjBBUV4AAAAASK5vzYskIy8kScUhKccvtTQ7lwnAhER5AQAAACApEzvJyAvLkgqLpdYmB1MBmIgoLwAAAAAkFz3JyAtJKiyWaWXkBYDMorwAAAAAkNwQW6UmFBZLlBcAMozyAgAAAEByJ1mwU5KswmKphWkjADKL8gIAAADIErt379Ytt9yir3/969q6deugx9XW1uozn/mMnn322cwGinZKvhxZPt/gxxSVSG3HZOyezGYBMKEN8V3IWY8++qieeOIJWZalGTNm6Oabb5bf73c7FgAAAOAI27a1efNm3XHHHQqHw1qzZo3mzZun6dOnDzju5z//uc4777zMh4pFpeDgoy4kxaeNGFs61hovMgAgA7Ji5EVjY6Mef/xxrVu3TlVVVbJtWzt27HA7FgAAAOCY2tpaTZs2TVOnTpXP59OCBQtUU1Mz4LjHH39cF198sQoLCzMfKtY59HoXkqzC3sKCdS8AZFBWlBdSvEHu6upST0+Purq6VFJCawsAAICJo7GxUeFwOPE8HA6rsbFxwDHPPfecli5d6kgmE+0ceqcRKT7yQqK8AJBRWTFtJBQK6WMf+5huuukm+f1+nXvuuTr33HMHHFddXa3q6mpJ0rp16xSJRJyOKp/P58p1Jyrut7O4387jnjuL++0s7rezuN9jnzFmwGuWZfV7/vDDD+tzn/ucPJ6T/z/I0f7b2efzyW/3yORPUmiI93Z3naYGSZPsbuVm6Gswm76+yZIcWZIjS3IjyZIV5UVbW5tqamr0wAMPKC8vT+vXr9fTTz+tyy67rN9xFRUVqqioSDyvr693OqoikYgr152ouN/O4n47j3vuLO63s7jfzuJ+p6a0tNTtCAOEw2E1NDQknjc0NAwYjfz666/rH/7hHyRJra2t2rVrlzwejy666KIB5xvtv50jkYi62o5J/sCQ7zV2vHQ5dugttWfoazCbvr7JkhxZkiNLckNlGez7c1aUFy+99JKmTJmSmLd38cUX67XXXhtQXgAAAADj1axZs3T48GHV1dUpFAppx44dWrlyZb9jHnjggX4fX3DBBUmLi7SJdr47LWQwgVzJ72faCICMyoryIhKJaN++fYrFYvL7/XrppZc0a9Yst2MBAAAAjvF6vVq+fLnWrl0r27a1ePFizZgxQ9u2bZMkx9a56CfaKetkC3ZallRYIrU0ORQKwESUFeXF7NmzNX/+fH3zm9+U1+vVqaee2m+IGwAAADARlJeXq7y8vN9rg5UWK1asyHyg4WyVKkmFxTKMvACQQVlRXkjStddeq2uvvdbtGAAAAAD6RDulwPDKCx09kvk8ACasrNkqFQAAAED2MN3dUvfxk2+VKskqLGHNCwAZRXkBAAAAYAAT64x/cJI1LyTFR160tcr09GQ2FIAJi/ICAAAAwACms6+8GMa0kaJiyRjpWEtmQwGYsCgvAAAAAAxgOjviHwxr2kjvdqpMHQGQIZQXAAAAAAYw0Xh5cbKtUiXFp41IUivbpQLIDMoLAAAAAAOYaO+0kWFtlVoSfw8jLwBkCOUFAAAAgAHsFKaNiGkjADKM8gIAAADAAImRF8NYsNMK5kr+gNRCeQEgMygvAAAAAAzQt+bFsLZKlaSiEkZeAMgYygsAAAAAA6Sy24gkqbBY5hjlBYDMoLwAAAAAMMC700YCw3tDYbHUwm4jADKD8gIAAADAAKazU/L7ZXm8wzreKixm2giAjKG8AAAAADCA6ewY/noXUnzkRVurTHd35kIBmLAoLwAAAAAMYKIdw1/vQpIKS+KPbS2ZCQRgQqO8AAAAADCAiXYOa5vUPlZhcfwDpo4AyADKCwAAAAADxKeNDL+8UFHvyIsWygsA6Ud5AQAAAGAAO9qZ4rSR+MgLw8gLABlAeQEAAABggPi0kdTLC6aNAMgEygsAAAAAA5hoh6xU1rwIBONlR2tTBlMBmKgoLwAAAAAMYDpT3G1EkgqLGHkBICMoLwAAAAAMYKKdUjCFBTslqahEpoWRFwDSj/ICAAAAQD+m+7jU3Z3amhdSfN0LRl4AyADKCwAAAAD9RTvjjylOG7EKSygvAGQE5QUAAACA/mLR+GMKC3ZKkoqKpfZjMsePpz8TgAmN8gIAAABAf9G+8iLFaSNFofgjO44ASDPKCwAAAAD9RTskSVaKC3Zaxb3lRXNjuhMBmOAoLwAAAAD0FxvlyIsWygsA6UV5AQAAAKC/vvIi1a1Si0skie1SAaQd5QUAAACAfkzfbiOpjryYVCR5PEwbAZB2lBcAAAAA+ouNcKtUj0cqLGHaCIC0o7wAAAAA0N9It0qVpKISpo0ASDvKCwAAAAD99U0b8QdSf29xiGkjANKO8gIAAABAf9GorGBefBpIiqyikMTICwBp5nM7AAAAAIDsYl3wQeXPPlMdI3lzUYl0rEWm+7gsX066owGYoCgvAAAAgCyxe/dubdmyRbZta8mSJaqsrOz3+ZqaGj3yyCOyLEter1fLli3TmWeemfYc1ulzlReJqKO+PvU3F4fij63NUmhyeoMBmLAoLwAAAIAsYNu2Nm/erDvuuEPhcFhr1qzRvHnzNH369MQxZ599tubNmyfLsnTgwAFt2LBB9913n4upB7KKQjJSfOoI5QWANGHNCwAAACAL1NbWatq0aZo6dap8Pp8WLFigmpqafscEg0FZliVJisViiY+zSt/ICxbtBJBGjLwAAAAAskBjY6PC4XDieTgc1r59+wYc99xzz+lf//Vf1dLSojVr1gx6vurqalVXV0uS1q1bp0gkklIen8+X8nskqccj1UvK7+lS3gjen84smUCW5MiSHFmSG0kWygsAAAAgCxhjBryWbGTFRRddpIsuukh79+7VI488ojvvvDPp+SoqKlRRUZF4Xp/i+hWRSCTl90iSsXsky6O2t98c2ZoZacySCWRJjizJkSW5obKUlpYmfZ1pIwAAAEAWCIfDamhoSDxvaGhQSUnJoMfPnTtXR44cUWtrqxPxhs3yeKXCYrZLBZBWlBcAAABAFpg1a5YOHz6suro6dXd3a8eOHZo3b16/Y44cOZIYofHGG2+ou7tbBQUFbsQdWlGJDOUFgDRi2ggAAACQBbxer5YvX661a9fKtm0tXrxYM2bM0LZt2yRJS5cu1bPPPqunn35aXq9Xfr9ft912W/Yu2tmUHcPTAYwPlBcAAABAligvL1d5eXm/15YuXZr4uLKyUpWVlU7HSplVVCLzt4GLjQLASDFtBAAAAEB6FYWkYy0yPT1uJwEwTlBeAAAAAEiv4pBkjNTa7HYSAOME5QUAAACAtLKKendJaWl0NwiAcYPyAgAAAEB6FYfij82UFwDSg/ICAAAAQHoVxcsLtksFkC6UFwAAAADSq7BYsixGXgBIG8oLAAAAAGlleb1SQRFrXgBIG8oLAAAAAOlXVMK0EQBpQ3kBAAAAIP2Kw0wbAZA2lBcAAAAA0s4qKpEYeQEgTSgvAAAAAKRfUYnU2ixj97idBMA4QHkBAAAAIP2KQ5KxpdYWt5MAGAcoLwAAAACknVUUin/A1BEAaUB5AQAAACD9ikrijyzaCSANKC8AAAAApF9xfOSFaaG8ADB6lBcAAAAA0i+/IP7Y2e5uDgDjAuUFAAAAgPTL8ccfo1F3cwAYFygvAAAAAKSd5fFIgaAU63Q7CoBxgPICAAAAQGYEglIs5nYKAOMA5QUAAACAzGDkBYA0obwAAAAAkBmBoEyMNS8AjB7lBQAAAIDMCAQlygsAaUB5AQAAACAzKC8ApAnlBQAAAIDMoLwAkCaUFwAAAAAywgrkUl4ASAvKCwAAAACZEQhQXgBIC8oLAAAAAJnByAsAaUJ5AQAAACAzAgGpKyZj97idBMAYR3kBAAAAIDMCufHHrpi7OQCMeZQXAAAAADIjEIw/xigvAIwO5QUAAACAzEiUF53u5gAw5lFeAAAAAMgIq6+8iLJoJ4DRobwAAAAAkBnB3vKii/ICwOhQXgAAAADIDD8jLwCkh8/tAAAAAADidu/erS1btsi2bS1ZskSVlZX9Pv/MM8/od7/7nSQpGAzqy1/+sk499VQXkg4TIy8ApAkjLwAAAIAsYNu2Nm/erG9961vasGGDtm/frrfffrvfMVOmTNFdd92le++9V9dcc41++tOfupR2mHq3SjWMvAAwSoy8AAAAAIZw//33D+s4n8+nG2+8ccTXqa2t1bRp0zR16lRJ0oIFC1RTU6Pp06cnjpkzZ07i49mzZ6uhoWHE13NEIBB/jFFeABgdygsAAABgCDt27NAnP/nJkx736KOPjqq8aGxsVDgcTjwPh8Pat2/foMc/8cQTOv/88wf9fHV1taqrqyVJ69atUyQSSSmPz+dL+T3vZSblq05Svs+r/FGcKx1Z0oUsyZElObIkN5IslBcAAADAEMLhsD796U+f9Ljt27eP6jrGmAGvWZaV9Ng9e/boySef1Pe+971Bz1dRUaGKiorE8/r6+pTyRCKRlN/zXsa2JctSe2ODOkdxrnRkSReyJEeW5MiS3FBZSktLk77OmhcAAADAEH784x8P67j77rtvVNcJh8P9poE0NDSopKRkwHEHDhzQT37yE61atUoFBQWjumamWR6P5A9IsU63owAY47Jm5EV7e7sefPBBvfXWW7IsSzfddJPOOOMMt2MBAAAAg3rnnXfk8Xg0efLkUZ9r1qxZOnz4sOrq6hQKhbRjxw6tXLmy3zH19fW699579bWvfW3Q/zuZdQJBKRZzOwWAMS5ryostW7bovPPO0+23367u7m7F+AYHAACALHPffffpox/9qObMmaMnn3xSDz30kDwej770pS/p8ssvH9W5vV6vli9frrVr18q2bS1evFgzZszQtm3bJElLly7Vr3/9a7W1temhhx5KvGfdunWj/n1lVCDIyAsAo5YV5UVHR4deeeUVrVixQlJ88Q6fLyuiAQAAAAl79uzR1772NUnxBTrvvPNO5efn65577hl1eSFJ5eXlKi8v7/fa0qVLEx/feOONo1oU1BWBoAy7jQAYpaxoCOrq6lRYWKhNmzbpwIEDmjlzppYtW6ZgMOh2NAAAACChu7tbPp9PjY2Namtr05lnnilJamlpcTlZFgsE2SoVwKhlRXnR09Oj/fv3a/ny5Zo9e7a2bNmirVu36rrrrut33Gi3e0qHbNpeZiLgfjuL++087rmzuN/O4n47i/vtjFNPPVW//e1vdfTo0cQIicbGRuXm5rqcLIsFglKUaSMARicryotwOKxwOKzZs2dLkubPn6+tW7cOOG602z2lQzZtLzMRcL+dxf12HvfcWdxvZ3G/ncX9Ts1IF7u88cYb9cgjj8jr9erzn/+8JOm1117TpZdems5440sgKLU0uZ0CwBiXFeVFcXGxwuGwDh06pNLSUr300kuaPn2627EAAAAASdIf//hHnX/++Zo2bZpuueWWfp+bP3++5s+f71Ky7GcFclnzAsCoZUV5IUnLly/Xxo0b1d3drSlTpujmm292OxIAAAAgSXr99df1m9/8Rvn5+SovL9f555+vOXPmyLIst6Nlv0CANS8AjFrWlBennnpq9m/zBAAAgAnpK1/5iiTpzTff1M6dO/WLX/xChw4d0gc+8AGdf/75Ou+881RYWOhyyiwVyKW8ADBqWVNeAAAAANnulFNO0SmnnKLKykp1dHRo9+7d2rlzp37+858rEono05/+tM477zy3Y2aXQFDqisnYPbI8XrfTABijKC8AAACAEcjLy9OCBQu0YMECSVJtba3LibJUIBh/7IpJwTx3swAYsygvAAAAgBS88sor2r9/v6LR/lMhrr76apcSZbm+8iJGeQFg5CgvAAAAgGH6p3/6J/3lL3/RmWeeKb/fn3idhTuHkCgvOiWVuBoFwNhFeQEAAAAM0zPPPKOqqiqFQiG3o4wZViAoI0lRFu0EMHIetwMAAAAAY0UkElFOTo7bMcaWYN+aF5QXAEaOkRcAAADAMN144436yU9+oksuuURFRUX9Pjd37lyXUmU5f295wcgLAKNAeQEAAAAM0xtvvKFdu3bplVde6bfmhST94z/+o0upshwjLwCkAeUFAAAAMEy/+MUv9M1vflPnnHOO21HGjkCuJMlEo2JZUwAjxZoXAAAAwDAFAgGmh6QqEIg/xhh5AWDkKC8AAACAYfrMZz6jhx9+WM3NzbJtu98vDKJ35AXTRgCMBtNGAAAAgGHqW9fiD3/4w4DPPfLII07HGRty/JJlsWAngFGhvAAAAACG6f7773c7wphjeTySPyDFOt2OAmAMo7wAAAAAhmny5MluRxibAkEpFnM7BYAxjDUvAAAAgCH88pe/HNZxv/rVrzKcZAwLBBl5AWBUGHkBAAAADOGxxx7T5ZdfLmPMkMc9/vjjuvbaax1KNcYEgjLsNgJgFCgvAAAAgCHEYjF9/etfP+lxOTk5DqQZowJBtkoFMCqUFwAAAMAQ2EUkDQJBKcq0EQAjx5oXAAAAADKLkRcARonyAgAAAEBGWYFcygsAo0J5AQAAACCzAgHKCwCjQnkBAAAAILMYeQFglCgvAAAAgBREo1E1NDQoGuWH8WELBKWumIzd43YSAGMUu40AAAAAJ/Hmm2+qurpaO3fu1NGjRxOvT5kyReedd56uuOIKnXLKKS4mzHKBYPyxKyYF89zNAmBMorwAAAAAhnDffffp7bff1oIFC/T1r39dZWVlys3NVWdnpw4ePKi9e/dq48aNmj59um699Va342anvvIiRnkBYGQoLwAAAIAhXHrppZo3b96A1ydNmqQ5c+Zozpw5+uQnP6kXXnjBhXRjRLCvvOiUVOJqFABjE2teAAAAAEM4sbjYt29f0mNqa2t1wQUXOBVpzLH8J4y8AIARoLwAAAAAhun73/9+0tfXrl3rcJIxpt/ICwBIHdNGAAAAgJOwbVuSZIxJ/OrzzjvvyOv1puU6u3fv1pYtW2TbtpYsWaLKysp+nz948KA2bdqk/fv367rrrtPHP/7xtFw34/pGXrBDC4ARorwAAAAATuKzn/1s4uPrrruu3+c8Ho8++clPjvoatm1r8+bNuuOOOxQOh7VmzRrNmzdP06dPTxwzadIkfelLX1JNTc2or+eovpEXXZQXAEaG8gIAAAA4ifvvv1/GGN111136+7//+8TrlmWpsLBQfr9/1Neora3VtGnTNHXqVEnSggULVFNT06+8KCoqUlFRkXbu3Dnq6zkqkCtJMtGoLJejABibKC8AAACAk5g8ebIkadOmTRm7RmNjo8LhcOJ5OBwedIHQ4aiurlZ1dbUkad26dYpEIim93+fzpfyewfT4PKqXNCnHq7wRnDOdWUaLLMmRJTmyJDeSLJQXAAAAwBB+9rOf6ROf+ISKi4sHPaa5uVm/+93v9MUvfnHE1zlxHY0+ljXycQoVFRWqqKhIPK+vr0/p/ZFIJOX3DMbE4tNF2hob1DGCc6Yzy2iRJTmyJEeW5IbKUlpamvR1ygsAAABgCKWlpVqzZo2mT5+us846S6WlpcrNzVVnZ6cOHz6svXv36tChQ7r66qtHdZ1wOKyGhobE84aGBpWUlIw2fnbI8UuWxYKdAEaM8gIAAAAYwhVXXKHFixfr+eef165du1RTU6OOjg7l5+frlFNO0RVXXKELLrhg1DuOzJo1S4cPH1ZdXZ1CoZB27NihlStXpul34S7L45H8AbZKBTBilBcAAADASfh8Ps2fP1/z58/P2DW8Xq+WL1+utWvXyrZtLV68WDNmzNC2bdskSUuXLlVzc7NWr16tzs5OWZalxx57TOvXr1deXl7GcqVNICjFYm6nADBGUV4AAAAAw/Twww/r0ksv1emnn56R85eXl6u8vLzfa0uXLk18XFxcrAcffDAj1864QJCRFwBGjPICAAAAGCZjjO655x4FAgFdeumluvTSSwddXA7vEQjKMPICwAhRXgAAAADD9KUvfUlf/OIXtWfPHv35z3/Wt7/9bU2ZMkULFy7UVVdd5Xa87MbICwCj4HE7AAAAADCWeDwenXPOObr55ptVVVWlgoIC/fM//7PbsbJfIFeKsdsIgJFh5AUAAACQgmg0queee07bt2/X3r17NXfuXK1YscLtWNkvEJBaGt1OAWCMorwAAAAAhmn9+vXatWuXZs6cqUsuuUQrVqxQYWGh27HGBCuQK8PICwAjRHkBAAAADNPMmTP1hS98QZFIxO0oY08gwLQRACNGeQEAAAAMU2VlpdsRxq5gLgt2AhgxFuwEAAAAkHmBXKmrS8bucTsJgDGI8gIAAABA5gWC8cdYzN0cAMYkygsAAAAAmRfMjT9GmToCIHWUFwAAAAAyLzHygvICQOooLwAAAABknNU38oIdRwCMAOUFAAAAgMxj2giAUaC8AAAAAJB5fdNGooy8AJA6ygsAAAAAmReIj7wwrHkBYAQoLwAAAABkHtNGAIwC5QUAAACAzEvsNsK0EQCpo7wAAAAAkHmJNS8YeQEgdZQXAAAAADLO8nolv19izQsAI0B5AQAAAMAZgVymjQAYEcoLAAAAAM4IBJk2AmBEKC8AAAAAOCOYKxNl5AWA1FFeAAAAAHBGIMiaFwBGhPICAAAAgDMCuUwbATAilBcAAAAAnBFkwU4AI0N5AQAAAMARViBIeQFgRCgvAAAAADgjyG4jAEaG8gIAAACAM4K5LNgJYEQoLwAAAAA4I5ArdXfLdB93OwmAMYbyAgAAAIAzAsH4I+teAEgR5QUAAAAAZwRz449RygsAqaG8AAAAAOCMQG95wboXAFJEeQEAAADAEVawd9oIO44ASBHlBQAAAABnJEZeMG0EQGp8bgcAAAAAELd7925t2bJFtm1ryZIlqqys7Pd5Y4y2bNmiXbt2KRAI6Oabb9bMmTNdSjsCfSMvmDYCIEWMvAAAAACygG3b2rx5s771rW9pw4YN2r59u95+++1+x+zatUtHjhzRxo0b9ZWvfEUPPfSQS2lHqHfkhWHaCIAUMfICAAAAyAK1tbWaNm2apk6dKklasGCBampqNH369MQxzz//vC677DJZlqUzzjhD7e3tahTzj6kAACAASURBVGpqUklJiVuxU8NuIxOCMUY63iV1xeK/jh+XurulnhN/2ZLdI8U6ZVqbpdYWqf2YlJsnFZbIKiyOf2yMJCPZtkzbMelYi9TaHD9/jj/+y++XcgJSTo7k96szf5Lsdw5J7W1SR7tk7HfDeX1SMC9+7mAwnq+zU+psl7qP9/+N+AOSP9g7YsiSuqLxr92uWDxTH4+3f4bu41JXl3S8S605ObJPLOt8Oe8eZ1nx38fxrvg9OvGclkfyeuN5vd74verpid/HE38/fdfvO86y4ve3uzv+nhO0BnP7Z8m08BR5Pnx12k5HeQEAAABkgcbGRoXD4cTzcDisffv2DTgmEon0O6axsTFpeVFdXa3q6mpJ0rp16/q9bzh8Pl/K7zkZOz9PRyXl+7zKT+HcmcgyUuM1i+mKyW5tkd3WKtN+TKajXXZ7m0xbq+zeX6ajQ+Z4TOrqkjneJRONykQ7ZKJRHY1FZbpiMl2x+A/jxpz8ou9hBfNkYp2SMRry3T6fLH9ApqtrYOEgqTVxQktWbn78h/o+x4/LRDuSXDs3XixYVu8NMfHfS1es/4E5flmBQLxc6NPTEz/2xCxer6xAUDHvCT9yGyN1H48fa9uJjPIHZOXk9D+nbUs93TJ9pY/XK8uXEy8pPJ7+57R73j3OSPL1Huvxvvv7kfSe30nG5cw8QyWf+0rSz43ka5fyAgAAAMgCJskPe9YJP3gM95g+FRUVqqioSDyvr69PKU8kEkn5PSdjen9ga2+oV2cK585ElpEaa1mM3SPV10n1R2Tq35Hq35GaGmXaWqRjrVJb76+TLaKalx8fseD3Sz5/fORAIBgfJREJKFBYpJhty8oJxI/xB979lZMjeXNk+XzvjibweOKPgYBUUCwVFMnKyYn/EN7WO7qiszP+w3ffr0kF8WPz8hNf98buiY9a6OpKjGIoCYXU1NUt5ebJ8gxcKcHYdnzHm1g0ni+YK+vEguO9968rJtlGCgQHPa5fFl9O4rjB/oxMd3e8ePD5Bv07nG5Of+32aPDvO0NlKS0tTfo65QUAAACQBcLhsBoaGhLPGxoaBoyoCIfD/f7Bn+yYbGZ5PPF1L1iwM+1MT4909Ih05G2ZI29Lh96SOfSmdPjN+A/2fbxeqTgsTSqMFwbvmyEVFMafTyqUlV8Qn1KRmx9/zC+IlwVD/NAuSUVp+sHY8vni+YrDJz9YkuXxSgFvvEjp5YtEZA2RxfJ44mVMXv7wzh/MSy3LcI718aN4qrhjAAAAQBaYNWuWDh8+rLq6OoVCIe3YsUMrV67sd8y8efP0X//1X7rkkku0b98+5eXljanyQlJ8/QC2Sh01c/y4uvbskv3cn2Ve2yO98Wp85EGfopBUdoqsyz4qlc6QNaVUikyVSkLxH7KBMYbyAgAAAMgCXq9Xy5cv19q1a2XbthYvXqwZM2Zo27ZtkqSlS5fq/PPP186dO7Vy5Ur5/X7dfPPNLqcegUAwPmQfKTPdx6W9u2We/7PM7ufU1Nken04x4zRZl31YmjFT1vumS9PKZOVNcjsukFaUFwAAAECWKC8vV3l5eb/Xli5dmvjYsix9+ctfdjpWegVzZRh5MWzGGOlArcz2P8o897TU0Sbl5ss6f74KL7tCx953CkUFJgTKCwAAAADOYeTFsJjDb8vsflbmr09JBw/Ed7k4f76sixdJc8+T5ctRMBJRW5YsHgpkGuUFAAAAAOcEcqVjLW6nyDqmvU16/RWZ1/bI/O/npCMH45847QxZn7tJ1kULGWGBCY3yAgAAAIBjrGBufMtOyBw9IvPnP8TLioMH4i96vdIZH5B1+VWyzr1IVmiyuyGBLEF5AQAAAMA5E3zaiDFG2v1X2X96XNq7S7I80plny5p3qazZc6VTz5AVCLgdE8g6WVVe2Lat1atXKxQKafXq1W7HAYCMq6jw6Ze/dDsFAAAOCuZO2K1STfsx2f/zfmnnX6RQRNbHr5d1SYWsUMTtaEDWy6ry4rHHHlNZWZk6OyduEwtgYnnmGY/bEQAAcFYgV4p1yhgjy7LcTuMY89rLsjdXSS1Nsj61TNYVn5Dl8bodCxgzsuZfzQ0NDdq5c6eWLFnidhQAAAAAmRIMSrYtHe9yO4lj7Or/Jfveb0tenzzf/JE8H76a4gJIUdaMvHj44Yd1ww03DDnqorq6WtXV1ZKkdevWKRJxfniVz+dz5boTFffbWdxvZ1RU+PqNuCgrK5UkLVxoq7q6261YEwJf487ifjuL+40xIxCMP8aikn/8r+1gP/mYzCMPSefPl2f5rbKCeW5HAsakrCgvXnjhBRUVFWnmzJl6+eWXBz2uoqJCFRUVief1LuxpHIlEXLnuRMX9dhb32xknrnFRVlaqgwcPJZ5z+zOLr3Fncb+dxf1OTWlpqdsRJq5gbvwx2ikVFLmbJcPsHX+U+dcHpXMvkucr/68sX1b8+AWMSVnxt+fVV1/V888/r127dqmrq0udnZ3auHGjVq5c6XY0AAAAAGlkBXJlJCk2vte5My9sl3n4x9JZ58rzVYoLYLSy4m/Q9ddfr+uvv16S9PLLL+s//uM/KC4ATAgLF9puRwAAwFl900ai43fHEbP/Ndn/f5U0a448K74tK8fvdiRgzMuaBTsBYCJyao2LqqoCR64DAMBJ9U0bGafbpZqOdtk/vUcqKpHna3fK6itrAIxK1pUXf/d3f6fVq1e7HQMAxpX16ykvAABZItg38mL8TRsxxsj8yyap8ag8//f/Iyt/ktuRgHEj68oLAAAAAONYID7ywozDNS/M9mqZmmdkffx6Waef5XYcYFyhvACAcaqqqkBlZaWJrVj7PmYKCQDAVeN02og5/LbML34qzTlb1kevcTsOMO5kxYKdAID0u/32Y7r99mOSBm7JCgCAawLjb9qIMUb2zzZKfr88X/6GLI/X7UjAuMPICwAAAADOyfFLlmdc7TZi/vqU9Pp/y7pmmazisNtxgHGJ8gIAJoBvfOOYI9dhSgoA4GQsy4pPHRkna16YaIfMrx+W3n+6rAVL3I4DjFuUFwAwAfRNH8k0djUBAAxLIJjxaSOm8ahM7d6MXkOSzH/+m9TSKM/1X5Xl4ccrIFP42wUAAADAWcFgxhfsNI//Wvb9azN7jXcOyfzhd7IWLJE1c05GrwVMdJQXAIBRYVcTAEDKArkymS4vWluk9mMyGRzhYT/ykJSTI+vqL2TsGgDi2G0EADAq7GoCAEhZMFeKdmT2Gh1t8cfmBmna9LSf3ry2R3rpeVmfWiarqCTt5wfQHyMvAAAAADgrmJvxaSNq713vqakhI6e3//PfpIIiWYv/R0bOD6A/ygsAQNo4tasJAGBsswLBzG+V2tEuSTIZKC/M/n3S3l2yrqiU5Q+k/fwABqK8AACkjVO7mkhsywoAY1ogmPmtUk+cNpJm9mP/JuXly/rQR9N+bgDJUV4AAMYktmUFgDEsmJvRkRemp0fq7F1TI83lRfeB16Xdz8q6/GOycvPSem4Ag6O8AAAAAOCsQK7UFZWx7cycv3fKiJT+aSPt//7PUiAoa8lVaT0vgKFRXgDASTA9IXuwLSsAjBPBoGSM1BXLzPnbT5jGmMbywtQdVvTP1bIWfUTWpMK0nRfAyVFeAMBJMD0he9x++zEdPHgosR1r38dOrrUBAEiDQG78MVM7jvStd1EUkpob03Zas+23ktcn64rKtJ0TwPBQXgAAAABwVjAYf4xmaNHO9t7youwUqbVJprt71Kc07W0yf3lSwcuWyioOjfp8AFJDeQEASTA9IfuxLSsAjF1WYuRFZsoL0zvywip7f3x6SmvT6M+5vVrqiinvymtGfS4AqaO8AIAkmJ6Q/Zz6s6CwAoAMCPaWF5nacaRvzYuyU+OPo1z3wtg9Mn96TDr9LOXMnDO6bABGhPICAIAhsOYJAGRAoHfaSIZGXigx8uKU+PPRrnvx0k7p6BFZl7PDCOAWygsAOAmmJwAAkGa900ZMxkZetMVHd4Snxq/TVD+q09lPPioVh2Sd/8F0pAMwAj63AwBAtmOqyMRTVVXQb8RF39on3/jGMb4eAGREW1ubNmzYoKNHj2ry5Mm67bbbNGnSpAHHbdq0STt37lRRUZGqqqpcSJomwcyueaH2NilvkjSpQPLlSM0jnzZijrwtvbxL1ieul+XjxyfALYy8AADgPVjzBIDTtm7dqrPPPlsbN27U2Wefra1btyY97kMf+pC+9a1vOZwuA3J7y4uO9oyc3nS0SfmTZFmWVBySmkY+bcQ8+Zjk88m67MNpTAggVZQXAAAAgMtqamq0aNEiSdKiRYtUU1OT9Li5c+cmHZEx5uTmSzl+qWX0u4Ak1TfyQpJKwjIjHHlhoh0yO/4oa96lsgpL0hgQQKoY9wQAwBCcWvOkqqqAkR3ABNbS0qKSkvgPxyUlJWptbR31Oaurq1VdXS1JWrdunSKRSErv9/l8Kb8nFfWhiHKi7SoaxjVSzVIf65Rv+vtVHImoeWqpumtfGdHvpWPb73Qs2qniT3xW/t73Z/q+pIIsyZElubGehfICAIAhOFUorF9PeQGMd3fffbeam5sHvH7ddddl5HoVFRWqqKhIPK+vT23RykgkkvJ7UtFTUKyedw7r+DCukWqWntYW2T6/6uvrZedNkmk4qqNHj8ankaSS8bHfSGXvV0t4mqze62f6vqSCLMmRJbmxkqW0tDTp65QXAAAAgAPuvPPOQT9XVFSkpqYmlZSUqKmpSYWFhQ4mc4dVHJJ5a39mTt5+rN+0ER3vim+fmj/87a/NgdelA7WyPvuVlEsPAOnHmhcAALikqqpAZWWlid1M+j6uqhr+P64BjA/z5s3TU089JUl66qmndOGFF7qcyAHFIal55AtpDsZ0xaTu41J+vLywisPxT6S4Xap55vdSjl/W/A+lOSGAkaC8AADAJexqAqBPZWWlXnzxRa1cuVIvvviiKisrJUmNjY364Q9/mDjuvvvu0x133KFDhw7pxhtv1BNPPOFW5NErDkmxTploR3rP294Wf+wtL5QoL4ZflJhop8xfn4ov1Jk3DhZIBcYBpo0AAAAALisoKNB3vvOdAa+HQiGtWbMm8fzWW291MlZmFYXij82N0rS89J23I15eJEqHkviigKa5QcOd/GFqnpGinWyPCmQRRl4AAJAFnNrVBACyhVV8QnmRTu2930/71rcoKpEsK6VpI+aZbdL7ZkizzkxvNgAjRnkBAEAWcHKqCGtqAMgKveWFSXd50Tvyom/BTsvnkwqKhl2SmDffkPa/JuuyD7NQJ5BFKC8AAJhg1q+nvACQBfpGXrSkt7ww713zQpKKwzJNDcN7/5+3Sb4cWR9cnNZcAEaH8gIAAACA46xgnhTIzcC0kf4jLyTFt0ttPnl5YY4fl/nr07LKPygrhW1VAWQe5QUAABMA27ICyEqZ2C61vU3yeKTcdxcBtUrC0nBGXrz4nNTRJuuDl6c3E4BRY7cRAAAmgNtvP5ZYV6OsrDSxPSsAuKo4lJk1L/Ly+69XURyW2o/JdMVk+QODvtXe8US8UJl7bnozARg1Rl4AAAAAcIVVFEr7mhfx8uI9o8pKwvHHIYoS09Ik7XlB1gcXy/J405sJwKhRXgAAMMGwLSuArNE7bcQYk7ZTmvZj/RfrlGQV95YXQ0wdMX/9k2Tbsj64JG1ZAKQP5QWAMY85+0BqnNqWtaLC2dmpTn4vcPr7zni+Ht/DJ7jikHS8S+poT985O9qlvPz+r/WOvDBNR5O+xRgjs+MJ6bQzZL1vevqyAEgbygsAYx7bPgLZ6ZlnnP1nhpPfC5z+vjOer8f38Amub7vUdK570X5s4E4hkalSQZHMnx6Xse2B73nzDengAVkLGHUBZCvKCwAAAACusIp6y4uWYewEMlwdbf23SZVk5fhlfepL0uv/LfPnPwx4i9nxR8mXI+vChenLASCtKC8AjEls+whkp2uuCSf9u3nNNeGMXM/J7wVOf98Zz9fjezgSekdepGvHEWPbUnv7gDUvJMn64GLpjA/I/OZnMq3N776nuUHmuadknXexrCTvA5AdLJPO1XEcduiQ89u8RSIR1dfXO37diYr77ayxer/H8raPY/Wej1Xcb2c5/XfTyetl4+8tnV/f4/le9iktLXX8mm5L9d/OTnzPNLGY7K99WtYnPy/PlZ8edRbT0S77ls/K+vRyeZZWDvz8oTdlf+8WWRctkmf5rTL79sr+yf8ndXbIc/v3Zc2cc9JrZNN/S8iSHFmSGytZBvv+7OxKWgAAAADQywoE4otrpmvNi/beBYnfu+ZF3/VKT5G1tFLm8d/IDubKPP1fUniKPLd9T1bZ+9OTAUBGMG0EwJjHto9Adlq4MMmieBnk5PcCp7/vjOfr8T0cKgrJtKSpvOjdtcTKzx/0EOt/XCeFp8g8+Z/S35XL8+0qigtgDGDkBYAxz6ltHwGkprq6W06OTnXye4HT33fG8/X4Hg4Vh9I/8iJv8PVTrEBAnpu/JXOgVtYlFbI8/P9cYCygvAAAAADgGqs4JPPqnvScrKMt/niShTetU2bKOmVmeq4JwBHUjAAAAADcUxySWpriO4WMkmkfXnkBYOyhvAAAAADgnqKQ1NP97pSP0egbeZFHeQGMN5QXAAAAAFxjFYfiH6Rj3Yv2Y1KOX5Y/MPpzAcgqlBcAAAAA3FOUxvKio51RF8A4RXkBAAAAwD29Iy9Mc8OoT2Xaj7HeBTBOUV4AAAAAcE/fyIuWdEwbaWPkBTBOUV4AAAAAcI2VkyNNKkjTmhdtjLwAxinKCwAAAADuKgrJpGXNizZZjLwAxiXKCwAAAADuKg6lceRFwejPAyDr+NwOAAAAAGBis4pDMq/uUc+mH8iafqqsGTOlcy+U5fEO+xymu1uKdUr5+RlMCsAtlBcAAAAAXGV96EqZzk7p4AGZ3X+VMUbWR6+RdfUXh/V+Y4zM7/89/qQ4nMGkANxCeQEAAADAVdaps+W9abUkycRiMg//g8wf/0NmycdlFZUM+V7TFZP52Y9lnnta1kWXybp4kQOJATiNNS8AAAAAZA0rEJBVeYPU3S3z+K+HPNY01sv+0RqZmmdkXf0FWV++XVaO36GkAJxEeQEAAAAgq1hTS2UtWCLz1OMyjUeTHmP++0XZ379NOnJQnhXfluejn5JlWQ4nBeAUygsAjqiqYuVvAAAwfNZVn5GMZP7z3/q9boyR/fhvZK//jpRfIM+375V17kUupQTgFMoLAI5Yv57yAgAADJ8VniLrsqUy2/+g7iMHZbqPy+x5QfaP75b595/JumBBvLh43wy3owJwAAt2AgAAAMhK1pWflvlztZrvvl12U4PU2S4FcmV95v+SteTjTBMBJhDKCwAZU1VV0G/ERVlZqSTpG984pttvP+ZWLAAAsk5bW5s2bNigo0ePavLkybrttts0adKkfsfU19frgQceUHNzsyzLUkVFha688kqXEjvDKg7L+sg1sp96XFb5fFnnL5DmnsuinMAERHkBIGNuv/3dkqKsrFQHDx5yOREAANlp69atOvvss1VZWamtW7dq69atuuGGG/od4/V69fnPf14zZ85UZ2enVq9erXPOOUfTp093KbUzPB//rCLLv676+nq3owBwEWteAAAAAC6rqanRokWLJEmLFi1STU3NgGNKSko0c+ZMSVJubq7KysrU2NjoaE4AcAsjLwA44hvfYJoIAACDaWlpUUlJiaR4SdHa2jrk8XV1ddq/f79OP/30QY+prq5WdXW1JGndunWKRCIpZfL5fCm/J1PIkhxZkiNLcmM9C+UFAEewxgUAYKK7++671dzcPOD16667LqXzRKNRVVVVadmyZcrLyxv0uIqKClVUVCSepzrtIhKJZM1UDbIkR5bkyJLcWMlSWlqa9HXKCwAAAMABd95556CfKyoqUlNTk0pKStTU1KTCwsKkx3V3d6uqqkoLFy7UxRdfnKmoAJB1WPMCAAAAcNm8efP01FNPSZKeeuopXXjhhQOOMcbowQcfVFlZma666iqnIwKAqygvAAAAAJdVVlbqxRdf1MqVK/Xiiy+qsrJSktTY2Kgf/vCHkqRXX31VTz/9tPbs2aNVq1Zp1apV2rlzp5uxAcAxTBsBAAAAXFZQUKDvfOc7A14PhUJas2aNJOnMM8/Ur371K6ejAUBWYOQFAAAAAADIapQXAAAAAAAgq1FeAAAAAACArEZ5AQAAAAAAshrlBQAAAAAAyGqUFwAAAAAAIKuxVSoAAAAwAZSWljrynkwhS3JkSY4syY3lLIy8AAAAADDA6tWr3Y6QQJbkyJIcWZIb61koLwAAAAAAQFajvAAAAAAAAFmN8gIAAADAABUVFW5HSCBLcmRJjizJjfUsljHGZCBLSurr6/XAAw+oublZlmWpoqJCV1555Unfd+jQIQfS9ReJRFRfX+/4dScq7rezrrtumn75yyNux5hQ+Bp3FvfbWdxvZ3G/U5NNi9YBAE4uK3Yb8Xq9+vznP6+ZM2eqs7NTq1ev1jnnnKPp06e7HQ2YUJ55hsFYAAAAALJPVpQXJSUlKikpkSTl5uaqrKxMjY2NlBcAAACAC3bv3q0tW7bItm0tWbJElZWVjl1706ZN2rlzp4qKilRVVSVJamtr04YNG3T06FFNnjxZt912myZNmpTxLIONEHc6T1dXl7773e+qu7tbPT09mj9/vq699lrX7osk2bat1atXKxQKafXq1a5mWbFihYLBoDwej7xer9atW+dKnvb2dj344IN66623ZFmWbrrpJpWWljqe49ChQ9qwYUPieV1dna699lotWrTIlT+jRx99VE888YQsy9KMGTN08803q6urK+UsWTFt5ER1dXX67ne/q6qqKuXl5fX7XHV1taqrqyVJ69atU1dXl+P5fD6furu7Hb/uRMX9zryKCl/SERcLF9qqrubeZxpf487ifjuL++0s7ndq/H6/2xGylm3buuWWW3THHXcoHA5rzZo1uuWWWxz7H4t79+5VMBjUAw88kCgv/uVf/kWTJk1SZWWltm7dqra2Nt1www0Zz9LU1KSmpqZ+I8RXrVqlP/3pT47mMcYoFospGAyqu7tb3/nOd7Rs2TI999xzrtwXKf4D6euvv564L279GUnx8uKHP/yhCgsLE6+5kef+++/XWWedpSVLlqi7u1uxWEy//e1vXbsvUvzv81e/+lX94Ac/0O9//3vHszQ2NurOO+/Uhg0b5Pf7tX79epWXl+vtt99OOUtWjLzoE41GVVVVpWXLlg0oLqT4oh4nLuzhxrxO5pM6i/udeb/85bsfl5WV6uDBd9eS4dZnHl/jzuJ+O4v77Szud2pY82JwtbW1mjZtmqZOnSpJWrBggWpqahwrL+bOnau6urp+r9XU1Oiuu+6SJC1atEh33XWXIz8ADjZC3Ok8lmUpGAxKknp6etTT0yPLsly7Lw0NDdq5c6euvvpqPfroo5Lc+zMajNN5Ojo69Morr2jFihWS4oWyz+dz/b689NJLmjZtmiZPnuxaFtu21dXVJa/Xq66uLpWUlOi3v/1tylmyprzo7u5WVVWVFi5cqIsvvtjtOAAAAMCE1NjYqHA4nHgeDoe1b98+FxNJLS0tiRKhpKREra2tjmeoq6vT/v37dfrpp7uSx7ZtffOb39SRI0f04Q9/WLNnz3btvjz88MO64YYb1NnZmXjN7T+jtWvXSpKuuOIKVVRUOJ6nrq5OhYWF2rRpkw4cOKCZM2dq2bJlrt+X7du365JLLpHkzp9RKBTSxz72Md10003y+/0699xzde65544oS1aUF8YYPfjggyorK9NVV13ldhxgwlq40HY7AgAAcFmyWeWWZbmQJHucbIS4Ezwej+655x61t7fr3nvv1ZtvvulKjhdeeEFFRUWaOXOmXn75ZVcyvNfdd9+tUCiklpYWff/733dlZFVPT4/279+v5cuXa/bs2dqyZYu2bt3qeI4TdXd364UXXtD111/vWoa2tjbV1NTogQceUF5entavX6+nn356ROca9tYCP/vZz/S3v/1tRBc5mVdffVVPP/209uzZo1WrVmnVqlXauXNnRq4FYHCscQEAAMLhsBoaGhLPGxoaEv+H1C1FRUVqamqSFF+H4sS1DTIt2QhxN/Pk5+dr7ty52r17tys5Xn31VT3//PNasWKF7rvvPu3Zs0cbN2509Z6EQiFJ8T+XCy+8ULW1tY7nCYfDCofDmj17tiRp/vz52r9/v6v3ZdeuXTrttNNUXFz8f9i79/CoynP//59nZUJOkMPMJIRwPoqIVk5S0FatqfbormjZUi27SGvVXatiPbSK/ixVcavo9lDbXbGlh23Vtnhqu/s12qo12oII0iIqWkUQDUmYTM4kmef3xwqBmITMJJlZk8n7dV1eM5lZs9adBXqZT+7nfiR58/d269atKioqUm5urnw+n+bPn6833nijT7VEHV60tbXpxhtv1OWXX65HH32003/Q+mv69Ol6+OGHddttt+nWW2/VrbfeqtmzZw/Y+QEAAABEZ/LkydqzZ48qKirU2tqq8vJyzZ0719Oa5s6dq2effVaS9Oyzz2revHkJuW5PHeKJriccDqu+vl6Su/PI1q1bNXr0aE/uy1e+8hX96Ec/0r333qtLL71UM2fO1Le//W3P/oyampo6lq80NTXp1Vdf1bhx4xJeT35+vgKBgN5/350ft3XrVo0ZM8az+yJ1XjIiefPvUTAY1Jtvvqnm5mZZa/v1dzem3UYikYheeeUVPf/889q0aZOmTp2qT37yk5o/f37HAJlEOvAXI5EYhpVY3O/E4n4nHvc8sbjficX9Tizud2wY2Hl4mzZt0rp16xSJRHTyySdr0aJFCbv2nXfeqW3btqm2tlZ5eXlavHix5s2bpzvuuEOVlZUKBoNasWJFQrZ43L59u6677jqNGzeuY+nMkiVLNHXq1ITW8+677+ree+9VJBKRtVYLFizQWWedpdraWk/uywH//Oc/kGR9jgAAIABJREFU9cQTT+jqq6/2rJYPP/xQt912myT3F+4nnHCCFi1a5Ek977zzjn70ox+ptbVVRUVFuuiii2St9eS+NDc368ILL9Q999zTsdTJqz+jhx9+WOXl5UpLS9OECRN0wQUXqKmpKeZa+rxV6nvvvae77rpLO3fu1LBhw3T88cdr8eLFHS07iUB4kfq434nF/U487nlicb8Ti/udWNzv2BBeAMDgEtPAzoaGBr300kt6/vnn9e6772r+/Plavny5gsGgnnzySd10000dqRcAAAAAAMBAiDq8uP3227VlyxYdeeSR+vSnP6158+YpPT294/2lS5fqa1/7WjxqBAAAAAAAQ1jU4cXUqVO1fPnyjkmlH+U4jn7yk58MWGEAAAAAAABSDOHF6aef3usxGRkZ/SoGAAAAAADgo6LeKhUAAAAAAMALhBcAAAAAACCpEV4AAAAAAICkRngBAAAAAACSGuEFAAAAAABIaoQXAAAAAAAgqRFeAAAAAACApEZ4AQAAAAAAkhrhBQAAAAAASGo+rwsAAAAAEH/vv/9+TMcHg0FVVlbGqZrYUEv3qKV71NK9wVJLSUlJt6/TeQEAAAAAAJIa4QUAAAAAAEhqhBcAAAAAACCpEV4AAAAAAICkRngBAAAAAACSGuEFAAAAAABIaoQXQIq5/fYRXpcAAAAAAAOK8AJIMWvWEF4AAAAASC2EFwAAAAAAIKkRXgAp4PbbR2j06BKNHl0iSR3PWUICAAAAIBX4vC4AQP9dfnmtLr+8VpIbXOze/b7HFQEAAADAwKHzAgAAAEAnkWee1L7vX+Z1GQDQgfACSDErVtR6XQIAAEgB+1/5m+yud7wuAwAkEV4AKefA8hEAAIC+MvM+ITlpsn971utSAEAS4QUAAACAjzAj8jRs1nzZvz0rG4l4XQ4AEF4AAAAA6CrrxNOkfZXSm//0upSo2S1/l21s8LoMAHFAeAEAAACgi4zjPiFlZMm+9JdOr9uNf1Xkwf+RtdabwnpgQ9WK3PMD2fKnvS4FQBwQXgAAAADowmRkysxeIPvyC7L7myVJdte/FFl7h+wzT0pvJFlHRjjkPu6r9LYOAHFBeAEAAACgW+bjJ0mNDdKrG2SbmxT58a1SznBp+AhFyh73urzO6mrcx5p93tYBIC58XhcAAAAAIElNP1rK8yvy0l9k/rFJ+nC3nMu+L7t9q+wfH5Hd+4FMYbHXVUqSbG3YfSS8AFISnRcAAAAAumWcNJn5n3Q7L14ok/ncl2WO/JjMyZ+VHEf26Se8LvGg2vbOi1C1t3UAiAvCCwAAAAA9MvNPkqyVphwp88Ul7mv5AZm5J8i+UJY8u3u0d16wbARITYQXAAAAAHpkxk2S8+3r5Vx0jUxa2sHXS0+XmhplX3jKw+oOUdceXjTUybbs97YWAAOO8AIAAADAYZmj58iMyO382oSp0pQZsk8/KRtp86awQ9gDAzslui+AFER4AQAAAKBPnNLTpcoPpVde8roUd+aFMe5zwgsg5RBeAAAAAOibY+dLo8Yq8tt1si0t3tZSG5YO7HxSw9BOINUQXgAAAADoE5OWJufsr0t7P5Ate8zbYupqpNHjJbFdKpCKCC8AAAAA9JmZMUs6dr7s7x+WDVV5UoONtEn1dTKjxknGkUKEF0CqIbwAAAAA0C/Ol8+T2lplf/tzbwqor3O3c83Nd/+p8SZEARA/hBcAAAAA+sUUjZL59JdkX/qz7FvbE19AbftOIyNypbwClo0AKYjwAgAAAEC/mc99Wcr3K/Lrn8ham9iL14bdGkbkSfl+lo0AKYjwAgAAAEC/mcwsmTO+Kr3zprTpxcRevO5g54XJK2C3ESAFEV4AAAAAGBDm4ye5W6c+9it3iGaC2APLRobnSnl+qS4s25a46wOIP8ILAAAAAAPCOGly/u0r0p73ZF96NnEXrnOXjbjhRYE7vDMcStz1AcQd4QUAAACAgTNrgTRusuwTD8q2tiTmmrVhKStHxpcuk1/gvsbSESClEF4AAAAAGDDGceR86Vyp8kPZvz6VmIvW1rg7jUjushFJYscRIKUQXgAAAAAYWDNnS1NmyD75sOz+5rhfztaFpRF57hd5bueFpfMCSCmEF0ASuv32EV6XAAAA0GfGGDlnfFWqqZb9/SPxv2BtjTvvQpJy891HtksFUgrhBZCE1qwhvAAAAIObmXaUzMJTZP/4iOz2V+N7sbqwTHt4YXzpbpDBshEgpRBeAAAAAIgLs+R8aWSJIvevObid6QCz1roDOw8sG5GkvAKWjQAphvACSBK33z5Co0eXaPToEknqeM4SEgAAMFiZzCw5518p1dcq8sCdspHIwF+ksUFqaz04sFNyh3bSeQGkFMILIElcfnmtdu9+X7t3vy9JHc8vv7zW48oAAAD6zoydKLN4ufSPl2XLHhv4C9S1d3QMP9h5YfIKCC+AFEN4AQAAACCuzEmflWYvkP3NOkVe/PPAnrw27F7j0M6L/AIpvC8+nR4APEF4ASShFSvotgAAAD3bvHmzLrnkEl188cV69NFHvS6nV8YYOeetkKYfLfvTOxX561Od3rfhkGx1Zd9OXtu180J5fqmtTarn/6mAVOHzugAAXbFUBAAA9CQSiWjt2rW69tprFQgE9N3vfldz587VmDFjvC7tsExGhpxvXavID2+SXXe3Iq0tMv5CRZ5/Stq6wT1m8XKZkz8vY0zU57V1bufFoTMvTF6BrCTVVHce5Alg0CK8AAAAAAaRHTt2qLi4WCNHjpQkLVy4UBs2bEj68EKSzLAMOf95jSL3rZb91Y/cgGFEnswpp8t+uFv2wf+R3n1LOvdCmfRh0Z209kB4cUhIke93H0P7pDETB/JbAOARwgsAAABgEKmurlYgEOj4OhAI6M033+xyXFlZmcrKyiRJq1evVjAYjOk6Pp8v5s9Ey668TfWP/q98YycqY+7xMj6fbCSi+ocfUP1DDyit4n3lf/cWpfmDvdZS27ZfDcMyVDj6YHjT2jpZVZKGR1qUNcDfQzzvS6yopXvU0r3BXgvhBQAAADCIWGu7vNbdMovS0lKVlpZ2fF1ZGdtMiWAwGPNnYnLyFyRJdaHQwddKvyQnUKzWtWtUef235Vy5WiYr+7C1RCo+kIaP6PS+bZ/TWbtrp+oH+HuI+32JAbV0j1q6N1hqKSkp6fZ1BnYCAAAAg0ggEFBVVVXH11VVVSooKPCwooFlZn1czoXflfa8p8h9N8u2thz2eFsb7jysU+7yFGXlsF0qkEIILwAAAIBBZPLkydqzZ48qKirU2tqq8vJyzZ071+uyBpQ5apbM0m9Jr22RXXd3t90mHerCnYZ1dsgrkK2pjl+RABKK8AIAAAAYRNLS0nTeeefpxhtv1GWXXaYFCxZo7NixXpc14JyFp8j82zmyL/1F9b/6cc8H1tbIdLejSF4BnReQrdl3+PDrwHHv75R94x+9Hxdpk938kmxTY+/H1obVVvlh78dVVSjyk9uiu7610X0/LfvdOvdV9XrsYMHMCwAAAGCQmT17tmbPnu11GXFnPr9Yqt6r+t/+XM6o8TIfm9f1oLquy0YkyeT5Zf/1egKqRKLYxgbpjX/Kbn9V8gflfPrfej7WWtknHpR94tcyi/5D5rNn9nxsXViRNSulurCcy1bJHDGz52N/u072/z0qTZwm55LrZXJGdH/cu28pcvf3VdXcLHPh1TIzju3+uOpKRW6/Vtr7gezGv8p86asyp50h43TuM7CRiOyG52V/93OpZb/MoqUyC0/pcpwk2be2K7LubmnPe5JxpJmz5ZxQqsjxn3JDF2MkGcmo83MZ9+tDnseybXG8EV4AAAAASErGGGnJ+Urb9S+1/uxOOdfdJVNwcKcVu79Zam6ShnfzA2S+23lhrU2qH8AQO9vYoMh9N0uvb5UiEfcHbGtlp82UGT+56/EtLbI//W/ZF5+R8vyyj/1Sdsax3R4ryd2ity4sFQQV+dFqOdfcLhMc2eW4yPP/zw0uZs6Wtr+qyK3fk3PpDTIHtuY9cL6tLyvy41uknBFKKypQ613fl1m+Qs68EzofF6p2g4u6sJwVq2Sf+5Ps79bJ7tgm5+xvSE6a1NYqVX6oyKO/lP71hjRukpQ+zF1O9dyf5Pz716Xi0VJbm9TaIvvU47JPPy4VBGS+frn0/k7Z8qcVuW+19t63uu9/CDEEHh3vTZiqtMt/0PdrfgThBQAAAICkZdKHKe/y76vq8mWKrF0jZ8X3ZZw0983asPvY3bKR3Hxpf7PU3ChlZieuYPRZT0GTfeEp6bUtbkfCzDlSyThFrv+WIo88IOfyH3T6jG2oV+ju78u+ulHm9K/InPw5RW74tvt359o17jDXQ8/9crns359zjz3uk4rcdLki994o5+r/ksnIPHjc61tlf3WfNGOWnG+tlN74hyL33qjIf10t5xvfkYZlSm2tsm9uk334fmnMBDkXr1TBqBLtveEy2Z/cqkhNtcysBW4g0dSgyE9ul2r2ybnsBpnJ06Xpx0hTZ8g+/IAi3zu/803I98ssu0Tm4ydLxsj+7S+yv1mnyOoru9wvc9LnZM5cKtP+996e/hVp2yvKrqlWfX2dZK1kJcm6z92b3/68/b0Dzw99r9Nn2t+LHPIZfeQ4f2EUf+rRI7wAAAAAkNR8o8fLLPmm7M/+W/YPv5H5wr+7b9TVSFL3My9y23dgCYcILzxiqyulf73hDk7NyJSGZcqk+2RD1dLeD2X37pFC1W7XQ32dtL9JzoXflTnm4PIgG4nI/vmP0qQj5Jy1rON188Wz3Y6JrRul9uNtU4Mid1yn/e+9LbPsEjkLT5EkOcsuUeSO62V/u05mycFQwIZDivzyh9L4KTKfPUvG55PzjSsUuev7ijxwp5zPfVlqbZEa6hVZu0YqHCXnm1fIpKVJR35MzmXfV+SuGxS56Tudv/GZs+V880qZzGw5OSPkXHqDO9PioftlH7r/4HHDMuRc8v+5wYXcTiPzqS/ITpspu+M1KS1NSvPJZGRIM+d0ClPMx0+WPXa+7N+ek1r2tx+bJjN2kszEaZ3KMWlp0tFzlRMMqjFJtkrtC8ILAAAAAEnPLPyUtG2z7OMPyk6cJnPUrEM6L7ruNmJy891fANeEpKKShNY6VNlQtezWjdI/X5F9a7sU6jossmPUZPowKVAk+YPuEo3hI9zlFr/5mZyZsw9212x/Vap4X+aLl3U6j/nkZ2SfftI9/qjZUiSiyA9vlna+pbyrblbdpCMPHjtjlswpX5R9+glFxk2SKSyWWlsVefoJqalBzrJLZXzuj8Zm5myZM/9D9jc/VWRT+cELDh8h5+KVMtnDD5538nQ5K++U3bHNDRnS0qTMLOmIY9znB44bliHngqtlN5VLTY2SL90NGsZNlike3eUemTETZMZM6PV+m8xsmRM/0+txqYLwAgAAAEDSM8ZI514ou/sdRe5ZJeebVx7c8aGbgZ3Ky3cfw+w4Ek+2rU32+T+p6sU/K/J2+4DU/IDMtJnSpCNkJh0hBYval/A0uV0CeQVSbkHXoZQvlyvyo9WyL/1Fpr1rIvLnP0jDc2XmHN/pWOPzyTnzPxS572bZ5/8k+9qr7tKSZZcq87hPqO4jHQZm0VLZ17bI/uwuHbpXhznrazKjx3U+9tQvyUyYKjXWS2k+yeeTRo+TOdDNc+ixwZHdzsfoclxamsy8T/R6HHpGeAEAAABgUDBZ2XKuuEmR/75BkftWS1OPct/oaeaF3KUBjOuMD/uPlxV5+AF3V4sp02W+dK67I8zoCX0bkjp7gTR+ittdM++TUm1I2vJ3mc+cIZM+rOvxsz4uTZkh+78/lqyV+fflchZ+qttTm2EZcq66RXp7e3sgkS5l58iMHt/1WGOkw+w4Am8QXgAAAAAYNEzOCDkrVinyw5uk17ZIjiNldTPTYniuu01kDZ0XA8W2tkg735Z9e7vsqxvd+180Ss5/fk/+Uz6vqqquy0RiYYyRc8ZXFbnzetnn/uTOK5GV+WT3SyOMMXIWn6fIf10tc9oiOaU9b50qSSY7R5o5p181wjuEFwAAAAAGFZOZJefilbJr75Ct/LDL8gNJ7syEEbntPwCjP2xDvewjD8j+7Vl32YckBYpkvnyezKc+L+NLH7jtaGccKx1xtOzvH3K/PmbeYZdlmInT5NzxK5nMrIG5PpIW4QUAAACAQcekD5O54CrZSKTng3ILZAkv+sVu26zIurukULXMCafKzPiYNGm6TEEgLtfr6L5o3wLUOemzvX+G4GJIILwAAAAAMGh113XRITefzos+sk0Nsr/9uexf/iAVj5Fz9X912YIzXszk6dKchdL770kzZiXkmkh+hBcAAAAAUpLJy5f9YJfXZQw6dvPfFPnfH0uhKpnSf5M541yZYRkJrcH5+nektrbDh1MYUggvAAAAAKSm9s4La+3AzWRIYTZUrciD/yNtKpdGj5fzzSvdLggPGF/7FqVAO/42AAAAAEhNuQVSa4vUWC9lD/e6mqRmt7+qyP/cKjU2yJzxVZlTz3ADBCBJ8LcRAAAAQGrKK3AfwyHCix5Ya2X/9DvZ3/1CGlki5zs3ypSM87osoAvCCwAAAAApyeTmy0pSTUgqHuN1OUnHNtQp8rO7pFdekplzvMzXLpbJzPa6LKBbSTP9ZPPmzbrkkkt08cUX69FHH/W6nB6tWpWW0OvdfvuIlLxWql/vzDPjs3UUAAAAYpDrdl6wXWpX9l9vKrLqMunVDTKLl8t880qCCyS1pAgvIpGI1q5dq+9973u644479MILL2jXruScCvyDHyQ2vFizJnE/cCfyWql+vZdeSuw0ZgAAAHQjN999DO/zto4kYq1VpOwxRW65SopE5Fxxs5xP/xsDTZH0kmLZyI4dO1RcXKyRI0dKkhYuXKgNGzZozBhauwAAAAD0Uc5wKS1NqiG8kCRbF3aXiWz5u3TsfDlf+7ZMTmJ/oQj0VVKEF9XV1QoEDrbZBwIBvfnmm12OKysrU1lZmSRp9erVCgaDCalv1aq0Th0Xo0eXSJKuvbZNK1e2DerrJfv35vP5+vXnnMjvr7TUp+efP9jMdOBan/hERGVlrQN6rXjp7/1G7LjnicX9Tizud2Jxv4GujONII/LcgZ1DnH19qyL3r5HqamT+/esyp3yRbgsMKsZaa70u4sUXX9SWLVt0wQUXSJKee+457dixQ+edd95hP/f+++8norxORo8u0e7dibtuIq+XjN9bMBhUZWVlwq43UBJ9LwfKQN5vRId7nljc78TificW9zs2JSUlXpeQcLH+v3My/Z3qTy1tqy6T8gqU9u3rPK9loEVTi7VW9okHZZ98SCocJeebV8iMm+xJLYlCLd0bLLX09N/npOi8CAQCqqqq6vi6qqpKBQUFHlYEAAAAICXk5g/pzgv716dkn/i1zMdPljnnApnMLK9LAvokKQZ2Tp48WXv27FFFRYVaW1tVXl6uuXPnel1Wt669duCXUhzOihW1KXmtVL/exz/enLBrAQAAoGcmL3/Izryw1ZWyjzwgHXG0zLJLCC4wqCVF50VaWprOO+883XjjjYpEIjr55JM1duxYr8vq1sqVbUpkp83llyfuB+5EXivVr/fb31b1fhAAAADiLzdfqq2RjUTcGRhDhLVWkV/+UGprk7P0W0Pqe0dqSorwQpJmz56t2bNne10GAAAAgFSSWyC1tUoNddLwXK+rSRj74jPS1o0yZ39DpmiU1+UA/ZY04QWSi93fLNXVSvW1Ul1Ytq5Wqg1JtTVSXViqr5Otr5OaG6WWFqm1RTJGMo6UmSll5cjk5kvBkTKjx0uTjpDJY44JAAAAEiw3330Mh4ZMeGFDVbIP3S9NmSFz8ue9LgcYEIQXQ4Bta3MDh3BIqg3JhmukuhqpNuwGEwcCiYZ6N6yor5Na9nd/MmPc/bKz2//JynYf09Pd99vapP3NUk217M63pJp96tjOZvwUmY+fKDP/JJkReYn41gEAADDEmbwC9/9Ha/ZJJeO8LifubFubIg/cKbW0yPnat1kugpRBeDFI2dYWKVzjdkOEQ7Jh9/HAPzYccv8DXRtyOyi64zhSzgh37+vhI6SiUTI509wwIme4NHyETE6u+97w3I7jjJMWfZ3NTdKud2Rf3yq76UXZh9bK/madzBfPlvP5xQN0NwAAAIAetHde2HBIxuNSEsH+5mfSa1tk/uNimZFDb0tgpC7CiyRhrXU7H2pr2v9p75DoeN6+ZCMcckOLhrruT5SRJeXmuf+RLh4tc8RMN3TIzZcZkd/+PM99zMqJexJrMjKlydNlJk+XPvdl2d07ZV94SmbMhLheFwAAAJDkzryQhsR2qZHyp2XLHpM55YtyTvi01+UAA4rwYoBZa905EA0NUmN9+z8N7nyI9vkRB2ZJ2Lpw+9KN9iUcba3dn/RAd0RuvszoCdKR7eFEbr47V2JE/sGvMzIT+v3GyoweJ7N4uddlAAAAYKjIzpF8vpQPL+zbr8v+4l5p+jEyZy3zuhxgwBFexCCydo2q91Wqbf9+KRJxw4bWVnc+RHOTO+thf/PhT2KMuyxjePtyDH9QZvzk9m6IfGlEbnuHRO7Br33pifkGAQAAgBRjjHF/0Vezz+tS4sbW7FPkhzdL+QE537xSxsePeUg9/K2ORfowmcwsKS1dSkuT0nwyaWnSsGHSsAwpI1Maltmx24aycmSyc9qHWua4gUV2TkwzIwAAAAD0U26BbDg1wwvb1qbIT26TGuvkXHqrzBDZUQVDD+FFDJyl31JBMKjKykqvSwEAAAAQrdx8aV9q/j+8ffx/pde3yiy7RGbMRK/LAeKGfXMAAAAApDSTV5CSMy+aXy6X/cMjMid8Ws7CU7wuB4grwgsAAAAAqW1EvlRbIxtp87qSAWOr9qrmzu9LYybKLDnf63KAuCO8AAAAAJDa8vLdgft1tV5XMmAiv7pPamuVc8FVMsMyvC4HiDvCCwAAAAApzeTmu09SZOmIfWu7tHWjcs5cKjOyxOtygIQgvAAAAACQ2nIL3McU2XEk8tivpBF5yvrcWV6XAiQM4QUAAACA1NbeeWFToPPCvvEP6bUtMp85U05WttflAAlDeAEAAAAgteW3d17sq/a2jn6y1rpdF3l+mZM+63U5QEIRXgAAAABIaSYzW8rKkfZVel1K/7y2RXrjnzKfO4shnRhyCC8AAAAApL6CgOy+Kq+r6LOOrgt/UOYTp3ldDpBwhBcAAAAAUp8/OLg7L95/T3r7dZnTFsmkp3tdDZBwhBcAAAAAUp4pGNzhhX3rNUmSOWq2x5UA3iC8AAAAAJD6CoJSOCTb2uJ1JX3z1nZpeK5UNMrrSgBPEF4AAAAASH0FAfdxkM69sG9vlyZPlzHG61IATxBeAAAAAEh5piDoPhmE4YWtC0sf7JaZdITXpQCeIbwAAAAAkPr8bnhhB+Pci7dflySZyUd6XAjgHcILAAAAAKmvY9nI4Asv7FuvS44jTZjidSmAZ3xeFwAAAACgsxdffFGPPPKIdu/erZtuukmTJ0/ueG/9+vV65pln5DiOli1bpmOPPdbDSgcPk5ktZeVI1YMxvHhNGjNRJiPT61IAz9B5AQAAACSZsWPH6jvf+Y6OPLLzMoFdu3apvLxca9as0TXXXKO1a9cqEol4VOUgVBCQHWQzL2xbm/TOmzKTp3tdCuApwgsAAAAgyYwZM0YlJSVdXt+wYYMWLlyo9PR0FRUVqbi4WDt27PCgwkHKHxx8y0Z2vys1N0mEFxjiWDYCAAAADBLV1dWaOnVqx9d+v1/V1dXdHltWVqaysjJJ0urVqxUMBmO6ls/ni/kz8TJQtYSLR6t51zv9Olei70vDhmdVKykwd4HSPnLdVPwzGgjU0r3BXgvhBQAAAOCBVatWKRQKdXn97LPP1rx587r9jLU26vOXlpaqtLS04+vKytg6DoLBYMyfiZeBqiWSNVw2VK29H+yR8aV7Wku0Iq++LOUVqNpJl/nIdVPxz2ggUEv3Bkst3XWdSYQXAAAAgCdWrlwZ82cCgYCqqg7ObKiurpbf7x/IslJb+3ap2lclFRZ7W0uU7FvbpUlHyBjjdSmAp5h5AQAAAAwSc+fOVXl5uVpaWlRRUaE9e/ZoyhS2z4yWKTgkvBgEbDgk7f2AYZ2A6LwAAAAAks7f//53PfDAAwqHw1q9erUmTJiga665RmPHjtWCBQu0YsUKOY6j5cuXy3H4fWTU2jsv7L5KDYo+hre3S5LMJMILgPACAAAASDLHHXecjjvuuG7fW7RokRYtWpTgilJEQcB9HCQ7jti3XpfSfNL4yV6XAniOmBYAAADAkGAys6WsHKl6kIQXe/dIhSNlhmV4XQrgOcILAAAAAENHQUB2kMy8UG2NlJvvdRVAUiC8AAAAADB0+IODZtmIwiGZEYQXgER4AQAAAGAIMQWDK7yg8wJwEV4AAAAAGDryA1I4JNvS4nUlh2VbWqSGesILoB3hBQAAAICho327VIWSfO5Fbch9JLwAJBFeAAAAABhCTEF7eJHsQztrayRJJjfP40KA5EB4AQAAAGDoaO+8sMk+9yLc3nnBwE5AEuEFAAAAgKGkIOA+Jnl4YcMsGwEORXgBAAAAYMgwmdlSVo5UndzhRUfnRW6Bt3UASYLwAgAAAMDQUhCQTfaZF+GQlJElk5HhdSVAUiC8AAAAADC0FASSftmIwjUSwzqBDoQXAAAAAIYUEyiS9u6RbWvzupQe2doQ8y6AQxBeAAAAABhSzIxjpYZ66c1/el1Kz8IhdhoBDkF4AQAAAGBomTlHGjZMdlO515X0LBySofMC6ODzugAAAABgMPjHP/4R1XGO42jGjBlxrgb9YTIypZlzZTe9KHv2+TJOcv1O17a1SXVhlo0AhyC8AAAAAKKwatUqFRYWylp72OPC4bB+8YtfJKgq9JWZs9DtvHhruzQ1ycKm+rBkLQM7gUMQXgAAAABRyMjI0D333NPrccuWLUtQN8WuAAAgAElEQVRANegvc8xcWV+67MsvyCRbeBEOSRLLRoBDJFd/FAAAAJCkrrjiiqiOu/zyy+NcCQaCycyWjprlLh2JRLwup7P28IKBncBBhBcAAABAFI4++uiojps5c2acK8FAMXOOl/ZVSu+86XUpndgD4QWdF0AHwgsAAAAgRn/961+1a9cuSdL777+v66+/XjfccIN2797tcWWIhfnYPCnNJ/tyku06QngBdEF4AQAAAMTooYce0vDhwyVJP//5zzV58mQdeeSRuv/++z2uDLEw2cOlIz8m+/ILvQ5iTahwjeTzSVnZXlcCJA3CCwAAACBG4XBY+fn52r9/v15//XUtWbJEZ511lt555x2vS0OMzJyFUlWFtPMtr0s5KByScvNljPG6EiBpEF4AAAAAMcrNzdUHH3ygzZs3a/LkyUpPT1dLS4vXZaEPzMfmS5Lsts0eV3KQrQ0xrBP4CLZKBQAAAGJ05pln6qqrrpLjOLrsssskSVu3btX48eM9rgyxMiNypUCR9N6/vC7loHBIyvN7XQWQVAgvAAAAgCg1NzcrIyNDJ510khYsWCBJysjIkCRNnTpVl156qZfloa/GTpJ9722vqzgoHJIZO8nrKoCkQngBAAAAROmiiy7SxIkTNWvWLM2ZM0fFxcUd7+Xl5XlYGfrDjJ0ou+Vvss1NMhmZntZirZVqa6Rc/j4BhyK8AAAAAKL04x//WK+99ppeeeUV3XLLLWpra9OsWbM0e/ZsHXXUUfL5+N/rwciMm+iGBrvekSZP97aYhjqprY1tUoGP4L+uAAAAQJR8Pp+OPvpoHX300Vq6dKkqKiq0adMm/eEPf9Bdd92ladOmadasWTruuOOUn88Pn4PG2MmSJPve2zJehxfhkPvIwE6gE8ILAAAAoI+Kior0mc98Rp/5zGe0f/9+bd26Va+88orS0tJ0yimneF0eouUPStnDk2NoZ3t4Yei8ADohvAAAAABiFIlEurzm8/k0Z84czZkzx4OK0B/GGGnsRNmd3g/ttAc6L3ILvC0ESDKEFwAAAECMlixZ0u3raWlpKigo0Pz587V48WJlZno7/BHRM+Mmyf7lj7JtbTJpad4V0hFeMLATOBThBQAAABCjZcuWacOGDfrSl76kQCCgyspKPf7445o9e7ZKSkr0yCOP6Gc/+5kuuOACr0tFtMZOklr2Sx/ulkrGeVdHuEZyHClnhHc1AEmI8AIAAACI0e9//3vdcsstys7OliSVlJRo8uTJuvrqq3X33Xdr3LhxuuqqqzyuErEwYyfKSrI735bxMryoDUkj8mQcx7sagCTEvxEAAABAjBoaGtTc3NzptebmZjU0NEiS8vPztX//fi9KQ18Vj5F86Z4P7bThEDuNAN2g8wIAAACI0Yknnqgf/OAH+uxnP6tgMKiqqir94Q9/0IknnihJ2rJli0pKSjyuErEwPp80erzsex4P7QyHJHYaAbogvAAAAABidO6556q4uFjl5eXat2+f8vPzddppp6m0tFSSdNRRR+mGG27wuErEyoydKLv5JVlr3R1IvBAOyYwk+AI+ivACAAAAiJHjODr11FN16qmndvv+sGHDElwRBsS4SdJfn5L2VUn+YMIvb611Z17QeQF0QXgBAAAAxMhaq6efflrl5eUKh8O67bbbtG3bNoVCIS1cuNDr8tBHB4Z26r1/eRJeqLlR2r+f8ALoBgM7AQAAgBg99NBD+vOf/6xTTjlFlZWVkqRAIKDHHnvM48rQL2MmSMbIvveWN9evDbuPw3O9uT6QxAgvAAAAgBg9++yzuuqqq3T88cd3zEYoKipSRUWFx5WhP0xmtlQ4StarHUeam9rryPLm+kASI7wAAAAAYhSJRJSZmdnptaampi6vYfAxYyd6t11qU6P7mEF4AXwU4QUAAAAQo1mzZunnP/+5WlpaJLkzMB566CHNmTPH48rQb8WjpcoK2dbWxF+7vfNCdF4AXRBeAAAAADFaunSpqqur9bWvfU0NDQ1aunSp9u7dq3POOcfr0tBfhcWSjUjVexN/7QOdF3TwAF2w2wgAAAAQo+zsbF155ZUKhUKqrKxUMBhUfj47RKQCEyx2dxyp/EAqGpXQa1uWjQA9IrwAAAAAohCJRLq8lpubq9zc3E7vOw7NzYNacKQkye79UCbR124+0HlBeAF8FOEFAAAAEIUlS5ZEddxDDz0U50oQVwV+Kc0nVX6Y+Gs3EV4APSG8AAAAAKJwzz33dDzftGmTXnrpJZ1xxhkKBoOqrKzUY489pvnz53tYIQaCcdKkQJG094PEX7y5SXIcyZee+GsDSY7wAgAAAIhCYWFhx/Mnn3xSq1evVk5OjiSppKREkyZN0ne/+12deuqp/b7WL37xC7388svy+XwaOXKkLrrooo5rrV+/Xs8884wcx9GyZct07LHH9vt6+IjgSFmvOi8ys2RMwhesAEmPBXkAAABAjBoaGtTc3Nzptf3796uhoWFAzn/MMcfo9ttv12233aZRo0Zp/fr1kqRdu3apvLxca9as0TXXXKO1a9d2O4sD/WMKR3qzbKS5kSUjQA/ovAAAAABidOKJJ2rVqlX6/Oc/r0AgoKqqKv3xj3/UiSeeOCDn/9jHPtbxfNq0aXrppZckSRs2bNDChQuVnp6uoqIiFRcXa8eOHZo2bdqAXBftCoul+lrZhnqZ7JyEXdY2NbLTCNADz8OLw7XEAQAAAMno3HPPVXFxscrLy7Vv3z7l5+frtNNOU2lp6YBf65lnntHChQslSdXV1Zo6dWrHe36/X9XV1d1+rqysTGVlZZKk1atXKxgMxnRdn88X82fiJdG1NE2cohpJ+W3NSg+OT1gt+yJtigwfoUCU5x/Kf0aHQy3dG+y1eB5eHHPMMfrKV76itLQ0/fKXv9T69et17rnnel0WAAAA0CPHcXTqqaf2a77FqlWrFAqFurx+9tlna968eZKk3/3ud0pLS9MnPvEJSZK1Nurzl5aWdgpTKisrY6rvwCDSZJDoWmyG+8vU0JvbZUb4E1ZLW21Y8qVHff6h/Gd0ONTSvcFSS0lJSbevex5e9NQSBwAAACSTLVu2dPp/1568+uqrOuaYY3o9buXKlYd9/y9/+YtefvllXXfddR0DHA8sUTmgurpafr+/p1Ogr4IjJUl274dK6OjMpkYpmJvIKwKDhufhxaEObYnrTn9b3wZCMrXaDAXc78Tifice9zyxuN+Jxf1OLO53/K1Zs0br1q3r9bg77rhDP/3pT/t1rc2bN+uxxx7TDTfcoIyMjI7X586dq7vuuktf+MIXtG/fPu3Zs0dTpkzp17XQlckZLmXnJH5oZ3OTTEZmYq8JDBIJCS/60hLXnf62vg2EZGq1GQq434nF/U487nlicb8Ti/udWNzv2PTUlnw4TU1NuvDCC3s9rrW1tS8ldbJ27Vq1trZq1apVkqSpU6fq/PPP19ixY7VgwQKtWLFCjuNo+fLlchw2EIyLYLFs5QeJvWYTu40APUlIeNGXljgAAAAgmVx//fVRHTcQ/z9799139/jeokWLtGjRon5fA70IjpTefzex12S3EaBHni8b6aklDgAAAEgmM2bM8LoEJJApHCn76t9lIxGZBHS32LY2qWU/nRdADzwPL3pqiQMAAAAAzwSLpdZWKVQt+RMwT6a50X0kvAC65Xl4cbiWOAAAAADwggmOlJXcoZ2JCC+amtxHBnYC3WK6DwAAAAB8VGGxJCVuaCedF8BhEV4AAAAAwEcFCiVjpL0J2i61vfPCMLAT6Jbny0YAAACAwaKtrU0bN27Upk2b9O6776q+vl45OTkaP368Zs2apXnz5iktLc3rMjEAjC9dKgi4y0YSoanBfaTzAugW4QUAAAAQhaeeekq/+93vNGbMGB155JGaM2eOMjMz1dTUpF27dunpp5/WunXrdMYZZ+jUU0/1ulwMhGCxB8tGmHkBdIfwAgAAAIjCnj17dPPNNys/P7/Le8cdd5wkad++fXriiScSXRrixARHym57JSHXsh0DO+m8ALpDeAEAAABEYenSpb0eU1BQENVxGCQKR0qhatn9zTLDMuJ7rSYGdgKHw8BOAAAAIEbLli3r9vWvf/3rCa4EcRV0dxxRVUX8r9Xc3nnBshGgW4QXAAAAQIza2tq6vNba2qpIJOJBNYgXExzpPknE0M4DnRfDCC+A7rBsBAAAAIjSddddJ2OMWlpadP3113d6r6qqStOmTfOoMsRFodt5Yfd+IBPvazU3ShmZMg6/Xwa6Q3gBAAAAROlTn/qUJGnHjh06+eSTO143xigvL08zZ870qjTEQ26+lD5Mqtob/2s1NTLvAjgMwgsAAAAgSieddJIkaerUqRo9erS3xSDujDGSv1C2KkHLRjJYMgL0hJ4kAAAAIAobN27seH644OLQ45ACAkUJ6bywzU10XgCHQecFAAAAEIUXXnhBDz74oE444QTNmDFDJSUlysrKUmNjo/bs2aNt27bp+eef1/jx4zV37lyvy8UAMYFC2ffejv+FWDYCHBbhBQAAABCFSy65RDt37tRTTz2le+65RxUVB7fPLC4u1qxZs3TppZdq7NixHlaJARcokmprZJubZTIy4ned5iZ3xgaAbhFeAAAAAFEaN26cli9fLklqbm5WfX29cnJylBHPH2rhrUCR+1hdIY2KYzDV1ChTNCp+5wcGOWZeAAAAADF69913lZGRIb/fT3CR4kywPbyoqjj8gf3FwE7gsOi8AAAAAGK0evVqNTc3a/r06ZoxY4ZmzJihiRMnurtTILX43fDCVlYorn+6zcy8AA6H8AIAAACI0X333acPP/xQr732mrZt26Y//elPqq2t1fTp03X11Vd7XR4GUn6BlOaLa+eFtdadeUHnBdAjwgsAAACgD0aOHKm2tja1traqtbVVmzdvVk1NjddlYYAZJ03yB+O7bGR/s2QtnRfAYRBeAAAAADG688479frrr8vv92vGjBk64YQT9I1vfENZWfzwmZICRbLxDC+aGt1HwgugRwzsBAAAAGL01ltvyXEcjR8/XuPHj9eECRMILlKYCRRKVXvjd4Hm9vAig79DQE/ovAAAAABidPfddysUCmnbtm3atm2bHnvsMe3fv19HHnmkLrjgAq/Lw0ALjJRqqmVb9sfn/O2dF4bOC6BHdF4AAAAAfZCfn6+SkhIVFxersLBQoVBIr7zyitdlIR4Che5jdWV8zn9g2QgDO4Ee0XkBAAAAxOiWW27R9u3blZWVpRkzZmjOnDn66le/qlGjRnldGuLABEbKSlLVh/G5QHOT+0jnBdAjwgsAAAAgRvPnz9eyZctUVFTkdSlIhKD752zjNPfCNjHzAugN4QUAAAAQo5NOOsnrEpBI+QHJcaTKOO04wm4jQK+YeQEAAAAAh2HS0qSCYByXjRBeAL0hvAAAAACA3gQK47ZsRE3tMy8Y2An0iPACAAAAAHphAkVSVRyXjfjSZXys6gd6QngBAAAAAL0JFEmhatmWloE/d3OjlEnXBXA4hBcAAAAA0JtAkWQjaotH90VTEzuNAL0gvAAAAACAXpiAu11qZO8HA35u29TIsE6gF4QXAAAAANCb9vCirWLPwJ+7mfAC6A3hBQAAAAD0xh+UjFFbHDov1NTIshGgF4QXAAAAANAL40uX8vxqq4hTeMHATuCwCC8AAAAAIBqBwjgtG2mSofMCOCzCCwAAAACIggmMjN+yETovgMMivAAAAACAaAQKFamqkI20Dex5m5sY2An0gvACAAAAAKLhL5Ta2qTQvgE7pW1pkdpaGdgJ9ILwAgAAAACiYAKF7pPqvQN30uZG95HOC+CwCC8AAAAAIBr+IkmSHcjwoqk9vMhg5gVwOIQXAAAAABANf9B9jEN4Yei8AA6L8AIAAAAAomCysmVyRkhVA7lspMl9JLwADovwAgAAAACilFZYHKdlI4QXwOEQXgAAAABAlJzCkQzsBDzg87oAAAAAAJ39+te/1saNG2WMUV5eni666CL5/X5J0vr16/XMM8/IcRwtW7ZMxx57rMfVDi1pwZHSP18ZsPNZBnYCUaHzAgAAAEgyp59+um677Tbdeuutmj17tn7zm99Iknbt2qXy8nKtWbNG11xzjdauXatIJOJxtUNLWuFIqaFetrFhYE7IzAsgKoQXAAAAQJLJzs7ueN7c3CxjjCRpw4YNWrhwodLT01VUVKTi4mLt2LHDqzKHpLTCYvfJQC0dYeYFEBWWjQAAAABJ6MEHH9Rzzz2n7OxsXX/99ZKk6upqTZ06teMYv9+v6urqbj9fVlamsrIySdLq1asVDAZjur7P54v5M/GSTLW0FZdIknJbm5UxADXVGanecRQsKekIqaKVTPeFWrpHLd3rSy2EFwAAAIAHVq1apVAo1OX1s88+W/PmzdOSJUu0ZMkSrV+/Xv/3f/+nxYsXy1ob9flLS0tVWlra8XVlZWVM9QWDwZg/Ey/JVEtBQaEkqeZfb8kZP63f54uE9kkZWaqqqor5s8l0X6ile9TSvcPVUlJS0u3rhBcAAACAB1auXBnVcSeccIJWr16txYsXKxAIdPoht7q6umOQJxLDKQhIaT6pumJgTtjUwLBOIArMvAAAAACSzJ49ezqeb9y4seM3kXPnzlV5eblaWlpUUVGhPXv2aMqUKV6VOSQZx5EKAlLVwPwG2zY2SFnZvR8IDHF0XgAAAABJ5le/+pX27NkjY4yCwaDOP/98SdLYsWO1YMECrVixQo7jaPny5XIcfh+ZcP5C2YEa2NnYIGXnDMy5gBRGeAEAAAAkme985zs9vrdo0SItWrQogdXgo4y/UPaNfwzMyRrqpRG5A3MuIIUR0wIAAABALAKFUqhKtq2t/+dqbJDJovMC6A3hBQAAAADEwl8oRSJSqPttamPSWM/MCyAKhBcAAAAAEAPjd7dL1UDMvWhskDIJL4DeEF4AAAAAQCwCRZLU76GdtqVFam2h8wKIAuEFAAAAAMTCH3Qfqyr6d56mBveRmRdArwgvAAAAACAGJiNTGj6i/8tGGuvdRzovgF4RXgAAAABArPyFstWV/TtHo9t5YbIJL4DeEF4AAAAAQKz8Rf3vvGg40HnBshGgN4QXAAAAABAjEyiUqipkre37SRoPzLyg8wLoDeEFAAAAAMTKH5SaGg/OregD28jATiBahBcAAAAAECPTvl1qv5aOMLATiBrhBQAAAADEyl/oPlb1Y2jngc6LTMILoDeEFwAAAAAQq/bwwlZX9P0cjfXSsGEyPt8AFQWkLsILAAAAAIjViDwpzSf1Z7vUpkbmXQBRIrwAAAAAgBgZx5EKAv0LLxrqmXcBRInwAgAAAAD6wh+U3df3gZ22sYHOCyBKhBcAAAAA0AemICjtq+r7CRrpvACiRXgBAAAAAH3RHl7YSKRvn29sILwAokR4AQAAAAB94Q9Kba1SbU3fPt/YIMOyESAqhBcAAAAA0AemIOg+6evQTpaNAFEjvAAAAACAvvC3hxd9GNppW1ul/c0M7ASiRHgBAAAAAH1RUChJsn3pvGhqcB/pvACiQngBAAAAAH0xfISUPkza14fwovFAeEHnBRANwgsAAAAA6ANjjLvjSF86L9rDC5OVNcBVAamJ8AIAAAAA+soflKXzAog7wgsAAAAA6CPT586LeveRmRdAVAgvAAAAAKCv/EEpVC3b1hbTxyydF0BMCC8AAAAAoK/8QclGpJrq2D5H5wUQE8ILAAAAAOgj075dasxLR+i8AGKSNOHF448/rsWLFyscDntdCgAAAABExx+UpNiHdjbWS750mfT0OBQFpJ6kCC8qKyu1detWBYNBr0sBAAAAgOgVBNzHvnResGQEiFpShBfr1q3TOeec4+6TDAAAAACDRVaOlJElxdx50cCSESAGnocXGzdulN/v14QJE7wuBQAAAABiYoyR/EHZ6r0xfc7SeQHExJeIi6xatUqhUKjL62effbbWr1+va6+9NqrzlJWVqaysTJK0evVqT5aZ+Hw+lrckEPc7sbjficc9Tyzud2JxvxOL+w14qCAo7auK7TON9YQXQAwSEl6sXLmy29d37typiooKXXHFFZKkqqoqXXXVVbr55puVn5/f5fjS0lKVlpZ2fF1ZGWNr1gAIBoOeXHeo4n4nFvc78bjnicX9Tizud2Jxv2NTUlLidQlIIcYflN39TmwfamyQcrv+zAOgewkJL3oybtw43X///R1f/+d//qduvvlm5ebmelgVAAAAAMSgICiFQ7KtLTK+KHcPaWyQYeYFEDXPZ14AAAAAwKDmD0rWxrZ0hGUjQEySKry499576boAAAAAMKgYf/u8mSh3HLGRNqmpkd1GgBgkVXgBAAAAAINOQaEkyVZHOXemqdF9pPMCiBrhBQAAAAD0R4ydF2pscB8JL4CoEV4AAAAAQD+YjEwpe7gUbedFY737uWyWjQDRIrwAAAAAgP7yB2Wj7bxooPMCiBXhBQAAAAD0V0FQqt4b3bFNB8ILOi+AaBFeAAAAAEA/mUChVBVdeGGZeQHEjPACAAAAAPqrsFhqqJOtr+v92PaZF8okvACiRXgBAAAAAP1kCke5T/bu6f3gRpaNALEivAAAAACA/iosliTZimjCi3opLU0aNizORQGpg/ACAAAASFKPP/64Fi9erHA43PHa+vXrdfHFF+uSSy7R5s2bPawOnbSHF9r7Qe/HNjZIWdkyxsS3JiCFEF4AAAAASaiyslJbt25VMBjseG3Xrl0qLy/XmjVrdM0112jt2rWKRCIeVokDTEamlOeXoum8aGhgyQgQI8ILAAAAIAmtW7dO55xzTqffzm/YsEELFy5Uenq6ioqKVFxcrB07dnhYJTopKpaNYuaFbaxnpxEgRoQXAAAAQJLZuHGj/H6/JkyY0On16upqBQKBjq/9fr+qq6sTXB16YgpHxbBshM4LIBY+rwsAAAAAhqJVq1YpFAp1ef3ss8/W+vXrde2113Z5z1ob9fnLyspUVlYmSVq9enWn5SfR8Pl8MX8mXgZLLXUTJqm+/GkFRgx3l5H0oKqlWWn+EuX383saLPcl0aile4O9FsILAAAAwAMrV67s9vWdO3eqoqJCV1xxhSSpqqpKV111lW6++WYFAgFVVVV1HFtdXS2/39/teUpLS1VaWtrxdWVlZUz1BYPBmD8TL4OllkhOniSpcvs/ZUaP7/EcbXW1aivx9ft7Giz3JdGopXuDpZaSkpJuX2fZCAAAAJBExo0bp/vvv1/33nuv7r33XgUCAd1yyy3Kz8/X3LlzVV5erpaWFlVUVGjPnj2aMmWK1yWjnSkc5T7pbWhnYz3LRoAY0XkBAAAADBJjx47VggULtGLFCjmOo+XLl8tx+H1k0ihywwu7d4962gTVWis1NkqZDOwEYkF4AQAAACSxe++9t9PXixYt0qJFizyqBodjcoZL2cMPP7SzuVGyESmb8AKIxf/f3r3FRlnnYRx/3nbaorgt7dAypcgKUi7Q3SYsTStCqhTRGE0Ma5p4SMBoCC2mUWICeqEmxiBq00YsVm+Q9E4Tq8EbNwTFA9lNbW2sRRFqhSItAy0CBaZlZv57Mdl6mlmhnc7/nbffz03n0Isn/5m+8+sz74GaFgAAAACSpTAg8/8OG7kwEvvJYSPAVaG8AAAAAIAkcYr+5HKpJ3/65fcAXDHKCwAAAABIlsJiaSgoEw7Hfdqc6I/dmHt9CkMB6Y/yAgAAAACSpahYikal4WD85weOx86L8ZdZqc0FpDnKCwAAAABIEqcwELsRjH/oiBnsl4rnyXESXY8EQDyUFwAAAACQLEWx8sKcSnDSzoHjcoo5ZAS4WpQXAAAAAJAseQVSdnbcPS/MyDnp/FmpeJ6FYEB6o7wAAAAAgCRxHEcqLI6/58XA8djvsOcFcNUoLwAAAAAgmQoDUvCP5YUZOBa7QXkBXDXKCwAAAABIIqeoWDo1KBON/vaJgeNSdo5UUGgnGJDGKC8AAAAAIJkKA1L4svTz8G8eNgP9UqBETgb/hgFXi78aAAAAAEgip6g4duPU707aOXBcToBDRoCJoLwAAAAAgGQqmitJMj/9OP6QCV2Shk9xpRFggigvAAAAACCZ/EVSYJ7Ml5//8tjJnyRxpRFgoigvAAAAACCJHMeRU3mbdPigzOmTkiRzoj/25FzKC2AiKC8AAAAAIMmciipJkvnP/tgDA/1SZqZUWGwxFZC+KC8AAAAAIMmc2XOk0iUy//5ExhiZgeNSYbEcn892NCAtUV4AAAAAwBRwKm+TBo9Lx3qlwX5O1glMAuUFAAAAAEwB5x8rJJ9P5vO9UnCAk3UCk0B5AQAAAABTwJl5nfT3cpnP/yVFo+x5AUwC5QUAAAAATJGMitukcFiS5BTPtxsGSGOUFwAAAAAwVf62TLp2Zux2oMRuFiCNUV4AAAAAwBRxsrLkrFgj/XWRnJwZtuMAaYvr9AAAAADAFHL+uU6O7RBAmqO8AAAAAIAp5GSwwzswWfwVAQAAAAAAV6O8AAAAAAAArkZ5AQAAAAAAXI3yAgAAAAAAuBrlBQAAAAAAcDXKCwAAAAAA4GqUFwAAAAAAwNUoLwAAAAAAgKtRXgAAAAAAAFejvAAAAAAAAK5GeQEAAAAAAFyN8gIAAAAAALiaY4wxtkMAAAAAAAAkwp4XV2nr1q22I0wrrHdqsd6px5qnFuudWqx3arHeSDY3vafIEh9Z4iNLfOmehfICAAAAAAC4GuUFAAAAAABwtcznn3/+edsh0s3ChQttR5hWWO/UYr1TjzVPLdY7tVjv1GK9kWxuek+RJT6yxEeW+NI5CyfsBAAAAAAArsZhIwAAAAAAwNUoLwAAAAAAgKv5bAdIF11dXdq1a5ei0aiqq6t133332Y7kaadPn1Zzc7N+/vlnOY6j1atX6+6777Ydy/Oi0ai2bt2qgoICV11KyYsuXLiglpYW9ff3y3Ec1dbWavHixbZjedaHH36offv2yXEcXX/99RFw3OgAAAg+SURBVKqrq1N2drbtWJ6yc+dOdXZ2Ki8vTw0NDZKkkZERNTY26tSpUyosLNSTTz6p6667znJSb4i33q2trero6JDP59OcOXNUV1enmTNnWk6KdGVz9nXT9iTRTJrqPGNjY3ruuecUDocViURUWVmpmpoaq9vZ38+NNrNs2rRJM2bMUEZGhjIzM/XSSy9ZyRNvvps7d27Kc5w4cUKNjY3j94PBoGpqalRVVWXlNYo3h42NjV19FoM/FYlEzOOPP24GBwfN5cuXzVNPPWX6+/ttx/K04eFh09vba4wx5uLFi6a+vp41T4E9e/aYpqYms23bNttRPG/Hjh1m7969xhhjLl++bEZGRiwn8q6hoSFTV1dnRkdHjTHGNDQ0mI8//thuKA/q6ekxvb29ZvPmzeOPtba2mra2NmOMMW1tbaa1tdVWPM+Jt95dXV0mHA4bY2Jrz3pjomzPvm7aniSaSVOdJxqNmkuXLhljYnPD008/bQ4dOmR1O/v7udFmlrq6OnP27NnfPGYjT7z5zvZnYSQSMY899pgJBoNWsiSawyaShcNGrsCRI0cUCAQ0Z84c+Xw+LV++XO3t7bZjeVp+fv742WevueYalZSUaHh42HIqbxsaGlJnZ6eqq6ttR/G8ixcv6ttvv9WqVaskST6fj29Hp1g0GtXY2JgikYjGxsaUn59vO5LnLFmy5A/fmLS3t6uqqkqSVFVVxWdnEsVb77KyMmVmZkqSFi9ezOcmJsz27Oum7UmimTTVeRzH0YwZMyRJkUhEkUhEjuNYW5d4c6PbtvmpzpNovrO9Lt3d3QoEAiosLLSWJd4cNpEsHDZyBYaHh+X3+8fv+/1+HT582GKi6SUYDKqvr0+LFi2yHcXT3n77bT388MO6dOmS7SieFwwGlZubq507d+ro0aNauHCh1q9fPz6UILkKCgp07733qra2VtnZ2SorK1NZWZntWNPC2bNnx4ui/Px8nTt3znKi6WPfvn1avny57RhIU26cfd2wPfn1TGojTzQa1ZYtWzQ4OKg777xTpaWl1tYl3txo+zV68cUXJUl33HGHVq9enfI8ieY72+vyxRdf6NZbb5Vk5zVKNIdNJAt7XlwBE+dqso7jWEgy/YRCITU0NGj9+vW69tprbcfxrI6ODuXl5bnqus9eFolE1NfXpzVr1ujll19WTk6O3n//fduxPGtkZETt7e1qbm7Wm2++qVAopE8//dR2LGDKvPfee8rMzNTKlSttR0GaYvb9IzfMpBkZGXrllVfU0tKi3t5eHTt2zEoON86NL7zwgrZv365nnnlGH330kQ4ePJjyDG6c78LhsDo6OlRZWWktQzLnMMqLK+D3+zU0NDR+f2hoiF2OUyAcDquhoUErV65URUWF7TiedujQIX355ZfatGmTmpqa9M033+i1116zHcuz/H6//H6/SktLJUmVlZXq6+uznMq7uru7VVRUpNzcXPl8PlVUVOj777+3HWtayMvL05kzZyRJZ86cUW5uruVE3vfJJ5+oo6ND9fX10/6fTUycG2dfm9uTeDOpzTwzZ87UkiVL1NXVZSVHornR5poUFBRIir0u5eXlOnLkSMrzJJrvbK7LV199pQULFmjWrFmS7LxvE81hE8lCeXEFbrzxRg0MDCgYDCocDuvAgQNatmyZ7VieZoxRS0uLSkpKdM8999iO43kPPvigWlpa1NzcrCeeeEI333yz6uvrbcfyrFmzZsnv9+vEiROSYhv1efPmWU7lXbNnz9bhw4c1OjoqY4y6u7tVUlJiO9a0sGzZMu3fv1+StH//fpWXl1tO5G1dXV364IMPtGXLFuXk5NiOgzTmxtnX1vYk0Uya6jznzp3ThQsXJMWuPPK/zzIb65JobrT1GoVCofHDV0KhkL7++mvNnz8/5XkSzXc2Pwt/fciIZOfvKNEcNpEsjom3Xxj+oLOzU7t371Y0GtXtt9+utWvX2o7kad99952effZZzZ8/f/ybowceeEBLly61nMz7enp6tGfPHi6VOsV+/PFHtbS0KBwOq6ioSHV1dVxCcgq98847OnDggDIzM3XDDTdo48aNysrKsh3LU5qamnTw4EGdP39eeXl5qqmpUXl5uRobG3X69GnNnj1bmzdv5n2eJPHWu62tTeFweHyNS0tLtWHDBstJka5szr5u2p4kmklLS0tTmufo0aNqbm5WNBqVMUa33HKL7r//fp0/f97qdvbXc6OtLCdPntSrr74qKXboxooVK7R27VoreeLNd8YYK+syOjqq2tpavf766+OHOtl6jeLNYaFQ6KqzUF4AAAAAAABX47ARAAAAAADgapQXAAAAAADA1SgvAAAAAACAq1FeAAAAAAAAV6O8AAAAAAAArkZ5AQAAAAAAXI3yAgAAAAAAuBrlBQAAAAAAcDXKCwDwiMHBQT3yyCP64YcfJEnDw8N69NFH1dPTYzkZAAAAMDmUFwDgEYFAQA899JB27Nih0dFRvfHGG6qqqtJNN91kOxoAAAAwKY4xxtgOAQBInu3btysYDMpxHG3btk1ZWVm2IwEAAACTwp4XAOAx1dXV6u/v11133UVxAQAAAE+gvAAADwmFQtq9e7dWrVqld999VyMjI7YjAQAAAJNGeQEAHrJr1y4tWLBAGzdu1NKlS/XWW2/ZjgQAAABMGuUFAHhEe3u7urq6tGHDBknSunXr1NfXp88++8xyMgAAAGByOGEnAAAAAABwNfa8AAAAAAAArkZ5AQAAAAAAXI3yAgAAAAAAuBrlBQAAAAAAcDXKCwAAAAAA4GqUFwAAAAAAwNUoLwAAAAAAgKtRXgAAAAAAAFf7L377S5xHQ5MQAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 1080x720 with 3 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"#plot trajectory\n",
|
|
"grid = plt.GridSpec(2, 3)\n",
|
|
"\n",
|
|
"plt.figure(figsize=(15,10))\n",
|
|
"\n",
|
|
"plt.subplot(grid[0:2, 0:2])\n",
|
|
"plt.plot(track[0,:],track[1,:],\"b+\")\n",
|
|
"plt.plot(x_sim[0,:],x_sim[1,:])\n",
|
|
"plt.axis(\"equal\")\n",
|
|
"plt.ylabel('y')\n",
|
|
"plt.xlabel('x')\n",
|
|
"\n",
|
|
"plt.subplot(grid[0, 2])\n",
|
|
"plt.plot(u_sim[0,:])\n",
|
|
"plt.ylabel('v(t) [m/s]')\n",
|
|
"\n",
|
|
"plt.subplot(grid[1, 2])\n",
|
|
"plt.plot(np.degrees(u_sim[1,:]))\n",
|
|
"plt.ylabel('w(t) [deg/s]')\n",
|
|
"\n",
|
|
"plt.tight_layout()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## OBSTACLE AVOIDANCE\n",
|
|
"see dccp paper for reference"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import dccp\n",
|
|
"track = compute_path_from_wp([0,3,4,6,10,13],\n",
|
|
" [0,0,2,4,3,3],0.25)\n",
|
|
"\n",
|
|
"obstacles=np.array([[4,6],[2,4]])\n",
|
|
"obstacle_radius=0.5\n",
|
|
"\n",
|
|
"def to_vehic_frame(pt,pos_x,pos_y,theta):\n",
|
|
" dx = pt[0] - pos_x\n",
|
|
" dy = pt[1] - pos_y\n",
|
|
"\n",
|
|
" return [dx * np.cos(-theta) - dy * np.sin(-theta),\n",
|
|
" dy * np.cos(-theta) + dx * np.sin(-theta)]\n",
|
|
" \n",
|
|
"# track = compute_path_from_wp([0,5,7.5,10,12,13,13,10],\n",
|
|
"# [0,0,2.5,2.5,0,0,5,10],0.5)\n",
|
|
"\n",
|
|
"sim_duration = 80 #time steps\n",
|
|
"opt_time=[]\n",
|
|
"\n",
|
|
"x_sim = np.zeros((N,sim_duration))\n",
|
|
"u_sim = np.zeros((M,sim_duration-1))\n",
|
|
"\n",
|
|
"MAX_SPEED = 1.25\n",
|
|
"MIN_SPEED = 0.75\n",
|
|
"MAX_STEER_SPEED = 1.57\n",
|
|
"\n",
|
|
"# Starting Condition\n",
|
|
"x0 = np.zeros(N)\n",
|
|
"x0[0] = 0\n",
|
|
"x0[1] = -0.25\n",
|
|
"x0[2] = np.radians(-0)\n",
|
|
"x_sim[:,0]=x0\n",
|
|
" \n",
|
|
"#starting guess\n",
|
|
"u_bar = np.zeros((M,T))\n",
|
|
"u_bar[0,:]=0.5*(MAX_SPEED+MIN_SPEED)\n",
|
|
"u_bar[1,:]=0.00\n",
|
|
"\n",
|
|
"for sim_time in range(sim_duration-1):\n",
|
|
" \n",
|
|
" iter_start=time.time()\n",
|
|
" \n",
|
|
" #compute coefficients\n",
|
|
" K=road_curve(x_sim[:,sim_time],track)\n",
|
|
" \n",
|
|
" #compute opstacles in ref frame\n",
|
|
" o_=[]\n",
|
|
" for j in range(2):\n",
|
|
" o_.append(to_vehic_frame(obstacles[:,j],x_sim[0,sim_time],x_sim[1,sim_time],x_sim[2,sim_time]) )\n",
|
|
" \n",
|
|
" # dynamics starting state w.r.t vehicle frame\n",
|
|
" x_bar=np.zeros((N,T+1))\n",
|
|
" \n",
|
|
" #prediction for linearization of costrains\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",
|
|
" A,B,C=get_linear_model(xt,ut)\n",
|
|
" xt_plus_one = np.squeeze(np.dot(A,xt)+np.dot(B,ut)+C)\n",
|
|
" x_bar[:,t]= xt_plus_one\n",
|
|
" \n",
|
|
" #CVXPY Linear MPC problem statement\n",
|
|
" cost = 0\n",
|
|
" constr = []\n",
|
|
" x = cp.Variable((N, T+1))\n",
|
|
" u = cp.Variable((M, T))\n",
|
|
" \n",
|
|
" #Prediction Horizon\n",
|
|
" for t in range(T):\n",
|
|
"\n",
|
|
" #cost += 30*cp.sum_squares(x[2,t]-np.arctan(df(x_bar[0,t],K))) # psi\n",
|
|
" cost += 50*cp.sum_squares(x[2,t]-np.arctan2(df(x_bar[0,t],K),x_bar[0,t])) # psi\n",
|
|
" cost += 20*cp.sum_squares(f(x_bar[0,t],K)-x[1,t]) # cte\n",
|
|
"\n",
|
|
" # Actuation rate of change\n",
|
|
" if t < (T - 1):\n",
|
|
" cost += cp.quad_form(u[:, t + 1] - u[:, t], 100*np.eye(M))\n",
|
|
" \n",
|
|
" # Actuation effort\n",
|
|
" cost += cp.quad_form( u[:, t],1*np.eye(M))\n",
|
|
" \n",
|
|
" # Kinrmatics Constrains (Linearized model)\n",
|
|
" A,B,C=get_linear_model(x_bar[:,t],u_bar[:,t])\n",
|
|
" constr += [x[:,t+1] == A@x[:,t] + B@u[:,t] + C.flatten()]\n",
|
|
" \n",
|
|
" # Obstacle Avoidance Contrains\n",
|
|
" for j in range(2):\n",
|
|
" constr += [ cp.norm(x[0:2,t]-o_[j],2) >= obstacle_radius ]\n",
|
|
"\n",
|
|
" # sums problem objectives and concatenates constraints.\n",
|
|
" constr += [x[:,0] == x_bar[:,0]] #<--watch out the start condition\n",
|
|
" constr += [u[0, :] <= MAX_SPEED]\n",
|
|
" constr += [u[0, :] >= MIN_SPEED]\n",
|
|
" constr += [cp.abs(u[1, :]) <= MAX_STEER_SPEED]\n",
|
|
" \n",
|
|
" # Solve\n",
|
|
" prob = cp.Problem(cp.Minimize(cost), constr)\n",
|
|
" solution = prob.solve(method=\"dccp\", verbose=False)\n",
|
|
" \n",
|
|
" #retrieved optimized U and assign to u_bar to linearize in next step\n",
|
|
" u_bar=np.vstack((np.array(u.value[0, :]).flatten(),\n",
|
|
" (np.array(u.value[1, :]).flatten())))\n",
|
|
" \n",
|
|
" u_sim[:,sim_time] = u_bar[:,0]\n",
|
|
" \n",
|
|
" # Measure elpased time to get results from cvxpy\n",
|
|
" opt_time.append(time.time()-iter_start)\n",
|
|
" \n",
|
|
" # move simulation to t+1\n",
|
|
" tspan = [0,dt]\n",
|
|
" x_sim[:,sim_time+1] = odeint(kinematics_model,\n",
|
|
" x_sim[:,sim_time],\n",
|
|
" tspan,\n",
|
|
" args=(u_bar[:,0],))[1]\n",
|
|
" \n",
|
|
"print(\"CVXPY Optimization Time: Avrg: {:.4f}s Max: {:.4f}s Min: {:.4f}s\".format(np.mean(opt_time),\n",
|
|
" np.max(opt_time),\n",
|
|
" np.min(opt_time))) "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#plot trajectory\n",
|
|
"grid = plt.GridSpec(2, 3)\n",
|
|
"\n",
|
|
"plt.figure(figsize=(15,10))\n",
|
|
"\n",
|
|
"ax=plt.subplot(grid[0:2, 0:2])\n",
|
|
"plt.plot(track[0,:],track[1,:],\"b+\")\n",
|
|
"plt.plot(x_sim[0,:],x_sim[1,:])\n",
|
|
"#obstacles\n",
|
|
"circle1=plt.Circle((obstacles[0,0], obstacles[1,0]), obstacle_radius, color='r')\n",
|
|
"circle2=plt.Circle((obstacles[0,1], obstacles[1,1]), obstacle_radius, color='r')\n",
|
|
"plt.axis(\"equal\")\n",
|
|
"plt.ylabel('y')\n",
|
|
"plt.xlabel('x')\n",
|
|
"\n",
|
|
"ax.add_artist(circle1)\n",
|
|
"ax.add_artist(circle2)\n",
|
|
"\n",
|
|
"plt.subplot(grid[0, 2])\n",
|
|
"plt.plot(u_sim[0,:])\n",
|
|
"plt.ylabel('v(t) [m/s]')\n",
|
|
"\n",
|
|
"plt.subplot(grid[1, 2])\n",
|
|
"plt.plot(np.degrees(u_sim[1,:]))\n",
|
|
"plt.ylabel('w(t) [deg/s]')\n",
|
|
"\n",
|
|
"\n",
|
|
"plt.tight_layout()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python [conda env:jupyter] *",
|
|
"language": "python",
|
|
"name": "conda-env-jupyter-py"
|
|
},
|
|
"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
|
|
}
|