{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Nyquist Stability\n",
"\n",
"TODO: Redraw diagrams using Python code.\n",
"TODO: Explanation for setting up animation function.\n",
"\n",
"By the end of this notebook, you should be able to:\n",
"\n",
"- [TODO Insert Learning Goals]"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# import libraries we need later\n",
"import control\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy as sp\n",
"\n",
"# https://stackoverflow.com/questions/49722298/drawing-a-nyquist-diagram-animation-on-python\n",
"from matplotlib import animation\n",
"from IPython.display import HTML\n",
"from matplotlib.patches import Arrow"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Nyquist plots and stability is an important topic in classical control, which is usually taught as a \"bonus\" section in undergraduate CHBE control courses. However, it is part of the regular curriculum in MECH/ECE undergrad control, so I consider it necessary to cover it here also in CHBE, at least for those who wish to have a strong background in control.\n",
"\n",
"To motivate the use of Nyquist, here are some of its strengths over other methods:\n",
"\n",
"- Computationally inexpensive to generate contours vs. root finding or solving Routh-Hurwitz arrays\n",
"- Contours and clockwise-encirclements are easy to visualize with fairly little practice\n",
"- Closed-loop stability can be determined regardless of multiple critical (phase cross-over) frequencies, or in the presence of RHP poles in L(s), two cases in which Bode may fail to provide the correct stability conclusions\n",
"\n",
"With that, let's look at the well-known, dreaded -1 point and see how it relates to closed-loop instability:\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# TODO replace diagram"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Above is an empty Nyquist plot with the blue line representing the imaginary axis, and the red arrow representing a vector from the origin to the point $-1$. Suppose that this vector represents a certain open-loop TF, $L(s)$, at some specific frequency $\\omega$. Using simple geometry, we can determine its magnitude and phase as $\\lvert L(j \\omega) \\rvert = 1 = 0dB$ and $\\angle L(j \\omega) = -180^o$, respectively.\n",
"\n",
"Recall from Bode that this combination implies $L$ is at its ultimate gain, or at the point of marginal stability. In other words, if any slight extra gain or delay is introduced to $L$, the closed-loop would become unstable. Equivalently, at this $-1$ point, the closed-loop TF $\\big( \\frac{L}{1+L} \\big)$ has a pair of poles located at $\\pm j \\omega_{180}$, which again implies marginal stability. In light of these realizations, as control engineers we strive to stay as far away from the $-1$ point as possible, so as to allow ourselves comfortably large gain and phase margins when designing controllers for our systems.\n",
"\n",
"However, at this point we still haven't justified mathematically why we're concerned with the point $-1$ (and not, for example, $-2$ or $-5$) and where the equation $Z = N + P$ comes from.\n",
"\n",
"To do this, let's discuss the following 2 concepts:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Contours \n",
"\n",
"Tracings or mappings of a rational expression on the complex plane, by cycling its argument (in this case frequency) between $- \\infty$ and $+\\infty$. For example, the expression $F(s) = \\frac{1}{s+1}$ has the spectrum $F(j \\omega) = \\frac{1-j \\omega}{1+\\omega^2}$ which can be evaluated at the following frequencies:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"| $\\omega$ | $F(j\\omega)$ |\n",
"|:-------------:|:---------------:|\n",
"| $-\\infty$ | $0+0^+j$ |\n",
"| $-1$ | $0.5+0.5j$ |\n",
"| $0$ | $1+0j$ |\n",
"| $1$ | $0.5-0.5j$ |\n",
"| $\\infty$ | $0+0^-j$ |\n",
"| | |\n",
"\n",
"*Note: Take a look at the Markdown code for the table above. The hidden `` in the last row is a [cool trick](https://stackoverflow.com/questions/36121672/set-table-column-width-via-markdown) to size the table.*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Knowing these points, the contour of F can be traced, and software such as MATLAB can do this relatively quickly:\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also do this easily in Python:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAFNCAYAAABST1gVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XecVNX5x/HPQ0dWpCOoFBUFbAgE\nRaOCgqJRsYsKYotRorErRuyJJTGmiRo0CmLBaCzYflgANcZCU6wIKiKCokhbQBD2+f1xZsMACztb\nZs6U7/v1mtfM3Htn53v2Muwz995zjrk7IiIiIpLdasQOICIiIiLlU9EmIiIikgNUtImIiIjkABVt\nIiIiIjlARZuIiIhIDlDRJiIiIpIDVLSJSF4zs7vN7OoMvl+xmW2fqfcTkcKhok1E0s7MZpvZt2bW\nIGnZWWY2Md3v7e7nuPuN5W2XyNhnM+t7mVlJoigrNrO5ZvYvM/vZBu9X5O6fV0d2EZFkKtpEJFNq\nARfEDlFF89y9CNgS2Bv4BHjdzA5K9xubWa10v4eIZDcVbSKSKX8ELjWzRhuuMLPhZvanDZY9Y2YX\nJh7vaWZTzWyZmT1qZmPM7HeJdaeZ2X82eK2b2Y6JxyOTtm1mZs+a2WIz+8HMXjezGmY2GmgDPJM4\ninb55hriwVx3vwa4F7h1E+/9CzObZmZLzewrM7tug5ynmtmXZrbQzK5OPtpnZteZ2eNm9qCZLQVO\nM7MeZvZmIv98M7vDzOps8N5DzGxm4nd1o5ntkHjN0sSRwTqISE5S0SYimTIZmAhcWsa6UcBJZlYD\nQnEFHAQ8kigyngJGA02Ax4BjK5nhEmAu0BxoCfyWUIMNAuYARyROb/6hAj/zCaBr8qnfJMuBU4FG\nwC+Ac83sKAAz6wzcCZwCtAK2ArbZ4PX9gccTr38IWAtcBDQDehJ+R0M2eE0/oBvhSODlwIjEe2wH\n7AqcVIG2iUgWUdEmIpl0DXC+mTVPXuju7wBLCEUIwABgort/Syg+agN/cfef3P1xYFIl3/8nQoHU\nNvGzXveqT8A8DzBCYbUed5/o7u+7e4m7TwceAQ5IrD4OeMbd/+Puqwm/mw2zvOnuTyVev9Ldp7j7\nW+6+xt1nA/9I+nmlbnX3pe7+IfAB8KK7f+7uS4AXgD2r2F4RiURFm4hkjLt/ADwLDC1j9ShgYOLx\nQMKRNYDWwNcbFFdfVjLCH4FZwItm9rmZlZWjorYhFFuLN1xhZnuZ2QQz+87MlgDnEI6SQWjXV6Xb\nuvsKYOEGP+Kr5CdmtlPi9O43iVOmNyX9vFLfJj1eWcbzopRbJiJZRUWbiGTatcAv2fhU4INAfzPb\nA+hEOCUKMB/Yxswsads2SY+XA1uUPjGzrTf1xu6+zN0vcfftgSOAi5M6EVT2iNvRwFR3X17GuoeB\nscB27r4VcDfhqByEdm2blLs+0HTDyBs8v4vQ+aGDuzcknN41RKQgqGgTkYxy91nAo8BvNlg+l3Da\nczTwb3dfmVj1JrAG+I2Z1TKzY4AeSS99D9jFzLqYWT3guk29t5kdbmY7JgrApYRrxNYmVn8LpDS+\nmgXbmNm1wFmE4qksWwI/uPuPZtYDODlp3ePAEWa2T+K6vespvwDbMpG72Mw6AuemkldE8oOKNhGJ\n4QagrAv3RwG7se7UKInrvY4BTgMWAScSLv4vXf9p4ue9DMwE1utJuoEOie2KCcXgne4+MbHuZmBY\nomdmWZ0lAFqbWXHi9ZMSWXu5+4ub2H4IcIOZLSNcs/avpNwfAucDYwhH3ZYBC4BVm8l/KaHwWwbc\nQyh+RaRAWNWvwRURqR5mtj/hNGk7dy/ZzHYjgbnuPixT2dLNzIoI18V1cPcvYucRkeyjI20ikhXM\nrDZh8N17N1ew5RMzO8LMtkgMF3Ib8D4wO24qEclWKtpEJDoz60Q4ytQK+EvkOJnUnzBkyDzCqdsB\n1TAEiYjkKZ0eFREREckBOtImIiIikgNUtImIiIjkgFqxA1S3Zs2aebt27dL+PsuXL6dBg7JGLMh/\nhdx2KOz2q+2F2XYo7PYXctuhsNufibZPmTLle3dvXv6WeVi0tWvXjsmTJ6f9fSZOnEivXr3S/j7Z\nqJDbDoXdfrW9V+wY0RRy+wu57VDY7c9E280s5Wn5dHpUREREJAeoaBMRERHJASraRERERHKAijYR\nERGRHKCiTURERCQHqGgTERERyQEq2kRERERygIo2ERERkRygok1EREQkB+TdjAgiIiLJFi6EsWNh\nzZr1l//859CpE3z/PTz55Mav690bdtwRFi6swz33bLz+4IOhbVuYMwfGjdt4/S9+Aa1bgzuYVU9b\npLCpaBMRkbyxbBm89Ra88Qb87GehcFqyBM44Y+Nt7747FG1z5sDZZ2+8/sEHQ9H21Vf1ueiijdc/\n/XQo2qZPL/v148eHou2xx+Cqq2DffdfdOnaEGjrXJRWkok1ERHJaSQlceCG8/noooEpKQkF0xRWh\naGvfHmbOhPr1139do0bhftddYe7cjX9u48bhvnPnpWWub9Ik3B90UNmvb9Ys3DdtCrvsAs89B6NG\nrXvtxx9DixbwzTew1VYb5xPZkIo2ERHJCWvWwHvvhaNob7wBDRrAffeFAu2dd0JxNGxYOJK1997Q\nsGF4nVk4YrYpderANttsbr1vdn39+pt//UEHhZt7KB7feCO0o3nzsP6SS+Bf/4KuXdc/Grf11pv+\nmVKYVLSJiEhWWrECttgiPL74YhgxApYvD8+32w4OO2zdtm++mf3XjZnBTjuFW7IzzwzteeMNuPNO\n+POfw9G/998P659/Htq0gc6ddUq10KloExGR6Nzhyy/XHUV74w2YMQMWLQpHsnbcEU4/fd1RqO22\nW//12V6wbc6BB4YbwOrVMHVquDYPwqnek08O1+U1agQ9e4b2H3YY7LlnvMwSh4o2ERHJuJ9+CqcI\nd9opnMYcPhzOPz+s23LLUJwceyysWhWKtiFD4ubNlDp1wqndUmYwefL6xewLL4QjjnvuGY5GDhsG\n++wTirlWreJll/RT0SYiIhkzZw7cdlu4Fm35cnjiCTj6aOjbNxRu++4bTg3WrBk7aXYovR5vxx1h\n8OCw7Icf1g1fMmMG3HVXOKUK0KULXHYZnHAC1NJf+LyjXSoiImm3fDmcd14YRgPCKb9f/AL22y88\n33nncJPylfZahXC0bckSePdd+M9/4J//hFNOCT1me/aMl1HSQ0WbiIikzfz54ZTdFluEnpNDhoTe\nkm3axE6WP+rUgR49wu3CC+G119YVbNdcE043n3NOuJfcpn4oIiJSrdzhlVegT59wzdqiReE032uv\nwV//qoItnWrUgF69wmN3mDYNLr88/M6vvjrM/iC5S0WbiIhUi5ISeOqpcCF9nz7w4YehUKhdO6zX\ncBWZZQbPPBPGsOvdG373u1C8Pfpo7GRSWTo9KiIi1eKjj0Kngu23D1NEDR4M9erFTiU/+1no8PHx\nx3DrreuGCpkxIxyN69gxbj5Jnb73iIhIpaxcGXp8XnFFeL7rruG06IwZ8KtfqWDLNp06wciR6wb3\nvfrqMGDvccfBlClRo0mKVLSJiEiFLF4MN90UJks/77wwQXvpEBQHHqihJnLF8OHw29/Cyy9D9+5w\n8MFh/lbJXiraREQkZc89F4q1q66Cbt3g1Vdh4kQVarmoefNwnducOeG06fTp8NJLYV1JSbhJdlHR\nJiIimzV7NnzxRZgEdI89whRKU6eGkfn33z+3p5CSMCPF5ZfDF1+EewgdGHbfPYyrt3atdnC2UNEm\nIiJl+vBDGDQojMZ/5507ArDttvDII5r3Mh/Vrw9FReFxvXqhGB80CAYO7MGdd4ZrGCUuFW0iIrKe\nL7+Eo44KHQuefBIuuAAuv/yT2LEkgw45JMwNO3YsNG26ml//et2k9hKPrkIQEZH1jB4N48fDtdeG\nSdybNoWJE1fHjiUZVqMGHHEEFBVNo0aNXhQXh+Vr14aps5Kn05LM0JE2ERFh1aowzhrAlVfCBx/A\nddeFgk0KmxkccECYKxbg978P17u99lrcXIVIRZuISIGbNQv22QcOOihM7F6zpqaakk074ohw/Vvp\nLAtr18ZOVDiiFm1m1s/MZpjZLDMbupntjjMzN7PumcwnIpLvxoyBrl1Dz8F//AMaNIidSLLdnnuG\n3sMnnhgG6D3kEPjmm9ipCkO0os3MagLDgUOBzsBJZta5jO22BH4DvJ3ZhCIi+Wv1ajj7bDjpJNht\nN3j3XTjyyNipJFdsuSU89BDcey9MnhzGepP0i3mkrQcwy90/d/fVwBigfxnb3Qj8Afgxk+FERPJZ\n7drw3Xfh+rWJE3U6VCrODM48M/Q27tEjLHvuuXWzY0j1M3eP88ZmxwH93P2sxPNBwF7ufl7SNnsC\nw9z9WDObCFzq7pPL+FlnA2cDtGzZstuYMWPSnr+4uJii0gFtCkwhtx0Ku/1qe2633R3GjWvJHnss\noVWrHykpCT0EU5EP7a+sQm47pN7+Tz8t4le/6s7uuy9m2LCPad58VQbSpVcm9n3v3r2nuHtql3+5\ne5QbcDxwb9LzQcDfk57XACYC7RLPJwLdy/u53bp180yYMGFCRt4nGxVy290Lu/1qe+5atsx90CB3\ncL/ggoq/PtfbXxWF3Hb3irV/9Gj3Bg3cmzZ1f/bZ9GXKlEzse2Cyp1g7xTw9OhfYLun5tsC8pOdb\nArsCE81sNrA3MFadEUREKua998KE4A89BNdfD3/6U+xEkq8GDoQpU8LMGYcfDtdcEztRfok5uO4k\noIOZtQe+BgYAJ5eudPclQLPS55s7PSoiImUbPz7MFdqkCbzyCvTqFTuR5Ludd4a33oJLLoH27WOn\nyS/RijZ3X2Nm5wHjgJrAfe7+oZndQDhUODZWNhGRfNGjB5xxRjjC1rx57DRSKOrVg+HD1z1/+OGw\n7Jhj4mXKB1HHaXP35919J3ffwd1/n1h2TVkFm7v30lE2EZHyTZoUBkBdvjxMAH7nnSrYJB73MDTI\nsceGadF+1FgQlaYZEURE8oQ7/OUvsO++MH06fPVV7EQiYWiQ//u/cLr0jjvC7BszZ8ZOlZtUtImI\n5IGVK8Opp4suCnNETpsGHTvGTiUS1KkDt90GzzwTxnXr2lWFW2XE7IggIiLV5OKL4emn4c9/hgsu\nCEc3RLLN4YeH2TfGj4cOHWKnyT0q2kRE8sDVV8MBB8CAAbGTiGzedtvB4MHh8eTJUFysXs2p0ulR\nEZEcVVIC//xnmDaodWsVbJJb3OE3vwmdZt7W7OIpUdEmIpKD3ENPvLPOgiefjJ1GpOLM4PHHoWVL\n6NcvDAItm6eiTUQkx7jD0KFhKI9LL4XjjoudSKRyWrcOgz4XFUHfvvDJJ7ETZTcVbSIiOeb3v4c/\n/AHOOSfcq9OB5LK2bUPhZga33x47TXZTRwQRkRwydy7ccgsMGhRGnFfBJvlgp53gv/+FNm1iJ8lu\nOtImIpJDtt02XLR9331QQ/+DSx7ZYQeoXRsWLIATT4TvvoudKPvoIy8ikgMeeSRcwwawyy5QS+dJ\nJE/NmgVjx8LBB8PixbHTZBcVbSIiWe7pp8Pp0H/9C9aujZ1GJL322Sf0iP7wQzjssDCOmwQq2kRE\nsthLL8EJJ0C3bmEKoJo1YycSSb9+/WDMGHjnHejfX5PMl1LRJiKSpf7zHzjqqDCH6AsvwJZbxk4k\nkjnHHAP33w/z58OiRbHTZAcVbSIiWerdd8OUPy++CE2axE4jknmDBoXPQatW4dKAQr88QEWbiEiW\nWbMm3J93HkybFkaMFylUdeqEKdsGDgwzgJSUxE4Uj4o2EZEssngxdO0aTo0C1K8fN49INqhRA3be\nGUaODANKFyoVbSIiWeS88+Cjj8J4VSKyzrXXwvHHw9VXw9SpsdPEoaJNRCRLPPIIPPRQ+OO0116x\n04hkFzO4++5wucApp8CKFbETZZ6KNhGRLDBnDpx7LvTsCVdeGTuNSHZq0gRGjQq9ST/9NHaazNOY\n2iIiWeCBB0LPuAcf1GwHIptz0EHw+eewxRaxk2SejrSJiGSBq64KQxtsv33sJCLZb4stwpec228P\nc5UWChVtIiIRffBBOM1jFibMFpHUfPYZ/Pa3YRgQ99hpMkNFm4hIJCtXhimqDj9cg4aKVNROO8Et\nt4Tp3e65J3aazFDRJiISyeWXw8cfw/DhmlNUpDJ+8xvo0wcuuqgwOiaoaBMRieCFF+COO+DCC6Fv\n39hpRHJTjRqhN2m9ejB4cP6fJlUfJRGRDPvuOzj9dNhlF7j55thpRHJb69YwejQ0bhyuDc1nKtpE\nRDKsQQM48UQ444xwhEBEquaww9Y9XrYMttwyXpZ00ulREZEMcg/DFfz1r7DHHrHTiOSXW24Jn6ul\nS2MnSQ8VbSIiGfLZZ2Ey+OnTYycRyU/77w9ffhk6KOQjFW0iIhkydCjMmhWm4hGR6rfPPnDFFaFz\nwqRJsdNUPxVtIiIZ8M478PjjcMklsO22sdOI5K+hQ6F583Cfb71JVbSJiKSZ+7o/JJdcEjuNSH5r\n2BCGDYM334QvvoidpnqpaBMRSbMJE8Jt2LD87dUmkk1+9atwDWm+zeWrok1EJM323x8eeCD8IRGR\n9KtbF1q1Cke5v/kmdprqo6JNRCSN3KFWLRg0KPwhEZHMOecc2HdfWL06dpLqoaJNRCRNVq8Ovdke\neyx2EpHC1L8/fP45jBgRO0n1UNEmIpImI0bAW2/pOjaRWA49FA44AG68EYqLY6epOhVtIiJpUFwc\n/lD06gWHHBI7jUhhMguzJCxYALffHjtN1aloExFJg9tvD38obrkl/yexFslme+8NRx8NjzwCa9fG\nTlM1UYs2M+tnZjPMbJaZDS1j/cVm9pGZTTezV8ysbYycIiIVsXgx/PGPcOyxsNdesdOIyF13wZQp\nULNm7CRVUyvWG5tZTWA40BeYC0wys7Hu/lHSZtOA7u6+wszOBf4AnJj5tCIiqWvUCJ58Etrqa6ZI\nVmjZMtyvWQMrV+budaYxj7T1AGa5++fuvhoYA/RP3sDdJ7j7isTTtwBN/iIiOaFPH+jQIXYKESm1\nahXsvDPccEPsJJUXs2jbBvgq6fncxLJNORN4Ia2JRESq6L774De/yZ9xoUTyRd260LVr+IyuXBk7\nTeWYR5pN1cyOBw5x97MSzwcBPdz9/DK2HQicBxzg7qvKWH82cDZAy5Ytu40ZMyat2QGKi4spKipK\n+/tko0JuOxR2+9X2zbfdHU4//WfUr7+Wu+6amqFkmaF9X5hth/xq/7Rpjbj44i5cccXH9Ov3bbnb\nZ6LtvXv3nuLu3VPa2N2j3ICewLik51cCV5axXR/gY6BFKj+3W7dungkTJkzIyPtko0Juu3tht19t\nL28bd3AfOTLtcTJO+75w5VP7S0rcO3Z079Ejte0z0XZgsqdYO8U8PToJ6GBm7c2sDjAAGJu8gZnt\nCfwDONLdF0TIKCKSsjvvhCZN4IQTYicRkbKYwZAh8M47MDUHD4ZH6z3q7mvM7DxgHFATuM/dPzSz\nGwhV51jgj0AR8JiFgY7muPuRsTKLiGzKvHmhx+iFF0L9+rHTiMimnHoq7LAD7LFH7CQVF61oA3D3\n54HnN1h2TdLjPhkPJSJSCWvXhj8Gv/pV7CQisjlbbQWHHRY7ReVoRgQRkWqw3Xbwz3/CjjvGTiIi\n5VmzBq68Mnxmc4mKNhGRKpo0CSZPjp1CRFJVqxa89hrceiuUlMROkzoVbSIiVTR0aJiyKtfnNRQp\nJEOGwMyZ8PLLsZOkTkWbiEgVzJ0L48fDWWfl/ryGIoXkuOOgcWN48MHYSVKnok1EpAqeeircH398\n3BwiUjF168KRR8LYsbkzg0nU3qMiIrnu3/+Gzp2hY8fYSUSkok4+OXRKWLIEmjePnaZ8KtpERCpp\nxQp47z0477zYSUSkMg4+ONxyhYo2EZFK2mILmD8ffvwxdhIRqYqPPoKdd87+61J1TZuISBXUrRsG\n6xSR3PTss7DLLvCf/8ROUj4VbSIilbB4MXTpAuPGxU4iIlXRqxfUqxeuT812KtpERCrh2WfD9Ww6\nyiaS24qK4JBD4Iknsn+gXRVtIiKV8PTT0Lo19OgRO4mIVNUxx8DXX8PUqbGTbJ6KNhGRCiopgQkT\noG9fqKH/RUVyXt++4X7ChLg5yqPeoyIiFfTBB7BwYbgWRkRyX6tW8MIL2X/kXEWbiEgFmcGJJ8KB\nB8ZOIiLVpV+/2AnKpwP7IiIVtNtuMGYMtGkTO4mIVJeFC+EPfwhjtmUrHWkTEamAkhL46ito2zZ2\nEhGpTiUlcMUV4b5z59hpyqYjbSIiFfDee9CuHTz2WOwkIlKdmjeHXXfN7s4IKtpERCqgdNT0vfaK\nm0NEqt8BB8B//5u947WpaBMRqYApU6BFC9huu9hJRKS6de8OxcXw6aexk5RNRZuISAVMmQLduoUe\npCKSX7p1g9q14bPPYicpmzoiiIikaMWK0LPsqKNiJxGRdNhlF1i2DOrWjZ2kbCraRERSZAajR4ch\nP0Qk/9Sokb0FG+j0qIhIyurXh5NPVtEmks+eeCJMIJ+NnRFUtImIpOjNN2Hy5NgpRCSdFiyAF1+E\nuXNjJ9mYTo+KiKToqqtg5cpQvIlIfurUKdx/8gnUqRM3y4Z0pE1EJEWffAIdO8ZOISLpVPoZ/+ST\nuDnKoqJNRCQFxcU1mT9fRZtIvmvRAho1UtEmIpKz5s2rD0CHDpGDiEhamUGfPlBUFDvJxnRNm4hI\nCr75ph4A7dtHDiIiaVc6t/DEiVFjbERH2kREUtC16yJef12nR0UkHhVtIiIpKCpay89/HsZqE5H8\n9swz4QvawoXZ1X1URZuISAreeKMpzz4bO4WIZEJJCcyYAd9/r6JNRCTnjBnThttvj51CRDKhVatw\n/8MP2TWnlYo2EZEU/PBDHVq2jJ1CRDKhtGjLttOjKfUeNbN6wOHAfkBrYCXwAfCcu3+YvngiItlh\nyZLaNG8eO4WIZEKzZuF+6dLacYNsoNyizcyuA44AJgJvAwuAesBOwC2Jgu4Sd5+evpgiIvGsWQPL\nl9eiadPYSUQkE+rXh8MPh2bNVsWOsp5UjrRNcvfrNrHudjNrAbSpvkgiItllyZJw36RJ3BwikjnP\nPAMTJ34LdIod5X/KLdrc/bkNl5lZDaDI3Ze6+wLC0TcRkbzUqBE88sibHHJIz9hRRKSApdwRwcwe\nNrOGZtYA+AiYYWaXpS+aiEh2qFkTtt56FY0bx04iIpnSvz9cf33n2DHWU5Heo53dfSlwFPA84ZTo\noLSkEhHJIl99BQ891IYvv4ydREQyZcmS0Gs8m1SkaKttZrUJRdvT7v4T4FV5czPrZ2YzzGyWmQ0t\nY31dM3s0sf5tM2tXlfcTEamMWbPg3nu354svYicRkUypXx9Wr86ukdEqkuYfwGygAfCambUFllb2\njc2sJjAcOBToDJxkZhsehzwTWOTuOwJ/Bm6t7PtVt3YjR8aOICIZsnp1uK+TXV+6RSSN6tSBn37K\n0aLN3f/m7tu4+2Hu7sAcoHcV3rsHMMvdP3f31cAYoP8G2/QHRiUePw4cZGZWhfesNu1GjSp/IxHJ\nC2vWhPtaKY1sKSL5oFYtWLs2K0qO/ym3aDOzgYneouvxYI2Z7WBmP6/Ee28DfJX0fG5iWZnbuPsa\nYAmgkZJEJKOWLYPWrVewKruGbBKRNKpVC7bYYk3sGOtJ5XtjU2CamU0BpgDfEQbX3RE4APge2Oh6\ntBSUVb5ueI1cKttgZmcDZwO0bNmSiRMnViJO+dqNHLn+EbbEQb/Zgwcz+7TT0vKe2ai4uDhtv+Nc\nUMjtL9S2f/55U+bN242PPprC2rXLYseJolD3PRR226Fw2//tt51Ztqx+VrU9lXHa/mpmdwAHAvsC\nuxOmsfoYGOTucyr53nOB7ZKebwvM28Q2c82sFrAV8EMZGUcAIwC6d+/uvXr1qmSkcvTqBaXXspmB\nh/qxXeJWKCZOnEjafsc5oJDbX6htX7Ei3Hfp0o299oqbJZZC3fdQ2G2Hwm1/06bw5ZfLs6rtKV2h\n4e5rgZcSt+oyCehgZu2Br4EBwMkbbDMWGAy8CRwHjE9cTycikjGl17Ktya4zJSKSRj/9BLVqlcSO\nsZ6UL6tNFFfnEw4q/e917n5kZd44cT3cecA4oCZwn7t/aGY3AJPdfSzwT2C0mc0iHGEbUJn3SofZ\ngwcX1NE1kUJW2mu0tBepiOS/1auhVq3sOk5Ukb5QTxGKqGeAaik93f15wkC9ycuuSXr8I3B8dbxX\ndZt92mkq2kQKRP364X7lyrg5RCRzVq6EunVz9Egb8KO7/y1tSUREslSDBuF++fK4OUQkc5Yvh3r1\n1saOsZ6KFG1/NbNrgReB/3V8d/ep1Z5KRCSLFBWF+2WF2XFUpCAtWwbbbJO7RdtuhLlGD2Td6VFP\nPBcRyVuNGoX7JUvi5hCRzFmyBHbaKbt6H1WkaDsa2D4xe4GISMFo2BDMnB9+yK7R0UUkPdxh0SIo\nKsquoq0ik2q9BzRKVxARkWxVowZsueUaFi6MnUREMqG4OPQebdTop9hR1lORI20tgU/MbBLrX9NW\nqSE/RERyyVZb/cR339WOHUNEMuC778L9VlvlbtF2bdpSiIhkuaZNVzF//haxY4hIBsyfH+6bNMmu\nCYdTLtrc/dV0BhERyWZNmqxm9uzYKUQkE0qLtqZNs+sy/nKLNjNbRhmTtBMmc3d3b1jtqUREskzz\n5qt44w0oKQnXuIlI/vr663DfrFmOFW3uvmUmgoiIZLOWLVexahUsWABbbx07jYik0+zZYVDthg2z\n65o2fV8UEUnB1lv/CKBTpCIFYPZsaNcOLMtG+VHRJiKSgtKi7YsvIgcRkbT74gto2zZ2io2paBMR\nSUHr1isxg08/jZ1ERNKppCR8znfeOXaSjaloExFJQd26JbRrB598EjuJiKTTV1/BypXQsWPsJBtT\n0SYikqKOHVW0ieS70s+4ijYRkRzWqRPMmBFOn4hIfvr443Cvok1EJIftums4bTJzZuwkIpIu770H\nzZpB8+axk2xMRZuISIq6dg33U6bEzSEi6TN1KnTvnn3DfYCKNhGRlHXuDHXrqmgTyVcrV8KHH0K3\nbrGTlE1Fm4hIimrXhj32UNFrnOy0AAAbgklEQVQmkq+mT4e1a1W0iYjkhW7dwukTdUYQyT+lX8hU\ntImI5IHu3WHZMg39IZKP3nkndEDYbrvYScqmok1EpAL23z/cv/pq3BwiUr3cYcKE8BnPxk4IoKJN\nRKRCdtgBtt02/OcuIvlj9myYMwd6946dZNNUtImIVIBZ+E994sTwzVxE8kPpF7FevaLG2CwVbSIi\nFdS7N3z3XRgaQETyw4QJ0KJFGNonW6loExGpoNLTJxMnRo0hItWk9Hq2Xr2y93o2UNEmIlJh7dpB\n27a6rk0kX8yaBV9/nd3Xs4GKNhGRSunTB15+GVavjp1ERKrqhRfC/YEHxs1RHhVtIiKVcNRRsHQp\njB8fO4mIVNUTT4Rr2XbaKXaSzVPRJiJSCX37wpZbwr//HTuJiFTFggXw+utw7LGxk5RPRZuISCXU\nrQuHHw5PPQVr1sROIyKV9dRTYVo6FW0iInnsmGPg++/hP/+JnUREKuuJJ8Kg2bvvHjtJ+VS0iYhU\n0qGHQv36OkUqkqsWLYJXXglH2bJ5qI9SKtpERCqpQQPo1y98Uy8piZ1GRCrqmWfC5Q3HHBM7SWpU\ntImIVMExx8C8efDOO7GTiEhFPfFEmEv4Zz+LnSQ1KtpERKrg8MOhdm0YMyZ2EhGpiEWLYNy48MWr\nRo5UQzkSU0QkOzVqFP7THzUKVqyInUZEUjVqFPz4I5x+euwkqVPRJiJSRUOGwOLFOtomkitKSuDO\nO6FnT+jSJXaa1KloExGpov32g112geHDw8TTIpLdxo+HmTPDF65coqJNRKSKzMJ//lOnwqRJsdOI\nSHnuvBOaNYPjjoudpGKiFG1m1sTMXjKzmYn7xmVs08XM3jSzD81supmdGCOriEgqBg6EoqJwtE1E\nstfcufD003DmmVCvXuw0FRPrSNtQ4BV37wC8kni+oRXAqe6+C9AP+IuZNcpgRhGRlDVsCIMGwaOP\nhlkSRCQ7jRgRLmP41a9iJ6m4WEVbf2BU4vEo4KgNN3D3T919ZuLxPGAB0DxjCUVEKujcc2HVKrj/\n/thJRKQsq1fDPffAYYdB+/ax01RcrKKtpbvPB0jct9jcxmbWA6gDfJaBbCIilbLbbqFTwl13aYYE\nkWz01FPwzTe51wGhlHmaujqZ2cvA1mWsugoY5e6NkrZd5O4bXdeWWNcKmAgMdve3NrHN2cDZAC1b\ntuw2JgP97ouLiykqKkr7+2SjQm47FHb71fby2z5+fAtuvLEzt9wynb32+iEDyTJD+74w2w751f4L\nL+zCggV1GT36bWrWLH/7TLS9d+/eU9y9e0obu3vGb8AMoFXicStgxia2awhMBY5P9Wd369bNM2HC\nhAkZeZ9sVMhtdy/s9qvt5Vu1yr1FC/eDD05vnkzTvi9c+dL+adPcwf3WW1N/TSbaDkz2FGucWKdH\nxwKDE48HA09vuIGZ1QGeBB5w98cymE1EpNLq1IFLL4UXX4SJE2OnEZFSv/0tNG4Mv/xl7CSVF6to\nuwXoa2Yzgb6J55hZdzO7N7HNCcD+wGlm9m7ilkPjFotIoTrvvDAJ9RVXaLBdkWwwcSK88AJceWUo\n3HJVrRhv6u4LgYPKWD4ZOCvx+EHgwQxHExGpsvr14brr4Kyz4Mknw9ykIhKHe/gCte224QtVLtOM\nCCIiaTB4MHTqFE7JrFkTO41I4XrySXjnHbj++vCFKpepaBMRSYNateCmm2DGDBg5MnYakcK0Zk34\n4tSpE5x6auw0VaeiTUQkTfr3h549w6nSlStjpxEpPPffH7443XRT+CKV61S0iYikiRnccgt8/TX8\n/e+x04gUlhUrwhemnj3DF6h8oKJNRCSN9t8/TJlz882waFHsNCKF4+9/h3nzwhcns9hpqoeKNhGR\nNLv5ZliyJPzxEJH0W7QofN5+8YvwxSlfqGgTEUmz3XeHgQPhb3+DuXNjpxHJf7fcEr4o3Xxz7CTV\nS0WbiEgG3HBDmET++utjJxHJb3Pnhi9IAwfCbrvFTlO9VLSJiGRAu3Zw7rlw330wZUrsNCL567LL\nwhekG26InaT6qWgTEcmQa6+F1q3hlFNCzzYRqV4PPwxjxsBVV4UvSvlGRZuISIY0bgyjRoVxoy69\nNHYakfzy5ZcwZAjss08YUDcfqWgTEcmgAw+ESy6Bu+6C556LnUYkP6xdG6aOW7sWRo/Oj4F0y6Ki\nTUQkw37/+9Cj9IwzYMGC2GlEct9tt8Grr4ax2bbfPnaa9FHRJiKSYXXrwkMPhSEJzjwT3GMnEsld\nU6fC1VfDsceGo235TEWbiEgEu+4Kt94Kzz4LI0bETiOSm1asCB17mjeHf/wjf2Y+2BQVbSIikZx/\nPhx8MFx0UeicICIVc/nl8MknMHIkNG0aO036qWgTEYmkRg24/36oXz8cLfjpp9iJRHLH88/D8OHh\nS0/fvrHTZIaKNhGRiFq3hnvuCQPuXndd7DQiuWHBgtCRZ9dd4aabYqfJHBVtIiKRHXNM+AN0883w\n+uux04hkN3f45S/DpPAPPQT16sVOlDkq2kREssBf/gLt28OgQaFXqYiU7d57YezYMCn87rvHTpNZ\nKtpERLLAllvCgw+Gya7PPz92GpHs9OmncOGF0KcPXHBB7DSZp6JNRCRL9OwJw4aFEd0feSR2GpHs\nsmoVDBwYxjkcOTJ05Ck0BdhkEZHsNWxYmDvx9NNh/PjYaUSyw5o1cNJJMGlS6LizzTaxE8Whok1E\nJIvUqhWu1+nQAY48Et58M3YikbhKSsKXmCefDNd+Hnts7ETxqGgTEckyTZvCSy9Bq1Zw6KEwbVrs\nRCJxuMOvfx2u9/zd7wrzOrZkKtpERLLQ1lvDK6/AVluFWRM+/jh2IpHMcofLLoO774YrroDf/jZ2\novhUtImIZKk2beDll6FmzdBb7vPPYycSyZwbboA//Skcabv55vyfVzQVKtpERLJYhw6hcPvxRzjo\noDAkiEi+u/32MEPIaafB3/6mgq2UijYRkSy3664wbhwsXBiOuC1YEDuRSPqMGAGXXALHHRd6ihbi\n0B6bol+FiEgO6N4dnnsO5swJ17gtWhQ7kUj1e+ghOOccOOyw8LhWrdiJsouKNhGRHLHffvDUU6FT\nwqGHwrJlsROJVJ8nn4TBg6FXL3j8cahTJ3ai7KOiTUQkhxx8MPzrXzB5chjHbeXK2IlEqu7FF2HA\ngHBE+emnoX792Imyk4o2EZEc078/PPAAvPpqGGh09erYiUQq7/XX4aijoFMneOGFMA+vlE1Fm4hI\nDjr5ZPjHP8IfuVNOCdP8iOSayZPhF7+Atm3D0bbGjWMnym66xE9EJEf98pdQXAwXXwxbbAH336+e\ndpI7PvgADjkEmjULw9q0aBE7UfZT0SYiksMuuigUbtdcA0VFcMcdGtNKst/MmWH4mnr1wswfhToB\nfEWpaBMRyXHDhoWepH/8IyxdCnfdFQo4kWz06qtw0kmwdi1MmADt28dOlDt0IF1EJMeZwa23hml/\nHn449MB7773YqUTWt3Zt+Dd64IGhs8GECaHzgaRORZuISB4wg6uvDqeali6FvfYKE227x04mAvPn\nQ9++cO21oePMlClhpg+pGBVtIiJ5pFevcJStd28491w48URYsiR2Kilk48bBHnvA22/DyJFhuBqd\nvq8cFW0iInmmefMw5dWtt8ITT0DXrjBpUuxUUmh++gmuvBL69YOWLcPwHoMHx06V26IUbWbWxMxe\nMrOZiftNjsxiZg3N7GszuyOTGUVEclmNGnD55WHg0jVrYN994S9/0elSyYw5c8JR31tuCUPTvPOO\nrl+rDrGOtA0FXnH3DsArieebciPwakZSiYjkmZ49Ydq0MAH3RReFked/+CF2KslnY8dCly4wfXro\nGDNihKalqi6xirb+wKjE41HAUWVtZGbdgJbAixnKJSKSd5o0CZNx//WvYQaFLl3gjTdip5J8s3p1\n+GLQvz+0awdTp4ahPaT6mEc4Vm5mi929UdLzRe7eeINtagDjgUHAQUB3dz9vEz/vbOBsgJYtW3Yb\nM2ZM2rKXKi4upqhAr6Qs5LZDYbdfbc/9ts+YsSU33NCZb76px5lnfsGAAXNSmkUhX9pfGYXcdkit\n/V9/XY8bb+zMjBkNOfrouZxzzmfUqZP75+Izse979+49xd27p7Sxu6flBrwMfFDGrT+weINtF5Xx\n+vOAyxOPTwPuSOV9u3Xr5pkwYcKEjLxPNirktrsXdvvV9vyweLH7CSe4g/vBB7t/+235r8mn9ldU\nIbfdvfz2P/qoe8OG7o0auT/xRGYyZUom9j0w2VOsrdI2I4K799nUOjP71sxauft8M2sFLChjs57A\nfmY2BCgC6phZsbtv7vo3EREpx1ZbwZgxcNBBcMEFYTiGhx8Ow4SIpGrlyjDv7d13h3EBx4wJp0Ul\nfWJd0zYWKO34Oxh4esMN3P0Ud2/j7u2AS4EHVLCJiFQPMzj77DB2VqNGoYC77rowar1IeT75BPbe\nOxRsl10WeimrYEu/WEXbLUBfM5sJ9E08x8y6m9m9kTKJiBSc3XcPY7ideipcf32YxHvevNipJJuN\nHh2mSps3D55/Hv7wB6hdO3aqwhClaHP3he5+kLt3SNz/kFg+2d3PKmP7kb6JTggiIlI1RUVhpPqR\nI8N4WnvsAf/3f7FTSbZZvhxOPz0U+N26wbvvwqGHxk5VWDQjgoiIAGG0+ilToFWr8Mf48MM1NIjA\nihU1ue026NABRo2Ca64Jc9xus03sZIVHRZuIiPxPx47hOrcbbwz3P/857L8/vP12E82mUGC+/z4U\naAMG7M1ll4UZDV57LZxGr5W2boyyOSraRERkPfXrw7BhMHt2mPrqiy9g6NDd6doVHn1UnRXy3dy5\nYZDctm1D8b7HHot5661wdO3nP4+drrCpaBMRkTI1aBCGBPnsM7j88k9YuRIGDAhHXO69F1atip1Q\nqtOnn8KZZ8L228Pf/w7HHQcffAA33vghe+0VO52AijYRESlHnTpw6KHf8OGH8Pjj0LBhmAR8++3h\n9tuhuDh2QqmKqVPh+OPDqfGHHw5DwcyaFa5f22WX2OkkmYo2ERFJSc2acOyxYYiQF1+EnXeGSy4J\np9Guuw4WLoydUFLlDq++CoccEnqCvvgiDB0aTonfcYfGXMtWKtpERKRCzKBvXxg/Ht58E/bbL1yc\n3rZtGCH/669jJ5RNKSmBZ56BffeFXr3CsB033wxz5sBNN0HLlrETyuaoaBMRkUrbe2946il4/304\n+mj429+gfXs46yyYOTN2Oim1Zg089FAYg+/II2H+fBg+PBxZGzo0TG0m2U9Fm4iIVNmuu4aR8mfO\nDNe7PfhgOH16wgkwbVrsdIXrxx/hrrtgp51g4MBwWnT06NDpYMiQ0FNYcoeKNhERqTbt24cjOF9+\nCVdcAePGQdeu0K9fGONLY71lxtKlcOut4dq0IUOgRQt4+mmYPj0Ub5p2KjepaBMRkWrXsuX610pN\nnQoHHBDG+Xr2WRVv6fLdd2GMvTZtwmnP3Xdfd+3hkUdCDf3Vz2nafSIikjZbbQVXXhmOvN1xR+ik\ncMQRsNtucNll4Xq4776LnTJ3rV0bOhMMHw4nnRQ6g9x0E/TpA5Mnh16hvXuHziOS+zQRhYiIpF39\n+vDrX4cxwB55BEaMCJ0WbrstrO/QIfRoLL117KhCoyzLlsFbb8F//xvmhX3rrbAMoHVrOOWUMAxL\nx45xc0p6qGgTEZGMqV0bTj013H78MRwNeuONUIQ88wyMHBm2a9IE9tlnXRHXvXthXjQ/Z074/ZTe\npk8Pw3aYhaOVAweu+x21batCN9+paBMRkSjq1QvXuJXOZ+keejUmFynPPhvW1a4dBoHdd991xVy+\njSm2Zk0oypLbP3duWNegAey1F1x1VWj73ntrmI5CpKJNRESyglkYJmTnneGMM8Ky779fdyrwjTfC\ndXF/+lNYt8MO659S7dQpty60X7o0nN4sbdvbb6+bEmzbbddv2+67Qy39xS54+icgIiJZq1mz0Ovx\nyCPD81WrYMqUdadUX3gBHnggrGvcGHr2XFfo/OxnsMUW8bInc9/4VOf774dTnTVqhKJs8OB12du0\niZ1YspGKNhERyRl164bTo/vsE567h8nNk4uh558P62rVCmPE7bsv9OhR+QLu/febsnRpxV/nHmYc\nKD1SWDq9V1FROL159dUh2157QcOGlcsmhUVFm4iI5Cyz0PO0Qwc47bSwbOHCMC5ZaRF3113w5z9X\n5V12q1LGNm1g//3XXY+322461SmVo382IiKSV5o2hcMPDzeA1athxgz46afK/bzJkyfTvXv3Sr22\nRYtwfZpIdVDRJiIiea1OnXB0q7KWLi2ma9fqyyNSWTnUz0ZERESkcKloExEREckBKtpEREREcoCK\nNhEREZEcoKJNREREJAeoaBMRERHJASraRERERHKAijYRERGRHKCiTURERCQHqGgTERERyQHm7rEz\nVCsz+w74MgNv1Qz4PgPvk40Kue1Q2O1X2wtXIbe/kNsOhd3+TLS9rbs3T2XDvCvaMsXMJrt75WYQ\nznGF3HYo7Par7YXZdijs9hdy26Gw259tbdfpUREREZEcoKJNREREJAeoaKu8EbEDRFTIbYfCbr/a\nXrgKuf2F3HYo7PZnVdt1TZuIiIhIDtCRNhEREZEcoKJtM8zseDP70MxKzGyTvUfMrJ+ZzTCzWWY2\nNGl5ezN728xmmtmjZlYnM8mrzsyamNlLiewvmVnjMrbpbWbvJt1+NLOjEutGmtkXSeu6ZL4VlZNK\n2xPbrU1q39ik5Tm73yHlfd/FzN5MfD6mm9mJSetybt9v6jOctL5uYl/OSuzbdknrrkwsn2Fmh2Qy\nd3VIoe0Xm9lHif38ipm1TVpX5mcgl6TQ/tPM7Lukdp6VtG5w4nMy08wGZzZ51aXQ9j8ntftTM1uc\ntC6n972Z3WdmC8zsg02sNzP7W+J3M93Muiati7ff3V23TdyATsDOwESg+ya2qQl8BmwP1AHeAzon\n1v0LGJB4fDdwbuw2VaDtfwCGJh4PBW4tZ/smwA/AFonnI4HjYrcjnW0HijexPGf3e6rtB3YCOiQe\ntwbmA41ycd9v7jOctM0Q4O7E4wHAo4nHnRPb1wXaJ35Ozdhtqua29076XJ9b2vbE8zI/A7lyS7H9\npwF3lPHaJsDnifvGiceNY7epOtu+wfbnA/fl0b7fH+gKfLCJ9YcBLwAG7A28nQ37XUfaNsPdP3b3\nGeVs1gOY5e6fu/tqYAzQ38wMOBB4PLHdKOCo9KWtdv0JmSG17McBL7j7irSmyoyKtv1/8mC/Qwrt\nd/dP3X1m4vE8YAGQ0uCQWajMz/AG2yT/Th4HDkrs6/7AGHdf5e5fALMSPy9XlNt2d5+Q9Ll+C9g2\nwxnTKZV9vymHAC+5+w/uvgh4CeiXppzpUNG2nwQ8kpFkGeDurxEONGxKf+ABD94CGplZKyLvdxVt\nVbcN8FXS87mJZU2Bxe6+ZoPluaKlu88HSNy3KGf7AWz8gf594rDyn82sbjpCpkmqba9nZpPN7K3S\n08Lk/n6HCu57M+tB+Kb+WdLiXNr3m/oMl7lNYt8uIezrVF6bzSqa/0zC0YdSZX0Gckmq7T828e/5\ncTPbroKvzVYp50+cEm8PjE9anOv7vjyb+v1E3e+1MvVG2crMXga2LmPVVe7+dCo/ooxlvpnlWWNz\nba/gz2kF7AaMS1p8JfAN4Y/5COAK4IbKJa1+1dT2Nu4+z8y2B8ab2fvA0jK2y6r9DtW+70cDg929\nJLE4q/d9GVL5rObs57wcKec3s4FAd+CApMUbfQbc/bOyXp+lUmn/M8Aj7r7KzM4hHHE9MMXXZrOK\n5B8APO7ua5OW5fq+L09WfuYLvmhz9z5V/BFzge2Snm8LzCPMVdbIzGolvpmXLs8am2u7mX1rZq3c\nfX7iD/OCzfyoE4An3f2npJ89P/FwlZndD1xaLaGrSXW0PXFaEHf/3MwmAnsC/ybL9ztUT/vNrCHw\nHDAscfqg9Gdn9b4vw6Y+w2VtM9fMagFbEU6tpPLabJZSfjPrQyjoD3D3VaXLN/EZyKU/3OW2390X\nJj29B7g16bW9NnjtxGpPmD4V+bc7APh18oI82Pfl2dTvJ+p+1+nRqpsEdLDQY7AO4R/3WA9XLE4g\nXOsFMBhI5chdthhLyAzlZ9/oWofEH/vSa7yOAsrsoZOlym27mTUuPe1nZs2AfYGP8mC/Q2rtrwM8\nSbjm47EN1uXavi/zM7zBNsm/k+OA8Yl9PRYYYKF3aXugA/BOhnJXh3LbbmZ7Av8AjnT3BUnLy/wM\nZCx59Uil/a2Snh4JfJx4PA44OPF7aAwczPpnG7JdKv/uMbOdCRfcv5m0LB/2fXnGAqcmepHuDSxJ\nfCGNu98z1eMhF2/A0YSqehXwLTAusbw18HzSdocBnxK+ZVyVtHx7wn/gs4DHgLqx21SBtjcFXgFm\nJu6bJJZ3B+5N2q4d8DVQY4PXjwfeJ/zBfhAoit2m6mw7sE+ife8l7s/Mh/1egfYPBH4C3k26dcnV\nfV/WZ5hwSvfIxON6iX05K7Fvt0967VWJ180ADo3dljS0/eXE/3+l+3lsYvkmPwO5dEuh/TcDHyba\nOQHomPTaMxL/JmYBp8duS3W3PfH8OuCWDV6X8/uecKBhfuL/sbmE6zXPAc5JrDdgeOJ38z5JI0jE\n3O+aEUFEREQkB+j0qIiIiEgOUNEmIiIikgNUtImIiIjkABVtIiIiIjlARZuIiIhIDlDRJiIiIpID\nVLSJSEEys7Vm9q6ZfWBmz5hZoxReU9/MXjWzmpvZ5uXEoJsiItVKRZuIFKqV7t7F3XclTEn16/Je\nQBhU8wlffw7GDY0GhlRHQBGRZCraRETCFD3blD4xs8vMbJKZTTez65O2O4XEtF5m1srMXks6Wrdf\nYpuxhKndRESqlYo2ESloiVOdB5GYd9HMDibMIdoD6AJ0M7P9E/Mzbu/usxMvPZkwtV0XYA/CFE+4\n+yKgrpk1zWhDRCTv1YodQEQkkvpm9i5h/twpwEuJ5QcnbtMSz4sIRdwsYHHS6ycB95lZbeApd383\nad0CwhzFC9OWXkQKjo60iUihWpk4StYWqMO6a9oMuDlxvVsXd9/R3f8JrCRMHA+Au78G7A98DYw2\ns1OTfna9xPYiItVGRZuIFDR3XwL8Brg0cdRsHHCGmRUBmNk2ZtYicdqzppnVSyxvCyxw93uAfwJd\nE8sN2BqYnfHGiEhe0+lRESl47j7NzN4DBrj7aDPrBLwZ6i+KgYGEU54vAj8HXgZ6AZeZ2U+JbUqP\ntHUD3nL3NZlthYjkO3P32BlERHKCme0JXOzugzazzV+Bse7+SuaSiUgh0OlREZEUufs0YMLmBtcF\nPlDBJiLpoCNtIiIiIjlAR9pEREREcoCKNhEREZEcoKJNREREJAeoaBMRERHJASraRERERHLA/wMd\nGjUlMcAZnwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Set up the plot parameters\n",
"plt.figure(figsize=(10,5))\n",
"plt.grid(True)\n",
"plt.title('Nyquist Diagram')\n",
"plt.xlabel('Re(s)')\n",
"plt.ylabel('Im(s)')\n",
"\n",
"# Nyquist plot\n",
"sys = control.tf(1,[1,1])\n",
"omega = np.logspace(-3,3,60)\n",
"real_arr, imag_arr, freq_arr = control.nyquist_plot(sys,omega)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And even create animations for it!"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def animate_Nyquist(sys,omega,xlims,ylims):\n",
" # Set up figure params\n",
" fig = plt.figure(figsize=(10,5))\n",
" ax = fig.add_subplot(1,1,1)\n",
" ax.grid(True)\n",
" ax.set_title('Nyquist Diagram')\n",
" ax.set_xlabel('Re(s)')\n",
" ax.set_ylabel('Im(s)')\n",
"\n",
" # Get the magnitude and phase of the system\n",
"\n",
" mag_tmp, phase_tmp, omega = sys.freqresp(omega)\n",
" mag = np.squeeze(mag_tmp)\n",
" phase = np.squeeze(phase_tmp)\n",
"\n",
" # Compute the primary curve\n",
" # https://www.youtube.com/watch?v=HbuGlYEoYEA\n",
" x = np.hstack([mag*sp.cos(phase),np.flip(mag*sp.cos(phase),axis=0)])\n",
" y = np.hstack([mag*sp.sin(phase),-mag*sp.sin(phase)]) # negative to get the mirror image\n",
"\n",
" xData, yData = [], []\n",
" primary_line, = ax.plot([], [], '-b', animated=True)\n",
" primary_arrow = ax.annotate('', xy = (0,0), \n",
" xytext = (0,0),\n",
" size=15,\n",
" arrowprops = {'arrowstyle': \"->\"})\n",
"\n",
" # Mark the -1 point\n",
" ax.plot([-1], [0], 'r+')\n",
"\n",
" def init():\n",
" ax.set_xlim(xlims)\n",
" ax.set_ylim(ylims)\n",
" return (primary_line,primary_arrow)\n",
"\n",
" def animate(i):\n",
" xData.append(x[i])\n",
" yData.append(y[i])\n",
"\n",
" # plot the actual line\n",
" primary_line.set_data(xData, yData)\n",
" primary_arrow.set_position((x[i],y[i]))\n",
"\n",
" # This is a crude way to handle the last data point\n",
" if i < len(x)-1:\n",
" primary_arrow.xy = (x[i+1],y[i+1])\n",
" else:\n",
" primary_arrow.xy = (x[i],y[i]) \n",
"\n",
" # Print out the current frame using \"\\r\", the 'carriage return' character, as our end character.\n",
" # This makes Python print the frame on the same line.\n",
" print(\"Frame: {:0d}\".format(i), end=\"\\r\")\n",
" return (primary_line,primary_arrow)\n",
" \n",
" return fig, animate, init, x"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Frame: 119\r"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sys = control.tf(1,[1,1])\n",
"omega = np.logspace(-3,3,60)\n",
"xlims = [-1.1, 1.1]\n",
"ylims = [-0.55, 0.55]\n",
"fig, animate, init, x = animate_Nyquist(sys,omega,xlims,ylims)\n",
"# blit=True will only draw the parts that have changed.\n",
"# The `interval` parameter is the delay between frames in milliseconds and it controls the speed of the animation.\n",
"anim = animation.FuncAnimation(fig, animate, init_func=init, frames=range(len(x)), interval=40, blit=True)\n",
"HTML(anim.to_jshtml())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can take any rational transfer function and apply the previous procedure to sketch its contour. Knowing this, let's move on to the next important concept, which explains the reason behind the $-1$ point as well as the $Z = N+P$ equation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cauchy's argument principle\n",
"\n",
"The contour of any rational expression will encircle the _origin_ once in a clockwise fashion, **for each additional RHP-zero it has compared to RHP-poles**. Mathematically, this is expressed as:\n",
"\n",
"$$ N = Z - P $$\n",
"\n",
"where $N$ is the number of clockwise encirclements the contour makes about the origin, $Z$ the number of RHP-zeroes in $F$, and $P$ the number of RHP-poles in $F$. This may look confusing, so let's go through several examples as exercise:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Case 1\n",
"$F(s)=\\frac{s-1}{s+1}$: $N=Z-P=1-0=1$, So we expect its contour to encircle the origin once clockwise."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Frame: 119\r"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"IOPub data rate exceeded.\n",
"The notebook server will temporarily stop sending output\n",
"to the client in order to avoid crashing it.\n",
"To change this limit, set the config variable\n",
"`--NotebookApp.iopub_data_rate_limit`.\n"
]
}
],
"source": [
"sys = control.tf([1,-1],[1,1])\n",
"omega = np.logspace(-3,3,60)\n",
"xlims = [-1.1, 1.1]\n",
"ylims = [-1, 1]\n",
"fig, animate, init, x = animate_Nyquist(sys,omega,xlims,ylims)\n",
"anim = animation.FuncAnimation(fig, animate, init_func=init, frames=range(len(x)), interval=40, blit=True)\n",
"HTML(anim.to_jshtml())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Case 2\n",
"\n",
"$F(s)=\\frac{(s-2)(s-3)}{(s-1)(s-4)}$: $N=Z-P=2-2=0$, so we expect no clockwise encirclements of the origin."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Frame: 119\r"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s = control.tf([1,0],[1])\n",
"sys = (s+1)/((s-2)*(s-3))\n",
"omega = np.logspace(-2,2,60)\n",
"xlims = [-2.5, 2]\n",
"ylims = [-2.5, 2.5]\n",
"fig, animate, init, x = animate_Nyquist(sys,omega,xlims,ylims)\n",
"anim = animation.FuncAnimation(fig, animate, init_func=init, frames=range(len(x)), interval=40, blit=True)\n",
"HTML(anim.to_jshtml())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So Cauchy's argument principle turns out to be pretty straightforward: it tells us **how many more RHP-zeroes than RHP-poles a TF has**. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Proof for Nyquist Stability\n",
"\n",
"But so far we've discussed encirclements about the **origin**; how does this translate to the well-known $-1$ point in Nyquist stability? Usually profound realizations (or proofs) are made by connecting a series of trivial statements, and this is no exception.\n",
"\n",
"Consider the following equivalent statements:\n",
"\n",
"1. The closed-loop is unstable if the TF $\\big( \\frac{L}{1+L} \\big)$ has at least one RHP-pole.\n",
"2. If the TF $\\big( \\frac{L}{1+L} \\big)$ has at least one RHP-pole, then the TF $(1+L)$ has at least one RHP-zero.\n",
"3. If the TF $(1+L)$ has at least one RHP-zero, then its contour will encircle the origin at least once clockwise, according to \\textbf{Cauchy's argument principle.}\n",
"4. The TF $L$ is equal to $(1+L)-1$. In other words, the contour of $L$ is the contour of $(1+L)$ shifted to the left by one unit, on the complex plane.\n",
"5. If the contour of $1+L$ encircles the origin at least once clockwise, \\textbf{then the contour of $L$ will encircle the point $-1$ at least once clockwise.}\n",
"\n",
"The final statement above completes the proof for Nyquist stability."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Margins\n",
"\n",
"As a last exercise, let's visually deduce the Gain and Phase Margins on a Nyquist plot, using our understanding from Bode and simple geometry:\n",
"\n",
"\n",
"\n",
"In the previous figure, the black solid curve represents the contour of $L(j\\omega)$. At the point **1**, $\\angle L(j\\omega)=-180^o$, therefore the corresponding frequency is $\\omega = \\omega_{180}$. The length of the vector spanning from the origin to the point **1** is $\\lvert L(j\\omega_{180})\\rvert$, which is the **inverse of the Gain Margin**.\n",
"\n",
"On the other hand, the point **2** is the intersection between the contour of $L(j\\omega)$ and the unit circle, so $\\lvert L(j\\omega)\\rvert = 1$ here, therefore the corresponding frequency is $\\omega=\\omega_{0dB}$. The **Phase Margin** is therefore the difference between the phase at this point and $-180^o$, i.e. $\\angle L(j\\omega_{0dB})-(-180^o)$.\n",
"\n",
"From the figure above, we can see that both the Gain and Phase Margins can be geometrically interpreted as ``distances\" from the $-1$ point, which we try to stay as far away from as possible."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"- Notes written by Yiting Tsai (Sept 2018)"
]
}
],
"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.6.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}