528 lines
60 KiB
Plaintext
528 lines
60 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# 1 System Modelling\n",
|
|
"\n",
|
|
"This notebook contains the theory on using the vehicle Kinematics Equations to derive the linearized state space model"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {
|
|
"ExecuteTime": {
|
|
"end_time": "2024-10-22T08:05:56.118290Z",
|
|
"start_time": "2024-10-22T08:05:46.550696Z"
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Matplotlib is building the font cache; this may take a moment.\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"from scipy.integrate import odeint\n",
|
|
"from scipy.interpolate import interp1d\n",
|
|
"import cvxpy as cp\n",
|
|
"\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"\n",
|
|
"plt.style.use(\"ggplot\")\n",
|
|
"\n",
|
|
"import time"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## kinematics model equations\n",
|
|
"\n",
|
|
"The variables of the model are:\n",
|
|
"\n",
|
|
"* $x$ coordinate of the robot\n",
|
|
"* $y$ coordinate of the robot\n",
|
|
"* $v$ velocity of the robot\n",
|
|
"* $\\theta$ heading of the robot\n",
|
|
"\n",
|
|
"The inputs of the model are:\n",
|
|
"\n",
|
|
"* $a$ acceleration of the robot\n",
|
|
"* $\\delta$ steering of the robot\n",
|
|
"\n",
|
|
"These are the differential equations f(x,u) of the model:\n",
|
|
"\n",
|
|
"$\\dot{x} = f(x,u)$\n",
|
|
"\n",
|
|
"* $\\dot{x} = v\\cos{\\theta}$ \n",
|
|
"* $\\dot{y} = v\\sin{\\theta}$\n",
|
|
"* $\\dot{v} = a$\n",
|
|
"* $\\dot{\\theta} = \\frac{v\\tan{\\delta}}{L}$\n",
|
|
"\n",
|
|
"Discretize with forward Euler Integration for time step dt:\n",
|
|
"\n",
|
|
"${x_{t+1}} = x_{t} + f(x,u)dt$\n",
|
|
"\n",
|
|
"* ${x_{t+1}} = x_{t} + v_t\\cos{\\theta}dt$\n",
|
|
"* ${y_{t+1}} = y_{t} + v_t\\sin{\\theta}dt$\n",
|
|
"* ${v_{t+1}} = v_{t} + a_tdt$\n",
|
|
"* ${\\theta_{t+1}} = \\theta_{t} + \\frac{v\\tan{\\delta}}{L} dt$\n",
|
|
"\n",
|
|
"----------------------\n",
|
|
"\n",
|
|
"The Model is **non-linear** and **time variant**, but the Numerical Optimizer requires a Linear sets of equations. To approximate the equivalent **LTI** State space model, the **Taylor's series expansion** is used around $\\bar{x}$ and $\\bar{u}$ (at each time step):\n",
|
|
"\n",
|
|
"$ f(x,u) \\approx f(\\bar{x},\\bar{u}) + \\frac{\\partial f(x,u)}{\\partial x}|_{x=\\bar{x},u=\\bar{u}}(x-\\bar{x}) + \\frac{\\partial f(x,u)}{\\partial u}|_{x=\\bar{x},u=\\bar{u}}(u-\\bar{u})$\n",
|
|
"\n",
|
|
"This can be rewritten usibg the State Space model form Ax+Bu :\n",
|
|
"\n",
|
|
"$ f(\\bar{x},\\bar{u}) + A|_{x=\\bar{x},u=\\bar{u}}(x-\\bar{x}) + B|_{x=\\bar{x},u=\\bar{u}}(u-\\bar{u})$\n",
|
|
"\n",
|
|
"Where:\n",
|
|
"\n",
|
|
"$\n",
|
|
"A =\n",
|
|
"\\quad\n",
|
|
"\\begin{bmatrix}\n",
|
|
"\\frac{\\partial f(x,u)}{\\partial x} & \\frac{\\partial f(x,u)}{\\partial y} & \\frac{\\partial f(x,u)}{\\partial v} & \\frac{\\partial f(x,u)}{\\partial \\theta} \\\\\n",
|
|
"\\end{bmatrix}\n",
|
|
"\\quad\n",
|
|
"=\n",
|
|
"\\displaystyle \\left[\\begin{matrix}0 & 0 & \\cos{\\left(\\theta \\right)} & - v \\sin{\\left(\\theta \\right)}\\\\0 & 0 & \\sin{\\left(\\theta \\right)} & v \\cos{\\left(\\theta \\right)}\\\\0 & 0 & 0 & 0\\\\0 & 0 & \\frac{\\tan{\\left(\\delta \\right)}}{L} & 0\\end{matrix}\\right]\n",
|
|
"$\n",
|
|
"\n",
|
|
"and\n",
|
|
"\n",
|
|
"$\n",
|
|
"B = \n",
|
|
"\\quad\n",
|
|
"\\begin{bmatrix}\n",
|
|
"\\frac{\\partial f(x,u)}{\\partial a} & \\frac{\\partial f(x,u)}{\\partial \\delta} \\\\\n",
|
|
"\\end{bmatrix}\n",
|
|
"\\quad\n",
|
|
"= \n",
|
|
"\\displaystyle \\left[\\begin{matrix}0 & 0\\\\0 & 0\\\\1 & 0\\\\0 & \\frac{v \\left(\\tan^{2}{\\left(\\delta \\right)} + 1\\right)}{L}\\end{matrix}\\right]\n",
|
|
"$\n",
|
|
"\n",
|
|
"are the *Jacobians*.\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"So the discretized model is given by:\n",
|
|
"\n",
|
|
"$ x_{t+1} = x_t + (f(\\bar{x},\\bar{u}) + A|_{x=\\bar{x}}(x_t-\\bar{x}) + B|_{u=\\bar{u}}(u_t-\\bar{u}) )dt $\n",
|
|
"\n",
|
|
"$ x_{t+1} = (I+dtA)x_t + dtBu_t +dt(f(\\bar{x},\\bar{u}) - A\\bar{x} - B\\bar{u}))$\n",
|
|
"\n",
|
|
"The LTI-equivalent kinematics model is:\n",
|
|
"\n",
|
|
"$ x_{t+1} = A'x_t + B' u_t + C' $\n",
|
|
"\n",
|
|
"with:\n",
|
|
"\n",
|
|
"$ A' = I+dtA|_{x=\\bar{x},u=\\bar{u}} $\n",
|
|
"\n",
|
|
"$ B' = dtB|_{x=\\bar{x},u=\\bar{u}} $\n",
|
|
"\n",
|
|
"$ C' = dt(f(\\bar{x},\\bar{u}) - A|_{x=\\bar{x},u=\\bar{u}}\\bar{x} - B|_{x=\\bar{x},u=\\bar{u}}\\bar{u}) $"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"-----------------\n",
|
|
"[About Taylor Series Expansion](https://courses.engr.illinois.edu/ece486/fa2017/documents/lecture_notes/state_space_p2.pdf):\n",
|
|
"\n",
|
|
"In order to linearize general nonlinear systems, we will use the Taylor Series expansion of functions.\n",
|
|
"\n",
|
|
"Typically it is possible to assume that the system is operating about some nominal\n",
|
|
"state solution $\\bar{x}$ (possibly requires a nominal input $\\bar{u}$) called **equilibrium point**.\n",
|
|
"\n",
|
|
"Recall that the Taylor Series expansion of f(x) around the\n",
|
|
"point $\\bar{x}$ is given by:\n",
|
|
"\n",
|
|
"$f(x)=f(\\bar{x}) + \\frac{df(x)}{dx}|_{x=\\bar{x}}(x-\\bar{x})$ + higher order terms...\n",
|
|
"\n",
|
|
"For x sufficiently close to $\\bar{x}$, these higher order terms will be very close to zero, and so we can drop them.\n",
|
|
"\n",
|
|
"The extension to functions of multiple states and inputs is very similar to the above procedure.Suppose the evolution of state x\n",
|
|
"is given by:\n",
|
|
"\n",
|
|
"$\\dot{x} = f(x1, x2, . . . , xn, u1, u2, . . . , um) = Ax+Bu$\n",
|
|
"\n",
|
|
"Where:\n",
|
|
"\n",
|
|
"$ A =\n",
|
|
"\\quad\n",
|
|
"\\begin{bmatrix}\n",
|
|
"\\frac{\\partial f(x,u)}{\\partial x1} & ... & \\frac{\\partial f(x,u)}{\\partial xn} \\\\\n",
|
|
"\\end{bmatrix}\n",
|
|
"\\quad\n",
|
|
"$ and $ B = \\quad\n",
|
|
"\\begin{bmatrix}\n",
|
|
"\\frac{\\partial f(x,u)}{\\partial u1} & ... & \\frac{\\partial f(x,u)}{\\partial um} \\\\\n",
|
|
"\\end{bmatrix}\n",
|
|
"\\quad $\n",
|
|
"\n",
|
|
"Then:\n",
|
|
"\n",
|
|
"$f(x,u)=f(\\bar{x},\\bar{u}) + \\frac{df(x,u)}{dx}|_{x=\\bar{x}}(x-\\bar{x}) + \\frac{df(x,u)}{du}|_{u=\\bar{u}}(u-\\bar{u}) = f(\\bar{x},\\bar{u}) + A_{x=\\bar{x}}(x-\\bar{x}) + B_{u=\\bar{u}}(u-\\bar{u})$\n",
|
|
"\n",
|
|
"-----------------"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## ODE Model\n",
|
|
"Motion Prediction: using scipy intergration"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {
|
|
"ExecuteTime": {
|
|
"end_time": "2024-10-22T09:21:29.100064Z",
|
|
"start_time": "2024-10-22T09:21:29.094806Z"
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Define process model\n",
|
|
"# This uses the continuous model\n",
|
|
"def kinematics_model(x, t, u):\n",
|
|
" \"\"\"\n",
|
|
" Returns the set of ODE of the vehicle model.\n",
|
|
" \"\"\"\n",
|
|
"\n",
|
|
" L = 0.3 # vehicle wheelbase\n",
|
|
" dxdt = x[2] * np.cos(x[3])\n",
|
|
" dydt = x[2] * np.sin(x[3])\n",
|
|
" dvdt = u[0]\n",
|
|
" dthetadt = x[2] * np.tan(u[1]) / L\n",
|
|
"\n",
|
|
" dqdt = [dxdt, dydt, dvdt, dthetadt]\n",
|
|
"\n",
|
|
" return dqdt\n",
|
|
"\n",
|
|
"\n",
|
|
"def predict(x0, u):\n",
|
|
" \"\"\" \"\"\"\n",
|
|
"\n",
|
|
" x_ = np.zeros((N, T + 1))\n",
|
|
"\n",
|
|
" x_[:, 0] = x0\n",
|
|
"\n",
|
|
" # solve ODE\n",
|
|
" for t in range(1, T + 1):\n",
|
|
"\n",
|
|
" tspan = [0, DT]\n",
|
|
" x_next = odeint(kinematics_model, x0, tspan, args=(u[:, t - 1],))\n",
|
|
"\n",
|
|
" x0 = x_next[1]\n",
|
|
" x_[:, t] = x_next[1]\n",
|
|
"\n",
|
|
" return x_"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 9,
|
|
"metadata": {
|
|
"ExecuteTime": {
|
|
"end_time": "2024-10-22T09:25:39.603454Z",
|
|
"start_time": "2024-10-22T09:25:39.598889Z"
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"CPU times: user 2.66 ms, sys: 55 μs, total: 2.71 ms\n",
|
|
"Wall time: 2.77 ms\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"%%time\n",
|
|
"\n",
|
|
"u_bar = np.zeros((M, T))\n",
|
|
"u_bar[0, :] = 0.2 # m/ss\n",
|
|
"u_bar[1, :] = np.radians(-np.pi / 4) # rad\n",
|
|
"\n",
|
|
"x0 = np.zeros(N)\n",
|
|
"x0[0] = 0\n",
|
|
"x0[1] = 1\n",
|
|
"x0[2] = 0\n",
|
|
"x0[3] = np.radians(0)\n",
|
|
"\n",
|
|
"x_bar = predict(x0, u_bar)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"metadata": {
|
|
"ExecuteTime": {
|
|
"end_time": "2024-10-22T09:25:42.569050Z",
|
|
"start_time": "2024-10-22T09:25:42.507079Z"
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": "<Figure size 640x480 with 2 Axes>",
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAEDCAYAAABTZPIVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCvElEQVR4nO3deVxVdf7H8de5LIIim4Kg4EKIW7lWttioOWqZM0ZjZo2VZWVhTVPjtNFv0sZyrKZlst21aXTMhSy13KYabbOpzAqX3EJFBEIWQdbz/f1xk0IWQYHLvbyfjwcPuWe59/Ph4uF9z/mecyxjjEFERERE3JrD1QWIiIiIyJlTqBMRERHxAAp1IiIiIh5AoU5ERETEAyjUiYiIiHgAhToRERERD6BQJyIiIuIBFOpEREREPIBCnYiIiIgH8HZ1Aa5y9OhRSktLT7lcWFgYGRkZjVCRa6g/9+XJvUHt+/P29iYkJKQRKnI/2s55dm+g/txdfW/nmm2oKy0tpaSkpMZlLMsqX9YT76am/tyXJ/cGnt9fY2nu2zlP7g3Un7triP50+FVERETEAyjUiYiIiHgAhToRERERD6BQJyIiIuIBFOpEREREGpkpKcac4kSmulKoExEREWlEJiuDsln3k/3a0/X6vM32kiYiIiIijc3s/Ab7lScgL4eCHzOwho6GoPq51qZCnYiIiEgDM8ZgNr6DWToPbBuiY2g3/VkybKverlOnUCciIiLSgExREeafszGffQiANXAwjhvuxLtdezh8uN5eR6FOREREpIGYzCPYLz4OB/aBw4F19c1Yw36D5aj/0xoU6kREREQagEneiv3qk5CfB62DcEy+D6vbOQ32egp1IiIiIvXIGINZl4RZ/joYGzrF4rjjQaw2YQ36ugp1IiIiIvXEFBViFvwD87/NAFgXDcP6/e1Yvi0a/LXd+jp1SUlJjBs3jgULFri6FBEREWnmTPph7Jl/dgY6Ly+s627HmviHRgl04MZ76nbv3s2GDRvo1KmTq0sRERGRZs58+wX2a09BQT4EBuO4/QGsrj0btQa33FNXWFjI888/z+TJk2nVqpWryxEREZFmyhiDvfpN7H886gx0Md1wPPxMowc6cNM9dXPmzKFfv3707t2bFStW1LhsSUkJJb+4t5plWfj7+5d/X5MT80+1nLtSf+7Lk3sDz+9PRDyDKSzAnv8cfPkJANYlI7CunYzl4+OSetwu1H300Ufs27ePmTNn1mr5pKQkli1bVv64S5cuzJo1i7Cw2p+BEhERUec63Yn6c1+e3Bt4fn8i4r5M2kHsF2fC4QPg7Y117WQcvxrp0prcKtRlZmayYMECEhMT8fX1rdU68fHxjB49uvzxiU/+GRkZlJaW1riuZVlERESQlpZWb7fwaErUn/vy5N6gbv15e3vX6UOaiMiZMl9vwZ77NBwvgOBQ5/i5s7q7uiz3CnV79+4lJyeHBx54oHyabdts376d9957j0WLFuE46QrNPj4++FSzG7S2fwyNMR75h/ME9ee+PLk38Pz+RMS9GNvGrFqCeWexc0JsTxy3348VFOLawn7iVqHunHPO4amnnqow7aWXXqJ9+/aMGTOmUqATERERqQ+mIB973jPw9RYArKFXYI27GcvbNePnquJWoc7f35+OHTtWmNaiRQtat25dabqIiIhIfTCpKc7xc0cOgbcP1oQEHBcPc3VZlbhVqBMRERFpTObLT7DnPQtFxyG0rfN2X527urqsKrl9qJs2bZqrSxAREREPY+wyzMpFmDVLnRO6nYNj8n1YrYNcW1gN3D7UiYiIiNQnk38Me85T8O2XAFi/HoM1diKWl5eLK6uZQp2ISCMpKytj6dKlbNq0iezsbEJCQhgyZAhXXXVV+YlexhiWLl3Kxo0bOXbsGF27dmXSpElER0e7uHqR5sEc3I/94uOQkQa+vlg33IVj4GBXl1UrCnUiIo1k5cqVrF+/nilTphAVFcXevXt58cUXadmyJaNGjSpfZvXq1SQkJBAZGcmKFSuYMWMGzz77bPndcESkYdifb8YseA6Ki6BNOI6Eh7A6xri6rFrTNUBERBrJrl27OPfcc+nfvz/h4eFccMEF9O7dmz179gDOvXRr1qwhPj6egQMH0rFjR6ZMmUJRURGbN292cfUinsuUlWEvW4B59QlnoOvZF8fDT7tVoAPtqRMRaTTdu3dn/fr1pKam0r59e/bv38/OnTu58cYbAUhPTyc7O5s+ffqUr+Pj40PPnj3ZuXMnw4cPr/J5dY/rqnlyb6D+6os5lot55UnM9q3O17vsdziuuh7L0bDj5xqiP4U6EZFGMmbMGAoKCrjnnntwOBzYts348eMZNGgQANnZ2QAEBVU8uy4oKIjMzMxqn1f3uK6ZJ/cG6u9MFO/ZSebMP2OOpGK18CP0nkdoeUnVH54aSn32p1AnItJIPv74YzZt2sQf/vAHoqOj2b9/PwsWLCg/YeKEkz+5n+pWabrHddU8uTdQf2fK/vQD7Nefh+JiCIvEMeUhcqI6k3P4cL2/VlUa4h7XCnUiIo3kjTfeYMyYMVx88cUAdOzYkYyMDN566y2GDBlCcHAwQPmZsSfk5uZW2nv3S7rHdc08uTdQf3V+vtJSzPIFmA1vOyecPQDHLX+CVgEu+TnWZ386UUJEpJEUFRVVuke1w+Eo36CHh4cTHBzMtm3byueXlpaSnJxMt27dGrVWEU9kcrOxn32kPNBZV4zDcdfDWK0CXFxZ/dCeOhGRRjJgwABWrFhB27ZtiYqKYv/+/axatYqhQ4cCzsMxo0aNIikpicjISCIiIkhKSqJFixbl4+5E5PSY/d9jvzQTsjKhhT+Om/+I1f9CV5dVrxTqREQayc0338ySJUuYM2cOOTk5hIaGMnz4cMaOHVu+zJgxYyguLmbOnDnk5+cTGxtLYmKirlEncgbsjzZi3ngRSksgooPz+nORnndBb4U6EZFG4u/vz8SJE5k4cWK1y1iWxbhx4xg3blzjFSbioUxpCebNuZj31zgn9Dkfx833YLVs5drCGohCnYiIiHgck3MU++VZsDsZAOu312FdMQ7L4bmnEyjUiYiIiEcxe3Zgv/w3yM4C/5Y4Jv0Jq895ri6rwSnUiYiIiMew//seZtGrUFYKkdHO8XMRHVxdVqNQqBMRERG3Z0pKMItfwWxa55zQ/0IcN92N5dfStYU1IoU6ERERcWvm6I/Oy5Xs2wWWhXXlBKzLx3rsfXGro1AnInKSG2+88bTWmz59Op07d67fYkSkRub7ZOf4udxsaNkKx61Tsc4e4OqyXEKhTkTkJIWFhfTr14/AwMBaLW/bNps2bcK27QauTEROMMZgPliDWTIHysqgQyfn+LnwSFeX5jIKdSIiVRg7diyxsbG1WrasrIxNmzY1cEUicoIpLsK88RLmk/8AYJ13CdaNd2G18HNxZa6lUCcicpLLL7+c4ODgWi/vcDjqvI6InB7zY4Zz/NwPu8FyYP3uRqwRVza78XNVUagTETlJTXd8qIplWXVeR0TqzuzYhv3KE3AsFwJa47j1z1g9+7q6rCZDoU5ERESaNGMMZuPbmKXzwbahYwyOOx7EatvO1aU1KQp1IiI1SE5Ornaew+GgZcuWdOjQAS8vr0asSqT5MEWFmIWzMVs+BMC6YAjW9VOwfFu4uLKmR6FORKQG06dPP+Uyfn5+jB49mquvvroRKhJpPkrTDlH2t/vhwF5wOLDGTcK6dLTGz1VDoU5EpAb3338/8+bNo3379lx88cUEBQWRnZ3NRx99RGpqKtdccw3bt29n+fLlBAQEcPnll7u6ZBGPYCdv5chrT0FeDrQOwjH5fqxuZ7u6rCZNoU5EpAZbt26lR48eTJkypcL0IUOGMHv2bHbs2MHkyZMB2Lhxo0KdyBkyxmDWJWGWvw7Ghs6xzvFzoWGuLq3Jc7i6ABGRpuzjjz/m4osvrnLeoEGD+OyzzwAYMGAAhw8fbszSRDyOKSrEvPokZtkCMDathv8Gr/tnKdDVkvbUiYjUoKioiNzc3Crn5eTkUFxcDDjH1elkCZHTZ9JTsV+cCYd+AC8vHONvI+Tam0lLS8MY4+ry3IJCnYhIDbp168aSJUuIjY2lffv25dMPHTrEm2++Sffu3QFIT0+nTZs2ripTxK2Zb77AnvMUFORDUAiOyffjiOulEyLqSKFORKQGEydO5JFHHuHee+8lOjqa4OBgsrOzOXDgAK1ateLGG28EICsri8GDB7u4WhH3YozBrFmKWfkvMAbO6o7j9vuxgvUB6XQo1ImI1CAqKoq///3vrFq1iu3bt3PkyBFat27N6NGjueKKKwgJCQGc94oVkdozhQXY856Frz4FwPrVZVjjb8Xy8XFtYW5MoU5E5BSCg4OZMGGCq8sQ8Rgm7aBz/NzhA+DtjXXd7TguGeHqstyeQp2ISC0UFBSwa9cu8vLy6NevHwEBAa4uScQtma2fYc97Bo4XQHAbHHc8gBXTzdVleQS3CnVJSUls2bKFQ4cO4evrS1xcHBMmTKgweFlEpL4tW7aMlStXlp/pOnPmTAICAnj00Ufp3bs3V155pWsLFHEDxrYxq/6NeeffzgmxPZ3j54JCXFuYB3Gr69QlJyczcuRIHnvsMR5++GFs22bGjBkUFha6ujQR8VBr165l2bJlDB06lAceeKDCvP79+/Pll1+6qDIR92EK8rFfeKw80FlDr8Dxp78q0NUzt9pTl5iYWOFxQkICt9xyC3v37qVnz54uqkpEPNl7773H6NGjmTBhArZtV5gXGRmpCw6LnIJJTcF+4XFITwVvH6wJCTguHubqsjySW4W6kxUUFADUOLalpKSEkpKS8seWZeHv71/+fU1OzPfU6+SoP/flyb1B0+ovPT2dPn36VDnP39+/fDskIpWZLz/GnvccFB2H0LbO23117urqsjyW24Y6YwwLFy6ke/fudOzYsdrlkpKSWLZsWfnjLl26MGvWLMLCan/LkYiIiDOqtalTf+7Lk3uDptFfy5YtycnJqXJeeno6gYGBjVyRSNNn7DLMykWYNUudE7qdg2PyfVitg1xbmIdz21A3d+5cUlJSePTRR2tcLj4+ntGjR5c/PvHJPyMjg9LS0hrXtSyLiIgIj71FifpzX57cG9StP29v7zp9SKurs88+m5UrV3Luuefi6+tbXl9ZWRnr16+vdi+eSHNl8o857w7xrXO8qfXrMVhjJ2LpNnoNzi1D3bx58/jiiy+YPn36KW/L4+Pjg081FzKs7R9DY4xH/uE8Qf25L0/uDZpGf9dccw0PPvgg9957L+effz7gHGe3f/9+MjMzueeee+r0fFlZWbzxxhts3bqV4uJiIiMjueOOO4iJiQGcPS9dupSNGzdy7NgxunbtyqRJk4iOjq733kTqmzm4z3n9uYw08PXFuuEuHAN1p5XG4lahzhjDvHnz2LJlC9OmTSM8PNzVJYmIh4uIiOCvf/0rCxcuZO3atQD897//pVevXtx11120bdu21s917Ngx/u///o9evXrx0EMPERgYyJEjR2jZsmX5MitXrmT16tUkJCQQGRnJihUrmDFjBs8++2z5eGCRpsj+fBNmwT+guAjahONIeAirY4yry2pW3CrUzZ07l82bN3Pffffh7+9PdnY24BzzcuKwiIhIfYuKiiIxMZGSkhLy8vIICAg4rW3OypUradOmDQkJCeXTfvnh1BjDmjVriI+PZ+DAgQBMmTKFW2+9lc2bNzN8+PAzb0aknpmyMkzS65i1Sc4JPfviuHUqVoDGmzY2twp169atA2DatGkVpickJDBkyJDGL0hEmhUfHx9CQ0NPe/3//e9/9OnTh6effprk5GRCQ0MZMWIEv/71rwHniRfZ2dkVxun5+PjQs2dPdu7cWW2o01n+VfPk3qBp9GfycjCvPoXZvtVZy2W/w3HV9ViOMx8/1xT6a0gN0Z9bhbo333zT1SWISDPwyzPma2Ps2LG1Wi49PZ3169dzxRVXEB8fz+7du5k/fz4+Pj4MHjy4/OhDUFDFMwSDgoLIzMys9nl1ln/NPLk3cF1/xXt2kDnzz5j0w1h+/oT+8S+0vKT+9ybr/as9twp1IiKNYenSpXVavrahzrZtzjrrLK677jrAGb4OHDjAunXrGDz458HkJ39yP9XJIjrLv2qe3Bu4tj/7k/exX58NJcUQFonjzkRyOnQipx4vxq3372e1PctfoU5E5CRLliwp//7w4cM8/vjjXHrppQwaNIjg4GCys7PZtGkT77//Pg899FCtnzckJISoqKgK06Kiovjss88ACA4OBiA7O5uQkJ9vn5Sbm1tp790v6Sz/mnlyb9C4/ZnSUsyy+ZiN7zgnnD0Axy1/glYBDVaD3r/ac6t7v4qINLYFCxYwePBg4uPjCQsLw8fHh7CwMK666ip+9atfMX/+/Fo/V7du3UhNTa0wLTU1tfwTeHh4OMHBwWzbtq18fmlpKcnJyXTr1q1+GhI5TSY3G/uZv5QHOuuKcTjuehirVfV3dZLGpVAnIlKD7du3VxuounXrxo4dO2r9XFdccQXff/89K1asIC0tjc2bN7Nx40ZGjhwJOA/HjBo1iqSkJLZs2UJKSgovvPACLVq0YNCgQfXSj8jpMPu+x55xL+z6Fvz8cSQ8hOPKCfVyQoTUHx1+FRGpgY+PD3v27OGcc86pNG/v3r14e9d+MxobG8vUqVNZtGgRy5cvJzw8nBtvvJFLLrmkfJkxY8ZQXFzMnDlzyM/PJzY2lsTERF2jTlzG/mgD5o2XoLQEIjo4rz8XqYthN0UKdSIiNTjvvPNYtmwZfn5+DBo0iICAAI4dO8bmzZtZtmxZnfegDRgwgAEDBlQ737Isxo0bx7hx4860dJEzYkpLMEvmYj5Y45zQ53wck+7F8m9Z84riMgp1IiI1uPHGGzly5Ajz589n/vz5eHl5UVZWBkCPHj248cYbXVyhSP0zOUexX/4b7N4OloX122uxRo3DcmjUVlOmUCciUgN/f38eeeQRtm7dynfffUdeXh6tW7emV69e9OnTx2MvjCrNl9mzA/ulv0FOFvi3wnHLvVi9z3N1WVILCnUiIrXQt29f+vbt6+oyRBqU/d/3MItehbJSiIzGMSURq117V5cltaRQJyIi0syZkhLM4lcwm5y346T/RThu+gOWn8bPuRMdHBcROcnUqVNJSUmp9fK2bTN16lQOHjzYgFWJNAxz9EfsJx90BjrLwrrqBhy3369A54a0p05E5CQHDhyguLi4wdcRcTWz6zvsV2ZBbja0DMBx61Sss/u7uiw5TQp1IiJVePLJJ6u99ZaIuzPGYD5Yg1kyB8rKIKqz8/pzYfV3c3lpfAp11TClpZCZRklpESYjHWMAywKH46d/vcDL8dO/3uDlBd7e4OWts+FE3NzgwYNPa73AwMB6rkSk/pniIswbL2E++Q8A1nmXYN14F1YLPxdXJmdKoa46eTmU/V8Caaezrrc3ePuAj6/zy9cXfFv84ssPy88PWvhBC3/w8wd/f/BrieXfCvxbQstWP30FQMtWuhWLSCNKSEhwdQkiDcL8mIH90kz4YTdYDqyxE7GGj9HOCA+hUFcdy4JWrXE4HNi2DcYGY8A2zu9t27nL2tiV1y0tdX4VHq/26U0dp+PfCloFQKvW0DoQKyAQAgKhdRC0DsIKDHZ+HxQCgcFYPr51bFhERDyZ2bEN+5Un4FguBLTGcdt9WD36uLosqUcKddWwgkPxfm4RkZGRHD58GGOqjlvGGGe4KyuDshLnvfFKS6Hkp++Li6GkyPlvcSGmuAiKiqDoOBQW/vTvcThegCk8Dsfz4XgBFOQ7v4p+CobH851fmUecr3tyHScX5t/KGfCCQrCCQ+GnLyukLZz4Cg6p15+ZiIg0PcYYzIa3McvmO3dIdIxxjp9rE+7q0qSeKdSdIcuyfjrc6g20OPXydXx+U1rqDHP5eXAsD/LzMMdynd/nZUNeLiYvx3nmUl4O5B51hsoTITDtYIXAVyH8ORyktgmnLLgNVmhbaBMObcKx2oZD23YQGoblrYHiIiLuyhQVYV6fjdnyIQDWBUOxrk/A8j313ytxPwp1TZzl7V1+iLV8Wg3LG2OcYS7nKOQcxWRnOW/1kp0FR3/EZP8IRzOdj22bsow0yEirOvhZDggJhbBIrPBICIv46d9ICI/E8vNvgI5FRKQ+mIw07BdnwsF94HBgXX0z1rDfaPycB1Oo8zCWZf10ckUAREZXGwCNXYaVm0MbyyZz1w7Mj0fgxwzMj+nOQ7w/HnEeMs7KhKxMzM5vnOv98kmCQ6FdB6x2HaBde6zIKIiIcu7t002fRURcxiR/hf3qU86jPK2DcEy+H6vb2a4uSxqYQl0zZTm8sELa0CIyEkdwWKUxg8YY5yHdjDRMZhqkH3Z+n34Y0lOdh3+znXsAKwU+X19n2IuMhvYdsTp0hPadoG07hT0RkQZkjMGsXYFZ8U/niXydu+K44wGs0DBXlyaNQKFOqmRZ1s8nWsT2qDTf5OfBkVTMkVQ4cgjSDmHSDjq/Ly6GA/swB/Y5lz2xkq8vRHbEat8RortgRXV2/huga3tJ05eamkpWVhbFxcUEBgbSvn17WrbUbZSk6TCFxzEL/oH54iMArIuHYf3+Dl0NoRlRqJPTYrVqDTHdsGK6VZhuysqch28PH8AcPgipKZjUHyD1gDPs/bAb88Nu+OQXYS+4jTPcRcdgdYqBjmc5D+Fq3Ie42K5du1i/fj1bt24lNze3wjyHw0Hnzp255JJLGDJkiAKeuJRJT8V+4XFITXFeBH/8LViDL9d2tJk5rVC3d+9eYmJi6rsW8QCWlxe0a+8cY9d3YPl0U1YGGWmQ+gPm4H7Mgf3OwbuZRyD7R8j+EfPN/34Oei0DoNNZWB1jsDp3hc5dFfSk0ezfv58FCxawfft2OnTowMCBA4mJiSEwMBBfX1+OHTvGkSNH+P7771m8eDFLliwhPj6e0aNH4+2tz8rSuOxt/8Oe85TzMlhBIThuvx8rtqeryxIXOK2tz4MPPkhsbCyXXXYZF154oTZickqWlxdEdICIDlj9Lyqfbo4XwMH9mIP7IGUvJmUPHEqBgmOw/WvM9q9/DnoBgdAlDqtzLFaXOOf3OnQrDeChhx5i0KBB3HDDDaf8AFtYWMjHH3/MypUrKSsr43e/+10jVSnNnbFtcv89F/uNl50Xxz+ruzPQBbdxdWniIqeVxhISEli7di2zZ8/m9ddfZ9iwYQwfPpw2bfSLJHVj+beErj2xuv78qdKUlDgP26bscR6u3b8bDu53XgX9m/9V3KMX3h7rrG4/HQruDh06OQOkyBn4+9//TmRkZK2W9fPz49JLL2XIkCFkZmY2cGUiTuZ4AWb+c+R89QkA1q8uwxp/K5aPri3anJ1WqBs8eDCDBw9m9+7dvPfee7zzzjusXLmSAQMGcPnll9OrV6/6rlOaEcvHx3notdNZcMkIAEzJTydf7P8e9n+P2bcL0g5BeiomPRU+ed8Z9Fr4OQNebA+srr2wg4NqfC2RqtQ20P2Sw+EgPFxX6JeGZ9IOOsfPpR0Ebx8c103G+mlbKc3bGR03jY2N5c477+SGG25gw4YNbNiwgUcffZSoqCguu+wyBg8ejK+vzrqRM2f5+FY6McPk58HeXZi9OzF7d8C+Xc5brP3isO2hZ72cJ2HE9sCKOxvieumQrdTJnXfeydSpU+ncuXOleSkpKTzxxBPMnj278QuTZsls/Qx77tPO20uGtCH84b+TFdSm2ltZSvNSL4PhvL29adGiRfnYuqKiIubMmcOKFSu45557iIuLq4+XEanAatUazhmAdc4AwDm+hNQUzO7tsDsZ830yZGWUn3FrNr7jXDGqM1bc2c4LcXY9G6u1Qp5ULyMjg9LS0irnlZSUkJGR0cgVSXNkbBvzzr8xq/7tnBDbE687HqBF955w+LBri5Mm44xC3Q8//MDatWvZvHkzpaWlXHDBBfzhD38gNjaWH374gVdffZXXXnuNJ598sr7qFamW5XA4A1tUZxjiPJU/zMviyMcfYL7/DrPzWzh84KcTM/Zj/rPKuWJUZ6zufbB69HbuyfPTpSmkdo4cOYK/v26XJw3LFBzDnvsMbPscAGvoFVjjbtb156SS0wp1H3/8MWvXrmXHjh0EBgYyevRoRowYQXBwcPkynTp14tprr+Wxxx6rr1pF6sw7PALHwMGY838FgMnNhu+/w+z8xhnyUlN+DnkbVoKXF3TuitWjrzPkxXR33n9XmpUPPviADz/8sPzxnDlzKoW34uJifvjhB3r21KUjpOGY1BTn+Ln0VPD2wbo+AcdFw1xdljRRp/XX6rnnnqNz587ccccdDBo0qNpLmoSFhXHJJZecUYEi9ckKDIYBF2MNuBgAk5eD2fEN7HCOwyMjDfbswOzZ4TzM4ecP3Xtj9eqH1as/VliEaxuQRlFcXFzhYsP5+fmUlJRUWMbHx4eLLrqIcePGNXZ50kyYLz7Gnv8cFB2H0LY4Eh7C6hTr6rKkCTutUDd9+nS6d+9+yuXatWtHQkLC6byESKOwWgdhnTcIzhsEgMk84gx3P51swbFc2PoZZutnzrNrwyN/CngDnGGvRQuX1i8NY8SIEYwY4TybcMqUKfzpT3+q8kQJkYZg7DLMW//CvLvMOaHbOTgm34fVWmfzS81OK9TVJtA1pLVr1/L222+TnZ1NVFQUEydOpEePyvcnFakrq20756UBLhnhPPHiwF7Mt19ikr+CPTsg/TAm/TDm/TXg4+sMdueci9X7XKw2upyFJ3rhhRca7LmTkpJYvHgxo0aNYuLEiYDzhuxLly5l48aNHDt2jK5duzJp0iSio6MbrA5pOkx+HvZrT8F3XwFgDR+D9buJuv6m1IrbDRb6+OOPWbBgAbfccgvdunVjw4YNPP744zzzzDO0bdvW1eWJB7EcDugU6zzcccU4590vdm5zhrxvvnCeWXviYsiLgPYdywMeXXvplmZuLDMz87S2J1lZWYSGhtZq2d27d7NhwwY6depUYfrKlStZvXo1CQkJREZGsmLFCmbMmMGzzz6rkzI8nDm4D/vFmc5hIL6+WDfchWPgYFeXJW7E7ULdqlWruPTSSxk2zDlQdOLEiXz99desW7eO6667rt5exxg4ftwiPx8KCiw88RJAlqX+6qYVxF0IcRdi4g0cPuAMeN99AXt3QUo6pKyBz7fi9chz9fGC1Wou752rerv77rv59a9/zeWXX05ERM3jKEtLS/n8889ZsWIFAwcOZOzYsad8/sLCQp5//nkmT57MihUryqcbY1izZg3x8fEMHOi8d/KUKVO49dZb2bx5M8OHDz+zxqTJsrf8F7PweSgugrbtnOPnoru4uixxM24V6kpLS9m7dy9XXnllhem9e/dm586dVa5TUlJSYYCzZVnln3Zr2pNy/LhFbOyJjbmnD45Xf6enPTCw6llvNNBLVuLZ792ePQ78/e1Gf92HH36YhQsX8t577xEbG0uvXr3o0qULQUFB+Pj4cOzYMY4cOcKuXbv4+uuvKSwsZNSoUYwePbpWzz9nzhz69etH7969K4S69PR0srOz6dOnT/k0Hx8fevbsyc6dO6sNdae7nfvlfE/cs+wOvZmyMuwVCzFrkwCwevbFcdufa3WRdHfo70yov7pzq1CXm5uLbdsEBVUcLBoUFER2dnaV6yQlJbFs2bLyx126dGHWrFmEhYXV+Fr5+WdcroicoXbt2tGqVeO/bo8ePfjb3/7GV199xfr163n33XcpLi6utFx4eDgjR45k+PDhhISE1Oq5P/roI/bt28fMmTMrzTuxHatqG1fTfWVPdzv3S6faI+nOmmpvZTnZ/DjrIYq+3gJA67E3EnRDQp3HzzXV/uqL+qs9twp1J1SVaqtLuvHx8RU+PZ9YrqarxIPzsM+ePQ7atWvHkSNHPPIWLJZlqT835cm9wc/95eamkZtbc3/e3t51Ci910a9fP/r160dpaSn79+/n6NGjFBcX07p1a6Kiomo9fu6EzMxMFixYQGJiYo23UDx5e3aq9/h0t3Mnlo2IiCAtLc3jfpeacm8mZQ9lLzwOP6ZDCz8cN93N8XMHcTw9vdbP0ZT7qw/q72e13c65VagLDAzE4XBU2iuXk5NT6ZPtCT4+Pvj4+FQ571Q/RH9/m1atnP966i+U+nNPntwb/Nxfbq5pEv15e3sTG3vm1wfbu3cvOTk5PPDAA+XTbNtm+/btvPfeezz77LOAc4/dL/f85ebmVruNgzPbzv1yuabws24ITa03+5P3Mf98AUqKITzSOX6uQ6fTrrGp9Vff1F/tuVWo8/b2JiYmhm3btnH++eeXT9+2bRvnnXeeCysTETm1c845h6eeeqrCtJdeeon27dszZswY2rVrR3BwMNu2baNLF+cg+dLSUpKTk/n973/vipKlHpnSUsyy+T/fh/qcc3Hcci9WywDXFiYew61CHcDo0aN5/vnniYmJIS4ujg0bNpCZmamzwkSkwRw+fJj169dz6NChSmPrLMviL3/5S62ex9/fn44dO1aY1qJFC1q3bl0+fdSoUSQlJREZGUlERARJSUm0aNGCQYMG1U8z4hImNxv7lSdg17cAWKOvwfrNtc5LJ4nUE7cLdRdddBF5eXksX76co0ePEh0dzYMPPthgY2pEpHlLSUkhMTGR0NBQ0tLS6NSpE3l5eWRlZdGmTRvatWtXr683ZswYiouLmTNnDvn5+cTGxpKYmKhr1Lkxs+977JdmwtFM8PPHcfM9WP0ucHVZ4oHcLtQBjBw5kpEjR7q6DBFpBhYvXkyfPn245557uO6667j99tuJiYnhyy+/5KWXXmL8+PFn9PzTpk2r8NiyLMaNG6d7ynoIe/N6zL9ehtISiOiAIyERKzLK1WWJh9J+XxGRGuzbt48hQ4aUn1F6YkBz//79+c1vfsOiRYtcWZ40Uaa0BPtfLzkvKFxaAn0H4njo7wp00qDcck+diEhjyc/PJyAgAIfDgZeXF/m/uIhlTExMhevDiQCY7Czsl//mvF+0ZWH99lqsUeM0fk4anEKdiEgNQkNDyc3NBZwXCU1OTqZ3796Ac7ydn5+fK8uTJsbs2YH90t8gJwv8WznPbu2tqzNI41CoExGpQbdu3di1axfnn38+gwYNYunSpWRnZ+Pt7c0HH3zAJZdc4uoSpQkwxmD+uxaz+FUoK4XIaBxTErHatXd1adKMKNSJiNTgqquu4ujRowBceeWVZGdns3nzZizL4sILL+T66693cYXiaqakGLPoFczm9c4JAy7CMfFuLD+dsSyNS6FORKQGERER5fdmdDgc3Hzzzdx8880urkqaCpOV6Rw/t2+Xc/xc/A1Yl13lsTehl6ZNozZFRGrw4osvkl7N/TgzMjJ48cUXG7kiaSrMrm+xZ9zjDHQtA3D84REcl/9OgU5cRqFORKQGH374YfmJEifLy8vjww8/bOSKxNWMMdj/WYX99P9BXg5Edcbx8NNYZ/d3dWnSzOnwq4jIaTp27Bg+Pj6uLkMakSkuwrzxIuaT9wGwzv8V1g13YrXQWdDiegp1IiInSU5OJjk5ufzxxo0b2bp1a4VliouL+fzzz4mK0sVkmwvzYzr2izMhZQ9YDqyxE7GGj9HhVmkyFOpERE7y3XffVbio8H/+858ql2vbti2TJk1qrLLEhcyObdivPAHHciEgEMdtf8bq0cfVZYlUoFAnInKSMWPGcNlll2GM4dZbbyUxMZEuXbpUWMbHx0cXHm4GjDGY9SsxyxeAbUPHs3AkPIjVJtzVpYlUolAnInISX19ffH19AZg9ezYhISF4e2tz2dyYoiLM689jtvwXAOvCoVgTErB8W7i4MpGqaSslIlKDsLAwAA4dOkRycjJ5eXlceumlBAcHk5WVRUBAQHkAFM9hMtKc4+cO7gOHA2vcLViXXqHxc9KkKdSJiNTAtm1eeeUVPvjgg/Jpffv2JTg4mFdffZUuXbpwzTXXuK5AqXfmu6+wX3sK8vOgdRCO2+/Hijvb1WWJnJKuUyciUoMVK1awefNmrr/+ev7+979XmNevX79KZ8WK+zLGYL+7HPu56c5A1yUOx8PPKNCJ29CeOhGRGnzwwQf87ne/Y/To0di2XWFeeHh4tXebEPdiCo9jL3gOvvgYAOviX2P9/nYsHx1aF/ehUCciUoOsrCzi4uKqnOfj40NhYWEjVyT1zaSnYr/wOKSmgJc31vhbsQZfpvFz4nYU6kREahAUFFTt3rjU1FRCQ0MbuSKpT+ab/2G/9nc4ng9BIThufwArtoeryxI5LQp1IiI16NevHytWrCg/OQLAsiwKCgp49913GTBggGsLlNNibBuzZinm7UVgDJzV3XlCRHAbV5cmctoU6kREajBu3Di++uor7rnnHnr16gXA4sWLOXDgAF5eXowdO9bFFUpdmeMF2POeha2fAjgPtY6/Fctb9/EV96ZQJyJSg+DgYGbOnMmbb77JV199hcPh4IcffqB///5cc801BAQEuLpEqQOTdtA5fi7tIHh7Y113O45LRri6LJF6oVAnInIKwcHB3Hbbba4uQ86Q2fop9txnoPA4BLdx3u6rS9UnwYi4I4U6ERHxaMa2sd9ejFn1b+eErj2d4+cCQ1xbmEg9U6gTETmFHTt2sHnzZjIyMiguLq4wz7Is/vKXv7ioMjkV+1ge9uwZmG2fA2AN+w3W2JuwdC9f8UD6rRYRqcH777/Pyy+/TEBAAJGRkfj4VBxMb4xxUWVyKubQDxx55QlMagr4+GJNSMBx0aWuLkukwSjUiYjU4O233+bCCy9kypQplQKdNF3mi4+x5z8LRYUQGuYcP9cp1tVliTQohToRkRpkZGRw0003KdC5CWOXYd76F+bdZQC06H0upTf9EQICXVuYSCNQqBMRqUGHDh3Iycmpl+dKSkpiy5YtHDp0CF9fX+Li4pgwYQLt27cvX8YYw9KlS9m4cSPHjh2ja9euTJo0iejo6HqpwZOZ/Dzs156C774CwBpxJWF3PkBaeoYOk0uz4HB1ASIiTdm1117LW2+9RVZW1hk/V3JyMiNHjuSxxx7j4YcfxrZtZsyYUeH+sStXrmT16tXcfPPNzJw5k+DgYGbMmMHx48fP+PU9mTmwD3vGvc5A5+uLdcuf8Bo3CctL+y6k+dBvu4jISWbNmlXhcUFBAXfffTedO3eudLFhy7K47777avW8iYmJFR4nJCRwyy23sHfvXnr27IkxhjVr1hAfH8/AgQMBmDJlCrfeeiubN29m+PDhZ9CV57K3/Bez8B9QXAxt2+FIeAgruouryxJpdAp1IiInSUlJqfDY4XAQGBhIVlZWveyxO6GgoACgPCimp6eTnZ1Nnz59ypfx8fGhZ8+e7Ny5s9pQV1JSQklJSfljy7Lw9/cv/74mJ+afarmmyJSVYS9fgFn3FgBWr344bv0zVkBr52M37q021J97a4j+FOpERE7ywgsvNPhrGGNYuHAh3bt3p2PHjgBkZ2cDEBQUVGHZoKAgMjMzq32upKQkli1bVv64S5cuzJo1i7CwsFrXExERUYfqXa8sJ5sfZz1I0dfO68+1vnoiQdffgeXlVWlZd+utrtSfe6vP/hTqRERqkJycTExMDH5+fpXmFRYWlh86rau5c+eSkpLCo48+WmneyZ/cTzXIPz4+ntGjR1daPyMjg9LS0hrXtSyLiIgI0tLS3OZkAvPDbspeeByyMqCFH46b7ub4uYM4np5eYTl37K0u1J97q0t/3t7etfqQ5jahLj09neXLl/Ptt9+SnZ1NaGgol1xyCVdddRXeujK4iDSQ6dOn89hjjxEbW/kaZ6mpqUyfPp0lS5bU6TnnzZvHF198wfTp02nTpk359ODgYMC5xy4k5OdbWOXm5lbae/dLPj4+1V5ypbZ/DI0xbvGH0/7kfcw/X4CSYgiPdI6f69CpxtrdpbfTpf7cW3325zZpKDU1FWMMt912GxERERw4cIBXXnmFwsJCbrjhBleXJyLNUGlpKQ5H7S8iYIxh3rx5bNmyhWnTphEeHl5hfnh4OMHBwWzbto0uXbqUv0ZycjK///3v67V2d2NKSzHL5mM2vuOccM65OG65F6tlQM0rijQjbhPq+vbtS9++fcsft2vXjtTUVNatW6dQJyL1qqCgoPwkBnDuOTt5TFtxcTEffvhh+d612pg7dy6bN2/mvvvuw9/fv3wMXcuWLfH19cWyLEaNGkVSUhKRkZFERESQlJREixYtGDRoUH205pZM7lHsV56AXd8BYI2+Bus312LVIVCLNAduE+qqUlBQUOnyAidrrmeF1Yb6c1+e3Bu4vr/Vq1dXOPHgySefrHbZ+Pj4Wj/vunXrAJg2bVqF6QkJCQwZMgSAMWPGUFxczJw5c8jPzyc2NpbExMTy7VZzY/btwn5xJmT/CH7+OG6+B6vfBa4uS6RJcttQl5aWxrvvvnvKvXTN8aywulJ/7suTewPX9denTx/8/PwwxvCvf/2Lyy67jLZt21ZYxsfHh44dO9bpJIk333zzlMtYlsW4ceMYN25cnev2NPbm9Zh/vQSlpRDRAUdCIlZklKvLEmmyXB7q3nzzzQqhqyozZ87krLPOKn+clZXF448/zoUXXsiwYcNqXLe5nRVWF+rPfXlyb9AwZ4XVRVxcHHFxcQAUFRUxbNgwQkND6/U1pHqmtATz79cwH77nnNB3oHMPnX9L1xYm0sS5PNRddtllXHzxxTUu88sNdlZWFtOnTycuLo7bbrvtlM/fnM4KO13qz315cm/QNPq7+uqrXfr6zY3JzsJ++W+wZwdYlnPs3BXjNH5OpBZcHuoCAwMJDAys1bInAl2XLl1ISEio01lnIiLStJnd27FfngU5WeDfynl2a+/zXF2WiNtweairraysLKZNm0bbtm254YYbyM3NLZ9Xl7PPRESkaTHGYD58D/Pv16CsFCKjcUxJxGrX3tWlibgVtwl127ZtIy0tjbS0NG6//fYK82oz+FhERJoeU1KMWfQKZvN654QBF+GYeDeWX/M821fkTLhNqBsyZEj5Kf8iIuL+TFamc/zcvl1gObDir8e67CqPvVSPSENzm1AnIiKew+z61jl+Li8HWgbguO3PWL36ubosEbemUCciIo3GGIP5z2rM0rlQVgZRnZ33bw3z7GsuijQGhToREWkUprgI88aLmE/eB8A6/1dYN9yJ1cLPxZWJeAaFOhERaXDmx3Tn7b5S9oDDgTX2Jqxf/1bj50TqkUKdiIg0KLP9a+xXn4RjuRAQ6Bw/16OPq8sS8TgKdSIi0iCMMZj1KzHLFoCxoeNZOBIexGoT7urSRDySQp2IiNQ7U1SIWfg85vNNAFgXDsWakIDl28LFlYl4LoU6ERGpVyYjDfvFx+HgfvDywrp6EtalV2j8nEgDU6gTEZF6Y777yjl+ruAYtA7Ccfv9WHFnu7oskWZBoU5ERM6YMQbz3nJM0j/BGOgSh+P2B7BC27q6NJFmQ6FORETOiCk8jr3gOfjiYwCsQcOxrpuM5ePr4spEmheFOhEROW3mSKpz/FxqCnh5Y117G9avRmr8nIgLKNSJiMhpMds+x57zNBzPh6AQ5+HW2B6uLkuk2VKoExGROjG2jVnzJubtxc7xc2d1dwa64FBXlybSrCnUiYhIrZnjBdjznoGtnwFgDbkc65pbsLx9XFuYiCjUiYhI7ZjDB7FffAzSDoG3N9Z1t+O4ZISryxKRnyjUiYjIKZmvPnXuoSs8DsFtnLf76hLn6rJE5BcU6kREpFrGtjHvLMasWuKcENcLx+T7sAJDXFuYiFSiUCciIlUyBcecZ7d+8z8ArGG/wRp7E5a3/nSINEX6nykiIpWYQynO8XPph8HHF+v6KTguHOrqskSkBgp1IiJSgfniI+z5z0FRIYSG4Uh4CKvTWa4uS0ROQaFOREQAMHYZ5q03MO8ud07o3hvHbX/Gah3k2sJEpFYU6kREmqC1a9fy9ttvk52dTVRUFBMnTqRHj4a7W4M5lof96pOQ/BUA1ogrsa66EcvLq8FeU0Tql8PVBYiISEUff/wxCxYs4KqrrmLWrFn06NGDxx9/nMzMzAZ5veK9uyibcY8z0Pn6Yt06FcfVNyvQibgZhToRkSZm1apVXHrppQwbNqx8L13btm1Zt25dvb+W/dmHpE+9CTKPQFgEjgefxHH+r+r9dUSk4enwq4hIE1JaWsrevXu58sorK0zv3bs3O3furHKdkpISSkpKyh9bloW/v3/599Wx167AXjrfudzZ/XHcMhUroPUZdtB0nOi9pp+BO1N/7q0h+lOoExFpQnJzc7Ftm6CgiicnBAUFkZ2dXeU6SUlJLFu2rPxxly5dmDVrFmFhYTW+VvGQkaSvWkLAb64haMLtHnu4NSIiwtUlNCj1597qsz+FOhGRJqiqT+/VfaKPj49n9OjRlZbLyMigtLS0+hdp0QqvGS8T3L0naWlpGGPOrOgmxrIsIiIiPLI3UH/uri79eXt7n/JDGijUiYg0KYGBgTgcjkp75XJycirtvTvBx8cHHx+fKued6o+FFRRSvpwn/uEEz+4N1J+7q8/+dKKEiEgT4u3tTUxMDNu2baswfdu2bXTr1s1FVYmIO9CeOhGRJmb06NE8//zzxMTEEBcXx4YNG8jMzGT48OGuLk1EmjCFOhGRJuaiiy4iLy+P5cuXc/ToUaKjo3nwwQdrNaZGRJovhToRkSZo5MiRjBw50tVliIgb0Zg6EREREQ/QbPfUeXvXvvW6LOuO1J/78uTeoHb9efrP4ExoO+fkyb2B+nN39bmds4wnnycsIiIi0kzo8GsNjh8/zv3338/x48ddXUqDUH/uy5N7A8/vrynx5J+1J/cG6s/dNUR/CnU1MMawb98+j73oofpzX57cG3h+f02JJ/+sPbk3UH/uriH6U6gTERER8QAKdSIiIiIeQKGuBj4+PowdO7baeyq6O/Xnvjy5N/D8/poST/5Ze3JvoP7cXUP0p7NfRURERDyA9tSJiIiIeACFOhEREREPoFAnIiIi4gEU6kREREQ8gGffUK0W1q5dy9tvv012djZRUVFMnDiRHj16VLt8cnIyCxcu5ODBg4SEhPDb3/6WESNGNGLFtZOUlMSWLVs4dOgQvr6+xMXFMWHCBNq3b1/tOt999x3Tp0+vNP2ZZ56hQ4cODVlunb355pssW7aswrSgoCBee+21atdxl/duypQpZGRkVJo+YsQIbrnllkrTm/r7lpyczNtvv82+ffs4evQoU6dO5fzzzy+fb4xh6dKlbNy4kWPHjtG1a1cmTZpEdHR0jc/76aefsmTJEo4cOUK7du249tprKzyvnFpdt39NVUP9jjUFtdmWu3N/69atY926deXbvKioKMaOHUu/fv0A9+6tKklJSSxevJhRo0YxceJEoJ57NM3YRx99ZMaPH282bNhgDhw4YObPn28mTJhgMjIyqlz+yJEjZsKECWb+/PnmwIEDZsOGDWb8+PHmk08+aeTKT23GjBnm/fffNykpKWbfvn1m5syZ5o477jDHjx+vdp1vv/3WXH311ebQoUPm6NGj5V9lZWWNWHntLFmyxNx7770V6szJyal2eXd673Jycir09fXXX5urr77afPvtt1Uu39Tfty+//NIsXrzYfPrpp+bqq682n332WYX5SUlJ5oYbbjCffvqp+eGHH8wzzzxjbrvtNlNQUFDtc+7cudNcc801ZsWKFebgwYNmxYoVZvz48WbXrl0N3Y7HqOv2rylriN+xpqI223J37u/zzz83X3zxhTl06JA5dOiQWbRokRk/frxJSUkxxrh3byf7/vvvTUJCgpk6daqZP39++fT67LFZH35dtWoVl156KcOGDSv/lNq2bVvWrVtX5fLr1q2jbdu2TJw4kaioKIYNG8bQoUN55513GrnyU0tMTGTIkCFER0fTuXNnEhISyMzMZO/evadcNygoiODg4PIvh6Np/po4HI4KdQYGBla7rDu9d4GBgRX6+vLLL2nXrh09e/ascb2m+r7169eP8ePHM3DgwErzjDGsWbOG+Ph4Bg4cSMeOHZkyZQpFRUVs3ry52udcvXo1vXv3Jj4+ng4dOhAfH8/ZZ5/N6tWrG7IVj1LX7V9T1hC/Y03Fqbbl7t7fueeeS//+/Wnfvj3t27fn2muvxc/Pj++//97te/ulwsJCnn/+eSZPnkyrVq3Kp9d3j01jq+8CpaWl7N27lz59+lSY3rt3b3bu3FnlOt9//z29e/euMK1v377s3buX0tLSBqu1PhQUFAAQEBBwymXvu+8+brvtNh599FG+/fbbhi7ttKWlpTF58mSmTJnCs88+y5EjR6pd1l3fu9LSUjZt2sTQoUOxLKvGZd3lfful9PR0srOzK/w/9PHxoWfPntX+PwTYtWtXpfezT58+7Nq1q8Fq9SSns/1zV6f7O9ZUnbwt96T+bNvmo48+oqioiLi4OI/qbc6cOfTr16/Sdqu+e2y2Y+pyc3OxbZugoKAK04OCgsjOzq5ynezs7CqXLysrIy8vj5CQkIYq94wYY1i4cCHdu3enY8eO1S4XEhLCbbfdRkxMDKWlpfz3v//lr3/9K4888sgp9xI1tq5duzJlyhTat29PdnY2K1as4OGHH+bpp5+mdevWlZZ31/duy5Yt5OfnM2TIkGqXcaf37WQn/q9V9d5kZmbWuF5wcHCFacHBwdX+35WKTmf7565O93esKapqW+4J/aWkpJCYmEhJSQl+fn5MnTqVqKio8lDjzr0BfPTRR+zbt4+ZM2dWmlff71+zDXUnVLX3o6Y9IifPMz/dkONUe1Fcae7cuaSkpPDoo4/WuNyJ3d8nxMXFkZmZyTvvvNPkwsGJQbQAHTt2JC4ujrvuuosPP/yQ0aNHV7mOO75377//Pn379iU0NLTaZdzpfatOde9NXRhjmvR72RTVdfvnzurjd8zVatqWu3N/7du358knnyQ/P5/PPvuMF154ocLJX+7cW2ZmJgsWLCAxMRFfX99ql6uvHpvt4dfAwEAcDkelT6U5OTmVEvMJVe0JyM3NxcvLq1aHNV1h3rx5fPHFFzzyyCO0adOmzuvHxcWRlpbWAJXVLz8/Pzp27Mjhw4ernO+O711GRgbbtm1j2LBhdV7XXd63E3vbqnpvqvt/eGK9uvzflYpOZ/vnrk73d6ypqW5b7gn9eXt7ExERwVlnncV1111H586dWbNmjUf0tnfvXnJycnjggQcYP34848ePJzk5mXfffZfx48eX91FfPTbbUOft7U1MTAzbtm2rMH3btm1069atynW6du1aafmvv/6amJgYvL2b1k5PYwxz587ls88+4y9/+Qvh4eGn9Tz79u2rdJirKSopKeHQoUPVHkZ1p/fuhPfff5+goCD69+9f53Xd5X0LDw8nODi4wntTWlpKcnJytf8PwRlav/nmmwrTtm3bRlxcXIPV6klOZ/vnrk73d6ypONW23N37q4oxhpKSEo/o7ZxzzuGpp57iiSeeKP8666yzGDRoEE888QTt2rWr1x6b5l+zRjJ69Gief/55YmJiiIuLY8OGDWRmZjJ8+HAAFi1aRFZWFnfeeSfgvE7Y2rVrWbhwIcOGDWPXrl385z//4e6773ZlG1WaO3cumzdv5r777sPf37/8U0DLli3LdwGf3N/q1asJCwsjOjq6fID+Z599xp/+9CdXtVGt119/nXPPPZe2bduSk5PD8uXLOX78OIMHDwbc+70D54DhDz74gMGDB+Pl5VVhnru9b4WFhRX2Gqanp7N//34CAgJo27Yto0aNIikpicjISCIiIkhKSqJFixYMGjSofJ3Zs2cTGhrKddddB8CoUaN45JFHeOuttzjvvPP4/PPP+eabb045xEB+dqrtnzupj9+xpupU23LLsty6v0WLFtGvXz/atGlDYWEhH330Ed999x2JiYlu3xuAv79/pbHsLVq0oHXr1uXT67PHZh3qLrroIvLy8li+fDlHjx4lOjqaBx98kLCwMACOHj1aYaBieHg4Dz74IAsXLmTt2rWEhIRw0003ccEFF7iqhWqduCzBtGnTKkxPSEgoH3R/cn+lpaX885//JCsrC19fX6Kjo3nggQdOa09RQ8vKyuK5554jNzeXwMBAunbtymOPPeYR7x3AN998Q2ZmJkOHDq00z93etz179lQYH/P6668DMHjwYKZMmcKYMWMoLi5mzpw55OfnExsbS2JiIv7+/uXrZGZmVhhz0q1bN/74xz/y73//myVLlhAREcEf//hHunbt2niNublTbf/cSX38jjVVtdmWu3N/OTk5zJ49m6NHj9KyZUs6depEYmJi+Vmi7txbbdVnj5ZxpxGHIiIiIlKlZjumTkRERMSTKNSJiIiIeACFOhEREREPoFAnIiIi4gEU6kREREQ8gEKdiIiIiAdQqBMRERHxAAp1IiIiIh5AoU5ERETEAyjUiYiIiHgAhToRERERD6BQJx6juLiY++67j7vuuouCgoLy6dnZ2dx6661MmzYN27ZdWKGIiEjDUagTj+Hr68s999xDbm4uL774IgC2bfOPf/wDgLvvvhuHQ7/yIiLimfQXTjxKZGQkkydPZsuWLaxZs4Zly5bx3XffcddddxESEuLq8kRERBqMt6sLEKlvF110EcnJyfzzn//Etm3i4+Pp3bu3q8sSERFpUNpTJx5p6NChlJWV4eXlxahRo1xdjoiISINTqBOPU1hYyOzZs4mMjMTX15eXX37Z1SWJiIg0OIU68TivvfYamZmZTJ06ldtvv53//e9/rFq1ytVliYiINCiFOvEoGzduZNOmTUyaNIno6GguuOACLrvsMv71r3+xe/duV5cnIiLSYBTqxGOkpKQwf/58Bg8ezJAhQ8qnX3/99XTq1IlnnnmG/Px81xUoIiLSgCxjjHF1ESIiIiJyZrSnTkRERMQDKNSJiIiIeACFOhEREREPoFAnIiIi4gEU6kREREQ8gEKdiIiIiAdQqBMRERHxAAp1IiIiIh5AoU5ERETEAyjUiYiIiHgAhToRERERD/D//WScEP2mn50AAAAASUVORK5CYII="
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"# plot trajectory\n",
|
|
"plt.subplot(2, 2, 1)\n",
|
|
"plt.plot(x_bar[0, :], x_bar[1, :])\n",
|
|
"plt.plot(np.linspace(0, 10, T + 1), np.zeros(T + 1), \"b-\")\n",
|
|
"plt.axis(\"equal\")\n",
|
|
"plt.ylabel(\"y\")\n",
|
|
"plt.xlabel(\"x\")\n",
|
|
"\n",
|
|
"plt.subplot(2, 2, 2)\n",
|
|
"plt.plot(np.degrees(x_bar[2, :]))\n",
|
|
"plt.ylabel(\"theta(t) [deg]\")\n",
|
|
"# plt.xlabel('time')\n",
|
|
"\n",
|
|
"\n",
|
|
"plt.tight_layout()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## State Space Linearized Model"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 15,
|
|
"metadata": {
|
|
"ExecuteTime": {
|
|
"end_time": "2024-10-22T09:30:50.489618Z",
|
|
"start_time": "2024-10-22T09:30:50.483390Z"
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"\"\"\"\n",
|
|
"Control problem statement.\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"N = 4 # number of state variables\n",
|
|
"M = 2 # number of control variables\n",
|
|
"T = 20 # Prediction Horizon\n",
|
|
"DT = 0.2 # discretization step"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 16,
|
|
"metadata": {
|
|
"ExecuteTime": {
|
|
"end_time": "2024-10-22T09:30:50.962772Z",
|
|
"start_time": "2024-10-22T09:30:50.957952Z"
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def get_linear_model(x_bar, u_bar):\n",
|
|
" \"\"\"\n",
|
|
" Computes the LTI approximated state space model x' = Ax + Bu + C\n",
|
|
" \"\"\"\n",
|
|
"\n",
|
|
" L = 0.3 # vehicle wheelbase\n",
|
|
"\n",
|
|
" x = x_bar[0]\n",
|
|
" y = x_bar[1]\n",
|
|
" v = x_bar[2]\n",
|
|
" theta = x_bar[3]\n",
|
|
"\n",
|
|
" a = u_bar[0]\n",
|
|
" delta = u_bar[1]\n",
|
|
"\n",
|
|
" A = np.zeros((N, N))\n",
|
|
" A[0, 2] = np.cos(theta)\n",
|
|
" A[0, 3] = -v * np.sin(theta)\n",
|
|
" A[1, 2] = np.sin(theta)\n",
|
|
" A[1, 3] = v * np.cos(theta)\n",
|
|
" A[3, 2] = v * np.tan(delta) / L\n",
|
|
" A_lin = np.eye(N) + DT * A\n",
|
|
"\n",
|
|
" B = np.zeros((N, M))\n",
|
|
" B[2, 0] = 1\n",
|
|
" B[3, 1] = v / (L * np.cos(delta) ** 2)\n",
|
|
" B_lin = DT * B\n",
|
|
"\n",
|
|
" f_xu = np.array(\n",
|
|
" [v * np.cos(theta), v * np.sin(theta), a, v * np.tan(delta) / L]\n",
|
|
" ).reshape(N, 1)\n",
|
|
" C_lin = DT * (\n",
|
|
" f_xu - np.dot(A, x_bar.reshape(N, 1)) - np.dot(B, u_bar.reshape(M, 1))\n",
|
|
" )\n",
|
|
"\n",
|
|
" return np.round(A_lin, 4), np.round(B_lin, 4), np.round(C_lin, 4)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 17,
|
|
"metadata": {
|
|
"ExecuteTime": {
|
|
"end_time": "2024-10-22T09:30:53.505536Z",
|
|
"start_time": "2024-10-22T09:30:53.502223Z"
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"CPU times: user 3.48 ms, sys: 2.16 ms, total: 5.64 ms\n",
|
|
"Wall time: 4.17 ms\n"
|
|
]
|
|
},
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:17: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n",
|
|
" A[0, 2] = np.cos(theta)\n",
|
|
"/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:18: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n",
|
|
" A[0, 3] = -v * np.sin(theta)\n",
|
|
"/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:19: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n",
|
|
" A[1, 2] = np.sin(theta)\n",
|
|
"/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:20: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n",
|
|
" A[1, 3] = v * np.cos(theta)\n",
|
|
"/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:21: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n",
|
|
" A[3, 2] = v * np.tan(delta) / L\n",
|
|
"/var/folders/hd/8kg_jtmd6svgg_sc384pbcdm0000gn/T/ipykernel_33582/46870782.py:26: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n",
|
|
" B[3, 1] = v / (L * np.cos(delta) ** 2)\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"%%time\n",
|
|
"\n",
|
|
"u_bar = np.zeros((M, T))\n",
|
|
"u_bar[0, :] = 0.2 # m/s\n",
|
|
"u_bar[1, :] = np.radians(-np.pi / 4) # rad\n",
|
|
"\n",
|
|
"x0 = np.zeros(N)\n",
|
|
"x0[0] = 0\n",
|
|
"x0[1] = 1\n",
|
|
"x0[2] = 0\n",
|
|
"x0[3] = np.radians(0)\n",
|
|
"\n",
|
|
"x_bar = np.zeros((N, T + 1))\n",
|
|
"x_bar[:, 0] = x0\n",
|
|
"\n",
|
|
"for t in range(1, T + 1):\n",
|
|
" xt = x_bar[:, t - 1].reshape(N, 1)\n",
|
|
" ut = u_bar[:, t - 1].reshape(M, 1)\n",
|
|
"\n",
|
|
" A, B, C = get_linear_model(xt, ut)\n",
|
|
"\n",
|
|
" xt_plus_one = np.dot(A, xt) + np.dot(B, ut) + C\n",
|
|
"\n",
|
|
" x_bar[:, t] = np.squeeze(xt_plus_one)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 18,
|
|
"metadata": {
|
|
"ExecuteTime": {
|
|
"end_time": "2024-10-22T09:30:54.635354Z",
|
|
"start_time": "2024-10-22T09:30:54.569920Z"
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": "<Figure size 640x480 with 2 Axes>",
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAEDCAYAAABTZPIVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/A0lEQVR4nO3deViVdf7/8ed92MSFTUVQQEVFs3KryWpstByXMWeMMrPdstSwpnGmLKPfpH2dHKtpmWx3bRrNXGjT0rRlss12S0ozVNwQCAEFWQ7n8/vjJEWCoAL3OTevx3VxwbnPfZ/zfnMfbl7nPp/7vi1jjEFERERE/JrL7gJERERE5OQp1ImIiIg4gEKdiIiIiAMo1ImIiIg4gEKdiIiIiAMo1ImIiIg4gEKdiIiIiAMo1ImIiIg4gEKdiIiIiAME2l2AXQ4cOIDb7a51vrZt25KTk9MIFdlD/fkvJ/cGde8vMDCQyMjIRqjI/2g75+zeQP35u/rezjXZUOd2uykvLz/mPJZlVc7rxKupqT//5eTewPn9NZamvp1zcm+g/vxdQ/Snj19FREREHEChTkRERMQBFOpEREREHEChTkRERMQBFOpEREREGpkpL8PUciDT8VKoExEREWlEJi+Hitl3kP/sQ/X6uE32lCYiIiIijc1s+RrP0/fDwQKKf8zBOn8khNfPuTYV6kREREQamDEGs/5VzLL54PFAfCLtZjxCjseqt/PUKdSJiIiINCBTWor5zxzMx+8CYPUfiOuamwls1x727au351GoExEREWkgJnc/nifug13bweXCuvR6rMF/xHLV/2ENCnUiIiIiDcCkf4nnmQeg6CC0Csc1cSpW99Mb7PkU6kRERETqkTEGszYNs+I5MB7o2BXXTdOwWrdt0OdVqBMRERGpJ6a0BLPw35hPNwBgnTsY68pJWMEhDf7cfn2eurS0NMaMGcPChQvtLkVERESaOJO9D8+s272BLiAA64pJWOP+3CiBDvx4T922bdtYt24dHTt2tLsUERERaeLMN5/hefZBKC6CsAhck+7E6tazUWvwyz11JSUlPPbYY0ycOJEWLVrYXY6IiIg0UcYYPKtexPPve72BrnMSrrsfbvRAB366p27u3Ln07duXXr16sXLlymPOW15eTvkvrq1mWRahoaGVPx/Lkftrm89fqT//5eTewPn9iYgzmJJiPAsehc8/BMD63TCssROwgoJsqcfvQt3777/P9u3bmTVrVp3mT0tLY/ny5ZW3O3fuzOzZs2nbtu5HoMTExBx3nf5E/fkvJ/cGzu9PRPyXydrjPf/cvl0QGIh1+URcvxtma01+Fepyc3NZuHAhqampBAcH12mZ5ORkRo4cWXn7yDv/nJwc3G73MZe1LIuYmBiysrLq7RIevkT9+S8n9wbH119gYOBxvUkTETlZ5qtP8Mz7Fxwuhogo7/i5Lj3sLsu/Ql1GRgYFBQXceeedldM8Hg/ffvstb7zxBosXL8b1qzM0BwUFEVTDbtC6/jM0xjjyH+cR6s9/Obk3cH5/IuJfjMeDeW0p5tUl3glde+KadAdWeKS9hf3Er0Ld6aefzoMPPlhl2pNPPkn79u0ZNWrUUYFOREREpD6Y4iI88x+GrzYCYA0agXXZeKxAe8bPVcevQl1oaCgJCQlVpoWEhNCqVaujpouIiIjUB7M3E88Ts2D/HggMwroqBddvB9td1lH8KtSJiIiINCbz+Yd45j8CpYchqo33cl+dutldVrX8PtRNnz7d7hJERETEYYynAvPyYszqZd4J3U/HNeF2rLAIW+s6Fr8PdSIiIiL1yRQdwjP3QfjmcwCs34/CGj0OKyDA5sqOTaFORKSRVFRUsGzZMt577z3y8/OJjIxk0KBBXHzxxZUHehljWLZsGevXr+fQoUN069aN8ePHEx8fb3P1Ik2D2b3De/65nCwIDsa6+mZcZw+yu6w6UagTEWkkL7/8Mm+++SaTJ08mLi6OjIwMnnjiCZo3b86IESMq51m1ahUpKSnExsaycuVKZs6cySOPPFJ5NRwRaRieTzZgFj4KZaXQOhpXyl1YCYl2l1VnOgeIiEgj2bp1K2eeeSb9+vUjOjqas88+m169evHDDz8A3r10q1evJjk5mf79+5OQkMDkyZMpLS1lw4YNNlcv4lymogLP8gWYZ+73BrpTeuO6+yG/CnSgPXUiIo2mR48evPnmm+zdu5f27duzY8cOtmzZwrXXXgtAdnY2+fn59O7du3KZoKAgevbsyZYtWxgyZEi1j6trXFfPyb2B+qsv5mAB5pkHMN9+5X2+4ZfgSr66wcfPNUR/CnUiIo1k1KhRFBcXM2XKFFwuFx6Ph7FjxzJgwAAA8vPzAQgPD6+yXHh4OLm5uTU+rq5xfWxO7g3U38ko++E7cmfdjsnehxXSjKgp99D8vOrfPDWU+uxPoU5EpJF88MEHvPfee/z5z38mPj6eHTt2sHDhwsoDJo749Tv32i6VpmtcV8/JvYH6O1meD9/G89wcKC+DtrG4Jt9FQVwnCvbtq/fnqk5DXONaoU5EpJE8//zzjBo1it/+9rcAJCQkkJOTw0svvcSgQYOIiIgAqDwy9ojCwsKj9t79kq5xfWxO7g3U33E/ntuNWb4As/5V74TTzsB1w9+gRUtbfo/12Z8OlBARaSSlpaVHXaPa5XJVbtCjo6OJiIhg06ZNlfe73W7S09Pp3r17o9Yq4kSmMB/Pw3+vDHTWiDG4brkbq0VLmyurH9pTJyLSSM444wxWrlxJmzZtiIuLY8eOHbz22mucf/75gPfjmBEjRpCWlkZsbCwxMTGkpaUREhJSOe5ORE6M2f49nidnwYFcCAnFdf1fsPqdY3dZ9UqhTkSkkVx//fUsXbqUuXPnUlBQQFRUFEOGDGH06NGV84waNYqysjLmzp1LUVERXbt2JTU1VeeoEzkJnvfXY55/Atzl0K4Drsl3YcU674TeCnUiIo0kNDSUcePGMW7cuBrnsSyLMWPGMGbMmMYrTMShjLscs3Qe5p3V3gm9z8J1/RSs5i3sLayBKNSJiIiI45iCA3iemg3b0gGw/nQF1oVjsFzOPZxAoU5EREQcxfzwHZ6n/gn5eRDaHNf4v2L1PsvushqcQp2IiIg4hud/b2AWPwMVboiN916/NaaD3WU1CoU6ERER8XumvByz5GnMe2u9E/qdg+u6W7GaNbe3sEakUCciIiJ+zRz40Xu6ku1bwbKwLroK6w+jHXtd3Joo1ImIiIjfMls343l6NhTmQ/MWuG68Deu0M+wuyxYKdSIiIuJ3jDGYd1Zjls6Figro0NE7fi461u7SbKNQJyIiIn7FlJdhnn8S88F6AKzfnId17S1YIc1srsxeCnUiIiLiN0xeDp4nZsHObWC5sEZfizXkoiY3fq46CnUiIiLiF8yWr/E8fT8cLICWrXDdeDtWzz52l+UzFOpERETEpxljMOtfwSxbAB4PJCTiumkaVpt2dpfmUxTqRERExGeZ0hLMc3MwH78LgHX2IKyrJ2MFh9hcme9RqBMRERGf5N6/l4p/3gG7MsDlwhozHuuCkRo/VwOFOhEREfE5nvQv2f/sg97xc63CcU28A6v7aXaX5dMU6kRERMRnGGMwa1ZiVv4HjAc6dfWOn4tqa3dpPk+hTkRERHyCKS3BLPw35tMNALQY8kdKLh4HgUH2FuYnFOpERETEdiZ7r/f8c3t2QkAArrETiLz8erKysjDG2F2eX1CoExEREVuZrz/DM/dBKC6C8EhcE+/AlXSqDog4Tgp1IiIiYgtjDGb1MszL/wVjILE7rpvuxIpobXdpfkmhTkRERBqdKSnGM/8R+OIjAKzfDcMaOwErSOPnTpRCnYiIiDQqk7XbO35u3y4IDMS6YhKu84baXZbfU6gTERGRRmO+2ohn3kNwuBgiWns/bk3sbndZjuBXoS4tLY2NGzeyZ88egoODSUpK4qqrrqJ9+/Z2lyYiIiLHYDwezGtLMa8u8U7o2hPXpDuwwiPtLcxB/CrUpaenM2zYMLp06UJFRQUvvPACM2fO5KGHHqJZs2Z2lyciIiLVMMVF3r1zmz4BwDr/Qqwx12Pp/HP1yq9CXWpqapXbKSkp3HDDDWRkZNCzZ0+bqhIREZGamL2Z3vFz+/dAYBDW1Sm4zh1sd1mO5Feh7teKi4sBaNmyZY3zlJeXU15eXnnbsixCQ0Mrfz6WI/c79Tw56s9/Obk3cH5/Ik2F+fwDPPMfhdLDENUGV8pdWB272l2WY/ltqDPGsGjRInr06EFCQkKN86WlpbF8+fLK2507d2b27Nm0bVv3a8jFxMScVK2+Tv35Lyf3Bs7vT8SpjKcC8/JizOpl3gndT8c1cSpWq3B7C3M4vw118+bNIzMzk3vvvfeY8yUnJzNy5MjK20fe+efk5OB2u4+5rGVZxMTEOPYSJerPfzm5Nzi+/gIDA4/rTZqINCxTdMh7dYhvPgfA+v0orNHjsAICbK7M+fwy1M2fP5/PPvuMGTNm0Lr1sc86HRQURFANJzKs6z9DY4wj/3Eeof78l5N7A2f2l5eXx/PPP8+XX35JWVkZsbGx3HTTTSQmJgLenpctW8b69es5dOgQ3bp1Y/z48cTHx9tcuUjtzO7t3vFzOVkQHIx19c24zh5kd1lNhl+FOmMM8+fPZ+PGjUyfPp3o6Gi7SxIRqbNDhw7x//7f/+PUU0/lrrvuIiwsjP3799O8efPKeV5++WVWrVpFSkoKsbGxrFy5kpkzZ/LII49UjgcW8UWeT97DLPw3lJVC62jv+LmERLvLalL8KtTNmzePDRs2MHXqVEJDQ8nPzwegefPmBAcH21uciEgtXn75ZVq3bk1KSkrltF++OTXGsHr1apKTk+nfvz8AkydP5sYbb2TDhg0MGTKk0WsWqY2pqMCsfA6zNs074ZTeuCbcjtUyzN7CmiC/CnVr164FYPr06VWmp6SkMGjQoMYvSESalL1795KXl0dZWRlhYWG0b9++yl622nz66af07t2bhx56iPT0dKKiohg6dCi///3vAcjOziY/P5/evXtXLhMUFETPnj3ZsmVLjaFOR/lXz8m9gW/0Zw4WYJ55EPPtl95ahl+C6+KrsVwnP37OF/prSA3Rn1+FuhdffNHuEkSkidm6dStvvvkmX375JYWFhVXuc7lcdOrUifPOO49BgwbVGvCys7N58803ufDCC0lOTmbbtm0sWLCAoKAgBg4cWPnpQ3h41SMEw8PDyc3NrfFxdZT/sTm5N7Cvv7IfviN31u2Y7H1YIc2ImnIPzc+r/73JWn9151ehTkSksezYsYOFCxfy7bff0qFDB/r3709iYiJhYWEEBwdz6NAh9u/fz/fff8+SJUtYunRp5dH2gYHVb1o9Hg9dunThiiuuALzha9euXaxdu5aBAwdWzvfrd+61HSyio/yr5+TewN7+PB++jee5OVBeBm1jcU2+i4K4ThTs21dvz6H197O6HuWvUCciUo277rqLAQMGcM0111QemVqTkpISPvjgA15++WUqKiq45JJLqp0vMjKSuLi4KtPi4uL4+OOPAYiIiAAgPz+fyMifr4dZWFh41N67X9JR/sfm5N6gcfszbjdmxULMule8E047A9cNf4MWLRusBq2/ulOoExGpxr/+9S9iY2PrNG+zZs244IILGDRo0DE/Ju3evTt79+6tMm3v3r2V78Cjo6OJiIhg06ZNdO7cGQC32016ejpXXnnlCXYiUj9MYT6eZx6ALV8DYF04ButPl9fL+DmpHy67CxAR8UV1DXS/5HK5jnmqpQsvvJDvv/+elStXkpWVxYYNG1i/fj3Dhg0DvB/HjBgxgrS0NDZu3EhmZiaPP/44ISEhDBgw4IR7ETlZZsf3eP7xV2+gCwnFddM0XBddpUDnY7SnTkSkFjfffDO33XYbnTp1Ouq+zMxM7r//fubMmVPr43Tt2pXbbruNxYsXs2LFCqKjo7n22ms577zzKucZNWoUZWVlzJ07l6KiIrp27UpqaqrOUSe28by/HvP8E+Auh3YdcE2+CytWJ8P2RQp1IiK1ONYBB+Xl5eTk5NT5sc444wzOOOOMGu+3LIsxY8YwZsyY465TpD4Zdzlm6TzMO6u9E3qfhev6KVjNW9hbmNRIoU5E5CTs379fe9HEcUzBATxP/RO2fQuW5R07N2IMlkujtnyZQp2ISDXeeecd3n333crbc+fOPSq8lZWVsXPnTnr27NnY5Yk0GPPDd95Al58HoS1w3fBXrF6/sbssqQOFOhGRapSVlVU52XBRUVGVqzaA91Qi5557rj4qFcfw/O8NzOJnoMINsfG4JqditWtvd1lSRwp1IiLVGDp0KEOHDgW811/929/+Vu2BEiJOYMrLMUuexrznvRwn/c7Fdd2fsZrV/TJ4Yj+FOhGRWjz++ON2lyDSYMyBH/E8OQu2b/WOn0u+Gmv4JY695qqTKdSJiFQjNzeXNm3aHPdyeXl5REVFNUBFIvXPbN2M5+nZUJgPzVviuvE2rNP62V2WnCAdxiIiUo1bb72VBQsWkJWVVeu8brebDz/8kNtvv5233nqrEaoTOTnGGDxvr8Lz0N3eQNehI67UfynQ+TntqauBcbshN4tydykmJxtjAMsClwWW66efAyAg4OfvAQEQEAgBAdptLeLn7r77bhYtWsQbb7xB165dOfXUU+ncuTPh4eEEBQVx6NAh9u/fz9atW/nqq68oKSlhxIgRjBw50u7SRY7JlJVinn8S86H3DYj1m/Owrr0FK6SZzZXJyVKoq8nBAir+Xwq1v0evhmVBYJD3KygIgoK9X8HBEBwCwc0gOASrWTPvz81CIeSn76HNsUKbQ2gLaN7i5+/NW2AFVn/BbhGpf6eccgr//Oc/+eKLL3jzzTd5/fXXKSsrO2q+6Ohohg0bxpAhQ4iMjLShUpG6Mz/meMfP7dwGlgtr9LVYQy7SjgiHUKiriWVBi1a4XC48FRWAAQMYDxgDHs/PX8ZTdVljoLzM+3W45qcwxzmdkFBo0dL71TIMq2UYtGgFrcKgVThWq3BoFQFh4RAW6Q2C+kMVOSl9+/alb9++uN1uduzYwYEDBygrK6NVq1bExcVp/Jz4DfPdJjxP3w+HCqFlK1wTpmKd0tvusqQeKdTVwIqIIvDRxcTGxrJv3z6MqTFqYTweqKjwntenosJ7fTx3OZSXg7vM+72sDMpLoawUU1oKpSVQVgIlh70/l5RASTHmcDGUFMPhn76Ki+BwkfeJSg97v/K8lyT6dUVHVRgYCGER3oAXEYUVEQXhUd6fI9tAVFs84WH19SsTcbTAwEC6du1qdxkix80Yg1n3Cmb5Au+OiIREXDdNw2rTzu7SpJ4p1NUDy+UCl8v7UWtd5j/OxzeeCm/AKzr001ch5tBB77utQ4VwsBBzqAAKC+BgARzM94ZBtxvycr1fVA19R37eAxDaHKLaQlRbrNZtISoaotp4/+DbtoNWEdrjJyLih0xpKea5OZiN3qujWGcPwrp6MlZwiM2VSUNQqPMDlivA+zFri1Y/T6tlGVNe5g15hQcgPw9TkOe95Et+Hib/Rzjw09fhIm9g3LMT9uysNvgRHAKto6FNO6zoWIiOxWrr/U7raKxAvYzE+fbt28ebb77Jnj17jhpbZ1kWf//7322qTKR6JicLzxOzYPd2cLmwxozHumCk3qQ7mP4bO5QVFAyt23q/qD4EWpZFu/Awsr7djPkxG5OX4/1o98ccTF425O73Br+yUti3C/btqgx6lYHP5YI27aBdB6x2HaBde+8lZWLjITxSGw9xhMzMTFJTU4mKiiIrK4uOHTty8OBB8vLyaN26Ne3a6WMs8S0m/Qs8zzwIRQehVTiuiXdgdT/N7rKkgSnUNXGu5i2w2sdDbFy1wc+4y71BL2c/JicLcrIw2fsg56evsjLI3gfZ+zBff+pd5sjCoS28jxsbB7EJWO0ToEMCRLZR2BO/smTJEnr37s2UKVO44oormDRpEomJiXz++ec8+eSTjB071u4SRYCfxs+tTcOseM57EF+nbrhuuhMrqq3dpUkjUKiTY7ICgyC6PUS3Pyr0GWO8H+nu34PZv/fn71m7IWe/96PdjC2YjC3e+Y8sGNoc2idgdegI7TtixXeCuM5YzVs0Ymcidbd9+3ZuuOGGyjcjRw6c6tevH3/84x9ZvHgxM2bMsLNEEUzJYcyixzCfbgDA+u3vsa6c5P3kRpoEhTo5YZZlQWRriGyN1aNXlftMeRns34vZt/vnj273ZsL+Pd4xfD98h/nhO++8RxZqHQ1xnbDiE7ESEiGhi/eADe3VE5sVFRXRsmVLXC4XAQEBFBUVVd6XmJjI8uXLbaxOBEz2Xu/4uT07ISAQa+wNWAP/oO1nE3NCoS4jI4PExMT6rkUcxAoK9ga0uE5Vpht3uTfsHTkwY89O2LX9p7F82fBjNuarjT8HvZatIKELVkIX7/dOXb0HbGhDJY0oKiqKwsJCAGJiYkhPT6dXL+8bmczMTJo105n4xT6erz/F8+yD3rMehEfimnQHVteedpclNjihUDdt2jS6du3K8OHDOeeccwjU0Y9SR1ZgEHTo6P3o9RdM0UHYvQOzazvs3o7JzIC9mXDoIKR/iUn/0jsfeINep25YP33RuRtWmM7kLw2ne/fubN26lbPOOosBAwawbNky8vPzCQwM5J133uG8886zu0RpgozHQ+EL8/A8/5T3pPddengDXURru0sTm5xQGktJSWHNmjXMmTOH5557jsGDBzNkyBBat9YLSU6M1aIVdD8dq/vpldNMeZl3b17mD7AzA7NzG+ze4Q1633yO+ebzKh/dWl16QGJ3rMQeEN9J40ik3lx88cUcOHAAgIsuuoj8/Hw2bNiAZVmcc845XH311TZXKE2NOVyMWfAoBV98CID1u+FYY2/EquP5UsWZLHOsSyXUYtu2bbzxxht8+OGHeDwezjjjDP7whz9w6qmn1meNDSInJ4fy8vJjzmNZVp2uKOGv/LE/U17u3aO3Yyts/x6z43vvgRm/rj8oGDp2oVWfsyiKTYAuPbzB0SH8cd0dj+PpLygoiLZtdWRfdZr6ds6pvZms3Xgev8+77QsMwnXFRKzzhtpdVr1z6vo7oiG2cyf1uWnXrl25+eabueaaa1i3bh3r1q3j3nvvJS4ujuHDhzNw4ECCg7W3ROqPFRTk/bi1czc43zvNHC6G7VsxR460zdjiPTfTtm85uO3bnxdun4DV9RTo2hMr6VSs1tH2NCF+54knnmD06NFERx/9msnJyWHZsmWkpKTYUJk0NebLj/HMf9h7wFlka6Lv/hd54a0dGXrk+NXLYLjAwEBCQkIqx9aVlpYyd+5cVq5cyZQpU0hKSqqPpxGplhXaHHr2werZB/jpdBP790LGd4Tu2UHRpk8haw/szfQegfu/Nd6PbVtHYyWd9tPHvqfpOohSo3fffZehQ4dWG+oOHjzIu+++q1AnDcp4PJjXXsC8+oJ3QreeBEy6k5AePWHfPnuLE59xUqFu586drFmzhg0bNuB2uzn77LP585//TNeuXdm5cyfPPPMMzz77LA888EB91StSK8uyIKYDVmwcUbGxlO7bh6cwH7Z9i9mWjvk+HXZu8x5p++Fb8OFb3pAX1dZ7xvUevbFO6Y0VqTGiUrtDhw4RpHFM0oBMcZF379xXGwG8l/q69HqNn5OjnFCo++CDD1izZg3fffcdYWFhjBw5kqFDhxIREVE5T8eOHbn88sv5xz/+UV+1ipwwq1U49D0bq+/ZAJiSYtj2HWbr15gt33hDXl4O5sO34cO3vSEvJg7rlF5Yp/SB7qdhNW9pZwvSyNLT00lPT6+8vX79er788ssq85SVlfHJJ58QFxfXyNVJU2H2ZnrHz2XvhaBgrKtScJ17gd1liY86oVD36KOP0qlTJ2666SYGDBhQ4ylN2rZtq0P9xSdZzZrDaf2wTusHeM/Ezg/fYb7bhPn2K8j8AbJ2Y7J2Y95eDZYLOnXFOrUf1ql9oXMSVkCAzV1IQ9q8eXOVkwq/9dZb1c7Xpk0bxo8f31hlSRNiPv8Az/xHofQwRLXFlTINq2NXu8sSH3ZCoW7GjBn06NGj1vnatWuncSbiF6xmoXBqX29gA0zRIdjyNebbrzDffeUdk7d9K2b7VsxrL0DzFt6PaU/zhjxdV9F5Ro0axfDhwzHGcOONN5Kamkrnzp2rzBMUFKQTD0u9M54KzEv/xbz+05uK7qfjmjjV+4mDyDGcUKirS6BrSGvWrOGVV14hPz+fuLg4xo0bxymnnGJrTeIsVouW0O8crH7nAGDycjHpX8DmL7wnQi4+BJ9/gPn8AwxgDb8E1yXX2lqz1K/g4ODKo/fnzJlDZGRkvZ9oPS0tjSVLljBixAjGjRsHeA/0WbZsGevXr+fQoUN069aN8ePHEx8fX6/PLb7JFB30Xh1i8xcAWENGYV0yTp8MSJ343aUgPvjgAxYuXMgNN9xA9+7dWbduHffddx8PP/wwbdq0sbs8cSgrqg3WgCEwYAjGUwE7tmE2f4HZ/DlkbIVfXQ5NnOXI+aH27NlDeno6Bw8e5IILLiAiIoK8vDxatmx53Kdv2rZtG+vWraNjx6pXV3n55ZdZtWoVKSkpxMbGsnLlSmbOnMkjjzxCaGhovfUkvsfs3u69fmtOFgQHY11zC67+A+0uS/yI34W61157jQsuuIDBgwcDMG7cOL766ivWrl3LFVdcUW/PYwwcPmxRVATFxdZR57Z1AstSfycmEGJ6eL8GX44pPgQBQVjFjXc92qay7nylN4/Hw9NPP80777xTOa1Pnz5ERETwzDPP0LlzZy677LI6P15JSQmPPfYYEydOZOXKlZXTjTGsXr2a5ORk+vfvD8DkyZO58cYb2bBhA0OGDKm3nsS3eDb+D7PoMSgrhdbRuCanYsV3rn1BkV/wq1DndrvJyMjgoosuqjK9V69ebNmypdplysvLq5xR3bKsyne7x7oo/OHDFl27xvx0K6bG+ZxB/fkvJ/cGP/zgIjTUY3cZrFy5kg0bNnD11VfTp08f/va3v1Xe17dvX955553jCnVz586lb9++9OrVq0qoy87OJj8/n969e1dOCwoKomfPnmzZsqXGUHei27lf3l/bfP7IH3ozFRV4Vi7CrEkDwOrZB9eE27FahtW6rD/0dzLU3/Hzq1BXWFiIx+MhPLzqYNHw8HDy8/OrXSYtLa3KEWydO3dm9uzZtV5uo6jopMsVkZPUrl07WrSwuwp45513uOSSSxg5ciQeT9WQGR0dTXZ2dp0f6/3332f79u3MmjXrqPuObMeq28bl5ubW+Jgnup37pZgY575B8NXeKgry+XH2XZT+dP65VqOvJfyalOMeP+er/dUX9Vd3fhXqjqgu1daUdJOTkxk5cuRR8+Xk5OB2u2t8DmO8ewnatWvH/v37HXkJFsuy1J+fcnJv8HN/hYVZFBYeu7/AwMAGv/ZrXl5ejVfGCQoKoqSkpE6Pk5uby8KFC0lNTT3mGLxfb89qW8cnup07Mm9MTAxZWVmOey35cm9m5zYqnpgFP2ZDSDNc193K4TMHcPg43iD4cn/1Qf39rK7bOb8KdWFhYbhcrqP2yhUUFBz1zvaIoKCgGs/2XtsvMTTUQ4sW3u9OfUGpP//k5N7g5/4KC41P9BceHl7j3ri9e/cSFRVVp8fJyMigoKCAO++8s3Kax+Ph22+/5Y033uCRRx4BvHvsIiMjK+cpLCyscRsHJ7ed++V8vvC7bgi+1pvnw7cx/3kcysugbYx3/FyHjidco6/1V9/UX935VagLDAwkMTGRTZs2cdZZZ1VO37RpE7/5zW9srExEnKxv376sXLmy8uAI8AbP4uJiXn/9dc4444w6Pc7pp5/Ogw8+WGXak08+Sfv27Rk1ahTt2rUjIiKCTZs2VZ4Tz+12k56ezpVXXlmvPUnjM243ZvkCzPpXvRNOPxPX+L96T6EkUg/8KtQBjBw5kscee4zExESSkpJYt24dubm5OipMRBrMmDFj+OKLL5gyZQqnnnoqAEuWLGHXrl0EBAQwevToOj1OaGgoCQkJVaaFhITQqlWryukjRowgLS2N2NhYYmJiSEtLIyQkhAEDBtRvU9KoTGE+nqdnw9bNAFgjL8P64+VYLpfNlYmT+F2oO/fcczl48CArVqzgwIEDxMfHM23atAYfUyMiTVdERASzZs3ixRdf5IsvvsDlcrFz50769evHZZddRsuW9benZdSoUZSVlTF37lyKioro2rUrqampOkedHzPbv8fz5Cw4kAvNQnFdP6XyOtQi9cnvQh3AsGHDGDZsmN1liEgTEhERwYQJE+r9cadPn17ltmVZjBkzhjFjxtT7c0nj82x4E/Pfp8BdDjEdcKXchRWrq4NIw/DLUCciIuLLjLscs3Qu5p3XvRN6n+UdPxfa3N7CxNEU6kRE6uC7775jw4YN5OTkUFZWVuU+y7L4+9//blNl4mtMfp53/Ny2b8GysP50OdaIMRo/Jw1OoU5EpBZvv/02Tz31FC1btiQ2Nvao04c4+XQLcnzMD9/hefKfUJAHoS1w3fBXrF46O4M0DoU6EZFavPLKK5xzzjlMnjy5xvPBiXj+9wZm8TNQ4YbYeO/559q1t7ssaUIU6kREapGTk8N1112nQCfVMuXlmCVPY95b651wxrm4xt2K1UxHLEvjUqgTEalFhw4dKCgosLsM8UEmLxfPU/+E7Vu94+eSr8EafrFjL0Ivvk2jNkVEanH55Zfz0ksvkZeXZ3cp4kPM1s14Zk7xBrrmLXH9+R5cf7hEgU5soz11IiLVmD17dpXbxcXF3HrrrXTq1Omokw1blsXUqVMbszyxkTEG8/YqzIvzoKIC4jp5zz/XNsbu0qSJU6gTEalGZmZmldsul4uwsDDy8vKO2mOnPTNNhykrxTz/BObDtwGwzvod1jU3Y4U0s7kyEYU6EZFqPf7445U/p6en07lz52ov1VVSUkJGRkZjliY2MT9m43liFmT+AJYLa/Q4rCGjFOrFZ2hMnYhILWbMmMGePXuqvW/v3r3MmDGjkSuSxma+24Rn5l+9ga5lGK4pM3ANvUiBTnyK9tSJiJwEt9uNS1cKcCxjDObNlzErFoLHAwldcKVMw2odbXdpIkdRqBMRqUZxcTHFxcWVt/Pz88nNza0yT1lZGe+++y4RERGNXJ00BlNainnuMczG/wFgnXM+1lUpWMEhNlcmUj2FOhGRaqxatYrly5dX3n7ggQdqnDc5ObkxSpJGZHKyvOPndm+HgACsS8djXXChPm4Vn6ZQJyJSjd69e9OsWTOMMfz3v/9l+PDhtGnTpso8QUFBJCQk0LNnT5uqlIZgNn+B59kHoeggtArHNekOrKTT7C5LpFYKdSIi1UhKSiIpKQmA0tJSBg8eTFRUlM1VSUMyxmDeWIFJex6MBzon4Zp0J1ZUm9oXFvEBCnUiIrW49NJL7S5BGpgpOYxn4aPw2QcAWAOGYF0xESso2ObKROpOoU5ERJo0k70Xz+P3wd5MCAjEunwC1u+Gafyc+B2FOhERabLMpk/wzH0IDhdBeBSum+7E6tLD7rJETohCnYiINDnG48GsXoZ5ZTEYA116eMfPRWjcpPgvhToREWlSzOFiPPMfhi8/BsAa9Aesy27ACgyytzCRk6RQJyIiTYbZtxvPE/dB1m4IDMS6YhKu84baXZZIvVCoExGRJsF8+RGeeQ9DyWGIaO293FfnJLvLEqk3CnUiIuJoxuPB88pizGtLvRO69fSeUDgs0t7CROqZQp2IiDiW59BBPHP+D7PpUwCsC0ZiXXo9VqD+/Ynz6FUtIiKOZPbsZP/TszF7d0FQMNZVKbjOvcDuskQajEKdiIg4jvnsfTwLHoXSEohq6x0/17Gr3WWJNCiFOhERcQzjqcC89Dzm9RUAhPQ6E/d1f4GWYfYWJtIIFOpERBpJWloaGzduZM+ePQQHB5OUlMRVV11F+/btK+cxxrBs2TLWr1/PoUOH6NatG+PHjyc+Pt7Gyv2DKTqI55kHIf0LAKyhF9H25jvJys7BGGNzdSINz2V3ASIiTUV6ejrDhg3jH//4B3fffTcej4eZM2dSUlJSOc/LL7/MqlWruP7665k1axYRERHMnDmTw4cP21i57zO7tuOZ+VdvoAsOxrrhbwSMGY8VoH0X0nQo1ImINJLU1FQGDRpEfHw8nTp1IiUlhdzcXDIyMgDvXrrVq1eTnJxM//79SUhIYPLkyZSWlrJhwwabq/ddno3/w/PPqZC7H9q0w3XnA7j6D7S7LJFGp7cwIiI2KS4uBqBly5YAZGdnk5+fT+/evSvnCQoKomfPnmzZsoUhQ4ZU+zjl5eWUl5dX3rYsi9DQ0Mqfj+XI/bXN54tMRQWeFQsxa18CwDq1L64bb8dq2cp72497qwv1598aoj+FOhERGxhjWLRoET169CAhIQGA/Px8AMLDw6vMGx4eTm5ubo2PlZaWxvLlyytvd+7cmdmzZ9O2bds61xMTE3Mc1duvoiCfH2dPo/SrTwBodek4wq++CSsg4Kh5/a2346X+/Ft99qdQJyJig3nz5pGZmcm999571H2/fude2yD/5ORkRo4cedTyOTk5uN3uYy5rWRYxMTFkZWX5zcEEZuc2Kh6/D/JyIKQZrutu5fCZAzicnV1lPn/s7XioP/92PP0FBgbW6U2a34S67OxsVqxYwTfffEN+fj5RUVGcd955XHzxxQTqzOAi4kfmz5/PZ599xowZM2jdunXl9IiICMC7xy4y8udLWBUWFh619+6XgoKCCAoKqva+uv4zNMb4xT9Oz4dvY/7zOJSXQXQsrpS7sDp0PGbt/tLbiVJ//q0++/ObNLR3716MMUyYMIGYmBh27drF008/TUlJCddcc43d5YmI1MoYw/z589m4cSPTp08nOjq6yv3R0dFERESwadMmOnfuDIDb7SY9PZ0rr7zSjpJ9hnG7McsXYNa/6p1w+pm4bvgrVvOW9hYm4kP8JtT16dOHPn36VN5u164de/fuZe3atQp1IuIX5s2bx4YNG5g6dSqhoaGVY+iaN29OcHAwlmUxYsQI0tLSiI2NJSYmhrS0NEJCQhgwYIC9xdvIFB7A8/T9sHUzANbIy7D+eDmWSydwEPklvwl11SkuLq48aqwmTfWosLpQf/7Lyb2Bc/tbu3YtANOnT68yPSUlhUGDBgEwatQoysrKmDt3LkVFRXTt2pXU1NTK7VZTY7ZvxfPELMj/EZqF4rp+Clbfs+0uS8Qn+W2oy8rK4vXXX691L11TPCrseKk//+Xk3sB5/b344ou1zmNZFmPGjGHMmDGNUJFv87y3FrP4KXC7IaYDrpRUrNg4u8sS8Vm2h7oXX3yxSuiqzqxZs+jSpUvl7by8PO677z7OOeccBg8efMxlm9pRYcdD/fkvJ/cGDXNUmPgP4y7HvPAs5t03vBP69PfuoQttbm9hIj7O9lA3fPhwfvvb3x5znl9usPPy8pgxYwZJSUlMmDCh1sdvSkeFnSj157+c3Bs4vz85msnPw/PUP+GH78CyvGPnLhyj8XMidWB7qAsLCyMsLKxO8x4JdJ07dyYlJQWX/shFRBzDbPsWz1OzoSAPQlt4j27t9Ru7yxLxG7aHurrKy8tj+vTptGnThmuuuYbCwsLK+46c20lERPyPMQbz7huYF56FCjfExuOanIrVrr3dpYn4Fb8JdZs2bSIrK4usrCwmTZpU5b66DD4WERHfY8rLMIufxmx40zvhjHNxjbsVq1nTPNpX5GT4TagbNGhQ5SH/IiLi/0xernf83PatYLmwkq/GGn6x405lI9JY/CbUiYiIc5it33jHzx0sgOYtcU24HevUvnaXJeLXFOpERKTRGGMwb63CLJsHFRUQ18l7/da2zjonoYgdFOpERKRRmLJSzPNPYD58GwDrrN9hXXMzVkgzmysTcQaFOhERaXDmx2zv5b4yf/COnxs9DmvIKI2fE6lHCnUiItKgzLdf4Xnmfjh0EFqGecfPndLb7rJEHEehTkREGoQxBvPmy5jlC8F4IKELrpRpWK2j7S5NxJEU6kREpN6Z0hLMoscwn7wHgHXO+VhXpWAFh9hcmYhzKdSJiEi9MjlZeJ64D3bvgIAArEvHY11wocbPiTQwhToREak3ZvMXeJ55AIoPQatwXJPuwEo6ze6yRJoEhToRETlpxhjMGysxaf/xjp/rnIRr0p1YUW3sLk2kyVCoExGRk2JKDuNZ+Ch89gEA1oAhWFdMxAoKtrkykaZFoU5ERE6Y2b/XO35ubyYEBGJdPgHrd8M0fk7EBgp1IiJyQsymT/DMfQgOF0F4FK6b7sTq0sPuskSaLIU6ERE5LsbjwaxehnllMRgDXXp4x89FRNldmkiTplAnIiJ1Zg4X45n/CHz5EQDWoD9gXXYDVmCQrXWJiEKdiIjUkdm32zt+Lms3BAZiXXkTrgFD7C5LRH6iUCciIrUyX36MZ/7DcLgYIlp7L/fVOcnuskTkFxTqRESkRsbjwby6BPPaUu+EpFNxTZyKFRZpb2EichSFOhERqZYpPuQ9uvXrTwGwBv8Ra/R1WIH61yHii/SXKSIiRzF7MvE88Q/I3gdBwVhXT8Z1zvl2lyUix6BQJyIiVZjP3sez4FEoLYGotrhS7sLq2MXuskSkFgp1IiICgPFUYF56HvP6Cu+EHr1wTbgdq1W4vYWJSJ0o1ImI+KA1a9bwyiuvkJ+fT1xcHOPGjeOUU05psOczhw7ieeYBSP8CAGvoRVgXX4sVENBgzyki9ctldwEiIlLVBx98wMKFC7n44ouZPXs2p5xyCvfddx+5ubkN8nxlGVupmDnFG+iCg7FuvA3Xpdcr0In4GYU6EREf89prr3HBBRcwePDgyr10bdq0Ye3atfX+XJ6P3yX7tusgdz+0jcE17QFcZ/2u3p9HRBqePn4VEfEhbrebjIwMLrrooirTe/XqxZYtW6pdpry8nPLy8srblmURGhpa+XNNPGtW4lm2wDvfaf1w3XAbVstWJ9mB7zjS+7F+B/5M/fm3huhPoU5ExIcUFhbi8XgID696cEJ4eDj5+fnVLpOWlsby5csrb3fu3JnZs2fTtm3bYz5X2aBhZL+2lJZ/vIzwqyY59uPWmJgYu0toUOrPv9Vnfwp1IiI+qLp37zW9o09OTmbkyJFHzZeTk4Pb7a75SUJaEDDzKSJ69CQrKwtjzMkV7WMsyyImJsaRvYH683fH019gYGCtb9JAoU5ExKeEhYXhcrmO2itXUFBw1N67I4KCgggKCqr2vtr+WVjhkZXzOfEfJzi7N1B//q4++9OBEiIiPiQwMJDExEQ2bdpUZfqmTZvo3r27TVWJiD/QnjoRER8zcuRIHnvsMRITE0lKSmLdunXk5uYyZMgQu0sTER+mUCci4mPOPfdcDh48yIoVKzhw4ADx8fFMmzatTmNqRKTpUqgTEfFBw4YNY9iwYXaXISJ+RGPqRERERBygye6pCwyse+vHM68/Un/+y8m9Qd36c/rv4GRoO+fl5N5A/fm7+tzOWcbJxwmLiIiINBH6+PUYDh8+zB133MHhw4ftLqVBqD//5eTewPn9+RIn/66d3BuoP3/XEP0p1B2DMYbt27c79qSH6s9/Obk3cH5/vsTJv2sn9wbqz981RH8KdSIiIiIOoFAnIiIi4gAKdccQFBTE6NGja7ymor9Tf/7Lyb2B8/vzJU7+XTu5N1B//q4h+tPRryIiIiIOoD11IiIiIg6gUCciIiLiAAp1IiIiIg6gUCciIiLiAM6+oFodrFmzhldeeYX8/Hzi4uIYN24cp5xySo3zp6ens2jRInbv3k1kZCR/+tOfGDp0aCNWXDdpaWls3LiRPXv2EBwcTFJSEldddRXt27evcZnNmzczY8aMo6Y//PDDdOjQoSHLPW4vvvgiy5cvrzItPDycZ599tsZl/GXdTZ48mZycnKOmDx06lBtuuOGo6b6+3tLT03nllVfYvn07Bw4c4LbbbuOss86qvN8Yw7Jly1i/fj2HDh2iW7dujB8/nvj4+GM+7kcffcTSpUvZv38/7dq14/LLL6/yuFK7493++aqGeo35grpsy/25v7Vr17J27drKbV5cXByjR4+mb9++gH/3Vp20tDSWLFnCiBEjGDduHFDPPZom7P333zdjx44169atM7t27TILFiwwV111lcnJyal2/v3795urrrrKLFiwwOzatcusW7fOjB071nz44YeNXHntZs6cad5++22TmZlptm/fbmbNmmVuuukmc/jw4RqX+eabb8yll15q9uzZYw4cOFD5VVFR0YiV183SpUvNX//61yp1FhQU1Di/P627goKCKn199dVX5tJLLzXffPNNtfP7+nr7/PPPzZIlS8xHH31kLr30UvPxxx9XuT8tLc1cc8015qOPPjI7d+40Dz/8sJkwYYIpLi6u8TG3bNliLrvsMrNy5Uqze/dus3LlSjN27FizdevWhm7HMY53++fLGuI15ivqsi335/4++eQT89lnn5k9e/aYPXv2mMWLF5uxY8eazMxMY4x/9/Zr33//vUlJSTG33XabWbBgQeX0+uyxSX/8+tprr3HBBRcwePDgynepbdq0Ye3atdXOv3btWtq0acO4ceOIi4tj8ODBnH/++bz66quNXHntUlNTGTRoEPHx8XTq1ImUlBRyc3PJyMioddnw8HAiIiIqv1wu33yZuFyuKnWGhYXVOK8/rbuwsLAqfX3++ee0a9eOnj17HnM5X11vffv2ZezYsfTv3/+o+4wxrF69muTkZPr3709CQgKTJ0+mtLSUDRs21PiYq1atolevXiQnJ9OhQweSk5M57bTTWLVqVUO24ijHu/3zZQ3xGvMVtW3L/b2/M888k379+tG+fXvat2/P5ZdfTrNmzfj+++/9vrdfKikp4bHHHmPixIm0aNGicnp99+gbW30buN1uMjIy6N27d5XpvXr1YsuWLdUu8/3339OrV68q0/r06UNGRgZut7vBaq0PxcXFALRs2bLWeadOncqECRO49957+eabbxq6tBOWlZXFxIkTmTx5Mo888gj79++vcV5/XXdut5v33nuP888/H8uyjjmvv6y3X8rOziY/P7/K32FQUBA9e/as8e8QYOvWrUetz969e7N169YGq9VJTmT7569O9DXmq369LXdSfx6Ph/fff5/S0lKSkpIc1dvcuXPp27fvUdut+u6xyY6pKywsxOPxEB4eXmV6eHg4+fn51S6Tn59f7fwVFRUcPHiQyMjIhir3pBhjWLRoET169CAhIaHG+SIjI5kwYQKJiYm43W7+97//8X//93/cc889te4lamzdunVj8uTJtG/fnvz8fFauXMndd9/NQw89RKtWrY6a31/X3caNGykqKmLQoEE1zuNP6+3XjvytVbducnNzj7lcRERElWkRERE1/u1KVSey/fNXJ/oa80XVbcud0F9mZiapqamUl5fTrFkzbrvtNuLi4ipDjT/3BvD++++zfft2Zs2addR99b3+mmyoO6K6vR/H2iPy6/vMTxfkqG0vip3mzZtHZmYm99577zHnO7L7+4ikpCRyc3N59dVXfS4cHBlEC5CQkEBSUhK33HIL7777LiNHjqx2GX9cd2+//TZ9+vQhKiqqxnn8ab3VpKZ1czyMMT69Ln3R8W7//Fl9vMbsdqxtuT/31759ex544AGKior4+OOPefzxx6sc/OXPveXm5rJw4UJSU1MJDg6ucb766rHJfvwaFhaGy+U66l1pQUHBUYn5iOr2BBQWFhIQEFCnjzXtMH/+fD777DPuueceWrdufdzLJyUlkZWV1QCV1a9mzZqRkJDAvn37qr3fH9ddTk4OmzZtYvDgwce9rL+styN726pbNzX9HR5Z7nj+dqWqE9n++asTfY35mpq25U7oLzAwkJiYGLp06cIVV1xBp06dWL16tSN6y8jIoKCggDvvvJOxY8cyduxY0tPTef311xk7dmxlH/XVY5MNdYGBgSQmJrJp06Yq0zdt2kT37t2rXaZbt25Hzf/VV1+RmJhIYKBv7fQ0xjBv3jw+/vhj/v73vxMdHX1Cj7N9+/ajPubyReXl5ezZs6fGj1H9ad0d8fbbbxMeHk6/fv2Oe1l/WW/R0dFERERUWTdut5v09PQa/w7BG1q//vrrKtM2bdpEUlJSg9XqJCey/fNXJ/oa8xW1bcv9vb/qGGMoLy93RG+nn346Dz74IPfff3/lV5cuXRgwYAD3338/7dq1q9ceffO/WSMZOXIkjz32GImJiSQlJbFu3Tpyc3MZMmQIAIsXLyYvL4+bb74Z8J4nbM2aNSxatIjBgwezdetW3nrrLW699VY726jWvHnz2LBhA1OnTiU0NLTyXUDz5s0rdwH/ur9Vq1bRtm1b4uPjKwfof/zxx/ztb3+zq40aPffcc5x55pm0adOGgoICVqxYweHDhxk4cCDg3+sOvAOG33nnHQYOHEhAQECV+/xtvZWUlFTZa5idnc2OHTto2bIlbdq0YcSIEaSlpREbG0tMTAxpaWmEhIQwYMCAymXmzJlDVFQUV1xxBQAjRozgnnvu4aWXXuI3v/kNn3zyCV9//XWtQwzkZ7Vt//xJfbzGfFVt23LLsvy6v8WLF9O3b19at25NSUkJ77//Pps3byY1NdXvewMIDQ09aix7SEgIrVq1qpxenz026VB37rnncvDgQVasWMGBAweIj49n2rRptG3bFoADBw5UGagYHR3NtGnTWLRoEWvWrCEyMpLrrruOs88+264WanTktATTp0+vMj0lJaVy0P2v+3O73fznP/8hLy+P4OBg4uPjufPOO09oT1FDy8vL49FHH6WwsJCwsDC6devGP/7xD0esO4Cvv/6a3Nxczj///KPu87f19sMPP1QZH/Pcc88BMHDgQCZPnsyoUaMoKytj7ty5FBUV0bVrV1JTUwkNDa1cJjc3t8qYk+7du/OXv/yFF154gaVLlxITE8Nf/vIXunXr1niN+bnatn/+pD5eY76qLttyf+6voKCAOXPmcODAAZo3b07Hjh1JTU2tPErUn3urq/rs0TL+NOJQRERERKrVZMfUiYiIiDiJQp2IiIiIAyjUiYiIiDiAQp2IiIiIAyjUiYiIiDiAQp2IiIiIAyjUiYiIiDiAQp2IiIiIAyjUiYiIiDiAQp2IiIiIAyjUiYiIiDiAQp04RllZGVOnTuWWW26huLi4cnp+fj433ngj06dPx+Px2FihiIhIw1GoE8cIDg5mypQpFBYW8sQTTwDg8Xj497//DcCtt96Ky6WXvIiIOJP+w4mjxMbGMnHiRDZu3Mjq1atZvnw5mzdv5pZbbiEyMtLu8kRERBpMoN0FiNS3c889l/T0dP7zn//g8XhITk6mV69edpclIiLSoLSnThzp/PPPp6KigoCAAEaMGGF3OSIiIg1OoU4cp6SkhDlz5hAbG0twcDBPPfWU3SWJiIg0OIU6cZxnn32W3NxcbrvtNiZNmsSnn37Ka6+9ZndZIiIiDUqhThxl/fr1vPfee4wfP574+HjOPvtshg8fzn//+1+2bdtmd3kiIiINRqFOHCMzM5MFCxYwcOBABg0aVDn96quvpmPHjjz88MMUFRXZV6CIiEgDsowxxu4iREREROTkaE+diIiIiAMo1ImIiIg4gEKdiIiIiAMo1ImIiIg4gEKdiIiIiAMo1ImIiIg4gEKdiIiIiAMo1ImIiIg4gEKdiIiIiAMo1ImIiIg4gEKdiIiIiAP8f322cq50LRmLAAAAAElFTkSuQmCC"
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"# plot trajectory\n",
|
|
"plt.subplot(2, 2, 1)\n",
|
|
"plt.plot(x_bar[0, :], x_bar[1, :])\n",
|
|
"plt.plot(np.linspace(0, 10, T + 1), np.zeros(T + 1), \"b-\")\n",
|
|
"plt.axis(\"equal\")\n",
|
|
"plt.ylabel(\"y\")\n",
|
|
"plt.xlabel(\"x\")\n",
|
|
"\n",
|
|
"plt.subplot(2, 2, 2)\n",
|
|
"plt.plot(np.degrees(x_bar[2, :]))\n",
|
|
"plt.ylabel(\"theta(t)\")\n",
|
|
"# plt.xlabel('time')\n",
|
|
"\n",
|
|
"\n",
|
|
"plt.tight_layout()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The results are the same as expected, so the linearized model is equivalent as expected."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"name": "python3",
|
|
"language": "python",
|
|
"display_name": "Python 3 (ipykernel)"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.8.5"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|