Commit 31b34c77 by Luz Andrea Silva Torres

### Mini project

parent ba0278b4

35.9 KB

This source diff could not be displayed because it is too large. You can view the blob instead.
 { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# MULTIDIMENSIONAL KALMAN FILTER" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 Vehicle location estimation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The aim of this notebook is to estimate the location of a vehicle in the XY plane using multidimensional Kalman filter. This is an implementation of a simple exampl proposed in an educational webpage.\n", "\n", "It is suppose that the vehicle has an onboard location sensor that reports X and Y coordinates of the system, this represents the observations with wich are corrected the projections.\n", "\n", "It is assumed a constant acceleration dynamics. \n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The state of the system is given by the position, velocity and acceleration of the car in the X and Y coordinates. The state extrapolated vehicle for time n+1 can be describe by the following system of equations:\n", "\n", "\\n", "\\begin{split}\n", "\\hat{x}_{n+1,n}=\\hat{x}_{n,n}+\\hat{\\dot{x}}_{n,n}\\Delta t+\\frac{1}{2}\\hat{\\ddot{x}}_{n,n}\\Delta t^2 \\\\\n", "\\hat{\\dot{x}}_{n+1,n} = \\hat{\\dot{x}}_{n,n}+ \\hat{\\ddot{x}}_{n,n}\\Delta t \\\\\n", "\\hat{\\ddot{x}}_{n+1,n} = \\hat{\\ddot{x}}_{n,n} \\\\\n", "\\hat{y}_{n+1,n}=\\hat{y}_{n,n}+\\hat{\\dot{y}}_{n,n}\\Delta t+\\frac{1}{2}\\hat{\\ddot{y}}_{n,n}\\Delta t^2 \\\\\n", "\\hat{\\dot{y}}_{n+1,n} = \\hat{\\dot{y}}_{n,n}+ \\hat{\\ddot{y}}_{n,n}\\Delta t \\\\\n", "\\hat{\\ddot{y}}_{n+1,n} = \\hat{\\ddot{y}}_{n,n} \\\\\n", "\\end{split}\n", "\\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Extended Kalman filter equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Kalman filter equations are given by:\n", "\n", "\\n", " \\left.\\begin{aligned}\n", " \\hat{x}_{n+1,n}= F\\hat{x}_{n,n}\\\\\n", " P_{n+1,n}= FP_{n,n}F^T+Q\n", " \\end{aligned}\n", " \\right\\}\n", " \\qquad \\text{Prediction equations}\n", "\\n", "\n", "where $F$ is the state transition matrix (comes from the dynamic model), $x$ is the state vector (what we want to predict), $Q$ is the process noise uncertainty and $P$ is the estimate uncertainty (covariance matrix).\n", "\n", "\\n", " \\left.\\begin{aligned}\n", " \\hat{x}_{n,n}= \\hat{x}_{n,n-1} + K_n(z_n- H\\hat{x}_{n,n-1}\\\\\n", " P_{n,n}= (I-K_nH)P_{n,n-1}(I-K_nH)^T+ K_n R_n K_n^T\\\\\n", " K_{n}= P_{n,n-1}H^T(HP_{n,n-1}H^T+R_n)^{-1}\n", " \\end{aligned}\n", " \\right\\}\n", " \\qquad \\text{Update (correction with observations)}\n", "\\n", "\n", "where $I$ is the indentity matrix, $K_n$ is the Kalman Gain, $P_{n,n}$ is the covariance update, $\\hat{x}_{n,n}$ is the state update, $z_n$ is the measurement equation that comes from the observations, $H$ is the observation matrix, that is a linear operator that allows to consider only the variables that are observed, and $R_n$ is the measurement uncertainty." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from numpy.linalg import inv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first step is to create the Kalman Filter arrays:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Creation of Kalman's filter arrays\n", "Xn = np.zeros([6,1])\n", "F = np.identity(6)\n", "Pn = np.zeros([6,6])\n", "Q = np.zeros([6,6])\n", "Kn = np.zeros([6,2])\n", "H = np.zeros([2,6])\n", "Rn = np.zeros([2,2])\n", "dt = 1\n", "# Variance of the process, \"random acceleration\":\n", "var_a = 0.15**2\n", "# Observation or measurement variance:\n", "var_mxy = 3**2\n", "\n", "#Observation matrix zn\n", "Zn = np.zeros([2,1])\n", "'observations'\n", "x_o = np.array([-393.66,-375.93,-351.04,-328.96,-299.35,-273.36,-245.89,-222.58,-198.03,-174.17,-146.32,-123.72,-103.47,-78.23,-52.63,-23.34,25.96,49.72,76.94,95.38,119.83,144.01,161.84,180.56,201.42,222.62,\t239.4,252.51,266.26,271.75,277.4,294.12,301.23,291.8,299.89,0])\n", "y_o = np.array([300.4,301.78,295.1,305.19,301.06,302.05,300,303.57,296.33,297.65,297.41,299.61,299.6,302.39,295.04,300.09,294.72,298.61,294.64,284.88,272.82,264.93,251.46,241.27,222.98,203.73,184.1,166.12,138.71,119.71,100.41,79.76,50.62,32.99,2.14,0])\n", "Z_o = np.zeros([2,36])\n", "Z_o[0,:] = x_o\n", "Z_o[1,:] = y_o" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the numerical application is assume that a vehicle is moving in a straight line in the X direction with a constant velocity. After traveling 400 meters the vehicle turns right, with a turning radius of 300 meters. During the turning maneuver, the vehicle experiences acceleration due to the circular motion (an angular acceleration). " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Y')" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAelUlEQVR4nO3deXhV9b3v8fc3A0kYQgKEEJIwCCiDCmJEkWMdqELpgNZjD7b1cm0r9lw9aseLem+H42OnU+u9p1Vb2trSXi1ah6qtVpDCqbMGBWWQEgYlEENAZkhCku/9Yy/iFhPGrL3WTj6v59nP3vu319rrQ5r6yZrN3REREQHIiDqAiIjEh0pBRERaqRRERKSVSkFERFqpFEREpFVW1AFORL9+/XzIkCFRxxARSStLlizZ6u5FbX2W1qUwZMgQKisro44hIpJWzOzt9j7T5iMREWmlUhARkVYqBRERaaVSEBGRVioFERFpFVopmFmumb1iZsvMbIWZfTcY72NmC8xsTfBcmDTPzWZWZWarzWxKWNlERKRtYa4pNAAXuftYYBww1czOAWYDC919BLAweI+ZjQZmAGOAqcDdZpYZYj4RETlEaOcpeOKa3HuCt9nBw4HpwAXB+FxgMfA/g/F57t4ArDezKmAC8GJHZ3t3Zz33v9zuYboi8WVG926Z9MjJomdOJj26ZdEzJ4sewaNnThYF3bPJzdbfU3J8Qj15LfhLfwkwHLjL3V82s2J3rwFw9xoz6x9MXgq8lDR7dTB26HfOAmYBDBo06Lhy1e6q56eLqo5rXpEoHc3tT8xgYO88TirqwbCinpxU1IOT+iWeB+TnkpFh4QeVtBVqKbh7MzDOzAqAR83s1MNM3tZv6of+L+Duc4A5ABUVFcd1h6Cx5QWs//7Hj2dWkUi1tDj7DzSzt6GJPQ1N7G1oZk9DE/sa339ft7uBdVv3sK5uL3+s3MjexubW+fOyMxnarwenl/Vm4rC+TBzWl/69ciP8F0ncpOQyF+6+w8wWk9hXUGtmJcFaQgmwJZisGihPmq0M2JyKfCLpIiPDWjcV9T/y5Lg7tbsaWFe3h7Vb9yae6/bylzdrmPfqRgCG9+/JucP6cu6wvpw9tC+FPbqF+4+QWLOwbsdpZkXAgaAQ8oD5wA+B84Ft7v4DM5sN9HH3b5rZGOB+EvsRBpLYCT0iWNtoU0VFhevaRyLHrrnFWbF5Jy+s3caLa7fx6ob32NfYjBmMGpDPucP6Mu30Es4oL8BMm5s6GzNb4u4VbX4WYimcTmJHciaJo5wedPd/N7O+wIPAIOAd4Ap3fy+Y51bgC0ATcJO7P3W4ZagURDrGgeYW3qjewQtV23hx3TYq395OY1MLQ/v14NJxpVx2RimD+naPOqZ0kEhKIRVUCiLh2F1/gKeWv8ujr23ipfXbcIeKwYVcNr6UT5w2kN7ds6OOKCdApSAix23zjv38aekmHnltE1Vb9tAtM4OLRvbniooyLjylv45mSkMqBRE5Ye7O8k27eOT1ap5Ytpmtexo5ubgn/+OC4Xzi9BKyMnXVnHShUhCRDnWguYUn36zh7kVrWV27m/I+eXz5/GFcPr5MJ86lAZWCiISipcVZ+NYWfraoimUbd9C/Vw7XnHcSnz17ED1y0vrGjp2aSkFEQuXuvLh2G3ctruL5qm30zsvm6klDuPrcodopHUMqBRFJmdfe2c7di9byzKpaCrpn880pI/mXs8rJ1A7p2DhcKWjPkIh0qPGDCvnVzAqevOE8Ti7uxS2Pvslldz/P0o07oo4mR0GlICKhGD0wnwdmncP/nTGOd3fWc9ndzzP74Td4b29j1NHkMFQKIhIaM2P6uFIWfu18vvRPQ3loSTUX/ngxv3/pbZpb0nfTdWemUhCR0PXKzebWj4/mqRvPY3RJPv/7T8uZftdzLHl7e9TR5BAqBRFJmRHFvbj/mrP56ZVnULe7gcvveYF/f2IlDU3tXvdSUkylICIpZWZ8cuxA/va1C5g5cTD3Pr+ey+95gXV1e448s4ROpSAikeiRk8V3p5/KL/9bBdXb9/OJnz7HQ0uqSefD5DsDlYKIROri0cU8deN5nFbam6//cRlfeWApu+sPRB2ry1IpiEjkSnrncf815/DVi0/m8WWb+cRPn2OZzmuIhEpBRGIhM8O4YfIIHrh2IgeaWrj8nheY8/e1tOjQ1ZRSKYhIrJw1pA9P3ngek0f153tPvsV//+2r7NyvzUmpolIQkdgp6N6Nn3/+TG679FReXLuVf77nBTbt2B91rC5BpSAisWRmXHXOYOZePYF3d9Vz2V3Ps3zTzqhjdXoqBRGJtXOH9+OhL59LVobxmV+8yKLVW6KO1KmpFEQk9k4Z0ItHr5vEkL49+NLcSua98k7UkTotlYKIpIXi/Fwe/PJEJg3vx+xH3uTHT6/WiW4hUCmISNromZPFr2dWMOOscn62qIqvPriMxqaWqGN1KrqJqoiklezMDL7/6dMoK8zjx/P/wbs76/n5VWfSO0+3/ewIoa0pmFm5mS0ys1VmtsLMbgzGv2Nmm8xsafCYljTPzWZWZWarzWxKWNlEJL2ZGddfNIKffGYslW+/x4w5L7Fjn27e0xHC3HzUBHzN3UcB5wDXmdno4LM73X1c8HgSIPhsBjAGmArcbWaZIeYTkTT36fFl/GrmWazdsoeZv3lV10zqAKGVgrvXuPtrwevdwCqg9DCzTAfmuXuDu68HqoAJYeUTkc7h/JOLuOtz41mxaSdf/G0l+xqboo6U1lKyo9nMhgBnAC8HQ9eb2Rtmdq+ZFQZjpcDGpNmqaaNEzGyWmVWaWWVdXV2YsUUkTVw8upg7/2UclW+/x7W/X0L9Ad2053iFXgpm1hN4GLjJ3XcB9wDDgHFADXDHwUnbmP1Dx5u5+xx3r3D3iqKionBCi0ja+eTYgfzw8tN5ds1Wrr//NQ4066ik4xFqKZhZNolCuM/dHwFw91p3b3b3FuCXvL+JqBooT5q9DNgcZj4R6VyuqCjntuljeGbVFm56YCnNusLqMQvz6CMDfg2scvefJI2XJE12GbA8eP04MMPMcsxsKDACeCWsfCLSOV01cQi3TBvJX96o4ZsPvaFLbx+jMM9TmARcBbxpZkuDsVuAK81sHIlNQxuAawHcfYWZPQisJHHk0nXurg2DInLMZn1kGPsbW7jzmX+Q1y2D26afSuLvVDmS0ErB3Z+j7f0ETx5mntuB28PKJCJdxw2Th7PvQBO/+K915GVncsu0USqGo6AzmkWkUzIzZk8dSX1jM798dj19euTwrxcMizpW7KkURKTTMjO+/ckxvLfvAD96+i1GlvTiwlP6Rx0r1nRBPBHp1DIyjB9dfjqjBuRzwx9eZ/3WvVFHijWVgoh0enndMvnFVWeSlWFc87tK9jTorOf2qBREpEso79Oduz47nvVb9/LVB5bqUNV2qBREpMs4d3g/bp02ivkra/np36qijhNLKgUR6VKunjSET48v5c5n/sH8Fe9GHSd2VAoi0qWYGd+77DROL+vNVx5Yypra3VFHihWVgoh0ObnZmfz882eS1y2TWb9fws79ug/DQSoFEemSBhbkcc/nz2Tje/u4cd7runheQKUgIl3WWUP68J1PjWHx6jp+smB11HFiQaUgIl3a584exJUTyrlr0VoWr94SdZzIqRREpEs7eCmMEf17MvvhN7v8/gWVgoh0ebnZmfz4irHU7Wngtj+vjDpOpFQKIiLA2PICvnz+STy0pJq/vVUbdZzIqBRERAI3TB7BKcW9EpuR9nXNzUgqBRGRQE5WJnd8Zizb9jby3SdWRB0nEioFEZEkp5b25roLh/PI65u65GUwVAoiIoe4/sLhjCrJ55ZHl7N9b2PUcVJKpSAicohuWRncccVYduxr5NuPd63NSCoFEZE2jB6Yzw2TR/D4ss089WZN1HFSRqUgItKOf71gGKeW5vO//rScbXsaoo6TEioFEZF2ZGdmcMcV49hVf4BvPdY1NiOpFEREDuOUAb246aMn85c3a/jzG5ujjhO60ErBzMrNbJGZrTKzFWZ2YzDex8wWmNma4LkwaZ6bzazKzFab2ZSwsomIHItrP3ISY8sL+PZjK9hd37lPagtzTaEJ+Jq7jwLOAa4zs9HAbGChu48AFgbvCT6bAYwBpgJ3m1lmiPlERI5KVmYGt00fw7a9jfziv9ZFHSdUoZWCu9e4+2vB693AKqAUmA7MDSabC1wavJ4OzHP3BndfD1QBE8LKJyJyLE4vK+BTYwfyq+fW8e7O+qjjhCYl+xTMbAhwBvAyUOzuNZAoDqB/MFkpsDFptupg7NDvmmVmlWZWWVdXF2puEZFk35hyCi0tcMf8zntDntBLwcx6Ag8DN7n7rsNN2sbYh+6P5+5z3L3C3SuKioo6KqaIyBGV9+nOzHMH89Br1ayqOdx/ztJXqKVgZtkkCuE+d38kGK41s5Lg8xLg4K2OqoHypNnLgM6/q19E0sr1F44gPzeb7z/1VtRRQhHm0UcG/BpY5e4/SfrocWBm8Hom8FjS+AwzyzGzocAI4JWw8omIHI/e3bP5t4uG8/d/1PHsms63CTvMNYVJwFXARWa2NHhMA34AXGxma4CLg/e4+wrgQWAl8FfgOndvDjGfiMhxuWriYMr75PG9J9+iueVDW7nTWlZYX+zuz9H2fgKAye3Mcztwe1iZREQ6Qk5WJt+YMpIb/vA6j76+iX8+syzqSB1GZzSLiByHT55ewtiy3twxfzX1BzrPRg2VgojIcTAzbpk2ipqd9fz6ufVRx+kwKgURkeN09kl9+eioYu5ZvLbTXEVVpSAicgJmf2wk+w80858L10QdpUOoFERETsDw/j2ZcVY59738Duvq9kQd54SpFERETtBNHz2ZnKwMfvjX9D+hTaUgInKCinrlcO35w3h6RS2vbngv6jgnRKUgItIBvnTeUPr1zEn7fQsqBRGRDtC9WxYzJw7m2TVbqdqyO+o4x02lICLSQT579iC6ZWXwm+c3RB3luKkUREQ6SN+eOUwfO5BHXtvEjn2NUcc5LioFEZEOdPWkoew/0My8VzceeeIYUimIiHSg0QPzOeekPvzuhQ00NbdEHeeYqRRERDrYFyYNZfPOeuavrI06yjFTKYiIdLDJo4op75PHvWl4oTyVgohIB8vMMGZOHELl29t5s3pn1HGOiUpBRCQEnzmrnB7dMvnN8+m1tqBSEBEJQX5uNldUlPPEG5vZsqs+6jhHTaUgIhKSmecOoanF+X8vvxN1lKOmUhARCcnQfj246JT+3PfS22lzy06VgohIiK6eNJRtext5YtnmqKMcFZWCiEiIJg3vy8nFPfnN8xtw96jjHJFKQUQkRGbG1ZOGsrJmFy+vj/+9FlQKIiIhu+yMUgq7Z6fF4amhlYKZ3WtmW8xsedLYd8xsk5ktDR7Tkj672cyqzGy1mU0JK5eISKrlZmdy5YRBLFhZy8b39kUd57DaLQUze9LMhpzAd/8WmNrG+J3uPi54PBksazQwAxgTzHO3mWWewLJFRGLlqomDMTPmvrAh6iiHdbg1hd8C883sVjPLPtYvdve/A0e7AW06MM/dG9x9PVAFTDjWZYqIxFVJ7zymnVbCA5Ub2dPQFHWcdrVbCu7+IHAGkA9UmtnXzeyrBx8nsMzrzeyNYPNSYTBWCiRffLw6GPsQM5tlZpVmVllXV3cCMUREUmvmxMHsrm9i/op3o47SriPtUzgA7AVygF6HPI7HPcAwYBxQA9wRjFsb07Z57Ja7z3H3CnevKCoqOs4YIiKpN35QIQPyc3k6xqWQ1d4HZjYV+AnwODDe3U9474i7t15c3Mx+Cfw5eFsNlCdNWgakx5keIiJHKSPDuHh0MX9cspH9jc3kdYvfrtPDrSncClzh7rM7ohAAzKwk6e1lwMEjkx4HZphZjpkNBUYAr3TEMkVE4mTKmAHUH2jh2TXx3Pzd7pqCu593Il9sZn8ALgD6mVk18G3gAjMbR2LT0Abg2mBZK8zsQWAl0ARc5+7pcaEQEZFjcPZJfcjPzWL+ylouGTMg6jgf0m4pnCh3v7KN4V8fZvrbgdvDyiMiEgfZmRlMHlXMwlW1NDW3kJUZr3OI45VGRKQLuGR0Mdv3HeDVDdujjvIhKgURkRQ7/5QicrIyYnkUkkpBRCTFunfL4rwR/ViwsjZ2V05VKYiIROCSMQPYtGM/KzbvijrKB6gUREQiMHlkfzKM2J3drFIQEYlA3545nDWkD0+vqD3yxCmkUhARicglYwawunY3G7bujTpKK5WCiEhELhldDMD8lfHZhKRSEBGJSHmf7owuyWd+jDYhqRRERCI0ZcwAlryznbrdDVFHAVQKIiKRumRMMe7wzKp4rC2oFEREIjRyQC8G9ekem7ObVQoiIhEyM6aMKeaFqm3srj8QdRyVgohI1C4ZM4DG5hYWr47+HgsqBRGRiI0fVEi/nt2YvzL6/QoqBRGRiGVmGB8dVcyit7bQ0BTt/cVUCiIiMTBlzAD2NDTx4tptkeZQKYiIxMDEYX3p0S0z8mshqRRERGIgNzuTC0b2Z8HKWlpaorvHgkpBRCQmLhldzNY9Dby+MbrbdKoURERi4sKR/cnOtEivhaRSEBGJifzcbCYO6xfpoakqBRGRGJkwpJD1W/eyvzGaQ1NVCiIiMVJW2B2ATTv2R7L80ErBzO41sy1mtjxprI+ZLTCzNcFzYdJnN5tZlZmtNrMpYeUSEYmz0sI8AKq374tk+WGuKfwWmHrI2GxgobuPABYG7zGz0cAMYEwwz91mlhliNhGRWCoLSqHTrSm4+9+B9w4Zng7MDV7PBS5NGp/n7g3uvh6oAiaElU1EJK7698olK8PYtL2TlUI7it29BiB47h+MlwIbk6arDsY+xMxmmVmlmVXW1UV/RUERkY6UmWGUFORS3UVKoT3Wxlibp/S5+xx3r3D3iqKiopBjiYikXllB9863+agdtWZWAhA8bwnGq4HypOnKgM0pziYiEgulhXldZvPR48DM4PVM4LGk8RlmlmNmQ4ERwCspziYiEgulBXnU7q6nsakl5csO85DUPwAvAqeYWbWZfRH4AXCxma0BLg7e4+4rgAeBlcBfgevcPdqLiouIRKSsMA93qNmZ+rWFrLC+2N2vbOejye1Mfztwe1h5RETSxfvnKuxncN8eKV12XHY0i4hIoKwgOKs5gv0KKgURkZgZ0DuXDIPqCI5AUimIiMRMt6wMivNzI7nUhUpBRCSGyiI6LFWlICISQ6UFeZGcwKZSEBGJodLCPGp21tPUnNpzFVQKIiIxVFbYneYWp3Z3Q0qXq1IQEYmh0oLgEtop3q+gUhARiaGobrajUhARiSGtKYiISKvc7Ez69cxJ+X0VVAoiIjFVWpj6w1JVCiIiMVWmUhARkYPKChJnNbe0tHkjylCoFEREYqqsMI/G5ha27knduQoqBRGRmGo9LDWFm5BUCiIiMVUa3FchlUcgqRRERGLq4JpCKs9VUCmIiMRUz5wsCrpns2lH6s5qVimIiMRYaUGeNh+JiEhCqm+2o1IQEYmx0oLubNqxH/fUnKugUhARibHSwjz2NTazfd+BlCxPpSAiEmNlKT4CKZJSMLMNZvammS01s8pgrI+ZLTCzNcFzYRTZRETi5OAltFN1X4Uo1xQudPdx7l4RvJ8NLHT3EcDC4L2ISJdWXpg4gS1VF8aL0+aj6cDc4PVc4NLoooiIxEN+XhY9c7JSdlhqVKXgwHwzW2Jms4KxYnevAQie+7c1o5nNMrNKM6usq6tLUVwRkWiYWUrPVchKyVI+bJK7bzaz/sACM3vraGd09znAHICKiorUXU9WRCQiqbyvQiRrCu6+OXjeAjwKTABqzawEIHjeEkU2EZG4KS3MY1Nn3dFsZj3MrNfB18AlwHLgcWBmMNlM4LFUZxMRiaPSgjx21Texqz78cxWi2HxUDDxqZgeXf7+7/9XMXgUeNLMvAu8AV0SQTUQkdsoOHoG0fT/5JdmhLivlpeDu64CxbYxvAyanOo+ISNwlX0J7VEl+qMuK0yGpIiLShlSewKZSEBGJuX49u5GTlZGSI5BUCiIiMWdmlBam5lwFlYKISBooK+yuNQUREUkoLUjNzXZUCiIiaaCsMI9texvZ19gU6nJUCiIiaeDgfRU2h7wJSaUgIpIG3j8sVaUgItLlHTyBTaUgIiL075VLdqaFfgSSSkFEJA1kZhglvcM/AkmlICKSJhI32wn3UhcqBRGRNJGKm+2oFERE0kRpYR61uxpoaGoObRkqBRGRNHHwsNSaHfWhLUOlICKSJlpvthPiJiSVgohImigrDP++CioFEZE0MaB3LhlGqIelqhRERNJEdmYGA/JzqdbmIxERAUK/2Y5KQUQkjZQVdtfmIxERSSgtyOPdXfU0NbeE8v0qBRGRNFJamEdzi/PurnDOVVApiIikkYOHpYa1CSl2pWBmU81stZlVmdnsqPOIiMRJ2DfbiVUpmFkmcBfwMWA0cKWZjY42lYhIfAwMSiGss5pjVQrABKDK3de5eyMwD5gecSYRkdjIzc6kqFdOl9l8VApsTHpfHYy1MrNZZlZpZpV1dXUpDSciEgcfP62Ekwf0CuW7s0L51uNnbYz5B964zwHmAFRUVHgb04uIdGrf+dSY0L47bmsK1UB50vsyYHNEWUREupy4lcKrwAgzG2pm3YAZwOMRZxIR6TJitfnI3ZvM7HrgaSATuNfdV0QcS0Sky4hVKQC4+5PAk1HnEBHpiuK2+UhERCKkUhARkVYqBRERaaVSEBGRVuaevud/mVkd8PYJfEU/YGsHxQlbOmWF9MqrrOFJp7zplBVOLO9gdy9q64O0LoUTZWaV7l4RdY6jkU5ZIb3yKmt40ilvOmWF8PJq85GIiLRSKYiISKuuXgpzog5wDNIpK6RXXmUNTzrlTaesEFLeLr1PQUREPqirrymIiEgSlYKIiLTqUqVgZl83MzezfkljN5tZlZmtNrMpSeNnmtmbwWf/aWZt3QAojIy3mdkbZrbUzOab2cC4Zg2W/R9m9laQ+VEzK4hrXjO7wsxWmFmLmVUc8lmssrbFzKYG+arMbHZUOZKZ2b1mtsXMlieN9TGzBWa2JnguTPqszZ9zirKWm9kiM1sV/B7cGNe8ZpZrZq+Y2bIg63dTltXdu8SDxM17niZxslu/YGw0sAzIAYYCa4HM4LNXgIkk7gb3FPCxFOXMT3p9A/DzuGYNln0JkBW8/iHww7jmBUYBpwCLgYqk8dhlbSN7ZpDrJKBbkHd0FFkOyfURYDywPGnsR8Ds4PXso/mdSFHWEmB88LoX8I8gU+zyBr9vPYPX2cDLwDmpyNqV1hTuBL7JB2/vOR2Y5+4N7r4eqAImmFkJif84v+iJn/jvgEtTEdLddyW97ZGUN3ZZg7zz3b0pePsSibvlxTKvu69y99VtfBS7rG2YAFS5+zp3bwTmBbkj5e5/B947ZHg6MDd4PZf3f2Zt/pxTkRPA3Wvc/bXg9W5gFYl7wMcuryfsCd5mBw9PRdYuUQpm9ilgk7svO+SjUmBj0vvqYKw0eH3oeEqY2e1mthH4HPCtYDiWWQ/xBRJ/TUN65D0oHbK2lzGOit29BhL/IQb6B+Ox+TeY2RDgDBJ/gccyr5llmtlSYAuwwN1TkjV2N9k5Xmb2DDCgjY9uBW4hsZnjQ7O1MeaHGe8Qh8vq7o+5+63ArWZ2M3A98O2ossKR8wbT3Ao0AfcdnK2dXJH+bNubrZ1Mof9sj0GcshyvWPwbzKwn8DBwk7vvOsxuokjzunszMC7YT/eomZ16mMk7LGunKQV3/2hb42Z2GoltbMuC//HLgNfMbAKJNi1PmrwM2ByMl7UxHmrWNtwP/IVEKUSSFY6c18xmAp8AJgebWSD+P9tkkf1sj0F7GeOo1sxK3L0m2AS3JRiP/N9gZtkkCuE+d38kGI5tXgB332Fmi4GppCJrKnaaxOkBbOD9Hc1j+ODOmXW8v4PxVRI7dg7uYJyWonwjkl7/G/BQXLMGy54KrASKDhmPZd5g+Yv54I7m2GZNypgV5BrK+zuax0SRpY1sQ/jgjub/4IM7Q390pJ9zinIaif1C/+eQ8djlBYqAguB1HvAsiT+8Qs8a+S9Uqh8klULw/lYSe+pXk3RkCVABLA8++xnB2d8pyPdwsNw3gCeA0rhmDZZdRWJb5tLg8fO45gUuI/EXVQNQCzwd16zt5J9G4oiZtSQ2h0WS45BMfwBqgAPBz/aLQF9gIbAmeO5zpJ9zirL+E4lNKm8k/b5Oi2Ne4HTg9SDrcuBbwXjoWXWZCxERadUljj4SEZGjo1IQEZFWKgUREWmlUhARkVYqBRERaaVSEOlAwZU415tZn+B9YfB+cNTZRI6GSkGkA7n7RuAe4AfB0A+AOe7+dnSpRI6ezlMQ6WDBpRSWAPcC1wBneOLKpiKx12mufSQSF+5+wMy+AfwVuESFIOlEm49EwvExEpd/ONyVLUViR6Ug0sHMbBxwMYkL6X0luJqlSFpQKYh0oOAezveQuFb/OySuavnjaFOJHD2VgkjHugZ4x90XBO/vBkaa2fkRZhI5ajr6SEREWmlNQUREWqkURESklUpBRERaqRRERKSVSkFERFqpFEREpJVKQUREWv1/AIxnhEOvgHYAAAAASUVORK5CYII=\n", "text/plain": [ "