Sympy_Laplace.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Laplace Transforms and Sympy\n",
    "\n",
    "Laplace transforms convert a function $f(t)$ in the time domain into function in the Laplace domain $F(s)$.\n",
    "\n",
    "$F(s) = \\mathcal{L}\\left(f(t)\\right) = \\int_0^\\infty f(t)e^{-s\\,t}dt$\n",
    "\n",
    "As an example of the Laplace transform, consider a constant $c$. The function $f(t) = c$ and the following expression is integrated. \n",
    "\n",
    "$\\mathcal{L}(c)=\\int_0^\\infty c \\, e^{-s\\,t} dt = -\\frac{c}{s}e^{-s\\,t} \\biggr\\rvert_0^\\infty = 0 - \\left(-\\frac{c}{s} \\right) = \\frac{c}{s}$\n",
    "\n",
    "Mathematicians have developed tables of commonly used Laplace transforms. <a href='https://apmonitor.com/pdc/index.php/Main/LaplaceTransforms'>Here</a> is summary table with a few of the entries that will be most common for analysis of linear differential equations in the <a href='https://apmonitor.com/pdc/'>Process Dynamics and Control course</a>. Notice that the derived value for a constant $c$ is the unit step function with $c=1$ where a signal output changes from 0 to 1 at $time=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sym\n",
    "from sympy.abc import s,t,x,y,z\n",
    "from sympy.integrals import laplace_transform\n",
    "from sympy.integrals import inverse_laplace_transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{5}{s^{2}}$"
      ],
      "text/plain": [
       "5/s**2"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Laplace transform (t->s)\n",
    "U = laplace_transform(5*t, t, s)\n",
    "U[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 5 t \\theta\\left(t\\right)$"
      ],
      "text/plain": [
       "5*t*Heaviside(t)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Inverse Laplace transform (s->t)\n",
    "inverse_laplace_transform(U[0],s,t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{5 s + 5}{\\left(s + 3\\right)^{2}}$"
      ],
      "text/plain": [
       "(5*s + 5)/(s + 3)**2"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Function\n",
    "F = 5*(s+1)/(s+3)**2\n",
    "F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{5}{s + 3} - \\frac{10}{\\left(s + 3\\right)^{2}}$"
      ],
      "text/plain": [
       "5/(s + 3) - 10/(s + 3)**2"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Partial fraction decomposition\n",
    "sym.apart(F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(s + 1\\right) \\left(s + 3\\right) \\left(s^{2} + 3 s + 1\\right)$"
      ],
      "text/plain": [
       "(s + 1)*(s + 3)*(s**2 + 3*s + 1)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# denominator of transfer function\n",
    "d1 = (s+1)*(s+3)*(s**2+3*s+1)\n",
    "d1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle s^{4} + 7 s^{3} + 16 s^{2} + 13 s + 3$"
      ],
      "text/plain": [
       "s**4 + 7*s**3 + 16*s**2 + 13*s + 3"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# expand polynomial\n",
    "d2=sym.expand(d1)\n",
    "d2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{-1: 1, -3: 1, -3/2 - sqrt(5)/2: 1, -3/2 + sqrt(5)/2: 1}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# find roots\n",
    "sym.roots(d2)"
   ]
  },
  {
   "attachments": {
    "input_laplace_function.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAApIAAAFxCAYAAADAsEjXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XtcVHX+P/DXmQsDKKCA3BQBTfOuhIViqJVgVK5mlunmpdW+9dW+m/J1/cVafdFvq1ttLrvrarteoiyV7Wtqu1ky7oagoCmJ611LEFQQr4yADHP7/QEzOgIqwwznzJnX8/HgUXPmzIf3mxk4bz+f8/l8BIvFYgERERERUSspxA6AiIiIiNwTC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyCAtJIiIiInIIC0kiIiIicggLSSIiIiJyiErsAOTAbDbjwoUL8PPzgyAIYodDRERELmSxWHDjxg1ERERAofDsPjkWkk5w4cIFREZGih0GERERtaOysjJ069ZN7DBExULSCfz8/AA0fKD8/f2d2rbBYEB2djaSk5OhVqud2raUME95YZ7ywjzlhXm2nU6nQ2RkpO3678lYSDqBdTjb39/fJYWkr68v/P39Zf8Lzzzlg3nKC/OUF+bpPLydjZNtiIiIiMhBLCSJiIiIyCEsJImIiIjIIbxHkoiIyAOZTCaYTCaxw3AZg8EAlUqFurq6VuepVquhVCpdFJm8sJAkIiLyIBaLBX5+fjhz5oysJ4tYLBaEhYWhrKzMoTw7deqEsLAwWf+MnIGFJBERkQeprKxE586d0aVLF3Ts2FG2hZLZbEZ1dTU6duzYqkXDLRYLamtrUVlZCQAIDw93VYiy4Fb3SC5btgwPP/ww/Pz8EBISggkTJuDkyZP3fN3mzZvRr18/aDQa9OvXD1u2bLF73mKxID09HREREfDx8cHo0aNx9OhRV6VBREQkCpPJBJ1Oh+DgYAQFBcHHxwfe3t6y/fLy8mr1a3x8fBAUFISQkBBcv35d1sP/zuBWheSuXbswd+5c7N27F1qtFkajEcnJyaipqWnxNQUFBZg8eTKmTZuGQ4cOYdq0aXjhhRewb98+2znvv/8+li9fjhUrVmD//v0ICwtDUlISbty40R5pERERtQuDwQAA8PLyEjkS6fP19QVw62dGzXOroe1vv/3W7vHHH3+MkJAQFBYWYuTIkc2+JiMjA0lJSUhLSwMApKWlYdeuXcjIyMDGjRthsViQkZGBRYsWYeLEiQCATz75BKGhodiwYQNeffVV1yZFRETUzuQ6nO1M/BndH7cqJO9UVVUFAAgMDGzxnIKCAsyfP9/u2NixY5GRkQEAKC4uRkVFBZKTk23PazQajBo1Cvn5+c0Wknq9Hnq93vZYp9MBaPhXi7P/5WJtT+7/Iro9z+2HK/Cn736CyWwROSrns8CC2holfn8qDwLk+0dKEID4AAFJHvS5lTPmKR8GgwEWS8PfVovFArPZLHJErtPWPM1mMywWCwwGQ5MZ3HL+jLSW2xaSFosFqampePTRRzFgwIAWz6uoqEBoaKjdsdDQUFRUVNietx6785yzZ8822+ayZcuwePHiJsezs7NtXeHOptVqXdKu1Gi1Wqw4qsCPOre666KVBKDupthBuNz5awoM3K5FB/nuwGbjSb+fnkDOeapUKoSFhQGALG/fys3NxYIFC7B3717bBJs789Tr9YiLi8Nnn32GIUOGtNhWfX09bt68idzcXBiNRrvnamtrnR+8m3LbQvL111/Hv//9b+zevfue597ZPW2xWJocu59zrNLS0pCammp7bN28PTk52SV7bWu1WiQlJcl+T1Rrnp+ePwjormPeEw/gkejOYofmVEajEQcO7MfQoQ9DpXLbX797WvyP4zh5sRpnvHrg/6X0ETscl/HE30/m6d7q6upQWloKAPDz85Pd8O2SJUuwaNEidOrUCRaLBYsWLcK3336LH374we68X/3qV3j33XeRnZ3dYlt1dXXw8fHByJEj4e3tbfecdSSS3LSQ/K//+i989dVXyM3NRbdu3e56blhYmK3X0aqystLWA2n9l1lFRYXdFP/bz7mTRqOBRqNpclytVrvsj48r25YStVoNY+OIdv+unZDQK0TcgJzMYDDg2ilg+ANdZP1+zh9jwGufF+Hz/efw6mO90cWv6e+LnHjS7yfzdG8mk8lWPAqC0KplcaQuPz8fp0+fxuTJk6FQKOyGs+/M86WXXsLChQtx8uRJ9O3bt9n2FAoFBEFo9vMg18+HI9zqE2SxWPD666/jyy+/xL/+9S/ExMTc8zXDhw9vMkyRnZ2NhIQEAEBMTAzCwsLszqmvr8euXbts51D7MjX+8quU8vqXsid5/MEuiOpowU2DGR/t+knscIhIBqKjo23zG6yGDBmC9PR0AMCmTZuQnJxs6z3MzMzEe++9h0OHDkEQBAiCgMzMTABAUFAQEhISsHHjxvZMQZbcqkdy7ty52LBhA7Zt2wY/Pz9bT2NAQAB8fHwAANOnT0fXrl2xbNkyAMAbb7yBkSNH4r333sP48eOxbds27Ny50zYkLggC5s2bh6VLl6JXr17o1asXli5dCl9fX0ydOlWcRD2c0dTQJalSsJB0V4IgICXSjI+OK7F+71m8ktgDYQHe934hEbU7i8WCm4b2XyvRR6106tB6bm4upkyZYns8efJkHDx4EN999x127twJoKFesHrkkUeQl5fntO/vqdyqkFy1ahUAYPTo0XbHP/74Y8ycORMAUFpaateFnZCQgE2bNuGtt97C22+/jZ49eyIrKwvx8fG2cxYuXIibN29izpw5uHbtGuLj45GdnQ0/Pz+X50RNGUyNPZIyGnLxRH0CLBga1QkHzl7Hn7/7Ef87oeVJcUQknpsGE/q9s6Pdv++xJWPh6+W8MqSkpAQRERG2xz4+PujQoYPdBKPbde3aFSUlJU77/p7KrQpJ61T+u8nJyWlybNKkSZg0aVKLrxEEAenp6bbucRKXddkfNYe23ZogAPOeeAAvrTuATftL8eqoHujW2TWrGhAR3bx5s8mkmLvx8fHh7GsncKtCkjyDoXFoW8mhbbcXHxOIhJ5ByP/pClb860f89rlBYodERHfwUStxbMlYUb5vaygUiiYdSrev5xgcHIxr167dd3tXr15Fly5dWhUDNcVCkiTnVo8kh7bl4L+TeyN/VQG+KDyH10b1RHRwB7FDIqLbCILg1CFmV+nSpQvKy8ttj3U6HYqLi22PY2NjcezYMbvXqNXqFvfKPnLkCGJjY10TrAfhlZokx8hZ27ISFxWIUb27wGS24K95Z8QOh4jc1OOPP47169cjLy8PR44cwYwZM+x2nBk7dmyTtaW7d++O4uJiFBUV4fLly3a70uXl5dntakeOYSFJkmPgrG3Z+c/RPQEAmwvP4Uq1/h5nExE1lZaWhpEjR+KZZ57BU089hQkTJqBnz56251966SUcO3YMJ0+etB372c9+hrFjx+Kxxx5Dly5dbMv9FBQUoKqq6q7zJ+j+SL8vmzyOdWibs7blIz4mEAO7BuDw+Sp8trcUb4zpJXZIRORm/P39kZWVZXdsxowZtv/v3LkzXn/9dSxfvhx/+ctfADRsIPLFF180WZB8+fLl+NWvfmVbOpAcxys1SY51+R9OtpEPQRAwO7FhA4H1e0tQJ8KadUQkf4sWLUJUVFSL90UCDXttDx48GPPnz2/HyOSLhSRJjpGTbWTpqYHhiAjwxuXqemwrOi92OEQkQwEBAfj1r39td+/knTQaDd566y32RjoJr9QkKRaL5dbQNifbyIpaqcDMEdEAgDV5xfe1LiwREUkbC0mSFGtvJMDJNnL04iPd0VGjwunKauScuiR2OERE1EYsJElSrPtsA4CKQ9uy4++txuSHIwEAa7gUEBGR2+OVmiSFPZLy9/KIaCgVAvb8eAVHL1SJHQ6RR+KtJfdmblzTmO6Oy/+QpBhv+8VlISlP3Tr7ImVAGP7x73KszSvG8slDxA6JyGN4eXlBoVDg0qVL0Gg00Gg0EAR5/q01m82or69HXV1dk+V/7sZisaC+vh6XLl2CQqGAl5eXC6N0fywkSVJuH9rm8j/y9UpiD/zj3+X46tAFLHyyD8ICvMUOicgjKBQKdO/eHYWFhVAoFLItIoGGgvDmzZvw8fFxKE9fX1907969VUWoJ2IhSZJya+kfQdZ/4Dzd4MhOeCQ6EN+XXMUnBSX4f0/2ETskIo+hVqtx9epVDB06VNZ/Zw0GA3JzczFy5Eio1epWvVapVEKlUsn65+MsLCRJUqxD2+yNlL9ZiTH4vuQqPt97Fq8/9gA6aPjniKg9qVSqVhdY7kSpVMJoNMLb21vWeYqN/bUkKdY1JNUcSpC9MX1DER3kC12dEV8cKBM7HCIicgCv1iQphsZ7JJVcjFz2lAoBsx5t2DZx3Z4S2z8iiIjIfbCQJEmxTrZRsUfSI0yKi0QnXzVKr9Yi+2iF2OEQEVEr8WpNkmK6bbINyZ+PlxIvxUcBANbsLhY5GiIiai0WkiQpBk628TjTE6LgpVSg8Ow1/FB6TexwiIioFVhIkqRYh7bV3B7RY4T4eWP8kAgA3DaRiMjduNXVOjc3F+PGjUNERAQEQcDWrVvvev7MmTMhCEKTr/79+9vOSU9Pb/J8WFiYq1OhFliHtrmrjWeZndgDAPDtkQqUXa0VORoiIrpfblVI1tTUYPDgwVixYsV9nf+HP/wB5eXltq+ysjIEBgbi+eeftzuvf//+ducdPnzYFeHTfeDQtmd6MMwPI3t3gdkCrOW9kkREbsOtVgBOSUlBSkrKfZ8fEBCAgIAA2+OtW7fi2rVrePnll+3OU6lU7IWUCA5te65XEmOQe+oS/nagDPOTeiPAhwsIExFJnVsVkm21du1ajBkzBlFRUXbHT58+jYiICGg0GsTHx2Pp0qXo0aNHi+3o9Xro9XrbY51OB6BhOyaDweDUmK3tObtdqbHmp2/8r1Ihz5w97f1sTZ7xUQF4MLQjTl6sxmcFxfiPxBhXhec0fD/lhXnKiyvzlPvPrjUEi8XilqsAC4KALVu2YMKECfd1fnl5OSIjI7Fhwwa88MILtuPffPMNamtr0bt3b1y8eBHvvvsuTpw4gaNHjyIoKKjZttLT07F48eImxzds2ABfX1/HEiIAQNEVAR+fUqKHnwVvDDCJHQ61s32VAjb8pESAlwXvxJqgYsc0EUlQbW0tpk6diqqqKvj7+4sdjqg8ppBctmwZPvzwQ1y4cAFeXl4tnldTU4OePXti4cKFSE1Nbfac5nokIyMjcfnyZad/oAwGA7RaLZKSkmS9V6g1T0P4ICz48hiGxXTG+l88LHZYTudp72dr89QbzXjsw1xcqq7H754bYJvNLVV8P+WFecqLK/PU6XQIDg5mIQkPGdq2WCxYt24dpk2bdtciEgA6dOiAgQMH4vTp0y2eo9FooNFomhxXq9Uu+6V0ZdtSYhEauqDUKqWs8/WU97O1earVwMwRMfhgx0msyy/Fc0O7QxCkP/GK76e8ME95cUWenvBzu18eMXC0a9cu/Pjjj5g1a9Y9z9Xr9Th+/DjCw8PbITK6k8HE5X883c/ju8NHrcSxch0KfroidjhERHQXblVIVldXo6ioCEVFRQCA4uJiFBUVobS0FACQlpaG6dOnN3nd2rVrER8fjwEDBjR5bsGCBdi1axeKi4uxb98+TJo0CTqdDjNmzHBtMtQs2zqSnLXtsTr5euH5od0AAKu5QDkRkaS51dX6wIEDiI2NRWxsLAAgNTUVsbGxeOeddwA0TKixFpVWVVVV2Lx5c4u9kefOncOUKVPw4IMPYuLEifDy8sLevXubzOym9mFsXEeSPZKe7RcjYiAIwHcnL+HHyhtih0NERC1wq3skR48ejbvNDcrMzGxyLCAgALW1Le+UsWnTJmeERk5iG9pmj6RHiw7ugOR+odhx9CLW7i7GsomDxA6JiIiawas1SYp1aFvNHkmP90rjtombfziPy9X6e5xNRERiYCFJkmI0cYtEahAX1RmDIzuh3mjG+oKzYodDRETNYCFJkmLgZBtqJAgCXmnc3Wb93rOoM3CBeiIiqeHVmiTFNrStZI8kAU/2D0PXTj64WlOPL384L3Y4RER0BxaSJCnGxsk2HNomoKFn+hePNvRKrtl9BmazW27ERUQkWywkSVKsy/+oObRNjSY/HAk/jQpnLtXgu5OVYodDRES34dWaJMVo5s42ZK+jRoUp8d0BcIFyIiKpYSFJkmLkFonUjJkJ0VApBOw9cxVHzleJHQ4RETViIUmSYuSsbWpGRCcfPD0oHAB7JYmIpIRXa5IU2xaJnLVNd7AuUP6Pf5fjwvWbIkdDREQAC0mSGA5tU0sGdA3AsB6BMJkt+CS/ROxwiIgILCRJYm5NtuFHk5qa/WhDr+SGfaWo1htFjoaIiHi1JkmxbpHIoW1qzuN9QtCjSwfc0BuRtb9M7HCIiDweC0mSFPZI0t0oFAJmNS5Qvm53se0fHkREJA5erUlSbs3aZo8kNe+5h7ohsIMXzl+/iW+PVogdDhGRR2MhSZLCyTZ0L95qJV4aFgUAWJ1XDIuF2yYSEYmFhSRJyq3lf/jRpJZNHx4FL5UCh8qu48DZa2KHQ0TksXi1JkkxNQ5tq9kjSXcR3FGDibFdAQCrc7lAORGRWFhIkqQYGoe2lSwk6R6sk260xy+i5HKNyNEQEXkmFpIkKbYeSQ5t0z30CvXD6Ae7wGIB1u0pFjscIiKPxKs1SQrXkaTWsG6b+MWBc7heWy9yNEREnsetCsnc3FyMGzcOEREREAQBW7duvev5OTk5EAShydeJEyfsztu8eTP69esHjUaDfv36YcuWLa5Mg+7CYObQNt2/hJ5B6Bvuj5sGEz7fVyp2OEREHsetCsmamhoMHjwYK1asaNXrTp48ifLycttXr169bM8VFBRg8uTJmDZtGg4dOoRp06bhhRdewL59+5wdPt0HDm1TawiCgFcSG+6VzMwvgd5oEjkiIiLPohI7gNZISUlBSkpKq18XEhKCTp06NftcRkYGkpKSkJaWBgBIS0vDrl27kJGRgY0bN7YpXmo969A2eyTpfj0zKALvfXsCF3V6fFV0Ac8PjRQ7JCIij+ER3T6xsbEIDw/HE088ge+++87uuYKCAiQnJ9sdGzt2LPLz89szRGpksC3/4xEfTXICL5UCMxMaeiXX7uYC5URE7cmteiRbKzw8HH/9618RFxcHvV6P9evX44knnkBOTg5GjhwJAKioqEBoaKjd60JDQ1FR0fLWa3q9Hnq93vZYp9MBAAwGAwwGg1NzsLbn7HalxpqfdWcbWEyyzNnT3s/2yvP5h8Lxp3+dxomKG/juRAUSHwhul+/L91NemKe8uDJPuf/sWkOwuOk/3wVBwJYtWzBhwoRWvW7cuHEQBAFfffUVAMDLywuffPIJpkyZYjvn888/x6xZs1BXV9dsG+np6Vi8eHGT4xs2bICvr2+r4iF7afuVqDUKSBtsRBh/lNQKm4sVyK1QoE+AGf/Zzyx2OEQkY7W1tZg6dSqqqqrg7+8vdjiiknWPZHOGDRuGzz77zPY4LCysSe9jZWVlk17K26WlpSE1NdX2WKfTITIyEsnJyU7/QBkMBmi1WiQlJUGtVju1bSmx5ikoVABMePyxUYgO6iB2WE7nae9ne+Y54GotkjJ240SVAg/EjUDvUD+Xf0++n/LCPOXFlXlaRyLJAwvJgwcPIjw83PZ4+PDh0Gq1mD9/vu1YdnY2EhISWmxDo9FAo9E0Oa5Wq132S+nKtqXE1NhB7u3lJet8PeX9bM88e4YGYGz/MHxzpAKZBWX44PnB7fJ9Ab6fcsM85cUVeXrCz+1+uVUhWV1djR9//NH2uLi4GEVFRQgMDET37t2RlpaG8+fP49NPPwXQMCM7Ojoa/fv3R319PT777DNs3rwZmzdvtrXxxhtvYOTIkXjvvfcwfvx4bNu2DTt37sTu3bvbPT+6dY8kFyQnR8xO7IFvjlRgW9EF/OrJBxHi5y12SEREsuZWU2MPHDiA2NhYxMbGAgBSU1MRGxuLd955BwBQXl6O0tJbixLX19djwYIFGDRoEBITE7F79258/fXXmDhxou2chIQEbNq0CR9//DEGDRqEzMxMZGVlIT4+vn2TI1gsgLFx1raKs7bJAXFRnfFQ906oN5nxaf5ZscMhIpI9t+qRHD169F2X9sjMzLR7vHDhQixcuPCe7U6aNAmTJk1qa3jURrdPj1CzR5Ic9EpiD/zn5z/gs31nMeexnvD1cqs/c0REboXdPiQZptsqSS5ITo5K7h+GyEAfXK81YHPhObHDISKSNRaSJBn2PZL8aJJjlAoBvxhxa4Fy67abRETkfLxak2Tc3iOpYo8ktcELQyPh761CyZVa/PP4RbHDISKSLRaSJBmm2zqOOLRNbdFBo8LU+CgAwJq8YpGjISKSLxaSJBnWEUiVQoAgsJCktpmZEA2VQsD3JVdxqOy62OEQEckSC0mSDGuPJHsjyRnCArzxs8ERAIDVeWdEjoaISJ5YSJJkWAtJTrQhZ5md2AMA8M2RCpy7VityNERE8sMrNkmGbWiba0iSk/SL8MeIB4JgMluQuadE7HCIiGSHhSRJhum2eySJnGX2ow29kpv2l0FXZxA5GiIieWEhSZJxa7INP5bkPKN6d8EDIR1RrTci6/syscMhIpIVXrFJMkwc2iYXUCgEzH60YYHyj/cUw3D7gqVERNQmLCRJMji0Ta4yIbYrgjt64UJVHbYfLhc7HCIi2WAhSZJhtjQUkCrO2iYn81YrMW1YNICGBcotFm6bSETkDLxik2SwR5Jc6aVh3aFRKXD4fBX2FV8VOxwiIllgIUmSwXskyZWCOmrwXFw3AMAaLlBOROQULCRJMjhrm1xtVuOkm53HK3HmUrXI0RARuT9esUkybu1swx5Jco2eXTpiTN8QAMDa3cUiR0NE5P5YSJJkmLnXNrUD67aJ/1d4Dldr6kWOhojIvbGQJMngXtvUHuJjAjGwawD0RjM+23tW7HCIiNwar9gkGSb2SFI7EAQBsxMb7pX8tKAEdQaTuAEREbkxFpIkGZxsQ+3lqYHhCA/wxuXqemwrOi92OEREbotXbJIMTrah9qJWKvDyiGgAXKCciKgt3KqQzM3Nxbhx4xAREQFBELB169a7nv/ll18iKSkJXbp0gb+/P4YPH44dO3bYnZOeng5BEOy+wsLCXJkGtYBD29SeXnykOzpqVDhdWY1dpy6JHQ4RkVtyq0KypqYGgwcPxooVK+7r/NzcXCQlJWH79u0oLCzEY489hnHjxuHgwYN25/Xv3x/l5eW2r8OHD7sifLoHMyfbUDvy91Zj8sORABp6JYmIqPVUYgfQGikpKUhJSbnv8zMyMuweL126FNu2bcPf//53xMbG2o6rVCr2QkoAt0ik9vbyiGhk5pdg94+XceyCDv0i/MUOiYjIrbhVIdlWZrMZN27cQGBgoN3x06dPIyIiAhqNBvHx8Vi6dCl69OjRYjt6vR56vd72WKfTAQAMBgMMBoNTY7a25+x2pcZgMNh6JBWCRbb5etL7eft/pSq0oxpP9gvF10cqsDr3R7z/3MBWvd5d8mwr5ikvzNN5bRMgWNz0LnNBELBlyxZMmDDhvl/zwQcf4Le//S2OHz+OkJCG3S2++eYb1NbWonfv3rh48SLeffddnDhxAkePHkVQUFCz7aSnp2Px4sVNjm/YsAG+vr6OJUTYXqbAjnMKPBpqxvM9zGKHQx7ibDWw/LAKSsGC/3nIhAAvsSMiIqmrra3F1KlTUVVVBX9/zx7J8JhCcuPGjZg9eza2bduGMWPGtHheTU0NevbsiYULFyI1NbXZc5rrkYyMjMTly5ed/oEyGAzQarVISkqCWq12attSYjAY8Mu1/8TO8wrMGN4dbz3VR+yQXMKT3k93ynPKmu9x4Ox1vJoYgwXJve77de6Wp6OYp7wwz7bT6XQIDg5mIQkPGdrOysrCrFmz8MUXX9y1iASADh06YODAgTh9+nSL52g0Gmg0mibH1Wq1y34pXdm2VJgbOyG9VErZ5+oJ7yfgPnm+MrInDqwvxMb9ZfjlmN7ooGndn0Z3ybOtmKe8MM+2tUkNZD89duPGjZg5cyY2bNiAp59++p7n6/V6HD9+HOHh4e0QHd3Our+IirO2qZ2N6RuK6CBf6OqM+OJAmdjhEBG5Dbe6YldXV6OoqAhFRUUAgOLiYhQVFaG0tBQAkJaWhunTp9vO37hxI6ZPn44PP/wQw4YNQ0VFBSoqKlBVVWU7Z8GCBdi1axeKi4uxb98+TJo0CTqdDjNmzGjf5MjWI8lZ29TelAoBsx5t2DZx3Z4SmMxueccPEVG7c6tC8sCBA4iNjbUt3ZOamorY2Fi88847AIDy8nJbUQkAf/nLX2A0GjF37lyEh4fbvt544w3bOefOncOUKVPw4IMPYuLEifDy8sLevXsRFRXVvsnRbcv/uNXHkmTiubhu6OSrRunVWmiPVYgdDhGRW3CreyRHjx59163MMjMz7R7n5OTcs81Nmza1MSpyFlshyS0SSQS+Xir8PL47/vzdT1idV4wnB/D2FiKie2HXD0mGmQuSk8hmDI+Gl1KBwrPX8EPpNbHDISKSPBaSJBm3eiT5sSRxhPh742dDIgAAa/LOiBwNEZH08YpNknFrr232SJJ4Zic2TLr59kgFyq7WihwNEZG0sZAkybD2SCo5tE0i6hPmj8RewTBibyqDAAAgAElEQVRbgLW7i8UOh4hI0lhIkmRYC0k1Z22TyF5J7AEA+NuBMlTd5J66REQt4RWbJMPMWdskEYm9gtEnzA+19SZs/L703i8gIvJQLCRJMji0TVIhCLcWKM/cU4J6o1nkiIiIpImFJEmG2dJQQKo5a5sk4GdDItDFT4MKXR2+PnxB7HCIiCSJV2ySDBPXkSQJ0aiUmDG8YYer1bnFd90MgYjIU7GQJMngzjYkNT+Pj4K3WoFj5ToU/HRF7HCIiCSHhSRJhpl7bZPEdO7ghefjIgEAq7lAORFRE7xik2RwaJukaNajMRAE4LuTl/Bj5Q2xwyEikhQWkiQZ3CKRpCg6uAOS+oYC4ALlRER34hWbJIPrSJJUvTKyYYHyzT+cx+VqvcjREBFJBwtJkgwObZNUDY3qjMHdAlBvNGN9wVmxwyEikgwWkiQZJk62IYkSBAGzG7dNXL/3LOoMJpEjIiKSBl6xSTKsQ9tqDm2TBKUMCEPXTj64WlOPL384L3Y4RESSwEKSJINbJJKUqZQKvDwiGgCwZvcZmM1coJyIiIUkSYbJ1iPJjyVJ0+SHI+GnUeHMpRp8d7JS7HCIiETHKzZJBmdtk9T5easxJb47AGBNHpcCIiJiIUmSYLFYYLI0FJAc2iYpm5EQDaVCQMGZKzh6QSd2OEREonKrQjI3Nxfjxo1DREQEBEHA1q1b7/maXbt2IS4uDt7e3ujRowc++uijJuesXLkSMTEx8Pb2RlxcHPLy8lwRPt3F7bebqTlrmySsaycfPD0wHACwbg+XAiIiz+ZWV+yamhoMHjwYK1asuK/zi4uL8dRTTyExMREHDx7Er3/9a/zyl7/E5s2bbedkZWVh3rx5WLRoEQ4ePIjExESkpKSgtLTUVWlQM4wms+3/ObRNUvdK41JA249U4DrXJyciD6YSO4DWSElJQUpKyn2f/9FHH6F79+7IyMgAAPTt2xcHDhzA7373Ozz33HMAgOXLl2PWrFmYPXs2ACAjIwM7duzAqlWrsGzZMucnQc0y3NYlyXUkSeoGdgtAfEwg9hVfxa4KBaaKHRARkUjcqpBsrYKCAiQnJ9sdGzt2LNauXQuDwQCLxYLCwkK8+eabduckJycjPz+/xXb1ej30+lvdEDpdw31SBoMBBoPBiRnA1p6z25WaOn297f8tZiMMBvNdznZfnvJ+ekKeLyd0x77iqyi4KOB69U106ih2RK7jCe8nwDzlxpV5yv1n1xqyLiQrKioQGhpqdyw0NBRGoxGXL19umOBhMjV7TkVFRYvtLlu2DIsXL25yPDs7G76+vs4J/g5ardYl7UrFDQNg/Thmf/stBJmPbsv9/bSSc55mCxDirURlnYClWTkYHS7/dSXl/H7ejnnKiyvyrK2tdXqb7krWhSTQsLXZ7SwWi+347f9/5zl3HrtdWloaUlNTbY91Oh0iIyORnJwMf39/Z4UOoOFfPVqtFklJSVCr1U5tW0rOXakGDuRDKQh4+umnxA7HZTzl/fSUPK8HlWDx16ew71oHLJ35KFQyXQPVU95P5ikvrszTOhJJMi8kw8LCmvQsVlZWQqVSISgoCBaLBUqlstlz7uylvJ1Go4FGo2lyXK1Wu+yX0pVtS4GgVAJomGgj5zyt5P5+Wsk9z0lxkfjdjpO4UFWHf566gmcGRYgdkkvJ/f20Yp7y4oo8PeHndr/k+c/nRsOHD2/SpZ2dnY2hQ4dCrVbDy8sLcXFxTc7RarVISEhoz1A9nrFxWxsV15AkN+KtVuLR0IbP7uq8YtsoBxGRp3CrQrK6uhpFRUUoKioC0LC8T1FRkW2pnrS0NEyfPt12/muvvYazZ88iNTUVx48fx7p167B27VosWLDAdk5qairWrFmDdevW4fjx45g/fz5KS0vx2muvtW9yHs7QuPwPl/4hd/NomBleKgUOlV1H4dlrYodDRNSu3Gpo+8CBA3jsscdsj633Kc6YMQOZmZkoLy+3W/8xJiYG27dvx/z58/HnP/8ZERER+OMf/2hb+gcAJk+ejCtXrmDJkiUoLy/HgAEDsH37dkRFRbVfYgST2doj6Vb/tiGCvxcwfnA4vig8j9V5ZzA0OlDskIiI2o1bFZKjR4++69BRZmZmk2OjRo3CDz/8cNd258yZgzlz5rQ1PGoDo5lD2+S+Xk6IwheF55F97CJKLtcgOriD2CEREbULdv+QJNgKSQ5tkxvqFdIRox/sAosFWLenWOxwiIjaDQtJkgTrFokc2iZ3Zd028YsD53C9tv4eZxMRyQOv2iQJ1h5JJYe2yU0l9AxC33B/3DSY8Pm+0nu/gIhIBlhIkiRYC0k1h7bJTQmCgFcSYwAAn+SXoN4oz20+iYhux0KSJME6tM0eSXJnzwyKQKi/BpU39Pjq0AWxwyEicjkWkiQJtgXJ2SNJbsxLpcCMhGgAwJq8M1ygnIhkj4UkSYJtaJuTbcjN/fyRKPh6KXGi4gZ2/3hZ7HCIiFyKV22SBE62IbkI8FXjhaGRABq2TSQikjMWkiQJRm6RSDLyixExUAhA7qlLOFlxQ+xwiIhchoUkSQKHtklOugf5Ymz/MADA2t1nRI6GiMh1eNUmSeDQNsnN7MalgLYevIDKG3UiR0NE5BosJEkSOLRNchMXFYjY7p1QbzJjfcFZscMhInIJFpIkCRzaJjmybpv42d6zuFlvEjkaIiLn41WbJIFD2yRHY/uHITLQB9dqDfi/H86JHQ4RkdOxkCRJ4ILkJEdKhYBfjGi4V3Ld7mKYzVygnIjkhYUkSQL32ia5emFoJPy9VSi+XIOdxy+KHQ4RkVOxkCRJ4F7bJFcdNCpMjY8CAKzhAuVEJDMsJEkSTI09kipOtiEZmpkQDZVCwPclV3Go7LrY4RAROQ2v2iQJBlshyR5Jkp+wAG+MGxwBAFidxwXKiUg+WEiSJHAdSZI76wLl3xypwLlrtSJHQ0TkHCwkSRI4tE1y1z8iAAk9g2AyW5C5p0TscIiInMItr9orV65ETEwMvL29ERcXh7y8vBbPHT16NARBaPL19NNP286ZOXNmk+eHDRvWHqlQIw5tkyewLlC+aX8ZdHUGkaMhImo7tysks7KyMG/ePCxatAgHDx5EYmIiUlJSUFpa2uz5X375JcrLy21fR44cgVKpxPPPP2933pNPPml33vbt29sjHWrEdSTJE4zq3QUPhHREtd6IrO/LxA6HiKjN3K6QXL58OWbNmoXZs2ejb9++yMjIQGRkJFatWtXs+YGBgQgLC7N9abVa+Pr6NikkNRqN3XmBgYHtkQ41Mpl5jyTJn0IhYPajDfdKfrynGIbGe4OJiNyVWxWS9fX1KCwsRHJyst3x5ORk5Ofn31cba9euxYsvvogOHTrYHc/JyUFISAh69+6NV155BZWVlU6Lm+7NYOI9kuQZJsR2RXBHL1yoqsP2w+Vih0NE1CYqsQNojcuXL8NkMiE0NNTueGhoKCoqKu75+u+//x5HjhzB2rVr7Y6npKTg+eefR1RUFIqLi/H222/j8ccfR2FhITQaTZN29Ho99Hq97bFOpwMAGAwGGAzOve/J2p6z25WaeqMJACBYzLLO1VPeT+bZMiWAqY9E4o//+glr8s4gpV8XCIK0e+L5fsoL83Re2wQIFovFbTZ/vXDhArp27Yr8/HwMHz7cdvw3v/kN1q9fjxMnTtz19a+++iry8/Nx+PDhu55XXl6OqKgobNq0CRMnTmzyfHp6OhYvXtzk+IYNG+Dr63uf2dDt1pxQ4PA1BV7oYcKIULf5SBI5pNoApBcqYbAI+K/+RjzgL3ZERNQatbW1mDp1KqqqquDv79m/wG7VIxkcHAylUtmk97GysrJJL+WdamtrsWnTJixZsuSe3yc8PBxRUVE4ffp0s8+npaUhNTXV9lin0yEyMhLJyclO/0AZDAZotVokJSVBrVY7tW0p+fJyIXDtCgb274enHokSOxyX8ZT3k3ne22HhGDbtP4djpnD88qlYF0XoHHw/5YV5tp11JJLcrJD08vJCXFwctFotnn32WdtxrVaL8ePH3/W1f/vb36DX6/HSSy/d8/tcuXIFZWVlCA8Pb/Z5jUbT7JC3Wq122S+lK9uWgsbVf6CReZ5Wcn8/rZhny14Z2ROb9p/Dv05eQtl1PXp06eii6JyH76e8MM+2tUkN3G5mQ2pqKtasWYN169bh+PHjmD9/PkpLS/Haa68BAKZPn460tLQmr1u7di0mTJiAoKAgu+PV1dVYsGABCgoKUFJSgpycHIwbNw7BwcF2xSq5ltHM5X/Is/Ts0hFj+obAYgHW7i4WOxwiIoe4VY8kAEyePBlXrlzBkiVLUF5ejgEDBmD79u2IimoYDi0tLYXijpm/p06dwu7du5Gdnd2kPaVSicOHD+PTTz/F9evXER4ejsceewxZWVnw8/Nrl5wItmVQuCA5eZLZiT2w83gl/q/wHP47+UEEdvASOyQiolZxu0ISAObMmYM5c+Y0+1xOTk6TY71790ZLc4p8fHywY8cOZ4ZHDuAWieSJ4mMCMaCrP46c1+GzvWfxyyd6iR0SEVGr8KpNkmAd2lZyaJs8iCAItm0TPy0oQZ3BJG5AREStxEKSJMG6ILmaQ9vkYZ4aGI7wAG9crq7HV0UXxA6HiKhVWEiSJHCLRPJUaqUCL4+IBgCs2X2mxdtwiIikiIUkSYKxsUdSyR5J8kAvPtIdHTUqnLpYjV2nLokdDhHRfWMhSZJgMFuHtvmRJM/j763G5IcjAQBr8rgUEBG5D161SRJMXEeSPNzLI6KhEIDdP17GsQvcNYOI3AMLSZIEY+M6khzaJk/VrbMvUgY27Ka1ZvcZkaMhIro/LCRJEowc2iayLQX090MXcFFXJ3I0RET3xqs2SQK3SCQChkR2wsPRnWEwWfBJfonY4RAR3RMLSZIEDm0TNZjd2Cv5+b5S1NYbRY6GiOjuWEiSJNiGtpX8SJJnG9M3FNFBvqi6acAXB86JHQ4R0V3xqk2SYNsikT2S5OGUCgG/eDQGALB2d7FtRQMiIiliIUmiM5ktsG7moWIhSYRJcd0Q4KNG6dVaaI9ViB0OEVGLWEiS6IyN2yMCgJqTbYjg66XCS8O6AwBWc4FyIpIwFpIkOuv2iACHtomsZgyPhpdSgcKz1/BD6TWxwyEiahYLSRLd7YWkiutIEgEAQvy98bMhEQCAteyVJCKJ4lWbRHf70DbvkSS6ZXZiw6Sbb46Uo+xqrcjREBE1xUKSRGedsS3AAgULSSKbPmH+SOwVDLMFWLeHvZJEJD0sJEl0Buti5KwhiZqwLlD+t/1lqLppEDkaIiJ7LCRJdNZ18lhIEjU1slcwHgz1Q029CRu/LxU7HCIiOywkSXSGxsk2HNUmakoQBMxqvFcyc08J6o3me7yCiKj9uGUhuXLlSsTExMDb2xtxcXHIy8tr8dzMzEwIgtDkq66uzuE2ybnYI0l0d+OHRKCLnwYVujp8ffiC2OEQEdm4XSGZlZWFefPmYdGiRTh48CASExORkpKC0tKWh3z8/f1RXl5u9+Xt7d2mNsl5eI8k0d1pVErMGB4FAFidWwyLhdsmEpE0uF0huXz5csyaNQuzZ89G3759kZGRgcjISKxatarF1wiCgLCwMLuvtrZJzmOdtc2hbaKW/Tw+Ct5qBY6V61Bw5orY4RARAXCzQrK+vh6FhYVITk62O56cnIz8/PwWX1ddXY2oqCh069YNzzzzDA4ePNjmNsl5TGb2SBLdS+cOXng+LhIAsIYLlBORRKjEDqA1Ll++DJPJhNDQULvjoaGhqKioaPY1ffr0QWZmJgYOHAidToc//OEPGDFiBA4dOoRevXo51KZer4der7c91ul0AACDwQCDwbnLc1jbc3a7UnJT35CbQpB3noBnvJ8A83SV6cO64bN9Z/GvE5U4fv4aHgjp2C7fl++nvDBP57VNblZIWgmCfdeVxWJpcsxq2LBhGDZsmO3xiBEj8NBDD+FPf/oT/vjHPzrU5rJly7B48eImx7Ozs+Hr63vfebSGVqt1SbtScPK6AEAJpSDvPG/HPOWlPfMc0EmBw9cUWJK1Gy/2bN8Z3Hw/5YV5Oq62ljtNWblVIRkcHAylUtmkp7CysrJJj2JLFAoFHn74YZw+fdrhNtPS0pCammp7rNPpEBkZieTkZPj7+7cmpXsyGAzQarVISkqCWq12attS0fH0ZeD4D1AqIOs8Ac94PwHm6Upd+l3D1LX7UXhVhd+/nIigjhqXf0++n/LCPNvOOhJJblZIenl5IS4uDlqtFs8++6ztuFarxfjx4++rDYvFgqKiIgwcONDhNjUaDTSapn+81Wq1y34pXdm22CyNt+oqIO88b8c85aU98xz+QBcM7haAQ+eqsPHABcxP6t0u3xfg+yk3zLNtbVIDt5psAwCpqalYs2YN1q1bh+PHj2P+/PkoLS3Fa6+9BgCYPn060tLSbOcvXrwYO3bswJkzZ1BUVIRZs2ahqKjIdv79tEmuZbROtnG7TyNR+xMEwbZt4vq9Z1FnMIkcERF5MrfqkQSAyZMn48qVK1iyZAnKy8sxYMAAbN++HVFRDWuslZaWQqG4VZFcv34d//Ef/4GKigoEBAQgNjYWubm5eOSRR+67TXItIxckJ2qVlAFh6NrJB+ev38SXP5zH1PjuYodERB7K7QpJAJgzZw7mzJnT7HM5OTl2j3//+9/j97//fZvaJNcy2rZI5CLLRPdDpVTg5RHRePfr41i7+wxefDgSCi7ESkQi4GAiiY472xC13uSHI+GnUeGnSzXIOVUpdjhE5KFYSJLouNc2Uev5easxpXFIe3UuFygnInGwkCTRGbhFIpFDZiZEQ6kQUHDmCo6crxI7HCLyQCwkSXQmDm0TOSSikw+eHhgOAFiTd0bkaIjIE7GQJNEZ2SNJ5LBXGpcC+se/y1FedVPkaIjI07CQJNEZTLxHkshRA7sFID4mEEazBZl7SsQOh4g8DAtJEp3JzKFtoraw9kpu+L4U1XqjyNEQkSdhIUmiM5g4tE3UFo/3CUGPLh1wo86IrP1lYodDRB6EhSSJzsgeSaI2USgEzHo0BgDw8Z5iGBsnsBERuRoLSRIdJ9sQtd3E2G7o7KvGuWs3sePoRbHDISIPwUKSRGfkZBuiNvPxUmLasCgAwOq8M7BYuOUoEbkeC0kSnZHrSBI5xbTh0fBSKVBUdh2FZ6+JHQ4ReQAWkiS6W0Pb7EEhaosufho8O6QrgIZeSSIiV2MhSaLj0DaR88xObJh0k33sIkou14gcDRHJHQtJEp21R5KFJFHb9Qr1w+gHu8BiAdbtKRY7HCKSORaSJDrr8j+ctU3kHNYFyr84cA7Xa+tFjoaI5IyFJImOQ9tEzpXQMwh9w/1x02DC5/tKxQ6HiGSMhSSJjj2SRM4lCAJmNy5Q/kl+CfRGk8gREZFcsZAk0bFHksj5xg2OQKi/BpU39Pj7oXKxwyEimWIhSaIzcLINkdN5qRSYkRANAFjDBcqJyEVYSJLoTBzaJnKJnz8SBV8vJU5U3MDuHy+LHQ4RyRALSRKdgUPbRC4R4KvGC0MjAQBr8rgUEBE5n1sWkitXrkRMTAy8vb0RFxeHvLy8Fs9dvXo1EhMT0blzZ3Tu3BljxozB999/b3fOzJkzIQiC3dewYcNcnQY1sm6RyB5JIuf7xYgYKARg16lLOHXxhtjhEJHMuF0hmZWVhXnz5mHRokU4ePAgEhMTkZKSgtLS5pe4yMnJwZQpU/Ddd9+hoKAA3bt3R3JyMs6fP2933pNPPony8nLb1/bt29sjHQJgst4j6XafRiLp6x7ki7H9wwA03CtJRORMbnfpXr58OWbNmoXZs2ejb9++yMjIQGRkJFatWtXs+Z9//jnmzJmDIUOGoE+fPli9ejXMZjP++c9/2p2n0WgQFhZm+woMDGyPdAi3DW2LHAeRXM1uXKB868ELqLxRJ3I0RCQnKrEDaI36+noUFhbizTfftDuenJyM/Pz8+2qjtrYWBoOhSaGYk5ODkJAQdOrUCaNGjcJvfvMbhISENNuGXq+HXq+3PdbpdAAAg8EAg8HQmpTuydqes9uVEoOpYY07hULeeQKe8X4CzFNqBkV0xJDIABSVVSFzdzHmj3mgVa93lzzbinnKiyvzlPvPrjUEixutCXHhwgV07doVe/bsQUJCgu340qVL8cknn+DkyZP3bGPu3LnYsWMHjhw5Am9vbwANw+UdO3ZEVFQUiouL8fbbb8NoNKKwsBAajaZJG+np6Vi8eHGT4xs2bICvr28bMvRMvzmoRGWdgP/qb8QD/mJHQyRPRVcEfHxKiQ4qC9IfMsGLQwBEDqutrcXUqVNRVVUFf3/PvnC5VY+klSDYz8qwWCxNjjXn/fffx8aNG5GTk2MrIgFg8uTJtv8fMGAAhg4diqioKHz99deYOHFik3bS0tKQmppqe6zT6RAZGYnk5GSnf6AMBgO0Wi2SkpKgVqud2rZU/O5EHlB3E0oBss4T8Iz3E2CeUjTWbIE2YzfOXbuJ2tCBmPBI5H2/1p3ybAvmKS+uzNM6EkluVkgGBwdDqVSioqLC7nhlZSVCQ0Pv+trf/e53WLp0KXbu3IlBgwbd9dzw8HBERUXh9OnTzT6v0Wia7alUq9Uu+6V0Zdtia5xrA4Ug7zxvxzzlxR3yVAOY9WgMFv/9GDILSjFteAwUrVwqwR3ydAbmKS+uyNMTfm73y60m23h5eSEuLg5ardbuuFartRvqvtMHH3yA//3f/8W3336LoUOH3vP7XLlyBWVlZQgPD29zzHRvhsblf7iOJJFrvTA0Ev7eKhRfrsE/T1SKHQ4RyYBbFZIAkJqaijVr1mDdunU4fvw45s+fj9LSUrz22msAgOnTpyMtLc12/vvvv4+33noL69atQ3R0NCoqKlBRUYHq6moAQHV1NRYsWICCggKUlJQgJycH48aNQ3BwMJ599llRcvQ0xsYuSa4jSeRaHTQqTI2PAgCs5lJAROQEbjW0DTTcz3jlyhUsWbIE5eXlGDBgALZv346oqIY/jqWlpVAobtXHK1euRH19PSZNmmTXzv/8z/8gPT0dSqUShw8fxqefforr168jPDwcjz32GLKysuDn59euuXkqI3skidrNzIRorMk7g++Lr+Lf565jULdOYodERG7M7QpJAJgzZw7mzJnT7HM5OTl2j0tKSu7alo+PD3bs2OGkyMgR1h5JFpJErhcW4I1xgyOw5eB5rM4rxp+mxIodEhG5Mbcb2ib5MZo4tE3UnmYnxgAAth8ux/nrN0WOhojcGQtJEp3RzKFtovbUPyIACT2DYDJbkLmnWOxwiMiNsZAkUZnNFtvyPywkidrPK43bJm76vgw36rhLBxE5hoUkicrQ2BsJcGibqD2N6t0FD4R0xA29EVn7y8QOh4jcFAtJEpXJfGuHTvZIErUfhULArEcb7pX8eE+JbfUEIqLWYCFJojKYWEgSieXZ2K4I6uCF89dvYvuRinu/gIjoDiwkSVS394JwaJuofXmrlZg2vGEN3jV5Z2CxWO7xCiIieywkSVSm23a1YSFJ1P6mDYuCRqXAv89V4fviq2KHQ0RuhoUkicpgXYycVSSRKII6ajDxoW4AgNV5XAqIiFqHhSSJytR4j6RayY8ikVisk27+eeIizlyqFjkaInInvHqTqKzL/7BHkkg8D4R0xBN9QmCxAOu4QDkRtQILSRKVdXtEFQtJIlHNblyg/P8Kz+FaTb3I0RCRu2AhSaKybo/IoW0icQ3rEYgBXf1RZzDjs71nxQ6HiNwEr94kKmuPJIe2icQlCAJmP9rQK/lJwVnUGUwiR0RE7oCFJInK2iPJoW0i8T09KBzhAd64XK3HV0UXxA6HiNwAC0kSldE2a5uFJJHY1EoFZiZEAwDW7OYC5UR0bywkSVRGriNJJClT4rujo0aFUxersePoRbHDISKJYyFJojKYrEPb/CgSSYG/t9rWK/l77SmYzeyVJKKW8epNorJukcihbSLpeCWxB/y8VTh58Qb+cbhc7HCISMJYSJKoDJy1TSQ5Ab5qvNK4rmTGzlMwNo4cEBHdiYUkicraI6niOpJEkvLyiGh08lXjzKUa/P3fFWKHQ0QS5ZZX75UrVyImJgbe3t6Ii4tDXl7eXc/fvHkz+vXrB41Gg379+mHLli12z1ssFqSnpyMiIgI+Pj4YPXo0jh496soUqJFtQXL2SBJJip+3Gq+O7AkA+NN3P4GdkkTUHLcrJLOysjBv3jwsWrQIBw8eRGJiIlJSUlBaWtrs+QUFBZg8eTKmTZuGQ4cOYdq0aXjhhRewb98+2znvv/8+li9fjhUrVmD//v0ICwtDUlISbty40V5peSwObRNJ14yEKAR39ELZtZvYd4m/o0TUlNsVksuXL8esWbMwe/Zs9O3bFxkZGYiMjMSqVauaPT8jIwNJSUlIS0tDnz59kJaWhieeeAIZGRkAGnojMzIysGjRIkycOBEDBgzAJ598gtraWmzYsKE9U/NIJuuC5JxsQyQ5vl4q/OfoBwAAO84poDeyW5KI7KnEDqA16uvrUVhYiDfffNPueHJyMvLz85t9TUFBAebPn293bOzYsbZCsri4GBUVFUhOTrY9r9FoMGrUKOTn5+PVV19t0qZer4der7c91ul0AACDwQCDweBYci34qug8/lGiwMGvj0Mpw/sIj1c09PpaM3P2z09qrPkxT3nwhDwnPxSOv+z6EZU36jE/6xC6B/mKHZLLmExmnD0r37+3Vp6WZ2JtHfyc/LGV8+98a7lVIXn58mWYTCaEhobaHQ8NDUVFRfM3g1dUVNz1fOt/mzvn7Nmzzba5bNkyLF68uMnx7Oxs+Po699P6tx8V2HdJAZSXObVdqam+WgkEAlqtVuxQ2gXzlBe55zkqWMAXN5TQnrgkdijtQP5/bxt4Tp47dv4Lvk6udmpra53boBtzq0LSShDsh0EtFkuTY609vzVtpqWlITU11fZYp+qVNycAAAzZSURBVNMhMjISycnJ8Pf3v68c7pfp0Hl0zD+M6OhoKGS6aLeXSoEJg0JxbH8ekpKSoFarxQ7JZQwGA7RaLfOUCU/J83F9PYyf/ROBETGy/TsEAGazGSUlJbL+ewt4Xp5jxzwOP19vp7ZtHYkkNyskg4ODoVQqm/Q+VlZWNulRtAoLC7vr+WFhYQAaeibDw8Pvq02NRgONRtPkuFqtdvrFZNzgrlCeP4SnUvrI+kJlMBhwDK75GUoR85QXT8hzdLjFI/4Obd9+hnnKhDVPP19vp+cp559ba7nVP0W8vLwQFxfXZBhJq9UiISGh2dcMHz68yfnZ2dm282NiYhAWFmZ3Tn19PXbt2tVim0RERETkZj2SAJCamopp06Zh6NChGD58OP7617+itLQUr732GgBg+vTp6Nq1K5YtWwYAeOONNzBy5Ei89957GD9+PLZt24adO3di9+7dABqGtOfNm4elS5eiV69e6NWrF5YuXQpfX19MnTpVtDyJiIiIpM7tCsnJkyfjypUrWLJkCcrLyzFgwABs374dUVFRAIDS0lK7ez4SEhKwadMmvPXWW3j77bfRs2dPZGVlIT4+3nbOwoULcfPmTcyZMwfXrl1DfHw8srOz4efn1+75EREREbkLtyskAWDOnDmYM2dOs8/l5OQ0OTZp0iRMmjSpxfYEQUB6ejrS09OdFCERERGR/LnVPZJEREREJB0sJImIiIjIISwkiYiIiMghLCSJiIiIyCEsJImIiIjIIW45a1tqLBYLANdsmWQwGFBbWwudTifrlfSZp7wwT3lhnvLCPNvOer23Xv89GQtJJ7hx4wYAIDIyUuRIiIiIqL3cuHEDAQEBYochKsHCcrrNzGYzLly4AD8/PwiC4NS2dTodIiMjUVZWBn9/f6e2LSXMU16Yp7wwT3lhnm1nsVhw48YNRERE2G2C4onYI+kECoUC3bp1c+n38Pf3l/UvvBXzlBfmKS/MU16YZ9t4ek+klWeX0URERETkMBaSREREROQQZTo3mJY8pVKJ0aNHQ6WS950IzFNemKe8ME95YZ7kLJxsQ0REREQO4dA2ERERETmEhSQREREROYSFJBERERE5hIUkERERETmEhaTErVy5EjExMfD29kZcXBzy8vLEDsnpcnNzMW7cOEREREAQBGzdulXskJxu2bJlePjhh+Hn54eQkBBMmDABJ0+eFDssp1u1ahUGDRpkWwB4+PDh+Oabb8QOy+WWLVsGQRAwb948sUNxqvT0dAiCYPcVFhYmdlgucf78ebz00ksICgqCr68vhgwZgsLCQrHDcqro6Ogm76cgCJg7d67YoTmV0WjEW2+9hZiYGPj4+KBHjx5YsmQJzGaz2KHJEgtJCcvKysK8efOwaNEiHDx4EImJiUhJSUFpaanYoTlVTU0NBg8ejBUrVogdisvs2rULc+fOxd69e6HVamE0GpGcnIyamhqxQ3Oqbt264be//S0OHDiAAwcO4PHHH8f48eNx9OhRsUNzmf379+Ovf/0rBg0aJHYoLtG/f3+Ul5fbvg4fPix2SE537do1jBgxAmq1Gt988w2OHTuGDz/8EJ06dRI7NKfav3+/3Xup1WoBAM8//7zIkTnXe++9h48++ggrVqzA8ePH8f777+ODDz7An/70J7FDkyUu/yNh8fHxeOihh7Bq1Srbsb59+2LChAlYtmyZiJG5jiAI2LJlCyZMmCB2KC516dIlhISEYNeuXRg5cqTY4bhUYGAgPvjgA8yaNUvsUJyuuroaDz30EFauXIl3330XQ4YMQUZGhthhOU16ejq2bt2KoqIisUNxqTfffBN79uyR5YjP3cybNw//+Mc/cPr0aQiCIHY4TvPMM88gNDQUa9eutR177rnn4Ovri/Xr14sYmTyxR1Ki6uvrUVhYiOTkZLvjycnJyM/PFykqcpaqqioADUWWXJlMJmzatAk1NTUYPny42OG4xNy5c/H0009jzJgxYofiMqdPn0ZERARiYmLw4osv4syZM/+/vbsJiWqN4zj+mzvpOL5kZGpajjS9+VqMDsSoICQtLFwUZUkvI7YJsszBsCwiyBdoIUSGOBFSiCiklRGllW9IC0uaEgutLG2RDS3MTBzBee7i0kDXy10M5/TY9PvAwTMH0e8BF//zzPGM7CTFtbW1wWw2Y8+ePYiIiIDJZMLVq1dlZ6lqbm4ODQ0NKCgo8KkhEgAyMjLw+PFjjIyMAABevHiBvr4+bN++XXKZb+Kj3hepL1++YH5+HpGRkT8dj4yMxMTEhKQqUoIQAjabDRkZGUhKSpKdo7jBwUFYLBbMzs4iODgYt27dQkJCguwsxTU1NWFgYADPnj2TnaKaLVu24MaNG9iwYQM+f/6M8vJypKWlYWhoCGFhYbLzFDM6Oora2lrYbDaUlZWhv78fx48fh06nw6FDh2TnqeL27duYnJxEfn6+7BTFlZaW4uvXr4iLi4NWq8X8/DwqKiqQl5cnO80ncZBc5P59pSiE8Lmrxz9NYWEhXr58ib6+Ptkpqti4cSMcDgcmJyfR0tICq9WKnp4enxomP378iKKiInR0dCAgIEB2jmqys7M9+8nJybBYLFi7di2uX78Om80msUxZbrcbZrMZlZWVAACTyYShoSHU1tb67CB57do1ZGdnIzo6WnaK4pqbm9HQ0IDGxkYkJibC4XDgxIkTiI6OhtVqlZ3nczhILlIrVqyAVqtdsProdDoXrFLS7+PYsWNoa2tDb28vVq9eLTtHFf7+/li3bh0AwGw24+nTp7h06RLq6uoklylnYGAATqcTqampnmPz8/Po7e1FTU0NXC4XtFqtxEJ1BAUFITk5GW/evJGdoqioqKgFFzrx8fFoaWmRVKSusbExPHr0CK2trbJTVHHy5EmcOnUK+/btA/DPRdDY2Biqqqo4SKqA90guUv7+/khNTfX8V90PDx8+RFpamqQq8pYQAoWFhWhtbUVnZyfWrFkjO+mXEULA5XLJzlBUVlYWBgcH4XA4PJvZbMb+/fvhcDh8cogEAJfLhdevXyMqKkp2iqLS09MXPI5rZGQEsbGxkorUVV9fj4iICOzYsUN2iipmZmbw118/jzdarZaP/1EJVyQXMZvNhoMHD8JsNsNiscBut2N8fBxHjhyRnaao6elpvH371vP6/fv3cDgcWL58OQwGg8Qy5Rw9ehSNjY24c+cOQkJCPCvNoaGh0Ov1kuuUU1ZWhuzsbMTExODbt29oampCd3c3Hjx4IDtNUSEhIQvubw0KCkJYWJhP3fdaUlKCnJwcGAwGOJ1OlJeXY2pqyudWdYqLi5GWlobKykrk5uaiv78fdrsddrtddpri3G436uvrYbVasWSJb44AOTk5qKiogMFgQGJiIp4/f47q6moUFBTITvNNgha1K1euiNjYWOHv7y9SUlJET0+P7CTFdXV1CQALNqvVKjtNMf91fgBEfX297DRFFRQUeP5ew8PDRVZWlujo6JCd9UtkZmaKoqIi2RmK2rt3r4iKihJ+fn4iOjpa7Nq1SwwNDcnOUsXdu3dFUlKS0Ol0Ii4uTtjtdtlJqmhvbxcAxPDwsOwU1UxNTYmioiJhMBhEQECAMBqN4syZM8LlcslO80l8jiQREREReYX3SBIRERGRVzhIEhEREZFXOEgSERERkVc4SBIRERGRVzhIEhEREZFXOEgSERERkVc4SBIRERGRVzhIEhEREZFXOEgSEf2PDx8+QKPRwOFwyE4hIlp0OEgS0R9Lo9H875afn4+YmBh8+vTJpz5Dm4hIKfyIRCL6Y01MTHj2m5ubce7cOQwPD3uO6fV6hIaGykgjIvotcEWSiP5YK1eu9GyhoaHQaDQLjv37re3u7m5oNBq0t7fDZDJBr9dj69atcDqduH//PuLj47F06VLk5eVhZmbG87uEELh48SKMRiP0ej02b96Mmzdvyjp1IiJFLJEdQET0Ozp//jxqamoQGBiI3Nxc5ObmQqfTobGxEdPT09i5cycuX76M0tJSAMDZs2fR2tqK2tparF+/Hr29vThw4ADCw8ORmZkp+WyIiLzDQZKIyAvl5eVIT08HABw+fBinT5/Gu3fvYDQaAQC7d+9GV1cXSktL8f37d1RXV6OzsxMWiwUAYDQa0dfXh7q6Og6SRPTb4iBJROSFTZs2efYjIyMRGBjoGSJ/HOvv7wcAvHr1CrOzs9i2bdtPP2Nubg4mk+nXBBMRqYCDJBGRF/z8/Dz7Go3mp9c/jrndbgDwfL137x5WrVr10/fpdDqVS4mI1MNBkohIZQkJCdDpdBgfH+fb2ETkUzhIEhGpLCQkBCUlJSguLobb7UZGRgampqbw5MkTBAcHw2q1yk4kIvIKB0kiol/gwoULiIiIQFVVFUZHR7Fs2TKkpKSgrKxMdhoRkdf+Bk3dgE+QJDEAAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example Solution\n",
    "\n",
    "Python SymPy computes symbolic solutions to many mathematical problems including Laplace transforms. A symbolic and numeric solution is created with the <a href='https://apmonitor.com/pdc/index.php/Main/TransferFunctions'>following example problem</a>.\n",
    "\n",
    "Compute the analytic and numeric system response to an input that includes a step and ramp function.\n",
    "\n",
    "![input_laplace_function.png](attachment:input_laplace_function.png)\n",
    "\n",
    "The system transfer function is a stable system with two poles (denominator roots) and one zero (numerator root):\n",
    "\n",
    "$G(s)=\\frac{5\\left(s + 1\\right)}{\\left(s + 3\\right)^2}$\n",
    "\n",
    "Create the step and ramp signals as three individual functions. Compute the system response to each of those three inputs and then sum the signals.\n",
    "\n",
    "**Symbolic Solution (Python Sympy)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y1\n",
      "-10*(-exp(3*t) + 6*exp(3)*log(exp(-t)) + 7*exp(3))*exp(-3*t)*Heaviside(t - 1)/9\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEGCAYAAAB1iW6ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxfElEQVR4nO3deXyU9bX48c+ZyQ4JCXtISAbZlFUWWeKKaIWCC7hSRbBSrte69XfvbbXeauvtfmuvWmuttQriAiraWop7cUNQCIuyiCIGSMIalgAhy8x8f388kxBjQrbv5JmZnPfrlVdmMjPPc+SFh2/Oc57zFWMMSimlop/H7QCUUkrZoQldKaVihCZ0pZSKEZrQlVIqRmhCV0qpGBHn1om7du1qfD6fW6dXSqmolJ+fv98Y062+11xL6D6fj9WrV7t1eqWUikoisr2h17TkopRSMUITulJKxQhN6EopFSNcq6ErpZQtVVVVFBYWUl5e7nYo1iQlJZGdnU18fHyTP6MJXSkV9QoLC0lNTcXn8yEibofTasYYSkpKKCwspE+fPk3+XKMlFxHpLSLLRGSziGwUkdvreY+IyEMislVEPhGRkc2MXymlWqy8vJwuXbrERDIHEBG6dOnS7N84mrJC9wP/YYxZIyKpQL6IvGmM2VTrPZOB/qGvscCfQt+VUqpNxEoyr9aS/55GE7oxZhewK/T4iIhsBrKA2gn9UuAp48ziXSki6SKSGfpsu3K4rIoFKwuo9AfdDiXiDeiZytRhvdwOQ6mY0awauoj4gBHAR3VeygJ21npeGPrZ1xK6iMwF5gLk5OQ0M9To8K8te/jdG58DEGMLBquMAY/A4F6d6NO1g9vhKBUTmpzQRaQjsBi4wxhTWvflej7yjZ0zjDGPAY8BjB49OiZ31qhema+463wyOyW7HE3k2nuknLN/s4xHlm3lf68c7nY4SoXVHXfcwfTp0znnnHN44IEHmDt3LikpKQBccMEFvPDCC2RkZLT6PE3qQxeReJxk/owx5qV63lII9K71PBsobnV0Uagq4Pw75fXo8vxkuqcmMWNMDi+vLWLngTK3w1EqbA4cOMDKlSs555xzAHjggQcoKzvxd37mzJk88sgjVs7V6ApdnMr8X4HNxpjfN/C2V4BbRGQhzsXQw+2xfg4QCDoJPc6j92w15qZz+/LsRzv407tf8stpQ90OR8WIn/1jI5uK6xYRWmdQrzTuvXjwSd9TUFDA1KlT2bBhAwC/+93vOHr0KL169WLSpEkAPPTQQxQXFzNhwgS6du3KsmXLuOSSSzj77LO5++67Wx1nU7LOmcBM4HwRWRf6+raI3CQiN4XesxTYBmwF/gLc3OrIolRVwCm5xHl1hd6Ynp2SuHJ0Ni+uLmTX4eNuh6NUWCxfvpxRo0YBcNttt9GrVy+WLVvGsmXLAMjIyKCiooKSkpJWn6spXS4fUH+NvPZ7DPD9VkcTA06s0DWhN8W/n9eXRat28ud3t/HTS06+AlKqKRpbSbe1Xbt20a1bvdNua3Tv3p3i4mK6dOnSqnNpXcAyv5ZcmiU7I4XpI7N47uMdHDhW6XY4SrVYXFwcweCJduXqm4KSk5MbvUGovLyc5OTWN1Fo1rHMH9AVenPNOfsUKvxBFq7a4XYoSrVYjx492Lt3LyUlJVRUVLBkyRIATjvtNLZu3VrzvtTUVI4cOVLz3BjD7t27sbHhjyZ0ywLBICLg0YTeZAN6pHJmvy48vWI7/oDekKWiU3x8PPfccw9jx45l6tSpnHrqqQBMmTKFd955p+Z9c+fOZfLkyUyYMAGA/Px8xo0bR1xc60dr6XAuy6qChngttzTbrPE+5i7I581Ne5g8NNPtcJRqkdtuu43bbrvtGz+/6667OHToEOnp6dx6663ceuutNa8tWLCAm2+200eimceyQNBoD3oLTDytB9kZycz7sMDtUJSy7v7772fHjvpLikOGDGHixIlWzqMJ3bKqQFBbFlvA6xGuH5/LR18dYPMuuz3ESrlt7NixDBs2rN7Xvve971k7jyZ0ywJBoxdEW+iq0b1JivcwX1fpSrWIJnTL/EFDnFf/WFsiPSWBaSOyeXltEQe1hVGpZtPMY5k/ENQVeivMzvNR4Q+yaPXOxt+slPoaTeiW+fWiaKsM7JnK+FO6sEBbGFUM2LVrF1OnTgVg3bp1LF26tOa1JUuWcO+991o9nyZ0y/wBQ7yWXFplVp6PokPHeWvzXrdDUapVfv/739dc9Kyb0KdMmcIrr7zytcmLraWZxzJtW2y9C07rTlZ6MvM+/MrtUJRqkp/85Cc8+OCDNc/vvvtuHnroIRYvXsykSZOorKzknnvuYdGiRZx++uksWrQIEeG8886ruaPUBr2xyDJ/UGvorRXn9TBzfC6/fvUzPttdyqk909wOSUWTV++E3Z/aPWbPoTD51w2+fOONNzJ9+nRuv/12gsEgCxcu5Pnnn2f+/PkkJiYCcN9997F69Woefvjhms+NHj2a999/n6uuuspKmLpCt8wfMNqHbsHVo3uTGKctjCo6+Hw+unTpwtq1a3njjTcYMWIEFRUVTZ6yaIuu0C1zLorqv5OtldEhgWkjsnh5bRE/mnQq6SkJboekosVJVtLhNGfOHObNm8fu3bv57ne/26ZTFqtp5rHMHwwSryUXK2bl+SivCrJolbYwqsg3bdo0XnvtNVatWsVFF13EgAEDKCgoqHm97pRFgM8//5whQ4ZYi0ETumX+gF4UteW0zDTG9unMgpXbazYOUSpSJSQkMGHCBK666iq8Xi8dOnSgb9++NaNzJ0yYwKZNm2ouigIsW7aMKVOmWItBE7pl/qC2Ldo0O89H4cHjvL15j9uhKHVSwWCQlStXcuONN9b87JZbbmHevHkAdO7cmVWrVrFu3Tquvvpq9uzZw/Hjxxk61N5+upp5LNMbi+y6cFAPenVK0imMKqJt2rSJfv36MXHiRPr371/z82nTpjW4ccWOHTu4//77rcahF0UtCwSDxGuXizVxXg/Xjc/lt69t4fM9RxjQI9XtkJT6hkGDBrFt27Z6X5szZ069Pz/jjDOsx6ErdMu0hm7fNWfkkKAtjKoRzl71saMl/z2a0C3zB41uEG1Z5w4JXHZ6L15aU8Thsiq3w1ERKCkpiZKSkphJ6sYYSkpKSEpKatbntORimV83uAiLWXk+nl9dyAv5O5lz9iluh6MiTHZ2NoWFhezbt8/tUKxJSkoiOzu7WZ/RhG6ZXhQNj8G9OnGGL4P5Kwq44cw++mesviY+Pp4+ffq4HYbrtDZgmT+gm0SHy+y8Puw8cJxln+kURqXqo5nHMn/Q4NWSS1h8a3APeqYlMX9FgduhKBWRNKFbFtBpi2ETH5rC+P4X+9m690jjH1CqndGEbpk/oF0u4XTNGb1DLYzb3Q5FqYijmccyZ5PoVqzQq8rh4HYoXgfb3oGyA7ZCiwldOiZy8bBeLF5TSGm5tjAqVZt2uVjWqg0u9m+FJyfBsVqtV0npcOF9MGIm6MofgBvO9LF4TSEvrC7kxrO0s0GpapohLHNuLGpBQi87AM9eBSYIFz8E1zwL1y6GHoPhH7c5iX7f5/YDjkJDsjoxKjeDp1YUENQpjErV0IRuUSBoMMaZP9Is/kp4/no4vNNJ5KNmwalToP8FMPufcOkjsP8LeOIi2LMxPMFHmdl5PraXlPHO59rCqFQ1TegW+YNBgObf9LL0P6HgfbjkYcgZ9/XXRGDEtTDnLYhLgvmXwN7NliKOXpOG9KRHWiLz9OKoUjU0oVtUvQlDs0ouezbCmvmQdxsMv7rh93XpC7OXgCcO5l8Mez9rZbTRLd7r4dqxubz3+T6+3HfU7XCUigia0C2qCoQSenNKLhsWg3ichN6Y6qQuHnj6cjjavssNM8bkkOD18JROYVQK0IRuVbNX6MY4Cb3PudDx5LuD1+jaH659AcpKYOG14K9oYbTRr1tqIlOHZfJifiFHtIVRqcYTuog8ISJ7RWRDA6+fJyKHRWRd6Ose+2FGB3/AqaE3uQ+9eA0cLIAhlzfvRJnDYdqfoPBjWPID5x+Gdmr2mT6OVQZ4Mb/Q7VCUcl1TVujzgEmNvOd9Y8zpoa/7Wh9WdPI3d4W+4SXwxMNpFzf/ZIOnwbl3wrpnYMXDzf98jBiWnc6InHSeWrFdWxhVu9doQjfGvAfo7YpNUF1y8TblBqBg0Eno/S+E5PSWnfDcH8Fpl8Cb98C2d1t2jBgwO8/HV/uP8e4XsTMLW6mWsFVDHy8i60XkVREZ3NCbRGSuiKwWkdWxNIi+WlWo5NKkPUV3rIAjxc0vt9Tm8cBlj0CX/vDiDXC4fZYdJg/JpFtqom5Rp9o9Gwl9DZBrjBkO/AH4W0NvNMY8ZowZbYwZ3a1bEy8CRpETK/QmJPQNiyEuGQY0Vs1qRGIqXP30iZuT2uFF0oQ4D9eOzeGdLfv4av8xt8NRyjWtTujGmFJjzNHQ46VAvIh0bXVkUaimbbGxkkvAD5v+BgMnQWLH1p+42wBnpV6UD6/+sPXHi0LfGZtDvFd0la7atVYndBHpKSISejwmdMyS1h43GjW5bXH3J07bYUsuhjZk0CVw5h2QPw/y59s7bpTonprElKFOC+PRCr/b4Sjliqa0LT4HrAAGikihiNwoIjeJyE2ht1wBbBCR9cBDwDUmVrbebqaqYBPbFqvnsWSebjeAiffAKROcUQKFq+0eOwrMyvNxtMLPYm1hVO1Uo+NzjTEzGnn9YaD99s3VcmKF3si/k3s2QHwHyLA8+tXjhSuegMfOhUUz4d/ehY7d7Z4jgo3IyWB473Tmryhg5rhcPLpzlGpn9E5Ri/yBJl4U3b0BegwKz3zzlM5w9TNw/GDoImml/XNEsNl5uWzbd4wPtu53OxSl2pwmdIuqpy2etG3RGNjzKfQYEr5AMofBpQ87rZHt7E7Sbw/NpGvHRObpxVHVDmlCt8jflLbF0iIoP+xsXBFOQ69wbjxa9zQsfzC854ogiXFevjM2h2Vb9lKgLYyqndGEblF1ySX+ZNMWd4dG4vQcGv6AzrsLBk+Ht34Km5eE/3wR4rqxOXhFeGqFzkpX7YsmdIsCTdngYs+nzvdwr9DB2RzjskcgayS89D3Y+XH4zxkBuqcl8e2hmbyweifHtIVRtSOa0C06cWPRyRL6RsjwOXd4toX4ZJixEDr2gGeubDe7Hc3K83Gkws9La7SFUbUfmtAtqmlbbKzkEs4LovXp2B2u/5uzhd2C6XBoR9ue3wUjc9IZlt2JeR8W0E5vi1DtkCZ0ixodn1tZBge+bPuEDs5vBTNfgqpj8NRlcGR328fQhkSEWeN9fKktjKod0YRuUaMbXOzdDCYIPV1I6ODU7b/zgpPMn/w2HC5yJ442MnV4Jl07Juh8F9VuaEK3qNG2xT2hDhc3VujVcsbCzJed/UjnfTumyy+JcV5mjMnh7c/2sqOkzO1wlAo7TegWVa/Q4xu6A3TPBkjoCOm5bRhVPXLGwvV/h7KD8OQU2P+Fu/GE0bVjc0MtjAVuh6JU2GlCt6hmhd5QyWX3BqfsEY5b/psrexTM+jtUlcHjF0DBcrcjCouenZKYNKQni7SFUbUDEZBZYsdJx+ca47QsulluqavXCJjzltMFs+Ay+OR5tyMKi9l5Po6U+3lpbWxfM1BKE7pF/pNNWzy8Eyra4Jb/5urcB258A3qPdW4+ev1uCFS5HZVVo3IzGJKVxlPawqhinCZ0i/wnu7Go5Evne7eBbRhREyVnwHUvwRlzYMXD8ORkOLTT7aisqW5h/GLvUT78sl3uvaLaCU3oFvmDQTxC/XO4S4ud72m92jaopopLgCn3O/PU934Gfz4bNrwUM5MaLx7ei84dEnQKo4ppmtAt8gdNw5tbHAkl9NQITejVhlzubIyR4YMXb4BF10HpLrejarWkeC8zxvTmrc172HlAWxhVbNKEbpE/EGy4B720GFK6QHxS2wbVEl36wo1vwYX3wda34I9j4aM/R31t/bpxuXhEWLBSpzCq2KQJ3SJ/0DR8l2hpceSvzmvzxsGZt8NNy6HXcHj1h/DIOPjsn1FbhsnslMykwT1Z+PEOyiq1hVHFHk3oFgWCpuE5LqVFkVs/P5mu/eD6V2DGIhAPLPwOPD7Rma8eGhccTWbl+Sgt9/O3tcVuh6KUdZrQLaoKmIYnLZbuis6EDs5c9YGT4N9XwMUPQlkJLLoW/jQe8udBxVG3I2yyM3wZDMpMY762MKoYpAndokAwWP8KvaocyvZDWlbbB2WTNw5GzYZb8mH64+CJg3/cDvef6uxdWpgf8eUYEWF2no8te46wYpu2MKrYogndIn+ggRr6kVCXSFpm2wYULt44GHYl3PQBfPcNOG0qrHsWHj8fHhwGb94DO1dBMOB2pPW65PReZKTE6xRGFXPi3A4gljTYtliT0KO05NIQEWfQV85YmPRr2LLU6V1f8UdnY+rkDOh7PpxyHuTkOd0zcpLdnNpIUryXa8bk8Od3v6TwYBnZGSluh6SUFZrQLfIHG2hbrLmpKMpLLieTnA6nf8f5KjsAX/4Ltr7ttD1uWOy8p0M3yB4DvU6HzOGQNRo6dHEl3OvG5fLYe9tYsHI7d00+zZUYlLJNE7pF/kADXS6loaFQqTFScmlMSmcYeoXzFQxCyRewYwVsXwFF+c5KHgPeBBgxE866A9Jz2jTErPRkvjWoBws/3skdEweQnOBt0/MrFQ6a0C0KNNSHXloMCamQlNb2QbnN43Hm13Qb6FxQBag4Ars/daY7rnkK1syHkbNg8m+d+nwbmZ3n49UNu/n7uiKuGdO2/6AoFQ56UdSiqoZq6KXFsVc/b43EVMjNg4sfgNvXO4l+9V9hye1t2iUzpk9nTu2ZqhtJq5ihCd2iBtsWS4tjp8PFtk5ZzlCwc38Ea5+Gf/28zU4tItxwpo/Pdh/ho68OtNl5lQoXTegWVQVMwxdFY/mCqA3n3QUjr4f3fwcf/6XNTnvp6VmkawujihGa0C0KBA3xde8UDfjh6B4tuTRGBKb8HwyY7MyN2bmqTU6bFO/lmjNyeH3jbooOHW+TcyoVLprQLap32uKxvWAC7afDpTW8cTD9MWeI2d9vdu6wbQPXjXMuiD6tUxhVlNOEbpE/aIiv2+XSHnrQbUpKg0sehP2fw7u/bpNTZmekcOGgHiz8eAflVZF5d6tSTaEJ3aJAsJ4aenUPupZcmq7fBTDiOlj+EBStaZNTzs7rw8GyKl5Zp1MYVfTShG5RVSD4zWmL1bv96Aq9eb71C+jYHf7+ffBXhv10405xWhif1BZGFcUaTegi8oSI7BWRDQ28LiLykIhsFZFPRGSk/TCjQ73z0EuLwJvo3D2pmi45Hab+H+zdBKvC3/UiIszK87F5VymrCg6G/XxKhUNTVujzgEkneX0y0D/0NRf4U+vDik71ti1W96BHwFCqqDNwslN+eec3cGx/2E932elZdErWFkYVvRpN6MaY94CT3XVxKfCUcawE0kWkXbZ0BIKG+Lp3ih7ZpeWW1vjWL6DyKCz7ZdhPlZzg5ZozevPaxt0UawujikI2auhZwM5azwtDP/sGEZkrIqtFZPW+ffssnDqy+IMG7ze6XIq0ZbE1up8KZ8yB/Cdhz6awn+66cbkYY3jmI21hVNHHRkKvr5ZQ71UlY8xjxpjRxpjR3bp1s3DqyOIPBomvXXIxRue42HDenZCYBq/fFfZZL707pzDxtB489/FObWFUUcdGQi8Eetd6ng20y96vQMDgrV1yKSuBQKWWXForpbMzGmDbO/DFm2E/3Q15Pg4cq+Qf69vlX2MVxWwk9FeA60PdLuOAw8aYXRaOG3WqgsGvj8+tualISy6tdsaNkNEH3v6ZM2M9jMb37cKAHh11CqOKOk1pW3wOWAEMFJFCEblRRG4SkZtCb1kKbAO2An8Bbg5btBHuG22L1VvPpWrJpdW88XD+f8OeDbDhxbCeqrqFcWNxKfnbtYVRRY+mdLnMMMZkGmPijTHZxpi/GmMeNcY8GnrdGGO+b4zpa4wZaoxZHf6wI48xhqq6OxYdDyUD7UG3Y/B06DnUGbEb5puNpo3IIi0pjnnawqiiiN4pakkw9Jv51+4UrU7oyRltH1As8nhg4k/h0HbInxfWU6UkxHHV6N68umE3uw+3zZAwpVpLE7ol/lBd92s3Fh0/5HxP6tT2AcWqfhMh9yx477dQcTSsp7p+vI+gtjCqKKIJ3RJ/wFmif23aYvkhSOwEHt2A2BoRuOCncGwffBTem5JzuqQw8dTuPPvRDir82sKoIp8mdEv8oZrL19oWjx+EZF2dW9f7DGcjjOV/OFHWCpPZeX0oOVbJkvXtsnFLRRlN6Jb4A07JJa5uyUXr5+Fx/t1QcRg+/ENYT3Nmvy70696R+Su0hVFFPk3olgRCK/Sv9aEfPwhJ6e4EFOt6DnW6XlY+CkfDN0ZCRJg1PpdPCg+zduehsJ1HKRs0oVtSVZ3QPXVq6MnprsTTLkz4MfiPwwf/F9bTTB+ZTWpiHPOWF4T1PEq1liZ0SwKB6oRet4auJZew6dofhs+AVY/D4aKwnaZDYhxXju7N0k93sadUWxhV5NKEbkl122JNycUYp4auJZfwOvdHYILw3v+G9TTXj88lYAzPfLQjrOdRqjU0oVtyossllNArj0GwSlfo4ZaRC6NmwdoFcOCrsJ3G17UDEwZqC6OKbJrQLfHXLbmUH3K+aw09/M7+T/DEwbu/CetpZuX52H+0gqWfagujikya0C2pKblUr9Cr7xLVFXr4pWXCmO/BJ4tg35awnebsfl05pVsHvTiqIpYmdEv8ddsWq2940Rp62zjzBxDfAZb9Imyn8HiEWeN9rC88zNodOoVRRR5N6JbU9KF/o+SiK/Q20aELjL8ZNv0diteF7TSXj8qmY2KcbiStIpImdEuqAnW6XGomLaa7E1B7NP77zm9E//p52E7RMTGOK0Zl889Pd7H3iLYwqsiiCd2SQN0bi7SG3vaSOsFZP4Ctb8L2D8N2muvH51IVMDyrLYwqwmhCt6S6y6WmbfH4QRAvJHR0Map2aMxc6NgT3r4vbBtKn9KtI+cN7MYzH+2g0h/e7fCUag5N6JZUXxSN99aqoSenO+NeVdtJSIFzfwg7VoR1Q+lZeT72Hang1Q3awqgihyZ0S6qnLX5tha7lFneMmAkZPmeVHqYNpc/t340+XTvoFnUqomhCt+TECr1WDV1bFt0RlwAT7oY9n8LGl8JyCo9HuH58Lmt3HGK9TmFUEUITuiWBuhtc6ArdXUOugO6DnY6XQFVYTnHFqGw6JHi1hVFFDE3ollTV3eBCR+e6y+Nxtqo7+FXYNpROTYrnilHZLPlkF/uPVoTlHEo1hyZ0S76xwYXuVuS+/heC72xnxkuYNpS+Ps9HZSDIc9rCqCKAJnRLqmpPWwwGofyw1tDdJgIX/MzZUHrFH8Nyir7dOnLOgG48/dH2mt/SlHKLJnRLAqH/meM9HmevS4yu0CNB9igYdCl8+FDYtqqbnZfLntIKXtuwOyzHV6qpNKFbUjMP3St623+kOf8eqDoetvG65w3oTm6XFG1hVK7ThG5JTduix6O3/Uearv1g9A2Q/yTs/8L64Z0WRh/52w/yaeFh68dXqqk0oVsSqF1D19G5kefcOyEuGd76aVgOf+XobFISvLpKV67ShG7J19oWdXRu5OnYDc7+AXy2BAqWWz98WlI8l4/M5h+fFFOiLYzKJZrQLQkEDR5xfv3WGnqEGnczpGXBG/8dlpEAs/JyqfQHWbhqp/VjK9UUmtAtqQqYE5tbVNfQteQSWeKTYeI9ULwGNiy2fvh+3VM5u39XFqzQFkblDk3olgSCwRM3FZUfcuq18UmuxqTqMfQqyBzu1NIry6wfftZ4H7tLy3lj4x7rx1aqMZrQLakKmDqTFtNdjUc1wOOBSb+G0kL48A/WDz/h1O707pys812UKzShWxIImq/vVqQXRCNXbh4MugyWPwCHi6we2hvaSPrjggNsKNIWRtW2NKFb4g8a4ry1auhaP49sF94HwQC8/TPrh75ydG+S4708taLA+rGVOhlN6Jb4A8E6kxZ1hR7RMnIh7xb4ZBHsXGX10J2S45k+Mou/rSvmwLFKq8dW6mSalNBFZJKIbBGRrSJyZz2vnycih0VkXejrHvuhRrZA0NSatKg19Khw1g+gYw947UfW2xhn5flCLYw6hVG1nUYTuoh4gT8Ck4FBwAwRGVTPW983xpwe+rrPcpwRrypYp21RV+iRLzHVKb0U5cP6Z60eekCPVPL6duHpFdtrtidUKtyaskIfA2w1xmwzxlQCC4FLwxtW9AkEQyUXfyVUHdMaerQYdjX0Hgtv3nvi/gFLZuf5KD5czpubtIVRtY2mJPQsoPatb4Whn9U1XkTWi8irIjK4vgOJyFwRWS0iq/ftC88oU7f4q9sWa277T3czHNVUIvDt/4WyEnjn11YPPfG0HmRnJOt8F9VmmpLQpZ6fmTrP1wC5xpjhwB+Av9V3IGPMY8aY0caY0d26dWtWoJHOX11D10mL0SdzuDON8ePHYM8ma4f1hjaS/uirA2wqLrV2XKUa0pSEXgj0rvU8Gyiu/QZjTKkx5mjo8VIgXkS6WosyCvira+g6aTE6nf8TSEqDpf8Fpu56peWuHp1DUrxHWxhVm2hKQl8F9BeRPiKSAFwDvFL7DSLSU0Qk9HhM6LgltoONZDVtizppMTqldIaJ98L2D+CT560dtlNKPNNGZPPy2iIOagujCrNGE7oxxg/cArwObAaeN8ZsFJGbROSm0NuuADaIyHrgIeAaYywuc6LAiZKLTlqMWiNnQdYoeONuqxdIZ+XlUuEPsmi1TmFU4dWkPnRjzFJjzABjTF9jzC9CP3vUGPNo6PHDxpjBxpjhxphxxpgPwxl0JHJW6B4oO+D8IKWzuwGp5vN4YOr/ORdI//Vza4c9tWca407pzAJtYVRhpneKWhIIhrpcykpAvJDYye2QVEtkDocxc2HV41C81tphZ+f1oejQcd7avNfaMZWqSxO6Jf6gId4bSugpnZ3VnopOE34MHbvDkh84814suOC07mSl6xRGFV6adSyp6UMvK4GULm6Ho1ojqRNc9Etnhf7xX6wcMs7rYeb4XFZsK+Gz3drCqMJDE7ol/mDQmbZYdgCStX4e9YZcDv0ugH/9DxwutHLIq0f3JjHOw/wPt1s5nlJ1aUK3xF89D7265KKimwhMud8puSz9oZVDZnRIYNqILF5eW8jhsiorx1SqNk3oltSUXI4f0JJLrMjwwYS7YMs/YfM/rBxyVp6P8qogi1brFEZlnyZ0S/zBIPGiNfSYM+5m6DHEuYO0vPU7EJ2WmcaYPp15asV2AsF2dauGagOa0C0JBA0pUgZBvyb0WOKNh0segqN7nImMFtyQ56Pw4HHe3qxTGJVdmtAt8QcNnYKhFZwm9NiSNcpZqec/CQUftPpwFw7qQa9OSczX+S7KMk3olvgDhtTgEeeJXhSNPRPudmrqr9wKVcdbdag4r4drx+WyfGsJn+85Yic+pdCEbo0/GCRVV+ixKyEFLn4IDmyDd37V6sPNGJNDQpxHbzRSVmlCt8QfMHQIhG4Y0RV6bDrlXBgxEz78AxStadWhOndI4NLhvXhpTRGHj2sLo7JDE7oFxhj8QUOHwCHnB7pCj13f+rmzsfTfvw/+ilYdalaej+NVAV7QKYzKEk3oFlR3n3UIlIInDhLT3A1IhU9yOlz8IOzdBO/9b6sONSSrE2f4MrSFUVmjCd2CqtBI1BT/IWd1LvXt2qdixoCLYPh34P3fQ/G6Vh1qdl4fdhwoY9lnOoVRtZ4mdAuqV1cp/kM6x6W9mPRL6NAtVHpp+U5E3xrcg55p2sKo7NCEboE/4CT05KrDWj9vL5Iz4OIHYM8GeO+3LT5MvNfDdeNyeP+L/Wzdqy2MqnU0oVvgDzoll6SqQ9rh0p4MnHyi9FKU3+LDnGhh1CmMqnU0oVvgD5VcEqsO6Qq9vZn0K6fr5eV/h6ryFh2iS8dELh7Wi8VrCikt1xZG1XKa0C3wBw1CkEQtubQ/yelw6cOwfwssa/k+pLPzfJRVBnhhtZ3Z66p90oRuQSBgSKUMjwloQm+P+k2EUTfAhw/D9hUtOsTQ7E6Mys3gqRUFBLWFUbWQJnQLqoJBOovOcWnXvvVzyMiFl/8NKlp2cXNWno/tJWW887m2MKqW0YRuQSBo6Ex1QtcVeruU2BEuexQO7YDXf9yiQ0we0pMeaYnM04ujqoU0oVtQFQiSLkedJ7pCb79yx8OZt8Oap2DLq83+eLzXw7Vjc3nv8318ue9oGAJUsU4TugWBoKlVctEVers24cfODkev3ArH9jf74zPG5JDg9fCUTmFULaAJ3YKqgCFDSy4KIC4Rpj/mbFf3ym1gmneBs1tqIlOHZfJifiFHtIVRNZMmdAuqV+hBTwIkdHQ7HOW2HoNh4r3O5tJrFzT747PyfByrDLA4X1sYVfNoQrfAHwySwRH8iRk6mEs5xt0Mfc6BV++Eki+b9dHhvdMZkZPO/BXbtYVRNYsmdAv8AWeF7k9KdzsUFSk8HqfrxRsHL82FgL9ZH5+d5+Or/cd474t9YQpQxSJN6BYEgoZ0OUogSTtcVC2dsmDqA1C0utmz0ycPyaRbaiLz9OKoagZN6BZUBYJ05ghBTeiqriHTYfgMZyLjjo+a/LGEOA/Xjs3hnS37+Gr/sTAGqGKJJnQLAkFDhhwhqLPQVX0m/xbSc+ClOVBe2uSPfWdsDvFe0Y2kVZNpQrfAH/CTzlGM3lSk6pOUBtP/AoeLYOl/Nvlj3VOTuGR4FgtX7WD/0dbtX6raB03oFngqDuMVg0nWHnTVgN5j4NwfwSeL4JMXmvyxmyf0pcIf5PH3vwpjcCpWaEK3wHP8gPNASy7qZM7+D+g9Dpb8AA40LUH37daRqcN6sWBFAQePtXyrO9U+NCmhi8gkEdkiIltF5M56XhcReSj0+iciMtJ+qJErruKQ86CDrtDVSXjj4PK/gHhg8RwINO1O0Fsm9ONYZYAnlusqXZ1cowldRLzAH4HJwCBghogMqvO2yUD/0Ndc4E+W44xocRXOCl00oavGpOc4e5EWrYZ3ftWkjwzsmcqkwT2Zt7yAw8d1HIBqWFwT3jMG2GqM2QYgIguBS4FNtd5zKfCUMcYAK0UkXUQyjTG7bAe8/bM17Fr5ou3Dtkrq3k8B8KZ0dTkSFRWGTIcv33b2Ij3lPOeO0kbccn4/Xtu4mzsXf8KQrE7hjzGGiAngDVbiDVbhNZV4glV4jT/0vQpP0I+YAB4TwGNCjwkixo/HBBETRAggxtmZTGp+FkSMAYIIBgwIQcCE3mtqHsOJO34FQ6e+YxiaN9n6f2tTEnoWsLPW80JgbBPekwV8LaGLyFycFTw5OTnNjRWA/dvWM67gjy36bDjtI4NOnXu5HYaKFpN/CztWOneR3rS80XLdkKxOXDy8F/9YX8yrG3a3UZDuS6GcrnKYLpSSIUfoLEfoxDHSpIw0jpEqx+nAcTpynA5STjIVdKCcZKkkkUqSqCRBAm7/Z3zDirLrIQwJXUwj0+BE5ErgImPMnNDzmcAYY8yttd7zT+BXxpgPQs/fBn5ojGlwK/TRo0eb1atXNzvgYCCA3x95v3Z6vXF445ry76NSIbvWw+MXQN+JMOO5RucAGWOoCsTQbBcThNIi5OB25OBXyKHtUFqIlBYhpcVwbC9SWf9ceINAYiokpmESU50NRuI7QkIKJj4F4lMgPhnikjBxyRCXAN5EZxqmNwG88eCJx3gTnGsbnjjwxDvfxRN6Xv3YA+LFiAc8XudnVD8WQJzv4qn/8Te+gyc+kbjElBb9sYlIvjFmdH2vNSUDFQK9az3PBopb8B4rPF4vCV5vOA6tVNvKHA4X3gev3QkfPwZj/+2kbxcREuKidPhb2QHnH7Ddn8CejbBvC+z/HKrKTrxHvJCa6YxM6DXcedyxO3ToDh26OaOpUzIgOQNJTHMSKhClfyJh0ZSEvgroLyJ9gCLgGuA7dd7zCnBLqL4+Fjgcjvq5UjFn7E2w7R14478hZzxkDnM7otYzBvZ9BtuXw86PYedHcLDgxOupvaD7qTBqNnTpB51Pgc59IC3bWS2rFmv0T88Y4xeRW4DXAS/whDFmo4jcFHr9UWAp8G1gK1AG3BC+kJWKISJw6SPw6Jnw4g0w912nfBBtjh+EL96ErW87/0AdDdX5O3SHnLEw6gbnN5LM4bpNYxg1WkMPl5bW0JWKSQUfwPyLYdjVMO1Rt6NpmmP7YePLsPkV2P4hBP1OWeSU85wv39mQ4dM9AixrbQ1dKRVuvrPgnB/Cu7+GPufC6TPcjqh+/gr4bAmsX+isxk0Aug6EvFth4BTIGuVcRFSu0ISuVKQ457+g4H34539A9mjo2t/tiE448BXkz4O1T0PZfkjLcpL4sKucLfdURNCErlSk8MbB5Y/Dn86EF26AOW9BfJK7MRWuhuUPwuZ/OK14AyfD6BvglPN1JR6BNKErFUnSejk19Gevgjfuhin3t30MxsBX78G7v3E6VRI7wVl3wBnfc1oKVcTShK5UpBlwEYy/BVY87IwFGHRp25274ANY9ksnkadmwkW/gpEznZt4VMTThK5UJJp4L+xYAX+/1Wn1y/CF93y7N8Bb98LWt6BjT2c0wchZ7pd8VLNoEUypSBSXAFc84Tx+4Qbwh2kW+uEiePnf4dGzoHCVc+fq7eucu1Y1mUcdTehKRaoMH1z2Ryhe46yebao8Bst+BX8YBRsWOx0rt6+HM293ZqCoqKQlF6Ui2WkXO+MBVj4CuXnO89YIBuGThfD2/8CRYhg8DS74GWTk2olXuUoTulKR7sL/cWai/O370HNoy+vpX70Hr9/tDMjqNRKufBJyxlkNVblLSy5KRbq4BCf5AjxzZZP3I61RmA9PX+6MFjh+EC7/K8x5W5N5DNKErlQ0yPDBjGfh2D74ywRntX0ywaDznmevhsfPh6I1TmnlllUw9Aq9KShGaclFqWjhOwu+9y949hpYMA3O/g9nc4zM4U5HSvlh2Pc5fPWOc4v+wQJISofzf+J0rWgveczTaYtKRZvyUnj5JtjyT+e5Jx6SM+DY3hPv8Z0NI2bCoEu0ayXG6LRFpWJJUppTfjmyx+kdL/wYjpVA137O5MOeQyG9d+PHUTFHE7pS0Sq1B5w21flSCr0oqpRSMUMTulJKxQhN6EopFSM0oSulVIzQhK6UUjFCE7pSSsUITehKKRUjNKErpVSMcO3WfxHZB2xv4ce7AvsthmNLpMYFkRubxtU8GlfzxGJcucaYbvW94FpCbw0RWd3QLAM3RWpcELmxaVzNo3E1T3uLS0suSikVIzShK6VUjIjWhP6Y2wE0IFLjgsiNTeNqHo2redpVXFFZQ1dKKfVN0bpCV0opVYcmdKWUihFRl9BFZJKIbBGRrSJyp9vxAIjIEyKyV0Q2uB1LbSLSW0SWichmEdkoIre7HROAiCSJyMcisj4U18/cjqk2EfGKyFoRWeJ2LNVEpEBEPhWRdSISMXs3iki6iLwoIp+F/p6Nj4CYBob+nKq/SkXkDrfjAhCRH4T+zm8QkedEJMnq8aOphi4iXuBz4EKgEFgFzDDGbHI5rnOAo8BTxpghbsZSm4hkApnGmDUikgrkA5dFwJ+XAB2MMUdFJB74ALjdGLPSzbiqicj/A0YDacaYiNgOSEQKgNHGmIi6SUZE5gPvG2MeF5EEIMUYc8jlsGqEckYRMNYY09IbGW3FkoXzd32QMea4iDwPLDXGzLN1jmhboY8BthpjthljKoGFwKUux4Qx5j3ggNtx1GWM2WWMWRN6fATYDGS5GxUYx9HQ0/jQV0SsLEQkG5gCPO52LJFORNKAc4C/AhhjKiMpmYdMBL50O5nXEgcki0gckAIU2zx4tCX0LGBnreeFRECCigYi4gNGAB+5HApQU9ZYB+wF3jTGRERcwAPAD4Ggy3HUZYA3RCRfROa6HUzIKcA+4MlQiepxEengdlB1XAM853YQAMaYIuB3wA5gF3DYGPOGzXNEW0KXen4WESu7SCYiHYHFwB3GmFK34wEwxgSMMacD2cAYEXG9VCUiU4G9xph8t2Opx5nGmJHAZOD7oTKf2+KAkcCfjDEjgGNARFzXAgiVgC4BXnA7FgARycCpKPQBegEdROQ6m+eItoReCPSu9Twby7+yxJpQjXox8Iwx5iW346kr9Cv6O8AkdyMB4EzgklC9eiFwvog87W5IDmNMcej7XuBlnPKj2wqBwlq/Xb2Ik+AjxWRgjTFmj9uBhFwAfGWM2WeMqQJeAvJsniDaEvoqoL+I9An963sN8IrLMUWs0MXHvwKbjTG/dzueaiLSTUTSQ4+Tcf6if+ZqUIAx5i5jTLYxxofzd+tfxhirK6iWEJEOoYvahEoa3wJc76gyxuwGdorIwNCPJgKuXnCvYwYRUm4J2QGME5GU0P+bE3Gua1kTZ/Ng4WaM8YvILcDrgBd4whiz0eWwEJHngPOAriJSCNxrjPmru1EBzopzJvBpqF4N8GNjzFL3QgIgE5gf6kDwAM8bYyKmRTAC9QBednIAccCzxpjX3A2pxq3AM6EF1jbgBpfjAUBEUnC64f7N7ViqGWM+EpEXgTWAH1iL5REAUdW2qJRSqmHRVnJRSinVAE3oSikVIzShK6VUjNCErpRSMUITulJKxQhN6CrmiUiXWpP3dotIUejxURF5xO34lLJF2xZVuyIiPwWOGmN+53YsStmmK3TVbonIedUzz0XkpyIyX0TeCM0eny4ivw3NIH8tNEIBERklIu+GhmS9HhpRrFRE0ISu1Al9cUbnXgo8DSwzxgwFjgNTQkn9D8AVxphRwBPAL9wKVqm6ourWf6XC7FVjTJWIfIozWqL69vpPAR8wEBgCvBm6Dd+LMwZVqYigCV2pEyoAjDFBEakyJy4wBXH+XxFgozHG9W3WlKqPllyUarotQLfqfTNFJF5EBrsck1I1NKEr1UShbQ+vAH4jIuuBdVieZ61Ua2jbolJKxQhdoSulVIzQhK6UUjFCE7pSSsUITehKKRUjNKErpVSM0ISulFIxQhO6UkrFiP8PEM1dxb6zXcYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import sympy as sym\n",
    "from sympy.abc import s,t,x,y,z\n",
    "import numpy as np\n",
    "from sympy.integrals import inverse_laplace_transform\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Define inputs\n",
    "# First step (up) starts at 1 sec\n",
    "U1 = 2/s*sym.exp(-s)\n",
    "# Ramp (down) starts at 3 sec\n",
    "U2 = -1/s**2*sym.exp(-3*s)\n",
    "# Ramp completes at 5 sec\n",
    "U3 = 1/s**2*sym.exp(-5*s)\n",
    "\n",
    "# Transfer function\n",
    "G = 5*(s+1)/(s+3)**2\n",
    "\n",
    "# Calculate responses\n",
    "Y1 = G * U1\n",
    "Y2 = G * U2\n",
    "Y3 = G * U3\n",
    "\n",
    "# Inverse Laplace Transform\n",
    "u1 = inverse_laplace_transform(U1,s,t)\n",
    "u2 = inverse_laplace_transform(U2,s,t)\n",
    "u3 = inverse_laplace_transform(U3,s,t)\n",
    "y1 = inverse_laplace_transform(Y1,s,t)\n",
    "y2 = inverse_laplace_transform(Y2,s,t)\n",
    "y3 = inverse_laplace_transform(Y3,s,t)\n",
    "print('y1')\n",
    "print(y1)\n",
    "\n",
    "# generate data for plot\n",
    "tm = np.linspace(0,8,100)\n",
    "us = np.zeros(len(tm))\n",
    "ys = np.zeros(len(tm))\n",
    "\n",
    "# substitute numeric values for u and y\n",
    "for u in [u1,u2,u3]:\n",
    "    for i in range(len(tm)):\n",
    "        us[i] += u.subs(t,tm[i])\n",
    "for y in [y1,y2,y3]:\n",
    "    for i in range(len(tm)):\n",
    "        ys[i] += y.subs(t,tm[i])\n",
    "\n",
    "# plot results\n",
    "plt.figure()\n",
    "plt.plot(tm,us,label='u(t)')\n",
    "plt.plot(tm,ys,label='y(t)')\n",
    "plt.legend()\n",
    "plt.xlabel('Time')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Numeric Solution (Python Gekko)**\n",
    "\n",
    "An alternative to a symbolic solution is to numerically compute the response in the time domain. The transfer function must first be translated into a differential equation.\n",
    "\n",
    "$G(s)=\\frac{Y(s)}{U(s)}=\\frac{5\\left(s + 1\\right)}{\\left(s + 3\\right)^2}$\n",
    "\n",
    "$Y(s)\\left(s + 3\\right)^2=5\\left(s + 1\\right)U(s)$\n",
    "\n",
    "$Y(s)\\left(s^2 +6s+9\\right)=\\left(5s + 5\\right)U(s)$\n",
    "\n",
    "$\\frac{dy^2(t)}{dt^2}+6\\frac{dy(t)}{dt}+9y(t)=5\\frac{du(t)}{dt}+5u(t)$\n",
    "\n",
    "There is additional information on solving differential equations with <a href='https://apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations'>Python GEKKO</a> or with <a href='https://apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations'>Python Scipy ODEINT</a>. The same solution is found with a analytic or numeric approach. The advantage of an symbolic (analytic) solution is that it is highly accurate and does not rely on numerical methods to approximate the solution. Also, the solution is in a compact form that can be used for further analysis. Symbolic solutions are limited to cases where the input function and system transfer function can be expressed in Laplace form. This may not be the case for inputs that come from data sources where there the input function has random variation. A symbolic solution with Laplace transforms is also not possible for systems that are nonlinear or complex while numeric solvers can handle many thousands or millions of equations with nonlinear relationships. The disadvantage of a numeric solution is that it is an approximation of the true solution with possible inaccuracies. Another disadvantage is that solvers may fail to converge although this is not typical on problems with an analytic solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEGCAYAAAB1iW6ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxiklEQVR4nO3deXxU9bn48c8zM9kJECAsISQDiqDsGFnCVUGwQkGt1AXcwGqprbjc9tfVXtt6r1281Vq01usKuIGKe3FrBVEWJWzKIogQIBD2LQGSzPL9/XEmIY2BbOfkzEye9+s1r2RmzpzzGOHhm+95vs9XjDEopZSKfR63A1BKKWUPTehKKRUnNKErpVSc0ISulFJxQhO6UkrFCZ9bF+7QoYPx+/1uXV4ppWLSihUr9htjMmt7z7WE7vf7KSgocOvySikVk0Rk26ne0ykXpZSKE5rQlVIqTmhCV0qpOOHaHLpSStklEAhQVFREWVmZ26HYJjk5mezsbBISEur9GU3oSqmYV1RURHp6On6/HxFxO5wmM8Zw4MABioqK6N69e70/V+eUi4h0E5EFIrJBRNaJyJ21HCMiMkNENovI5yIyuIHxK6VUo5WVldG+ffu4SOYAIkL79u0b/BtHfUboQeAnxpiVIpIOrBCRD4wx66sdMw7oGXkMBf4e+aqUUs0iXpJ5pcb899SZ0I0xxUBx5PsSEdkAdAWqJ/TLgdnG6sW7TETaikiXyGdbpOIjJ5i7fAfhsLYnPpVEn4epI7rTKkln/pSyQ4P+JomIHxgEfFrjra7AjmrPiyKv/VtCF5FpwDSAnJycBoYaW+atKOKhf34FQJwNHGxT2Yp/+kU93Q1EqThR74QuIq2AecBdxpijNd+u5SPfGJoaYx4HHgfIy8uL66FreTCMCGz9w3i3Q4laNzz1Kc8t284PLjyDBK9W0Kr4dddddzFx4kQuuOACHnroIaZNm0ZqaioAY8aM4eWXXyYjI6PJ16nX3yIRScBK5s8bY16t5ZAioFu159nAriZHF8MCIUOCR5PU6UwZ7mf30TLeW7fb7VCUcszBgwdZtmwZF1xwAQAPPfQQx48fr3r/hhtu4NFHH7XlWnWO0MWamX8K2GCMefAUh70JTBeROVg3Q4+05PlzgEAoTIJX51pOZ1TvjuS0S2XWkkIm9M9yOxwVJ3731jrW76o5idA052S15jeX9jntMYWFhUyYMIG1a9cC8Oc//5nS0lKysrIYO3YsADNmzGDXrl2MGjWKDh06sGDBAi677DLOP/987r777ibHWZ8h5AjgBuAiEVkdeXxbRG4VkVsjx8wHtgCbgSeAHzU5shgXDIXx6TTCaXk9wo3Dc1leeIi1O4+4HY5Sjli8eDHnnnsuAHfccQdZWVksWLCABQsWAJCRkUF5eTkHDhxo8rXqU+XyCbXPkVc/xgC3NTmaOFIRMjovXA9X5XXjgfc3MWtJIf971QC3w1FxoK6RdHMrLi4mM7PWbrdVOnbsyK5du2jfvn2TrqUZxyFBnXKplzYpCUwc3JU31uzi4LEKt8NRqtF8Ph/hcLjqeeWioJSUlDoXCJWVlZGSktLkGDShO8SaQ9cfb31MyfdTEQwzZ/l2t0NRqtE6derE3r17OXDgAOXl5bz99tsAnH322WzevLnquPT0dEpKSqqeG2PYvXs3dmz4oxnHIYGwwacj9Ho5q1M6+We059ml2wiGwnV/QKkolJCQwD333MPQoUOZMGECvXv3BmD8+PEsXLiw6rhp06Yxbtw4Ro0aBcCKFSsYNmwYPl/TF9jpEj2HBENhEnWEXm9T8/1Me3YFH6zfw7h+XdwOR6lGueOOO7jjjju+8fovf/lLDh8+TNu2bbn99tu5/fbbq9579tln+dGP7Kkj0YzjkEBIR+gNMfrsTmRnpDBzSaHboShluwceeIDt22ufUuzbty+jR4+25Tqa0B2ic+gN4/UINwzL5dOtB9lQbG8NsVJuGzp0KP3796/1ve9///u2XUczjkMCobCuFG2ga87rRnKCh1k6SleqUTTjOCQYMiT4dMqlIdqmJnLFoK68vnonh7SEUakG04TukEAojE9H6A02Jd9PWSDM3IIddR+slPo3mnEcEtCVoo3Su3NrhvVox7NLtxHSXvIqxhUXFzNhwgQAVq9ezfz586vee/vtt/nNb35j6/U04zhEm3M13tR8PzsPn+CfG/a4HYpSTfLggw9W3fSsmdDHjx/Pm2+++W+dF5tKE7pDgmEdoTfWmLM7kdUmmZmLC90ORal6+a//+i/++te/Vj2/++67mTFjBvPmzWPs2LFUVFRwzz33MHfuXAYOHMjcuXMREUaOHFm1otQOurDIIRXBsNahN5LP6+GG4X7+9O6XbNxdQq/O6W6HpGLJO7+A3V/Ye87O/WDcH0/59s0338zEiRO58847CYfDzJkzh5deeolZs2aRlJQEwL333ktBQQGPPPJI1efy8vL4+OOPufrqq20JU4eQDgmGtWyxKSad140kn4dZSwvdDkWpOvn9ftq3b8+qVat4//33GTRoEOXl5fXusmgXHaE7JKBli02SkZbIdwZ25bWVO/n5Jb1pk5rgdkgqVpxmJO2kW265hZkzZ7J7926+973vNWuXxUo6hHSIli023ZR8PycCIV7SEkYVA6644greffddli9fziWXXMJZZ51FYWFh1fs1uywCbNq0ib59+9oWg2YchwRCYRJ9+uNtinOyWjOkeztmLS3UEkYV9RITExk1ahRXX301Xq+XtLQ0zjjjjKrWuaNGjWL9+vVVN0UBFixYwPjx9m0krxnHIcGQwefRKZemmprvp+jQCT78cq/boSh1WuFwmGXLlnHzzTdXvTZ9+nRmzpwJQLt27Vi+fDmrV6/mmmuuYc+ePZw4cYJ+/frZFoMmdAcYY7Rs0SbfOqcTXdoka38XFdXWr1/PmWeeyejRo+nZs2fV61dcccUpN67Yvn07DzzwgK1x6E1RBwRC1vSALixqOp/Xw/XDcvnf9zby1Z4SenbSEkYVfc455xy2bNlS63u33HJLra+fd955tsehQ0gHBCP7CuoI3R6Th+SQ6PNor3R1WtZe9fGjMf89mnEcEAha/yN8mtBt0S4tkcsHZPHqyp0cORFwOxwVhZKTkzlw4EDcJHVjDAcOHCA5OblBn9MpFwcEIiP0RJ1ysc2UfD8vryji5YId3HJ+D7fDUVEmOzuboqIi9u3b53YotklOTiY7O7tBn9GE7oBAZKNjHaHbp2/XNuTlZjB76TZuGtEdr1YQqWoSEhLo3r2722G4TjOOA4JVN0X1x2unqSP8bD94nIUbtYRRqdpoxnFARajypqiOIu10SZ/OdG6drDdHlToFTegO0BG6MxK8Hq4flsPHX+1n895St8NRKupoxnFA1Ry6zvPabtKQHBK9HmZrF0alvkETugMqE3qCXb1cghVQstuec8W4Dq2SuHRAFq+sKOJomZYwKlWdJnQHVK0UbWq3xfJSWPII/HUAPNAL5l4PezfYEGFsm5rv53hFiJcLitwORamoogndAcGqssVGTrmEgrDwT/BQX3j/bmh/Boy4E75eCI8Oh3nfh4O1LzNuCfplt2FwTlueXVpIWLswKlVFE7oDTla5NPLHu3ImLPw95AyHm/8JU9+Gi++Fuz63EvuGt+D/RsLOlbbFHGumjuhO4YHjfLQpfhaSKNVUmtAdEGxKc65wCJb+DbrmwaQXoFu1Bj6p7eDi38Ftn0JKG3j2Cij+3KaoY8u4vp3pmJ7EM1rCqFQVTegOCDRlhL5xvjWdkn87yCn+QcjIhSlvQWIrmH057FnfhGhjU4LXw3VDc1m0aR9f79MSRqVAE7ojAuEmjNCXPAJtc+HsS09/XIYfprwJviSYdSns29jwa8W4a4fmkOAVnl26ze1QlIoKmtAdEAg2coS+4zPYsQyG3wYeb93Htz/DGql7vPD8lXDsQCOijV2Z6UlM6J/FywU7KNESRqXqTugi8rSI7BWRtad4f6SIHBGR1ZHHPfaHGVsq+6E3uDnXkochuS0MvK7+n+nQEya/CCV74OUpEGpZiW1qvp9jFSHmrdASRqXqk3FmAmPrOOZjY8zAyOPepocV2xq1Y9HBLVb1St73IKlVwy7Y9Vy4bAYUfgzv3d2wz8a4Ad3aMrBbW2Yt3aYljKrFqzOhG2MWAQebIZa4UXVTtCELi5Y+Ch4fDP1B4y46YBIMnw6f/R+sfLZx54hRU/P9bN1/jEVfaQmjatnsmkMfLiJrROQdEelzqoNEZJqIFIhIQTw1oq+pqmyxvkv/K47B6ueh/9WQ3rnxFx7zO+gxCv7xY2s+voX4dr8uZKYn6UbSqsWzI6GvBHKNMQOAh4HXT3WgMeZxY0yeMSYvMzPThktHp4qGNufa/QUEjtdd2VIXrw+ufBpad4WXbrTm1VuARJ+Ha4fksGDjPrbuP+Z2OEq5pskJ3Rhz1BhTGvl+PpAgIh2aHFkMa3D73OI11tcuA5t+8dR2cM1zUHakRd0kvW5oDj6PaBdG1aI1OaGLSGcRawWMiAyJnLNl1c/VEAiF8Qj13yZt12pI69i06ZbqOveFyx6G7UtbzE3Sjq2TGd+/C68UFFFaHnQ7HKVcUZ+yxReBpUAvESkSkZtF5FYRuTVyyJXAWhFZA8wAJpl42Xq7kQLhcMNq0IvXQJcBp14Z2hj9roRht1k3SdfMse+8UWxKvp+S8iCvrtQSRtUy1blJtDFmch3vPwI8YltEcSAQNPVP6IETsO9L6DXO/kAuvhd2fw5v3QmZvSBrkP3XiCKDurVlQHYbZi0p5IZhuYid/0AqFQN0pagDguFw/Vvn7lkPJmSN0O3m9cGVz0BaJrx4LRwttv8aUUREmJLv5+t9x/hk8363w1Gq2WlCd0Ag1IApl+LV1lcnEjpAq0yYPMe6STpnMlQcd+Y6UWJ8/y50aJXIzMWFboeiVLPThO6AQMiQUN8bosVrrOX+bXOcC6hzX/juk9bN1zd+BHF8iyPJ5+XaITl8uHEv2w5oCaNqWTShOyAQCtd/UZETN0Rr0/vbVi/1da/Bwj86ey2XXTcsF69oF0bV8mhCd0AwZOq3qChYAXvXOzfdUlP+HVbjr4/+CCtnN881XdCpdTJj+3ZmbsEOjmkJo2pBNKE7oKK+c+j7voRQBWQNdDwmwPotYMJDcOYYq/Jl/ZvNc10X3DTCT0lZkNdW7XQ7FKWajSZ0BwTrm9DtXCFaX75EuPpZyD4P5t0MXy9ovms3o8E5GfTt2ppZSwpp4csiVAuiCd0BwbCpX+vc4tWQmA4Z3R2P6d8kpsK1c6F9T5hzHRQVNO/1m4GIMDW/O1/tLWXJ1y164bJqQTShO6AiGK7f5hbFa6BLf2hIm127pGTADa9aZY3PTozL7owT+nehXVoiM7ULo2ohNKE7IBg2JNaV0ENB2L22+W6I1ia9s7WFXVp7mP0d2LLQvVgckJzgZfKQbvxrwx52HIzv+nulQBO6IwKheqwUPfAVBE+4m9DBqn+/6V1r0+nnr4Iv/+FuPDa7PtIC4NllWsKo4p8mdAcEQvXo5VJ1Q9TlhA6Q3gmmvg2d+8HcG2D1C25HZJsubVIY26czcz7bzvEKLWFU8U0TugOspf91jNCL14AvxboxGQ1S28GNb4D/P+D1H1ptd0PxkQCn5Ps5Whbk9VW73A5FKUdpQndAvcoWi9dYS/K9dTa8bD5J6XD9PBh6Kyx9BJ6bCMdjfzvZ8/wZnNNFSxhV/NOE7oBAyOCrq3Jl/yboeHbzBNQQ3gQY9ye4/FHYvgwevxB2rnQ7qiaxShj9bNxTwrItsf8PlFKnogndAYFQmETfaaZcQgE4ts/a+zNaDboObnoHwiF4cgx8eJ/VqiBGXTYwi4zUBGYu2ep2KEo5RhO6AwKh8OlH6KWRzZvt2nLOKdnnwg+XQP+rYdH98MQoKP7c7agaJTnBy6QhOXywfg9Fh7SEUcUnTegOCIbM6csWS3ZbX9O7NE9ATZHSFq54DCa9CKV7raT+7i9jcm79+mG5AFrCqOKWJnQHVITCp19YdDRSbRHtI/Tqen8bbvsUBl4Lnz4GMwbB0r/F1DRM17YpXNKnM3OX76AsEHI7HKVspwndAcFwfUfoWc0TkF1S28FlD8MPPrb2J33vV/C382D5k9beqDFgSr6fw8cDvLFauzCq+KMJ3WbhsCEUrmNhUUkxeHyQ2r75ArNT575ww2tw3SuQ0g7+8RP4Sx9r44xj0d0Ia2j3dvTunM4zi7WEUcUfTeg2C4TDAHUk9N3QqrM7TbnsIgI9L4bvfwhT51vteBf+AR7oBXOvt1oIROF0TGUJ45e7S/hsa+zdB1DqdGI4o0SnQMga9Z12pWhJcWzNn5+OCPhHWO14b/sMhkyz6tfnXAsP9oa37oJN70XVlMzlA7vSJiVBuzCquBNFyxTjQzBkjdBPW7ZYUgwdomTJv50ye8HY31t7l27+F6x5ET5/CVY8Y7U56HEh9BgJOcOtvjEerythpiR6mTSkG09+vJVdh0+Q1TbFlTiUspsmdJtVjdBPt0l0STF0v6CZInKBNwF6jbUewXIo/MQapW9613qAtbFHzlDoNQ7OvgxadWzWEG8YlssTi7bw3LJt/Gxs72a9tlJO0SkXmwUiI/SEU20SXXEcyo7Ez5RLXXxJcOZo+Pb9cNfn8J/rYOKT1mKlQ4XWDdUHesHMCbDqOWimG5XZGalcfE4nXvxsu5YwqrihCd1mwao59FP8aEtjaFGRE9pkQ/+rYMKDML3AWol6/v+zfmt54zb46P5mC2VKvp9DxwO8uUa7MKr4oAndZhWVc+inuil6tNj62lJG6KcjAp36wEV3W8l9wGRY+HtY/WKzXH54j/b06pSuXRhV3NCEbrNgpGzxlCtFSyoTeowtKnKaCFw6w7q38OZ02PJRM1xSmJLvZ92uoxRsO+T49ZRymiZ0mwWC1kjvlJtEV60S1RH6N/gS4epnrU0/5t4Aezc4fsnvDMqidbKPmYsLHb+WUk7ThG6zkwuLTjHlUlJslfAlt2nGqGJISlu47mVISIEXroGKY45eLjXRxzXndePddbspPhI9tfJKNYYmdJsFgnWsFK1cVCR1bFHXkrXtBlc9A4e3We0EHHbjcD9hY3h+2XbHr6WUkzSh2ywYrqPKpWQ3tNb58zrl5sPgG62Ojru/cPRS3dqlMrq3ljCq2KcJ3WZ1VrnE07J/p435HaRkWO0DIlNZTrlphJ8Dxyp4+/NiR6+jlJM0odusqg69tqX/xlgj9JZag95Qqe3gkvtgZwGseNrRS+Wf0Z6eHVtpCaOKaXUmdBF5WkT2isjaU7wvIjJDRDaLyOciMtj+MGNH1UrR2vYULT8KgeM6Qm+I/tdYpYz/vPdkhZADRIQb8/18sfMIK7drCaOKTfUZoc8Exp7m/XFAz8hjGvD3pocVuwKna85VtahIR+j1JgLj/wLBE/D+rx291MRBXUlP9jFziW5Rp2JTnQndGLMIOF3j6MuB2cayDGgrIi02Y1U256p1YVGJJvRG6XAm5N8OX7wMO1c6dpm0JB9X53XjnS+K2XO0zLHrKOUUO+bQuwI7qj0virz2DSIyTUQKRKRg3759Nlw6+gRPd1NUFxU13oi7rB2ePrjH0QZeNw7PJWQMz+tG0ioG2ZHQayvnqPVvnDHmcWNMnjEmLzMz04ZLR5+qOfTTjtA1oTdYcmu48OdQ+DFs/qdjl8ltn8ZFvTrywmfbKQ9qCaOKLXYk9CKgW7Xn2UCLbV932h2LSoohqQ0kpjVzVHHi3Jsgozt88BsIO5dsp+T72V9awT+0hFHFGDsS+pvAjZFql2HAEWNMi/2bEDzdnqJag940vkQYfQ/sXQdr5jh2mfN7duCMzDRm6RZ1KsbUp2zxRWAp0EtEikTkZhG5VURujRwyH9gCbAaeAH7kWLQxoHKEfso59NZ6Q7RJ+lwBWYNhwX2O7VNa2YVxTdERVmkJo4oh9alymWyM6WKMSTDGZBtjnjLGPGaMeSzyvjHG3GaMOcMY088YU+B82NHr5I5FtY3QdVFRk4nAxffC0Z3w6f85dpmJg7NJT/LpRtIqpuhKUZsFQmG8HsFTcwu6cDiS0HXKpcm6nw9nXgyLH7K283NAqyQfV+ZlM/+LYvaWaAmjig2a0G0WDJnab4gePwDhgI7Q7XLRr+HEIVj6qGOXmDLcTzBseOFT7cKoYoMmdJtVhMKnmG7RRUW2yhoIZ19mdWM8frp1b43n75DGyLMyeW7ZdiqCzjYHU8oOmtBtFgwZEnynmD8HTeh2GnU3VJTCJ39x7BJWCWM576xtsYVbKoZoQrdZIBTGV3P+HHRRkRM69ob+V8NnTzjWuOuCnpn06JCmN0dVTNCEbrNAyJx+lWirTs0bULwb+Qvr3sSiPztyeo9HuHF4Lqu2H2bNjsOOXEMpu2hCt1kgFD71KtHUDtbiGGWfdj1g0PWwYiYccqb/ynfPzSYt0asLjVTU04Rus2A4jK+2EXrpPh2dO+WCn4F4YNH9jpw+PTmBq/K68dbnu9hXUu7INZSygyZ0m1UETzHlUnYEkts0f0AtQZuucN7NsPpF2L/ZkUvcODyXQMjw4mdawqiilyZ0mwXDp5hyKT9idQxUzviP/wRfEiz8gyOn75HZigvOyuS5ZduqVgMrFW00odvMmkOvbYR+FJI0oTumVUcYeiusnQd71jlyiZvy/ewtKeedtc5thadUU2hCt1kgZGovWyw/qiN0p424w/pHc8HvHTn9hWdl4m+fyszFWx05v1JNpQndZsFQmMSaC4uM0RF6c0jJgPzp8OXbsHOF7ae3Shj9rNx+mC+KnOkho1RTaEK3Wa0j9MBxMCEdoTeHYT+ElHbw4X2OnP7KvGxSE7260EhFJU3oNqt1Dr3sqPVVR+jOS0q3bpB+/S8oXGz76VsnJ/Ddwdm8tWYX+0u1hFFFF03oNqs1oZdHErqWLTaPId+3eub863eObCg9JT+XilCYOVrCqKKMJnSbBcO1tM/VEXrzSkixNpTe8Slses/205/ZMZ3ze3bguWXbtYRRRRVN6DYLBGtZKVoeuYGmc+jNZ9D1VluAD//b2lzEZlPz/ew+WsZ767SEUUUPTeg2C4RrWSmqI/Tm502w2uvuWWvVpttsZK+O5LRL1f4uKqpoQrdZrc25qubQNaE3qz4ToXM/WPA/EKyw9dTeSBfG5YWHWLtTSxhVdNCEbrNgbe1zdYTuDo8HLroHDhXCqtm2n/6qvG6kJGgXRhU9NKHbrCIUxlfrCF0gsZUrMbVoPS+GnOHw0f1QcczWU7dJSWDi4K68sWYXB4/Z+xuAUo2hCd1mwdr2FK1cJVrbXqPKWSIw5ndQuseRDaWn5vupCIa1C6OKCpphbBQKG8KG2uvQdf7cPTlDofcEWPxXqy+9jXp2SmfEme15ftk2glrCqFymCd1GlTXJ35hy0T4u7hvzW6sFgwObYEwZ7mfXkTI+WL/H9nMr1RCa0G1UmdATdYQefTr0hHOnQMHTcOBrW089+uxOZGek8IzeHFUu04Ruo2DIWmb+zRG67lYUFS78BXiT4F/32nrayhLGz7YeZP2uo7aeW6mG0IRuo8oReq1z6Drl4r70TpB/O6x/HYoKbD311XndSE7waAmjcpUmdBsFwtYIvdZeLjrlEh3yp0NaR3j/17Y27mqbmsgVg7J5ffVODmkJo3KJJnQbBWsboRujI/RokpQOo34F25fC+jdsPfWU/FzKg2HmLN9h63mVqi9N6DY6WeVS7ccaOAHhoI7Qo8ngG6FjH/jgHgiU2Xba3p1bM7xHe57TEkblEk3oNgpEboomVp9yKddl/1HH44VL7oPD2+DTx2w99ZR8PzsPn+CfG/bael6l6kMTuo2qRujVV4RW9nHRKpfocsYoOGssLPqzrYuNxpzdka5tU/TmqHKFJnQbVY7QE6pvEq0j9Oj1rf+B4AlYYN/+oz6vhxuG57J0ywG+3K0ljKp5aUK3UVXZYvVNost0c4uo1aEnnHcLrJwFe9bZdtpr8rqR5PMwa8k2286pVH1oQrdRUEfosefCn1v/b979hW1ljBlpiVwxqCuvrSri8HEtYVTNp14JXUTGishGEdksIr+o5f2RInJERFZHHvfYH2r0OzmHXn2ErptbRLXUdnDRr2HrIlvLGKfk+ykLhHmpQEsYVfOpM6GLiBf4GzAOOAeYLCLn1HLox8aYgZGHvWurY0StK0V1hB79zr0JOvW1FhtVHLfllGd3ac2Q7u2YvXQbobB9C5iUOp36jNCHAJuNMVuMMRXAHOByZ8OKTVU3Rb01q1x0c4uo5vXBuPvhyA745C+2nfamfD9Fh07wrw3ahVE1j/ok9K5A9d8biyKv1TRcRNaIyDsi0qe2E4nINBEpEJGCffvs7UsdDYLhWtrnVq4S1c0topt/BPS7yuqZfnCrLae8+JxOZLVJZtbSQlvOp1Rd6pNlpJbXav4OuRLINcYMAB4GXq/tRMaYx40xecaYvMzMzAYFGgsqgrW0z9U+LrHj4nvB44P37rbldD6vh+uH57J48wE27Smx5ZxKnU59EnoR0K3a82xgV/UDjDFHjTGlke/nAwki0sG2KGNEMFxL+1zt4xI7WmfBhT+Fjf+Arz6w5ZSTzssh0addGFXzqE9CXw70FJHuIpIITALerH6AiHQWEYl8PyRy3gN2Bxvtar0pWnZER+ixZNiPoH1PmP//rD48TdQuLZHLB2Tx6sqdHDkRsCFApU6tzoRujAkC04H3gA3AS8aYdSJyq4jcGjnsSmCtiKwBZgCTjLGxN2mMqLop6qlR5aIj9NjhS4LxD8ChQvj4QVtOOSXfz4lAiJe1hFE5rF536owx840xZxljzjDG3Bd57TFjzGOR7x8xxvQxxgwwxgwzxixxMuhoVTVC99WoQ9cRemzpcSH0uxoWPwT7Nzf5dH27tuE8f4aWMCrHaemFjYK1NefSEXpsuuQ+8KXA/J/YsoJ0Sr6f7QePs3CjdmFUztGEbqOTdeiREboxOkKPVa06wuj/gi0LYe28Jp/ukj6d6dw6mZl6c1Q5SBO6jQKhMD6PELk/DMEyCAd0hB6r8r4HWYPgvV+dbLLWSAleD9cPy+Hjr/azea+WMCpnaEK3UTBsalklio7QY5XHCxP+Asf2wT9/2+TTTR6SQ6JXuzAq52hCt1FFMPzNGnSAJN3cImZlDbJKGQuehm1Lm3Sq9q2SuHRAFvNWFnG0TEsYlf00odsoGA5/c5Uo6Ag91o38JbTJgbfugGB5k041Nd/P8YoQrxQU2RScUidpQrdRIGhqjNAj8646hx7bklpZUy/7NzW5eVe/7Dacm5vB7KWFhLWEUdlME7qNAuGwzqHHq55jrOZdHz8A+zY26VRT8v0UHjjOR5vir0GdcpcmdBsFQkZ7ocezS/4AiWnw5h0Q6azZGOP6dqZT6ySe0RJGZTNN6DYKRsoWq+gIPb60yoRLfg87lsHyJxt9mgSvh+uG5rJo0z6+3ldqY4CqpdOEbqNAKFzLCF0gMd21mJTNBkyGM8dYZYyHCht9msoSxmeXagmjso8mdBtZUy41RuhJ6bq5RTwRgQkPgXisqZdGtgXITE9iQv8uvFywgxItYVQ20Uxjo1pH6Dp/Hn/adoNv3QtbP4KVsxt9min5fo5VhJi3QksYlT00odsoGKpRtqi90OPX4KngP9/aWPrIzkadYkC3tgzKacuspdu0hFHZQhO6jSp0hN5yeDxw2QwIB+GtOxs99TI138/W/cdY9JWWMKqm04Ruo2Btdeg6Qo9f7XrAmN/C5g9g1XONOsW4vl3ITE/SLeqULTSh2ygQrHFTVEfo8e+871tTL+/+Eg43fEeiRJ+H64bmsGDjPrbuP+ZAgKol0YRuo0A4jE9H6C2LxwOXPwIYeHN6o6Zerh2aQ4JXmL200PbwVMuiCd1GwZD59+ZcOkJvGTL88K3/tjbDKHiqwR/vmJ7M+H5deKWgiNLyoO3hqZZDE7qNAtVXigbKIFShI/SW4tyb4IyL4P174ODWBn98Sr6fkvIgr67UEkbVeJrQbRQIGRJ8kR+p9nFpWUTgsoetTTFe/yGEQw36+KCcDAZkt2HmEu3CqBpPE7qNAqEwCZUj9Ko+Lrq5RYvRJhvG3Q/bl8KShxv88akj/GzZd4xPNu93IDjVEmhCt1Gweh269kJvmQZMgrMvhQ//B3Z/0aCPfrtfFzq0StQSRtVomtBtFAiZk1Uu2mmxZRKBCX+FlAx49QcN2uEoyefl2qG5fLhxL9sOaAmjajhN6DYxxhAIh0msrEPXOfSWK629Vcq4d501Um+A64bm4BVhtnZhVI2gCd0mobDBGHSErixnXQLnTrXm0gsX1/tjnVonM65fF14q2MExLWFUDaQJ3SbBSGWCT0foqtK37rNq1F/7AZw4XO+PTc33U1IW5NVVjWv6pVouTeg2qQhZW5JVLSzau97a2EITesuV1Aq++yQc3QXzf1rvjw3OaUu/rm2YvaQQ08imX6pl0oRuk2AoMkL3CAQrYMNbcPYE3dyipcvOgwt/Dl+8BF+8Uq+PiAhT8v18tbeUJV8fcDhAFU8029gkEBmhJ/g88PWHVi/0vt91OSoVFc7/CWQPgbd/DIe31+sjE/p3oX1aIs8sLnQ2NhVXNKHbpCqhezywdp5VttZjpLtBqejg9cHEx8GE4LVb67WKNDnBy+QhOfzryz3sOHi8GYJU8UATuk0CkSmXZMpg43w453LwJrgclYoa7bpbq0i3LYbFD9XrI9cPy8Uj2oVR1Z8mdJsEIyP0rH2fQEWpTreobxp4LfSZCB/eB0UFdR7euU0yY/t2Zu7yHRyv0BJGVTdN6DapHKFn73wHWnWC3BEuR6SijghM+Au0zoJ5N0N5SZ0fmZrv52hZkNdX7WqGAFWs04Ruk0AoTCuOk7l7IfS5wuq6p1RNKW2t+fTD22H+z+o8PC83gz5ZrZm5ZKuWMKo61Suhi8hYEdkoIptF5Be1vC8iMiPy/uciMtj+UKNbMBxmjGcl3lC5Treo08vNhwt+CmteqLOUsbKEcdOeUpZu0RJGdXp1JnQR8QJ/A8YB5wCTReScGoeNA3pGHtOAv9scZ9SrCBou9S6lLK0rZJ/ndjgq2l3ws0gp43/CocLTHnrZgCzapWkXRlU3Xz2OGQJsNsZsARCROcDlwPpqx1wOzDbW74TLRKStiHQxxhTbHfCOzV+wc/lbdp+2yQ6VlnGx53P2d7+FLiJ1f0C1bF6ftYr0sfNh3i1w0zunrIpKTvAy6bxuPPbR1zyxaMu/b0SuGsaE8YYr8IQDeE0ATziAxwTxhIN4TAAxITwmFPkaREw48ghZXwmffI0wGHPyNUzVc+taxnoNU/UegGBo12MQZ+eNsv0/rz4JvStQfTvzImBoPY7pCvxbQheRaVgjeHJychoaKwB7NxUwbOOfGvVZpwXwIgMmuR2GihUZuXDpQ/DKTbDwDzD6nlMeev2wXGYuKeS++RuaL76oYWjFCTrIEdpRQlsppQ3HrK9yjFacsB5yglaUkSLlpFJGChUkSwXJWI8kAiRIw3aScsrSQzeCSwm9tuFAzbsz9TkGY8zjwOMAeXl5jbrD0/eiazg8ZFxjPuq4xOQUOrfSHYpUA/SdCFsWwMcPQvcLoceFtR6W1TaFFb++mLJAdCQkWwXL8RzdgefwVryHC/Ec3YmnZBeekp1I6W48x/YiwbJTftwkpGGS0jGJraxHQntMQiokpGISUsCbjPElE/QlEfQmYXyJ4E0CTwLGmwheH8aTaBUyeBLA48N4PCC+yGtejHhBvCCeSMGDWN+LByMeq4Ip8tx6D+srcvI9rNeMCAPSnMkT9UnoRUC3as+zgZo1VPU5xhZJyakkJac6cWql3DH2j7B9Gbw6DX64GNI61HpYSqKXlMQYrp4Kh+HAZtj9udW8bu8G6+uhbfzb+M+bBG26QuuukDscWnWEtExIi3xNybCqhVIyILkN4vHWOqJsieqT0JcDPUWkO7ATmARcW+OYN4Hpkfn1ocARJ+bPlYpLiWlw5dPwxEXWBtPXvmSN6mJdean1D9W2xbBzBexadbKttMcH7XtC1iDofw206wEZ3a0VtWmZ8fHf74I6E7oxJigi04H3AC/wtDFmnYjcGnn/MWA+8G1gM3AcuMm5kJWKQ537Wf3T3/kpLP0b5E93O6KGMwZ2rYRN78GWj2BnAYSDVvLu1Af6XQVdz4WsgVYy9yW6HXHcEbcWK+Tl5ZmCgrqXPyvVYhgDc6+3EuLN71nJL9qFQ1D4sdUu+sv5ULLLmi/uMtC6H9D9Aug2DBJ1mtQuIrLCGJNX23v1mXJRSjUHEbjsYauU8ZXvwQ8WQXKU3mTftxFWvwCfv2QlcV8KnDkaet9jbb+X2s7tCFskTehKRZPUdnDlU/DMt+Gtu6y59WiZTw4FYP0b8OljULTcqvo4cwyM/T30vERH4VFAE7pS0SZnGIz6FXz439a0xblT3Y3n2AFY8TQsfwpKiqHdGdZ8f/+rrQoUFTU0oSsVjf7jx9bc9Ds/t1oEdKrZbaMZlO6DJTOsRB44BmeMhktnWKNy3VoxKmlCVyoaeTww8Qn4+wh4eSpMW2CVNzaH0r2w+K9WIg+VQ98rrW30OvZunuurRtOErlS0atURvvsEzP4OzP8pfOdRZ69XXgpLH4ElD0PghFUffv5PoMOZzl5X2UYTulLRrMdIq9XuovvBfz4MnGz/NUIBWDkLFv4Jju21tk+86B5N5DFIE7pS0W7kL2DbEvjHT6DrYMjsZd+5v/oA3vsV7N9k7bI1+UXIrrXEWcUAvbOhVLTzeK1Wu4mp8MI1cPxg08+5dwM8OxGev9JazTnpBZj6D03mMU4TulKxoHUXK+ke3QVzb4BgRePOc2QnvDEd/p5vLc2/5Pfwo0+h9/joqXdXjaYJXalY0W0IXP4IbPsE/vHjqg0T6uXEYfjgN/DwYPh8Lgy9FW5fBcNv054qcUTn0JWKJf2vtua7F/2vNZeef/vpj9+/GT573FqmX1FqVa6M+pW1uYaKO5rQlYo1IyM3Md//Nax73apKOecyyPBD2RHYtwn2bbDe+/pf1qYNfSdayb9zP5eDV07SbotKxaLACfj0/2Dda1C82nottQMc33/ymPQukPc9q3WALtGPG9ptUal4k5AC/3GX9Ti4FTa8aXVAbH8mZPa2VnW2zY1sl6ZaCk3oSsW6dt1hxJ1uR6GigFa5KKVUnNCErpRScUITulJKxQlN6EopFSc0oSulVJzQhK6UUnFCE7pSSsUJTehKKRUnXFv6LyL7gG2N/HgHYH+dRzW/aI0Lojc2jathNK6Gice4co0xmbW94VpCbwoRKThVLwM3RWtcEL2xaVwNo3E1TEuLS6dclFIqTmhCV0qpOBGrCf1xtwM4hWiNC6I3No2rYTSuhmlRccXkHLpSSqlvitURulJKqRo0oSulVJyIuYQuImNFZKOIbBaRX7gdD4CIPC0ie0VkrduxVCci3URkgYhsEJF1IhIVuyCISLKIfCYiayJx/c7tmKoTEa+IrBKRt92OpZKIFIrIFyKyWkSiZu9GEWkrIq+IyJeRP2fDoyCmXpGfU+XjqIjc5XZcACLyn5E/82tF5EURSbb1/LE0hy4iXmATcDFQBCwHJhtj1rsc1wVAKTDbGNPXzViqE5EuQBdjzEoRSQdWAN+Jgp+XAGnGmFIRSQA+Ae40xixzM65KIvJjIA9obYyZ4HY8YCV0IM8YE1WLZERkFvCxMeZJEUkEUo0xh10Oq0okZ+wEhhpjGruQ0a5YumL9WT/HGHNCRF4C5htjZtp1jVgboQ8BNhtjthhjKoA5wOUux4QxZhFw0O04ajLGFBtjVka+LwE2AF3djQqMpTTyNCHyiIqRhYhkA+OBJ92OJdqJSGvgAuApAGNMRTQl84jRwNduJ/NqfECKiPiAVGCXnSePtYTeFdhR7XkRUZCgYoGI+IFBwKcuhwJUTWusBvYCHxhjoiIu4CHgZ0DY5ThqMsD7IrJCRKa5HUxED2Af8ExkiupJEUlzO6gaJgEvuh0EgDFmJ/BnYDtQDBwxxrxv5zViLaFLLa9FxcgumolIK2AecJcx5qjb8QAYY0LGmIFANjBERFyfqhKRCcBeY8wKt2OpxQhjzGBgHHBbZJrPbT5gMPB3Y8wg4BgQFfe1ACJTQJcBL7sdC4CIZGDNKHQHsoA0EbnezmvEWkIvArpVe56Nzb+yxJvIHPU84HljzKtux1NT5Ff0hcBYdyMBYARwWWS+eg5wkYg8525IFmPMrsjXvcBrWNOPbisCiqr9dvUKVoKPFuOAlcaYPW4HEjEG2GqM2WeMCQCvAvl2XiDWEvpyoKeIdI/86zsJeNPlmKJW5ObjU8AGY8yDbsdTSUQyRaRt5PsUrD/oX7oaFGCM+aUxJtsY48f6s/WhMcbWEVRjiEha5KY2kSmNbwGuV1QZY3YDO0SkV+Sl0YCrN9xrmEyUTLdEbAeGiUhq5O/maKz7Wrbx2XkypxljgiIyHXgP8AJPG2PWuRwWIvIiMBLoICJFwG+MMU+5GxVgjThvAL6IzFcD/MoYM9+9kADoAsyKVCB4gJeMMVFTIhiFOgGvWTkAH/CCMeZdd0OqcjvwfGSAtQW4yeV4ABCRVKxquB+4HUslY8ynIvIKsBIIAquwuQVATJUtKqWUOrVYm3JRSil1CprQlVIqTmhCV0qpOKEJXSml4oQmdKWUihOa0FXcE5H21Trv7RaRnZHvS0XkUbfjU8ouWraoWhQR+S1Qaoz5s9uxKGU3HaGrFktERlb2PBeR34rILBF5P9J7fKKI3B/pQf5upIUCInKuiHwUaZL1XqRFsVJRQRO6UiedgdU693LgOWCBMaYfcAIYH0nqDwNXGmPOBZ4G7nMrWKVqiqml/0o57B1jTEBEvsBqLVG5vP4LwA/0AvoCH0SW4Xux2qAqFRU0oSt1UjmAMSYsIgFz8gZTGOvvigDrjDGub7OmVG10ykWp+tsIZFbumykiCSLSx+WYlKqiCV2peopse3gl8CcRWQOsxuZ+1ko1hZYtKqVUnNARulJKxQlN6EopFSc0oSulVJzQhK6UUnFCE7pSSsUJTehKKRUnNKErpVSc+P/cg11fHWGf7AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from gekko import GEKKO\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Create GEKKO model\n",
    "m = GEKKO(remote=False)\n",
    "\n",
    "# Time points for simulation\n",
    "nt = 81\n",
    "m.time = np.linspace(0,8,nt)\n",
    "\n",
    "# Define input\n",
    "# First step (up) starts at 1 sec\n",
    "# Ramp (down) starts at 3 sec\n",
    "# Ramp completes at 5 sec\n",
    "ut = np.zeros(nt)\n",
    "ut[11:31] = 2.0\n",
    "for i in range(31,51):\n",
    "    ut[i] = ut[i-1] - 0.1\n",
    "\n",
    "# Define model\n",
    "u = m.Param(value=ut)\n",
    "ud = m.Var()\n",
    "y = m.Var()\n",
    "dydt = m.Var()\n",
    "m.Equation(ud==u)\n",
    "m.Equation(dydt==y.dt())\n",
    "m.Equation(dydt.dt() + 6*y.dt() + 9*y==5*ud.dt()+5*u)\n",
    "\n",
    "# Simulation options\n",
    "m.options.IMODE=7\n",
    "m.options.NODES=4\n",
    "m.solve(disp=False)\n",
    "\n",
    "# plot results\n",
    "plt.figure()\n",
    "plt.plot(m.time,u.value,label='u(t)')\n",
    "plt.plot(m.time,y.value,label='y(t)')\n",
    "plt.legend()\n",
    "plt.xlabel('Time')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}