Introduction to Simulation-based Optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lifting line theory: Background and numerical example\n",
    "====================================\n",
    "\n",
    "In this notebook, we'll start to look at simulation-based design optimization. For the purposes of this section, we'll look at some simple aerodynamic optimization problems using lifting-line theory.\n",
    "\n",
    "A good review of lifting line theory can be found here: \n",
    "\n",
    "https://en.wikipedia.org/wiki/Lifting-line_theory\n",
    "\n",
    "In lifting line theory, the bound circulation can be written as a function of the span-wise coordinate $y$ as follows:\n",
    "\n",
    "$\\Gamma(\\theta) = 2 b V_{\\infty} \\sum_{n=1}^{\\infty} A_{n} sin(n\\theta)$\n",
    "\n",
    "here $y = \\frac{b}{2} \\cos(\\theta)$. \n",
    "\n",
    "The Fourier coefficients $A_{n}$ are treated as the unknowns. In computations, the coefficients are truncated for $n = 1, \\ldots, N$ so that we solve for $A_{1}, A_{2}, \\ldots, A_{N}$.\n",
    "\n",
    "The following system of $N$ equations can be obtained to find these coefficients:\n",
    "\n",
    "$\\begin{bmatrix} a_{11} & a_{12} & \\ldots & a_{1N} \\\\ a_{21} & a_{22} & \\ldots & a_{2N} \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ a_{N1} & a_{N2} & \\ldots & a_{NN} \\end{bmatrix} \\begin{bmatrix} A_{1} \\\\ A_{2} \\\\ \\vdots \\\\ A_{N} \\end{bmatrix} = \\begin{bmatrix} b_{1} \\\\ b_{2} \\\\ \\vdots \\\\ b_{N} \\end{bmatrix}$\n",
    "\n",
    "where the entries $a_{nm}$ and $b_{n}$ are given by\n",
    "\n",
    "$a_{mn} = \\sin(n\\theta_{m})\\sin(\\theta_{m}) + c_{l,\\alpha} \\frac{c(\\theta_{m})}{4b} \\sin(\\theta_{m})$\n",
    "\n",
    "$b_{m} = \\frac{c(\\theta_{m})}{4b} c_{l,\\alpha} (\\alpha_{\\infty} + x_{m})$\n",
    "\n",
    "Note that here $x_{m}$ is our design variable: The twist of the wing at the $m-th$ collocation point.\n",
    "\n",
    "Output functions of interest:\n",
    "---------------------------------\n",
    "\n",
    "We're not interested in determining just the coefficients $A_{n}$, we want to know what the lift, drag, and perhaps root bending moment are on the wing as a function of the design variables.\n",
    "\n",
    "The coefficient of lift produced by the wing is:\n",
    "\n",
    "$C_{L} = \\pi AR A_{1}$\n",
    "\n",
    "The coefficient of drag produced by the wing is:\n",
    "\n",
    "$C_{D} = \\pi AR \\sum_{n=1}^{N} n A_{n}^{2}$\n",
    "\n",
    "The root bending moment coefficient defined is defined as\n",
    "\n",
    "$C_{b} = \\frac{M_{b}}{\\frac{1}{2}\\rho V_{\\infty}^2 S b} = AR \\sum_{n=1}^{N} c_{n} A_{n}$,\n",
    "\n",
    "where $c_{n}$ are given as follows $c_{2} = 0$ and all other values are given by:\n",
    "\n",
    "$c_{n} = \\frac{(2\\sin(n\\pi/2) - \\sin(\\pi n))}{4 - n^2}$\n",
    "\n",
    "Optimization problem formulations:\n",
    "------------------------------------------\n",
    "\n",
    "Now we want to solve an optimization problem. We want to design the wing so that it has minimum $C_{D}$ while meeting a requirement for lift $C_{L} = C_{L,\\text{target}}$. This constrained optimization problem can be expressed as:\n",
    "\n",
    "$\\min_{x} \\qquad C_{D} \\qquad \\text{such that} \\qquad C_{L} = C_{L,\\text{target}}$\n",
    "\n",
    "This will give the classical lifting line result that the minimum induced drag is produced by an ellipitical lift distribution.\n",
    "\n",
    "However, this results in a lift distribution which can overload the wing tip or a heavier structure. Therefore, we can consider an additional constraint where we seek to reduce the bending moment at the root of the wing by limiting $C_{b} \\le C_{b,\\text{target}}$.\n",
    "\n",
    "$\\min_{x} \\qquad C_{D} \\qquad \\text{such that} \\qquad  C_{L} = C_{L,\\text{target}} \\qquad \\text{and} \\qquad C_{b} \\le C_{b,\\text{target}}$\n",
    "\n",
    "This will produce a more triangular lift distribution that shifts lift inboard.\n",
    "\n",
    "Computational cost\n",
    "-----------------------\n",
    "\n",
    "These problems are solved with scipy.minimize. Note that as $N$ increases, the computational cost grows rapidly. That's because the lifting line theory code is invoked $N+1$ times for a gradient calculation for the objective and each constraint. This computational cost adds up quickly as the number of design variables increase. In addition, a larger design problem takes longer to solve, making these larger design problems even more computationally expensive.\n",
    "\n",
    "To address this issue, we need a fast way to compute the gradient!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4lOW9//H3dyaZ7AnZWRIS9n01rIJ1q6K2qK0LWK07rUuP9tieWrW1tb/WtnZv8bTU2larWHFFxa2KPYpsYSdAIBAgISEr2bdJ5v79MZM0QiADTPLM8n1dl9eVTB4yH4fwyT33cz/PLcYYlFJKBReb1QGUUkr5npa7UkoFIS13pZQKQlruSikVhLTclVIqCGm5K6VUENJyV0qpIKTlrpRSQUjLXSmlglCYVU+ckpJisrOzrXp6pZQKSJs2bao0xqT2dpxl5Z6dnU1ubq5VT6+UUgFJRA55c5xOyyilVBDScldKqSCk5a6UUkFIy10ppYKQlrtSSgUhLXellApCWu5KKRWELFvnrlR/qG9xUlrbQmVDK9WNbVQ3tlHZ0EZds5OmtnaanS6a29ppdnbgbD9xy8kwuxDtsBPlCCM63E6Uw05CVDjJsQ6SYtz/pcRGMHhAFLER+s9J+Q/9aVQBr7bZSUF5PXvLGigob6CouokjNc0UH2umttnZbzkGRIczZEAUGYlRZCZGMzItllHpcYxKjyU+MrzfcigFWu4qwJTXtbCtuJZtRTVsP1JL/tE6yupaT3p8ZLiNwQOiSImNINkz0k6OcZAQ7SDG4R6JR3lG5OF2G9LtzxrA2eGiua2DZmcHTW3u/2qa2qhqbKO6ofOdQCtHapqpaXJS0+Qkr6TuhBwD4yMZMzCOyRkJTM4YwJSMBNLiI33/AinloeWu/JYxhn3lDaw7UMW6A1VsPlTD0bqWE46LDLe5R8lpcYxMiyUrOZqMxGgyEqNIjnEgIj18d99nrWxoo/hYE8XHmjlc3cS+snr2lbvfTRyta+FoXQv/3lvR9WcGJUQyPSuR2cOTmTM8iRGpsf2SVYUGMebEecb+kJOTY/TeMup4R2tbWJ1fzsf7Klh/oJqqxrbPfD0uIoxJ3Ua/4wfHk5EYjd3mv6XY4TIUH2tiV0kdW4tr2F5Uy44jtTS0tn/muJRYB7OGJ3PeqBQuGJOmI3vVIxHZZIzJ6fU4LXdlJZfLsK24hg/3lPPB7nJ2lX52SiMtLoI5I5KZPTyZGdmJDE+JxebHRe4tl8twoLKBDYXHWHegirUHqqio/+z00sQh8Vw4Jo0Lx6UzJSNBR/UK0HJXfswYw9aiGt7cXsqqHaWU1v5nqiUq3M65I1O4YGwqc0ekkJ0cHRKlZoyhsLKRNfur+GhPOWv2V9LidHV9fciAKK6YPIgrJg1ishZ9SPNpuYvIAuC3gB14yhjz0+O+fgvwBHDE89AfjDFPnep7armHnoLyBl7aVMwb20o4UtPc9fighEguGZ/OBWPTmD08mchwu4Up/UOLs4O1+6v4cE857+8q+8y5hozEKL44ZTDXnJPBiNRYC1MqK/is3EXEDuwFPg8UAxuBxcaYXd2OuQXIMcbc621ALffQUN/i5M3tpazILWLz4ZquxwfGR3L5pEFcMXkQ0zIHBMVUS19xuQybDx/reqdT3m36JicrketyMrl88iBdZx8ifFnuc4AfGGMu9Xz+XQBjzOPdjrkFLXfVTV5JLX//9CArt5V0TS/EOOx8ccpgvjQ9g5ysRC30M+ByGTYcrOaVzcW8ub2UprYOAKIddq6cOpib52YzdmC8xSlVX/K23L35VT8EKOr2eTEwq4fjviwi5+Ee5X/TGFPUwzEqiDk7XLyXV8bfPi1k48FjXY/PGpbEdTmZXDZpINEOHV2eDZtNmD3cfYL50S9O4K0dpbyUW8yGg9Us31DE8g1FzBqWxK3nZnPxuHTC7HqHkVDlq39pbwDLjTGtIvI14O/AhccfJCJLgCUAQ4cO9dFTK6s1trazfMNh/vJJYdfJ0diIMK45J4OvzsliuM4L94mYiDCuy8nkupxMCsobeGbtQV7eVMz6wmrWF1YzZEAUd8wfxqIZQ4ly6HmMUOOTaZnjjrcD1caYhFN9X52WCXw1TW387dOD/O3Tg9Q0uS/zH54awy1zs/nS9AydA7ZAfYuTlzYV88zaQxRWNgKQFOPgtnOzuWlONglRehuEQOfLOfcw3FMtF+FeDbMRuMEYk9ftmEHGmFLPx1cD3zHGzD7V99VyD1zHGtv44//t59m1h7rmfKcPHcDd54/kwrFpOpfuB1wuw/u7y3jyo/1sK3KfyI6NCOPmuVksmT+ChGgt+UDl66WQlwO/wb0U8mljzI9F5DEg1xizUkQeBxYC7UA1cJcxZs+pvqeWe+Cpb3Hyl08Keerjwq6rK+ePSuGeC0Yya1iSrr32Q8YY1u6v4smP9vNJQSUAcZFhLJk/nFvnDdN3VwFIL2JSPtPi7OCZtQf534/2c8wz/XLe6FQe+PxopmQOsDac8tqWw8f45Xt7u0o+KcbB3eeP4KY5WUSE6Zx8oNByV2fNGMOqHUd5/O3dFB9zX3Q0IzuRb10yhlnDky1Op87Up/sr+cW7+V3XHWQlR/Pdy8Zx6YR0ffcVALTc1VnZXlzDj97c1bWkcezAOB68bCyfG52qBRAEjDF8uKecx9/eQ0F5AwCzhyfxvS+MZ8LgU66FUBbTcldnpLqxjZ++vZsXc4sBSI5x8MAlY7h+RqZf33lRnZn2DhfPbzjMr9/fy7EmJyKweOZQvnPpWD3p6qe03NVpMcawYlMxj6/azbEmJw67jVvnZXPPBSN1F6EQUNvk5Pcf7uNvnx6k3WVIiXXw8BXjuGrqEH2n5me03JXX9pXV8/BrO9lQWA3A3BHJ/OiqiXpTqhCkPwv+T8td9aq9w8Uf/72f336wD2eHe7T2yBXjuXLqYB2thTBjDC9tKuYnne/iwmw88PnR3DF/uE7N+QEtd3VKe8vq+daKbWwvrgXc86wPLtB5VvUfxxrb+PGq3by0yX3+ZdrQAfzi2ik6ireYlrvqUXuHi2UfH+A37++jrcPFkAFR/OzLk5k3KsXqaMpPrd5TzoOvbKesrpWIMBvfumQMt80bpqN4i2i5qxMcqWnm/he2dC1vXDxzKA9dPpY4PWGqelHb5OSxN3fx8mb3KH7O8GR+ff1UBiboPq/9TctdfcaqHaU8+PJ26lraSYuL4Ilrp/C50alWx1IB5l+7yvjOy9upamwjMTqcn18zhc+PT7c6Vkjxttz1Zs9Brrmtg+++sp27n9tMXUs7F49L4537z9NiV2fk4vHpvH3/fOaPSuFYk5M7n8nle6/tpMXZYXU0dRwt9yBWWNnIVUvXsHxDEY4wGz9cOIE/fzWHpBiH1dFUAEuLi+Tvt87kkSvGEW4Xnl13iC//76cUVTdZHU11o+UepN7NO8rC339Cflk9w1NjeP2ec7l5brYucVQ+YbMJd8wfzqt3n0tWcjR5JXV84fefsHpPudXRlIeWe5Bp73Dx07f38LVnN1Hf2s5lEwfy+j3nMm6Q7qupfG/ikARW3juPi8elUdvs5Na/beRX7++lw2XNuTz1H1ruQaS2ycktf93IH/+9H7tNePjycTz5lem6Gkb1qYSocJbdlMO3Lx2DTeB3H+zjjr9vpL7FaXW0kKblHiQOVDRw9ZNr+KSgkpRYB8/fMYs7zxuu0zCqX9hswj0XjOSZ22YxIDqc1fkVOg9vMS33ILCmoJKrlq7hQGUj4wbF8/q98/R+68oS80al8Po95zIyLZa9ZQ1cuXRN131qVP/Scg9wz68/zFef3kBdSzufH5/OS1+fw5ABUVbHUiEsKzmGV+6ey3mjU6lubOMrT63ruoWB6j9a7gHKGMMv38vnoVd30OEy3HX+CP504znE6J6Yyg/ER4bz9M053HpuNs4Ow7dWbOMPH+7DqosmQ5E2QQBydrh4+NUdvJhbjN0m/PiqiSyaOdTqWEp9RpjdxqNfnEB2cgw/eCOPX7y3l7K6Vn6wcILel6YfaLkHmKa2du55bjOr8yuIDLex9IbpXDROL/9W/uvmudmkxkVw/z+38uy6Q1TUt/KbRVOJDNdNufuSTssEkNpmJ195aj2r8ytIjA5n+Z2ztdhVQLh80iCevW0mcZFhvJN3lJuf3kBDa7vVsYKalnuAqGpoZfGydWw5XMOQAVG8fNdcpg1NtDqWUl6bNTyZl74+l/T4CNYXVnPjU+upbdK18H1Fyz0AlNe1sGjZOnaV1jEsJYYVX5/DcN0wQQWgMQPjePFr7hVdW4tqWPzndVQ1tFodKyhpufu5IzXNXPentewrb2B0eiz//NpsButSRxXAspLdA5RhKTHsKq3j+mXrKKtrsTpW0NFy92MlNc0sWraWg1VNTBwSzwtL5pAWp5sjqMA3eEAU//zabMakx1FQ3sDiZeso14L3KS13P3W0toXFf15HUXUzUzISeO6O2XqrXhVU0uIiWb5kNuMGxXOgspEbnlpPpU7R+IyWux8qr2vhhj+v45BnxP7MbbNIiNKbf6ngkxTj4B+3z2R0eiwF5Q3c+NR6qhvbrI4VFLTc/UxlQys3PLW+6z4xz942i4RoLXYVvJJjI3jujtmMSI1hz9F6bnxqPTVNWvBny6tyF5EFIpIvIgUi8uApjvuyiBgR6XV/P3Wi+hYnNz+9gQLPydN/3D6TRJ2KUSEgNS6C5XfO7jrJetvfNtLUpuvgz0av5S4idmApcBkwHlgsIuN7OC4OuA9Y7+uQoaDF2cGSZzaRV1JHVnI0/7hjFsmxEVbHUqrfpMVH8twdsxicEMnmwzXc/dxmnB0uq2MFLG9G7jOBAmPMAWNMG/ACcGUPx/0I+Bmgp7xPU4fLcP8LW1l7oIrUuAievW2WropRIWnwgCieuX0WidHhfJRfwbdXbMOluzqdEW/KfQhQ1O3zYs9jXURkOpBpjHnrVN9IRJaISK6I5FZUVJx22GBkjOGR13bwTt5R4iLDeOa2mQxNjrY6llKWGZkWy99unUm0w85rW0v40Vu79G6SZ+CsT6iKiA34FfBAb8caY5YZY3KMMTmpqaln+9RB4Q8fFrB8QxERYTaevmWG7nWqFDAlcwDLbsrBYbfx1zUHeXrNQasjBRxvyv0IkNnt8wzPY53igInARyJyEJgNrNSTqr17fesRfvn+XkTg94unMSM7yepISvmNeaNS+MV1UwD4f2/t4t28oxYnCizelPtGYJSIDBMRB7AIWNn5RWNMrTEmxRiTbYzJBtYBC40xuX2SOEhsKKzm2yu2A/DIFeO5ZMJAixMp5X8WThnMty8dgzFw3wtb2F5cY3WkgNFruRtj2oF7gXeB3cCLxpg8EXlMRBb2dcBgVFjZyJJnc2nrcHHznCxuOzfb6khK+a27zx/Btedk0OJ0cfvfcyk+pptue0OsOlGRk5NjcnNDb3Bf2+zkas9m1heOTWPZTecQZtdryZQ6lbZ2F7f8dQOf7q9i7MA4Xrl7LtGO0NxrSEQ2GWN6nfbWVulH7iWPWzhQ2cjYgXH8bvE0LXalvOAIs/G/N57D8BT3VazfXrFdV9D0QpulH/3q/XxW51cwIDqcZTflEKubWSvltYSocJZ91f3v5q0dpTz50X6rI/k1Lfd+8tb2Upau3o9N4A+Lp+tadqXOwMi0WH5z/VRE4Bfv5fPhnjKrI/ktLfd+kH+0nm+t2AbAQ5ePY96oFIsTKRW4Lh6fzn9fPNq9gmb5VgorG62O5Je03PtYY2s7dz23iWZnB1dPG8Lt84ZZHUmpgHfPBSNZMGEg9a3t3PPcZlqcHVZH8jta7n3IGMNDr+7gQEUjo9Nj+fHVExERq2MpFfBsNuHn104mKzmaXaV1/PCNXVZH8jta7n3o+Q2HeX1rCdEOO09+ZXrILt1Sqi/ER4az9IbpOMJsLN9wmFe3FFsdya9oufeRnUdqu0YTj39pEiPT4ixOpFTwmTgkgR98cQIAD72yk31l9RYn8h9a7n2gqa2d/1q+hbZ2FzfMGsqVU4f0/oeUUmdk8cxMrpo6mGZnB99YvkXn3z203PvAj97cxYHKRsakx/H9L5ywr4lSyodEhB9fPYns5Gj2HK3niXfzrY7kF7TcfeydnUdZvqEIR5iN3y6eSmS43epISgW9mIgwfrNoGnab8JdPCvl4n+4XoeXuQ2V1LTz4ivtOj9+9bCxjB+q92ZXqL1MzB/DNi0cB8MCL26huDO1NtrXcfcQYw7dWbKOmycnnRqdyy9xsqyMpFXLuOn8kM7ITKa9v5cGXQ/v+M1ruPvL8hsN8vK+SxOhwnrh2sq5nV8oCdpvw6+unEhcRxnu7ynh9a4nVkSyj5e4DRdVN/OSt3QD86KqJurm1UhbKSIzmkS+MA+DRlXmU17VYnMgaWu5nyRjDg69sp7Gtg8snDeQLkwdbHUmpkHddTiafG51KbbOTh17dGZLTM1ruZ+n5DYdZU1BFUoyDx66caHUcpRTu5ZE//fIk4iLC+Nfu0Jye0XI/CyU1zV3TMY9dOYGU2AiLEymlOg1KiOJ7nutMHl2ZR2VDq8WJ+peW+1l4dGUejW0dXDohXadjlPJD1+ZkMH9UCrXNTn7sGYiFCi33M/Ru3lHe31VGbEQYP1yo0zFK+SMR4cdXTSIizMarW47wyb5KqyP1Gy33M9DQ2s6jr+cB8O1LxzAwQVfHKOWvhiZHc5/n4qZHXtsRMvee0XI/A798L5+jdS1MyUjgxtlZVsdRSvXizvnDGZMex8GqJpauLrA6Tr/Qcj9NeSW1/P3Tg9htwk++NAm7TS9WUsrfhdtt/ORL7unTP/57P/srGixO1Pe03E+DMYYfrMzDZeDmOdlMGJxgdSSllJfOyUri+pxMnB2Gx97YFfRr37XcT8PKbSVsPHiMlFgH939+lNVxlFKn6dsLxhAXGca/91bw4Z5yq+P0KS13LzW2tvP4qj0A/M+lY4mPDLc4kVLqdKXERvDNi0cD8Nibu2htD96Tq1ruXnryowKO1rUwOSOBa87JsDqOUuoM3TQni1FpsRyqauIvnxRaHafPaLl7oai6iT//n/uH4AcLJ2DTk6hKBaxwu41HPfuu/uHDAsrrg/PGYl6Vu4gsEJF8ESkQkQd7+PrXRWSHiGwVkU9EJKj2lvvFe/m0dbi4aupgpg9NtDqOUuoszRuVwsXj0mlq6+A3/9pndZw+0Wu5i4gdWApcBowHFvdQ3s8bYyYZY6YCPwd+5fOkFtlRXMvrW0twhNn41qVjrI6jlPKRBy8bi90m/HNjEQXl9VbH8TlvRu4zgQJjzAFjTBvwAnBl9wOMMXXdPo0BgmKNkTGGn6xy34/i1rnZZCRGW5xIKeUrI9NiuX5GJh0uw8/eCb5Ntb0p9yFAUbfPiz2PfYaI3CMi+3GP3P/LN/Gs9VF+BWsPVJEQFc7d54+0Oo5Sysfuv3gU0Q477+8qY0NhtdVxfMpnJ1SNMUuNMSOA7wCP9HSMiCwRkVwRya2o8O/dyTtchp++7V76+I0LR5IQrUsflQo2aXGRLDlvOACPv707qC5s8qbcjwCZ3T7P8Dx2Mi8AV/X0BWPMMmNMjjEmJzU11fuUFnhzewn5ZfUMTojkpjl6/xilgtWd84eTHONgy+EaVucHz4VN3pT7RmCUiAwTEQewCFjZ/QAR6X655hVAQJ9+bu9w8VvPGfRvXDSKiDC7xYmUUn0lJiKMu84fAcCv3t8bNKP3XsvdGNMO3Au8C+wGXjTG5InIYyKy0HPYvSKSJyJbgf8Gbu6zxP3gta0lHKhsZGhStF6wpFQIuHF2FmlxEew8Use7eWVWx/GJMG8OMsasAlYd99j3u318n49zWcbZ4eJ3H7hH7fddNIpwu17npVSwiwy3c88FI3l0ZR6/fn8vl4xPD/iLFbW5jvPypmIOVzcxPDWGq6adsChIKRWkFs3MZHBCJPll9by1o9TqOGdNy72b9g4XT360H3CP2vVe7UqFjogwO/dc6F7yvHR1QcDPvWu5d/PWjlIOVzeRnRytG14rFYKuOSeD9PgI9hyt54Pdgb1yRsvdw+UyXdtv3XX+CB21KxWCIsLs3Dnfve79DwE+etdy9/jX7jL2ljUwKCGSq6fpChmlQtUNs4aSGB3O1qIa1h6osjrOGdNyx30PmaWeufYl5w3HEaYvi1KhKtoRxm3nDgMI6M20tcWA9YXVbCuqISnGwaIZQ62Oo5Sy2FfnZBMbEcaagip2Hqm1Os4Z0XIHnvrYvRHHTbOziHLo1ahKhbqE6HCun+G+60qg7tYU8uVeWNnIB3vKcITZuHG23kNGKeV2y9xsbAJvbCvhaG3g7dYU8uX+1zWFGANXTx1CalyE1XGUUn4iMymayyYOot1leGbtQavjnLaQLveapjZW5BYDcPv8YRanUUr5m85eeG79YZra2i1Oc3pCutyXbyii2dnB/FEpjE6PszqOUsrPTB+ayLShA6htdvLy5lPd6dz/hGy5d7gMz60/BNC17EkppY7X2Q//WHsooC5qCtly/7+9FRQfayYzKYrPjfbvjUOUUta5dMJAUmIjyC+rJ/fQMavjeC1ky/3Zde5R+1dmZQX8rT2VUn3HEWZjkWdZ5LNrD1mcxnshWe5F1U2szi/HYbdxrW7GoZTqxeJZQ7EJvL2zlMqGVqvjeCUky/35DYcxBq6YPIjkWF3+qJQ6tSEDorhwbDrODsM/NxZZHccrIVfuzg5X1/LHr8zSWw0opbxz42x3X/xzYxEul/+fWA25cv8ov4LKhlZGpsVyTlai1XGUUgFi/qhUBidEcri6ifWF1VbH6VXIlfuKXPdbqmvPyUBET6Qqpbxjtwlf9pyj6+wRfxZS5V5R38qHe8qx24Srp+v+qEqp03ONp9xX7SylrsVpcZpTC6lyf23LEdpdhgvGpJEWF2l1HKVUgMlKjmH28CRanC7e2u7fm2iHTLkbY1ixyf1W6rocXf6olDoz1+W417y/6OdTMyFT7rtL69lb1kBSjIMLxqZZHUcpFaAumziIaIedLYdrOFTVaHWckwqZcn99m/umP1dMGkS4PWT+t5VSPhblsHPphIEArNxaYnGakwuJlnO5DG94/hKunDrY4jRKqUC30NMjr2094rc3EwuJcs89dIyS2haGDIhi+lBd266UOjvzRqaQFONgf0Uju0rrrI7To5Ao99e3uqdkFk4drDcJU0qdtXC7jSsmDQJg5Tb/nJoJ+nJv73Cxaod7ydLCKTolo5Tyjc4p3je2lvjl1IxX5S4iC0QkX0QKROTBHr7+3yKyS0S2i8gHIuI3O02vL6zmWJOTEakxjBsUb3UcpVSQmD40kfT4CEpqW9heXGt1nBP0Wu4iYgeWApcB44HFIjL+uMO2ADnGmMnAS8DPfR30TL290z1qXzBxoMVJlFLBxGYTFnhWzby986jFaU7kzch9JlBgjDlgjGkDXgCu7H6AMWa1MabJ8+k6wC+uEnK5DO/mlQHutalKKeVLCzy98s7OUr+bmvGm3IcA3S/FKvY8djK3A2+fTShf2Xz4GBX1rWQkRjFhsE7JKKV8a0Z2IkkxDg5WNZFfVm91nM/w6QlVEbkRyAGeOMnXl4hIrojkVlRU+PKpe9T5VmnBhIF6B0illM+F2W1cMj4dgHf8bGrGm3I/AmR2+zzD89hniMjFwMPAQmNMj/tQGWOWGWNyjDE5qal9uym1MYb3dnnKXefblVJ9pLNfOqeA/YU35b4RGCUiw0TEASwCVnY/QESmAX/CXezlvo95+vZXNFBU3UxSjINpeuGSUqqPzB6eTLTDzu7SOkprm62O06XXcjfGtAP3Au8Cu4EXjTF5IvKYiCz0HPYEEAusEJGtIrLyJN+u33y4x/075vzRqdj1wiWlVB+JDLdz7sgUAFbv6fvpZm+FeXOQMWYVsOq4x77f7eOLfZzrrH2w213uegdIpVRfu3BsGu/vKuPDPeXc4Cd7MwflFaq1zU5yDx3DbhPOG923c/tKKXXBGPcgck1BJS3ODovTuAVluX+8r4IOlyEnK5GEqHCr4yilgtzAhEgmDI6n2dnBugNVVscBgrTc/53vnvfSKRmlVH+50NM3H+X7x7x70JW7MYY1BZUAzB+VYnEapVSomOc5qfrp/kqLk7gFXbkXVjZSUttCcoyDcQP1qlSlVP+YNjSRaIedvWUNlNe1WB0n+Mr9E8+ofe7IFL13u1Kq3zjCbMwalgTAGj8YvQdfue9zv6jzRiZbnEQpFWo617t/vE/L3afaO1ys9Zyp7nyRlVKqv8zznOdbU1Bp+V0ig6rcd5XWUd/STlZyNBmJ0VbHUUqFmDHpcaTEOiira6WwstHSLEFV7hsKqwGYmZ1kcRKlVCgSEWZ4+mfjwWpLswRluc8YpuWulLJGZ7lvKDxmaY6gKXdjDLmH3C/mLC13pZRFZnr6Z8NBa69UDZpy31/RQHVjG2lxEQxN0vl2pZQ1xg2KJzYijKLqZo7WWrfePWjKfX23KRnddUkpZRW7TTgny72HxAYL592Dptw3H6oBYEaWbsyhlLLWjGx3D23Scj97W4vc8+2665JSympTM909tLW41rIMQVHudS1O9lc04rDbGDsozuo4SqkQNzkzAYDdJXW0tltzf/egKPftRe7fjuMGxxMRZrc4jVIq1MVHhjMiNYa2Dhe7S+styRAU5b6t2D3fPi1zgMVJlFLKrXNqZltRjSXPHxTlvtXz4k3xvBVSSimrTfX00VYt9zO3w3PSYkqGjtyVUv5himcmYXuxlvsZqW5s42hdC9EOO9nJMVbHUUopAEanx2G3CYWVjTS39f9J1YAv992ldQCMHRinm3MopfxGZLidkamxuAzkl/X/SdWAL/ddJe5yHz9Yt9RTSvmXzl7q7Kn+FPjl7hm5jx+kJ1OVUv5l/CBPuZf2/8VMgV/uOnJXSvmpcZ5yt2Kte0CXe1u7i/0VDYi4d0BRSil/Ms5zxfye0rp+33YvoMv9UFUj7S5DZmI0UQ69MlUp5V+SYyNIjnHQ2NbB0br+vf1vQJf7vvIGAEamxVqcRCmlejbC00/7yhr69Xm9KncRWSAi+SJSICIP9vD180Rks4i0i8g1vo/ZswJPuY/ScldK+anOwWdnX/WXXsvR41YmAAAJyElEQVRdROzAUuAyYDywWETGH3fYYeAW4HlfBzyVzhdrhJa7UspPjUz1lHtF/5Z7mBfHzAQKjDEHAETkBeBKYFfnAcaYg56vufog40kV6LSMUsrP+e3IHRgCFHX7vNjz2GkTkSUikisiuRUVFWfyLbq4XIYDlZ6Re6qWu1LKP3WW+34/LHefMcYsM8bkGGNyUlNTz+p7lde30uJ0kRTjICEq3EcJlVLKtwbGR+IIs1HV2EZDa3u/Pa835X4EyOz2eYbnMUsdrm4CYGhStMVJlFLq5Gw26eqpw1VN/fe8XhyzERglIsNExAEsAlb2bazeHapqBCArWctdKeXfsjrLvbqx356z13I3xrQD9wLvAruBF40xeSLymIgsBBCRGSJSDFwL/ElE8voyNOjIXSkVOIZ6BqGH+nHk7s1qGYwxq4BVxz32/W4fb8Q9XdNvOl8kLXellL/rHLkfqvavaRm/pCN3pVSg6By5F2m59660thmAIYlRFidRSqlTGzLAXe4lNc399pwBWe7ODhfl9a2IQHp8pNVxlFLqlAYmuHuqtLal3+4OGZDlXl7fijGQGhtBuD0g/xeUUiEkPjKMGIedprYO6lr6Z617QDbjUc+UzKAEHbUrpfyfiHSN3o/W9s+tfwOy3Es9L86gBJ1vV0oFhs6+Kqntn3n3gCz3zt98A3XkrpQKEF3z7jU6cj+pivpWQE+mKqUCR3p8BACVDa398nwBWe6VDW0AJMc4LE6ilFLeSYpxl3uVlvvJVTW6X5zkWC13pVRgSPH0VWVjW788X2CWe+fIPTbC4iRKKeWdZB25967zxdFpGaVUoOicaegcnPa1gCt3Y0zX2xqdllFKBYrOvqrWaZmeNbV10NbuIiLMRrTDq5taKqWU5RKj3eV+rKmtX25BEHDlXu+5dFe31lNKBZJwu41ohx2Xgca2jj5/voAr97oWJwDxWu5KqQATH+nurbpmZ58/V8CVe72n3OMidUpGKRVYOnurvh9uHhZw5V7X7H5ROn8DKqVUoOicceicgehLgVfuOi2jlApQ8Z6Ru07L9KDz7UxshE7LKKUCS6xnxkGnZXrQ4nSfZY4Kt1ucRCmlTk9UuLtyO3usLwVsuUeGB1x0pVSI6xyUarn3oFlH7kqpABXp6a1mp6vPnyvgyr3F86JEarkrpQJMpI7cT06nZZRSgaqr3Nu13E/QOXKP0JG7UirAdA5KW/T2Aydqd7nLPdwuFidRSqnTE2Z3V267S28cdoIOz4tiEy13pVRgsXt6y+Uvd4UUkQUiki8iBSLyYA9fjxCRf3q+vl5Esn0dtFPni2K3abkrpQKLZ+DeNUjtS72Wu4jYgaXAZcB4YLGIjD/usNuBY8aYkcCvgZ/5OminzhfFriN3pVSA6Zxx6Oj7lZBejdxnAgXGmAPGmDbgBeDK4465Evi75+OXgItE+qZ9O18Um47clVIBpnPGwV+mZYYARd0+L/Y81uMxxph2oBZI9kXA43W+KGFa7kqpANNZ7n4xLeNLIrJERHJFJLeiouKMvkdsRBhJMQ4iwnQppFIqsESE2UmKcRAT0ff95c2tFY8Amd0+z/A81tMxxSISBiQAVcd/I2PMMmAZQE5Ozhn96vrd4mln8seUUspyCyYOZMHEgf3yXN6M3DcCo0RkmIg4gEXAyuOOWQnc7Pn4GuBD0x87wCqllOpRryN3Y0y7iNwLvAvYgaeNMXki8hiQa4xZCfwFeFZECoBq3L8AlFJKWcSrHS+MMauAVcc99v1uH7cA1/o2mlJKqTMVcFeoKqWU6p2Wu1JKBSEtd6WUCkJa7kopFYS03JVSKgiJVcvRRaQCOHTcwylApQVxToe/Z/T3fOD/Gf09H2hGX/D3fNBzxixjTGpvf9Cycu+JiOQaY3KsznEq/p7R3/OB/2f093ygGX3B3/PB2WXUaRmllApCWu5KKRWE/K3cl1kdwAv+ntHf84H/Z/T3fKAZfcHf88FZZPSrOXellFK+4W8jd6WUUj7gd+UuIlNFZJ2IbPVs7DHT6kzHE5FviMgeEckTkZ9bnedkROQBETEikmJ1luOJyBOe13C7iLwqIgOszgS9bwZvNRHJFJHVIrLL8/N3n9WZeiIidhHZIiJvWp2lJyIyQERe8vwM7haROVZn6k5Evun5+90pIstFJPJ0v4fflTvwc+CHxpipwPc9n/sNEbkA956xU4wxE4BfWBypRyKSCVwCHLY6y0m8D0w0xkwG9gLftTiPt5vBW60deMAYMx6YDdzjhxkB7gN2Wx3iFH4LvGOMGQtMwY+yisgQ4L+AHGPMRNy3Wj/t26j7Y7kbIN7zcQJQYmGWntwF/NQY0wpgjCm3OM/J/Br4H9yvp98xxrzn2W8XYB3uHb6s5s1m8JYyxpQaYzZ7Pq7HXUrH72lsKRHJAK4AnrI6S09EJAE4D/c+FBhj2owxNdamOkEYEOXZ2S6aM+hBfyz3+4EnRKQI96jY8hHdcUYD80VkvYj8W0RmWB3oeCJyJXDEGLPN6ixeug142+oQeLcZvN8QkWxgGrDe2iQn+A3ugYXL6iAnMQyoAP7qmTp6SkRirA7VyRhzBHf3HQZKgVpjzHun+3282qzD10TkX0BPGwk+DFwEfNMY87KIXIf7t+vFfpQvDEjC/ZZ4BvCiiAzv720Fe8n4EO4pGUudKqMx5nXPMQ/jnmp4rj+zBToRiQVeBu43xtRZnaeTiHwBKDfGbBKR863OcxJhwHTgG8aY9SLyW+BB4HvWxnITkUTc7xiHATXAChG50Rjzj9P5PpaUuzHmpGUtIs/gnq8DWIEFb+16yXcX8IqnzDeIiAv3/R8q+isfnDyjiEzC/UOxTUTAPd2xWURmGmOO9mPEU76OACJyC/AF4CI/2XPXm83gLSci4biL/TljzCtW5znOucBCEbkciATiReQfxpgbLc7VXTFQbIzpfMfzEu5y9xcXA4XGmAoAEXkFmAucVrn747RMCfA5z8cXAvsszNKT14ALAERkNODAj24+ZIzZYYxJM8ZkG2Oycf8gT+/vYu+NiCzA/dZ9oTGmyeo8Ht5sBm8pcf/G/guw2xjzK6vzHM8Y811jTIbnZ28R8KGfFTuefwtFIjLG89BFwC4LIx3vMDBbRKI9f98XcQYnfC0ZuffiTuC3nhMJLcASi/Mc72ngaRHZCbQBN/vJqDPQ/AGIAN73vMNYZ4z5upWBTrYZvJWZenAucBOwQ0S2eh57yLPPsfLeN4DnPL/EDwC3Wpyni2eq6CVgM+4pyy2cwZWqeoWqUkoFIX+cllFKKXWWtNyVUioIabkrpVQQ0nJXSqkgpOWulFJBSMtdKaWCkJa7UkoFIS13pZQKQv8f49PEcBoigKQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VGXax/HvPakkhCSEAJKEptSEJhHEtmsHV4OKBSyra1vffW0sFtQVFRuW1dVdt1hfyyoiIiKgWNYuKKEFAgQwtFBCqGmkTOZ+/5iBzbLBDDDJmZncn+vy2szJIfNbCD9OnvOc5xFVxRhjTHhxOR3AGGNM4Fm5G2NMGLJyN8aYMGTlbowxYcjK3RhjwpCVuzHGhCErd2OMCUNW7sYYE4as3I0xJgxFOvXG7dq1065duzr19sYYE5IWLFiwXVVTGzvPsXLv2rUrubm5Tr29McaEJBFZ7895NixjjDFhyMrdGGPCkJW7McaEIb/KXUSGi0iBiKwRkfENfP4ZEVns+2+ViOwOfFRjjDH+avSGqohEAM8DZwJFwHwRmaGqy/edo6pj651/MzCoCbIaY4zxkz9X7kOANapaqKo1wGRg5M+cPwZ4OxDhjDHGHB5/yj0N2FjvdZHv2H8RkS5AN+BfRx7NGGPM4Qr0PPfRwFRVrWvokyJyA3ADQOfOnQP81sb4R1XZsqeKdTsqWLe9kvU7Kli3o4JdlbW+E/59blJcFF3bxdMlJY6uKd7/TUtqhYg4E94YP/lT7puAjHqv033HGjIa+N+DfSFVfQF4ASA7O9s2bzXNZk9lLd+u2c5Xq7bx1aoSikurD/trtU+I4ZSeqZzSM5WTj2lHcnx0AJMaExj+lPt8oIeIdMNb6qOByw48SUR6A8nA3IAmNOYwVVS7mb54E+8v3MTCDbvw1LucSI6L4ujU1nRJiadbuzi6pMTTrnUM+y7IBe8F/PbyatbvqGTd9grW76zkp23lbCurZuqCIqYuKEIEBqQnceGxaVx4bDqtYxx76NuY/9Dod6KqukXkJmAOEAG8oqr5IjIRyFXVGb5TRwOTVdWuyI2jCkvKeWPeeqbmFlFW7QYg0iUM6ZbML3q25xc9U+lzVMJhDa2oKgXFZXxVUMLXq0uYv3YXizfuZvHG3Tz+0UpGDU7nyuO70KNDQqD/bxlzSMSpLs7OzlZbW8YE0qINu3j601V8s3r7/mODuyRz5fFdOL1PexJiowL+nhXVbj5fuY03567nx3U79x8/4egUfn9mT7K7tg34e5qWTUQWqGp2o+dZuZtQt3FnJY9/vJKZeVsAiIl0cf7ANK4c1oWstMRmy7Fyaymvz13P9EWbqKzxzikYntmRu0b0plu7+GbLYcKblbsJe3sqa/nLF6t57fv11NR5iIl0ce1J3bjhlO4kxTl3k7O0qpaXvi7khW8Kqar1EOkSrji+C7ec3oO2dvPVHCErdxPWZuZt5r7py/ZPX7xwUBrjzu5FWlIrh5P929Y9VTz9aQHvLihCFdrERjJxZBYjB3ayqZTmsFm5m7C0u7KGCR/kM2PJZgCGdmvLfef2bdbhl0O1Ykspj85esf9ewIisjjx8fhYprWMcTmZCkZW7CTtfFmzjrvfyKC6tplVUBPf8qg9XDO0cElfBqsqU3I1M/HA5FTV1tGsdzWMX9ufMvh2cjmZCjJW7CRs1bg8PzVzOG/O8G9AM7pLMHy8eQNcQvEm5cWcld0xdwrxC78ya0cdl8EBOJrFREQ4nM6HC33K39dxNUCspq+ayF+fxxrz1REe4uGt4b6b8dlhIFjtARts43rrueCac25eYSBeT52/k0hfmUVxa5XQ0E2as3E3QyivaTc5fviV3/S46tonl3RuH8T+/PJoIV/APw/wcl0u45qRuTPvdCaQltWLJxt2c++dvWbB+l9PRTBixcjdB6f1FRVz897ls2VPF4C7JzLj5RAZkJDkdK6AyOyUy46YTOb57W0rKqhnzwjzemb/B6VgmTFi5m6Ciqjw5ZyVj31lCtdvDmCEZvHX9UNonxDodrUmktI7hjWuHcvUJXamp83DXe0t5ZNZybBUPc6RslSMTNOo8yn0fLOOtHzYQ4RIeyMkMmdkwRyIqwsUDOZn0PaoN905fyovfrGVXZS2TLuxHZIRdf5nDY+VugkKN28PYKYuZlbeFmEgXz192LGe0sGmClxyXQYfEWG58YwFTFxRRureW58YMspk05rDYZYFxXGWNm+tfz2VW3hYSYiJ5/ZohLa7Y9/lFz1TevG4obWIj+WR5Mb95dT7lvpUtjTkUVu7GUXv21nLlyz/y1aoSUuKjefuG4xnaPcXpWI4a3CWZKTcOIzUhhrmFO7jsxXnsqqhxOpYJMVbuxjHl1W6ufvVHFqzfRafEWKbcOCyolxFoTr07tuG9G08go20r8or2cNWrP1JaVet0LBNCrNyNI6pq67jutfks2rCbtKRWTLlxGEentnY6VlDpnBLHu789gc5t48gr2sNvXp1PhQ3RGD9ZuZtmV+2u47dvLGBe4U7aJ8Tw1vVDSU+OczpWUOqYGMs/rxtKp8RYFqzfxfWv51JV2+D+88b8Byt306xq6zzc/Nai/WPsb10/lC4pobmUQHPJaBvHP68/ntSEGL7/aQc3vrmAarcVvPl5Vu6m2Xg8yu3vLuGT5cW0iY3kjWuHckx722vUH93axfPWdUNpGx/NlwUl3Pr2Yuo89qCTOTi/yl1EhotIgYisEZHxBznnEhFZLiL5IvJWYGOacPDYRyv4YPFmWsdE8vq1Q+nbqY3TkUJKjw4JvH7NENrERvJx/lYmfphvT7Kag2q03EUkAngeGAH0BcaISN8DzukB3A2cqKqZwG1NkNWEsNfnruPFb9YS6RL+ceVgBobZOjHNJSstkRd/nU10hIvX5q7n5W/XOh3JBCl/rtyHAGtUtVBVa4DJwMgDzrkeeF5VdwGo6rbAxjSh7LPlxTwwIx+ASaP6c+Ix7RxOFNqGdk/hqUsGAPDwrBXMXrrF4UQmGPlT7mnAxnqvi3zH6usJ9BSR70RknogMD1RAE9qWbNzNzW8vwqMw9oyeXDQ43elIYSFnQCfuGt4bgNveWcyC9TsdTmSCTaBuqEYCPYBfAmOAF0Xkv37uFpEbRCRXRHJLSkoC9NYmWG3cWcm1r81nb20dFw1O55bTj3E6Uli58RfduXxoZ2rcHq57LZfCknKnI5kg4k+5bwIy6r1O9x2rrwiYoaq1qroWWIW37P+Dqr6gqtmqmp2amnq4mU0IqKh2c91ruWwvr+GkY9rx2IX9wn51x+YmIjyYk8lpvduzq7KW617PtadYzX7+lPt8oIeIdBORaGA0MOOAc6bjvWpHRNrhHaYpDGBOE0I8HmXclCUUFJdxdGo8f73iWKJs6domERnh4s9jBtG7YwKFJRWMnbwYj02RNPhR7qrqBm4C5gArgCmqmi8iE0Ukx3faHGCHiCwHvgDuUNUdTRXaBLfnv1jDx/lbSYiJ5MVfZ9MmNsrpSGEtPiaSF67MJikuis9XbuOZz1Y5HckEAXFqnmx2drbm5uY68t6m6Xy2vJjr3/D+ub58VTan9W6ZS/c64dvV2/n1Kz/gUfjb5ccyot9RTkcyTUBEFqhqdmPn2c/KJmDWbCtn7DuLUYXbz+plxd7MTurRjnvO6QPAuHeXsHJrqcOJjJOs3E1AlFXVcsMbuZRVuzmnX0d+98ujnY7UIl17UjfOH9iJypo6bnh9AXsq7QZrS2Xlbo6YqjL+vaUUllTQq0MCT140wGbGOEREmDSqP1lpbdiws5Lbpy6xJQpaKCt3c8TemLeeWUu3EB8dwd+vHEx8jG3N66TYqAj+dvlgEmIj+XR5sS1R0EJZuZsjsrRoDw/PXAHAY6P6062dLd8bDDLaxvHkRd4lCiZ9tJKFG3Y5nMg0Nyt3c9hKq2r537cWUlPn4fKhnckZ0MnpSKae4VkduebEbrg9ys1vLWJ3pe3D2pJYuZvDoqrcNTWPDTsr6XtUG+47t2/jv8g0u/EjejMgI4lNu/dy+7s2/t6SWLmbw/L63PV8tGwrrWMi+evlxxIbFeF0JNOA6EgXfxkziDaxkXy2YhsvfWPj7y2Flbs5ZCu2lPLILO84+6RR/ehq4+xBLaNtHH+8ZCAAT8xZybJNexxOZJqDlbs5JFW1ddw2eTE1dR7GDMng3P42zh4KzuzbgauGdaG2Trll8iL21tgerOHOyt0ckkkfraSguIzu7eJtnD3E3H1OH3q0b01hSQWPzF7udBzTxKzcjd++LNjG/32/jkiX8KfRA4mLtvnsoSQ2KoJnRw8iOsLFm/M28NnyYqcjmSZk5W78sqO8mtvfzQNg7Jk96Z9ue6CGor6d2nDn8F4A3PleHtvKqhxOZJqKlbtplKpy13t5bC+vZki3ttz4C1s3JpRdc2I3TjqmHTsrarj93Txb/z1MWbmbRr0zfyOfrdhGQmwkz1w6kAiXrRsTylwu4Y+XDCApLoqvV5Xw5g/rnY5kmoCVu/lZG3dW8tBM7823h8/PIi2plcOJTCB0aBPLYxf0A+Cx2StZt73C4UQm0KzczUF5PMqdU/OoqKljRFZHW14gzIzodxQ5Azqxt7aOO6Yuoc6GZ8KKlbs5qNfnrmNu4Q5S4qN5+PwsW8Y3DE0cmUlqQgzz1+3i1e/s6dVwYuVuGrR2ewWTPl4JwCMXZJHSOsbhRKYpJMVF8/go7/DME3MKWLOtzOFEJlCs3M1/qfMo46YspqrWw/kDOzE8y/biDGen9e7AxYPTqXF7GDdlCe46j9ORTAD4Ve4iMlxECkRkjYiMb+DzV4tIiYgs9v13XeCjmuby0jeFLNywmw5tYngwJ8vpOKYZ3HdeXzolxrKkaA9//+onp+OYAGi03EUkAngeGAH0BcaISEPPnb+jqgN9/70U4JymmRSWlPP0p6sAmHRhfxLjohxOZJpDm9goHr+oPwDPfb6G1cU2PBPq/LlyHwKsUdVCVa0BJgMjmzaWcYLH490LtdrtYdSx6Zzau73TkUwzOrlHKqOPy6CmzsNd7+XZ7JkQ50+5pwEb670u8h070CgRyRORqSKSEZB0pln984f1/LhuJ+1ax3DfuX2cjmMccPc5fejQJoaFG3bz2vfrnI5jjkCgbqh+CHRV1f7Ap8BrDZ0kIjeISK6I5JaUlATorU0gFO2qZNJH3tkxD43MJCku2uFExgmJraJ4+Hzv7Jkn5xSwcWelw4nM4fKn3DcB9a/E033H9lPVHapa7Xv5EjC4oS+kqi+oaraqZqemph5OXtMEVJV731+2/2GlEf1sdkxLdmbfDpzne7jp7mlLbWu+EOVPuc8HeohINxGJBkYDM+qfICL12yAHWBG4iKapTVu4ia9WlZDYKooHR2Y6HccEgQfO60tyXBTfrtnOlNyNjf8CE3QaLXdVdQM3AXPwlvYUVc0XkYkikuM77RYRyReRJcAtwNVNFdgE1vbyaib61o6579y+tE+IdTiRCQYprWN4IMf7D/3Ds1awrdSWBg41fo25q+psVe2pqker6iO+YxNUdYbv47tVNVNVB6jqqaq6silDm8B5aOZy9uyt5eQe7Rh1bEP3yU1LlTOgE6f1bk9ZlZsHP7Sdm0KNPaHagn1ZsI0PFm8mNsrFoxf0s7VjzH8QESaOzKRVVASzlm7h8xW2c1MosXJvoSpr3Pxh+jIAxp7Rk4y2cQ4nMsEoPTmOcWf1BGDCB/lUVLsdTmT8ZeXeQj372WqKdu2lz1FtuOakbk7HMUHs6hO6kpXWhk279+5/etkEPyv3Fih/8x5e+nYtIjDpwn5ERdi3gTm4yAgXky7sj0vg1e/WsrRoj9ORjB/sb3ULU+dR7p62lDqPctWwrgzIsI2uTeOy0hK55sRueBTGT8uzlSNDgJV7C/PG3HXkFe3hqMRYbj+7l9NxTAgZe2ZP0pJakb+5lP+zpQmCnpV7C1JcWsVTn3jHTB/MyaR1TKTDiUwoiY+JZKLvIbdnPl3Flj17HU5kfo6Vewvy0MzllFe7OaNPB87K7Oh0HBOCTu/TgbMzO1BRU8dEm/se1KzcW4ivV5UwM28LraIieCCnoeX4jfHP/edlEhcdwUfLtvJFwTan45iDsHJvAapq67jvA++c9lvP6EF6ss1pN4evU1Irxp6xb+77Mqpq6xxOZBpi5d4C/PXLn1i/o5KeHVpzrc1pNwFw9Yld6d0xgY079/KXf61xOo5pgJV7mCssKefvX3r3xHz4fJvTbgIjKsLFIxd499f9x9c/sWZbucOJzIHsb3oYU1UmfJBPTZ2HiwenM6RbW6cjmTAyuEtbRh+XQW2dMuGDZbbue5Cxcg9js5du5ds120lsFcX4Eb2djmPC0F3De5McF8X3P+3gw7wtTscx9Vi5h6nyajcP+dZpv3N4L1JaxzicyISj5Pho7hruvXB42DfV1gQHK/cw9dznq9laWsWA9ERGH9fZ6TgmjF2SncHAjCS2lVXzJ1tYLGhYuYehVcVlvOJbGOyh87OIcNk67abpuFzCw+dnIQKvfr+OlVtLnY5ksHIPO6rKfdOX4fYolw3pTP90WxjMNL2stESuGNqFOo8yYXq+3VwNAlbuYeaDxZv5Ye1O2sZHc4ctDGaa0e1n9SIlPpof1+3k/UWbnI7T4vlV7iIyXEQKRGSNiIz/mfNGiYiKSHbgIhp/lVXV8sjsFQCMH96bpLhohxOZliQx7t+zsh6dvZLSqlqHE7VsjZa7iEQAzwMjgL7AGBH5r8VJRCQBuBX4IdAhjX+e/Ww1JWXVDOqcxEWD052OY1qgUcemM7hLMtvLq/nTp6udjtOi+XPlPgRYo6qFqloDTAZGNnDeQ8DjQFUA8xk/rSou49Xv13lvoo7MwmU3UY0DXC7hwZxMXAKvzV1HwdYypyO1WP6Uexqwsd7rIt+x/UTkWCBDVWcFMJvxk6py/wf51PluomalJTodybRgWWmJXL7v5qo9ueqYI76hKiIu4GlgnB/n3iAiuSKSW1JScqRvbXxm5m1hbuEOkuOi7CaqCQrjzupJclwUP6zdyYwlm52O0yL5U+6bgIx6r9N9x/ZJALKAL0VkHXA8MKOhm6qq+oKqZqtqdmpq6uGnNvtVVLt5ZJb3JuqddhPVBImkuH8/ufro7BX25KoD/Cn3+UAPEekmItHAaGDGvk+q6h5VbaeqXVW1KzAPyFHV3CZJbP7Dn/+1hq2lVfRPT+SS7IzGf4ExzeSS7AwGpCdSXFrNnz+3m6vNrdFyV1U3cBMwB1gBTFHVfBGZKCI5TR3QHNxPJeW8/G0hABNH2pOoJri4XMLEkd4nV1/+dq0tC9zM/BpzV9XZqtpTVY9W1Ud8xyao6owGzv2lXbU3PVXlwQ+XU1unXOpb28OYYDMgI4lLszNwe5QHP7QnV5uTPaEaoj5dXszXq0pIiI3kjuF2E9UErzvO7kWb2Ei+Wb2dT5YXOx2nxbByD0FVtXVM9C3nO+7MnrSz5XxNEEtpHcO4s7wXIA/NXG57rjYTK/cQ9I+vCinatZfeHRO44vguTscxplGXD+1M744JFO3ay9+/+snpOC2ClXuI2bizkr9+6d2Q+IGcTCJtT1QTAiIjXDyYkwnA3778iY07Kx1OFP6sGULMw7OWU+32kDOgE8d3T3E6jjF+G9o9hZwBnah2e3h41nKn44Q9K/cQ8s3qEubkFxMXHcE95/RxOo4xh+yec/oQFx3BnPxivlltT6k3JSv3EFHj9vDAjHwAbj6tBx0TYx1OZMyh65gYy82n9QDggRn51Lg9DicKX1buIeK179fxU0kF3drFc81JXZ2OY8xhu+akrnRrF89PJRW89v06p+OELSv3ELCttIpnfY9vTzivLzGREQ4nMubwxURGMOE875YQz36+mm2ltkp4U7ByDwGTPl5JebWbM/p04NRe7Z2OY8wRO7VXe87o057yajeTPl7pdJywZOUe5Bas38m0hZuIjnQx4dz/2gDLmJB137l9iY50MW3hJhas3+l0nLBj5R7E6jzK/b6bqL89pTudU+IcTmRM4HRJieeGk7sDMMG32YwJHCv3IDZ5/gaWbSolLakVv/vlMU7HMSbgfnfq0XRKjCV/cymT529wOk5YsXIPUrsqanhyTgEA9/6qD62i7SaqCT9x0ZHc+yvvcOOTcwrYVVHjcKLwYeUepJ76pIDdlbWceEwKI7I6Oh3HmCZzTr+ODOuewu7KWv74aYHTccKGlXsQWrZpD2/9uIFIl/DAeZmI2CYcJnyJCA+OzCTCJbz1wwbyN+9xOlJYsHIPMqrem6iqcPUJXenRIcHpSMY0uZ4dErhqWFc86n1y1Tb1OHJW7kHm/UWbWLB+F+1ax3DrGT2cjmNMs7ntzB60ax3N/HW7mL54k9NxQp6VexAprarl0dneBzruHtGbhNgohxMZ03zaxEZx1/DeADw6eyVlVbUOJwptVu5B5NnPVrO9vJrBXZK5YFCa03GMaXajjk1nUOckSsqqec635IY5PH6Vu4gMF5ECEVkjIuMb+PyNIrJURBaLyLciYo9SHqKCrWX83/frcAlMHJmJy2U3UU3L43IJE3OyEIFXv1vH6uIypyOFrEbLXUQigOeBEUBfYEwD5f2WqvZT1YHAE8DTAU8axrw3UZdR51EuH9qFzE6JTkcyxjH90hO5bEhn3J59kwvs5urh8OfKfQiwRlULVbUGmAyMrH+CqpbWexkP2J/GIfgwbwvzCnfSNj6acWf1dDqOMY67/axeJMVF8f1PO5i9dKvTcUKSP+WeBmys97rId+w/iMj/ishPeK/cb2noC4nIDSKSKyK5JSW2CwtARbWbR2etAODOs3uRFBftcCJjnJccH80dZ/cCvFtLVlS7HU4UegJ2Q1VVn1fVo4G7gD8c5JwXVDVbVbNTU1MD9dYh7bnPV7O1tIoB6Ylckp3hdBxjgsbo4zqTldaGLXuq+MsXa5yOE3L8KfdNQP3WSfcdO5jJwPlHEqqlWF1cxsvfrkUEJo7MspuoxtQT4RIeGpkFwEvfFPJTSbnDiUKLP+U+H+ghIt1EJBoYDcyof4KI1H/a5leAzWFqhKoy4YN83B5lzJDODMhIcjqSMUFnUOdkLs3OoLZO7cnVQ9RouauqG7gJmAOsAKaoar6ITBSRHN9pN4lIvogsBn4PXNVkicPEh3lbmFu4g+S4KO44q5fTcYwJWncO70Viqyi+Wb2dj5bZzVV/RfpzkqrOBmYfcGxCvY9vDXCusFZe7ebhmcsBuGt4b5Lj7SaqMQeT0jqGO87uxR+mL+Ohmcv5Rc9U4mP8qq4WzZ5QdcCzn61iW1k1AzOS7CaqMX4YM6Qz/dIS2bKnij//y26u+sPKvZmt3FrKK9+tQwQespuoxvglwiU8dL73ydWXvim0J1f9YOXejDwe5Q/ve59EvfL4LvRLtydRjfHXwIwkRh/nfXL1vg+W2c3VRli5N6P3FhaR61vOd5zdRDXmkN15di/axkczr3AnHyze7HScoGbl3kx2VdTw2Efe5Xz/8Ks+JLay5XyNOVTJ8dGMH+FdFvjhWSvYs9eWBT4YK/dm8sScAnZW1HB897aMHNjJ6TjGhKyLjk0nu0sy28ureebTVU7HCVpW7s1g4YZdTJ6/gagI4eHzs2xPVGOOgMt3czXCJbw+dx1Li2zP1YZYuTcxd52He6YtRRWuO7k7x7S3PVGNOVJ9jmrDb07w7rl67/Sl1Hns5uqBrNyb2KvfrWPl1jIy2rbiltNsT1RjAmXsmT05KjGWvKI9vDlvvdNxgo6VexPatHsvT/vGBCfmZNEqOsLhRMaEj/iYSB7IyQTgyTkFFJdWOZwouFi5N6EHZuSzt7aOc/p15NTe7Z2OY0zYOatvB87o057yajcTfUt6GC8r9ybySf5WPl1eTOuYSO4/L9PpOMaEJRHhgZxMWkVFMCtvC18WbHM6UtCwcm8C5dVu7p+RD8C4s3rSoU2sw4mMCV/pyXHcdob3ftZ9Hyxjb02dw4mCg5V7E3hqTgFb9nh3V/r1sK5OxzEm7F1zUjd6d0xg4869/Olzm/sOVu4Bt3jjbl6bu44Il/Dohf2IsIXBjGlyUREuJo3q71tYbC35m23uu5V7ANXWeRj/Xp5vTns3MjvZwmDGNJeBGUlcNawrdR7lnmk2993KPYBe+mbt/jntt53e0+k4xrQ4t5/di6MSY1lStIfX565zOo6jrNwDZN32Cv70mXes75Hz+9mcdmMc0Domkgd9c9+fmlPApt17HU7kHCv3AFBVxk/Lo9rt4YJBaZzSM9XpSMa0WGdldmREVkcqauq49/2lLXbdd7/KXUSGi0iBiKwRkfENfP73IrJcRPJE5HMR6RL4qMFr8vyNzCvcSdv4aO47t6/TcYxp8R4cmUmb2Ei+LChpseu+N1ruIhIBPA+MAPoCY0TkwAZbBGSran9gKvBEoIMGq+LSKh6dvQKA+8/rS1vb7NoYx7VPiOXeX/UB4MEP89lRXu1woubnz5X7EGCNqhaqag0wGRhZ/wRV/UJVK30v5wHpgY0ZnFSV+6Yvo6zKzWm925MzwNZpNyZYXJKdwQlHp7CrspaHWuDSBP6Uexqwsd7rIt+xg7kW+OhIQoWKWUu38IlviQFbp92Y4CIiPHZhP2KjXExfvJnPVxQ7HalZBfSGqohcAWQDTx7k8zeISK6I5JaUlATyrZvdjvJqJnzgXWJg/IjedEpq5XAiY8yBuqTEc7tvv+J73l/aorbl86fcNwEZ9V6n+479BxE5A7gXyFHVBge4VPUFVc1W1ezU1NCeUXL/jHx2VtQwrHsKlw3p7HQcY8xB/ObEbgzqnERxaTWPzlrhdJxm40+5zwd6iEg3EYkGRgMz6p8gIoOAf+At9rBflu3jZVuZmbeFVlERPD6qPy5bYsCYoBXhEp68qD/RkS7eyd3I16tCe9TAX42Wu6q6gZuAOcAKYIqq5ovIRBHJ8Z32JNAaeFdEFovIjIN8uZC3q6KG+z5YBsBdw3vROSXO4UTGmMYc0z5h/8qRd09bSllV+A/PRPpzkqrOBmYfcGxCvY/PCHCuoHX/jHxKyqoZ0rWtrfhoTAi54eTufLR0K0s37eHR2St57MJ+TkdqUvaE6iH4aOkWZizZTKuoCJ7TQDjoAAAOb0lEQVS82IZjjAklkREu/njJAKIjXLz94wa+CvPhGSt3P20vr+be6d7hmLvP6U2XlHiHExljDlXPDgmMPdO7qN9dU/PCevaMlbsf9j2stLOihhOOTuGKoS1qdQVjwsr1J3djYEYSW0urwvrhJit3P0xfvImPlm0lPjqCJy6y4RhjQllkhIunLh5ATKSLqQuK+HR5eD7cZOXeiM279+5/WOn+8zJJT7bZMcaEumPat+bO4b0BuHtaXliuPWPl/jM8HuX2d5dQVuXmjD4duDi7RSyZY0yL8JsTujKsewrby2u4e1r4LQ1s5f4z/u/7dXz/0w5S4qN57MJ+tnaMMWHE5RKeumQACTGRfLK8mKkLipyOFFBW7gexqriMSR+vBODRC/uRmhDjcCJjTKClJbXift/OTQ9+uJyNOysb+RWhw8q9AdXuOm6dvJgat4dLstM5O7Oj05GMMU1k1LFpjMjqSHm1m7HvLA6bjbWt3Bvwx09WsWJLKV1S4rj/vEyn4xhjmpCI8OgF/WifEEPu+l38/aufnI4UEFbuB/h+zXZe/KaQCJfwzKUDiY/xa4UGY0wIS46P5qmLBwDwzKerWLJxt8OJjpyVez27KmoYO2UxqnDzacdwbOdkpyMZY5rJKT1T+c2JXXF7lFsnL6Ki2u10pCNi5e6jqtz5Xh7FpdUM7pLMTace43QkY0wzu2t4b3p3TGDdjkoemJHvdJwjYuXu888fNvDp8mISYiP506UDiYyw3xpjWprYqAieGzOImEgX7y4o4sMlm52OdNiswfBOe9y3xsQjF/Qjo609hWpMS9WzQwJ/OLcv4N2aL1SnR7b4ct9bU8dNby2k2u3hosHp5Azo5HQkY4zDrhjamTP7dqCsys0tkxdRW+dxOtIha/HlPnHmclYVl9M9NZ6JI23aozHGOz3yiVH9OSoxlkUbdvPMp6ucjnTIWnS5z8rbwts/biA60sWfxwwiLtqmPRpjvJLjo/nTpQNxCfztq5/4ZnVobe7RYst9/Y4Kxr+XB8C95/Qhs1Oiw4mMMcFmaPcUbjm9B6ow9p0lbCurcjqS31pkuVe767jprUWUVbsZntmRXw+zzTeMMQ27+bQeHN+9LdvLq0NqeQK/yl1EhotIgYisEZHxDXz+FBFZKCJuEbko8DEDa9JHK1m6aQ/pya14/KL+ttqjMeagIlzCs6MHkRIfzXdrdvDXL9Y4HckvjZa7iEQAzwMjgL7AGBHpe8BpG4CrgbcCHTDQ5uRv5dXv1hHpEv48ZhCJraKcjmSMCXId2sTy9KUDAXjms1X8ULjD4USN8+fKfQiwRlULVbUGmAyMrH+Cqq5T1TwgqOcLbdxZyR3vLgFg/IjeDLLlBYwxfvpFz1R+98uj8SjcMnlR0O/e5E+5pwEb670u8h07ZCJyg4jkikhuSUnz3nmucXu46a2FlPp2Vbr2pG7N+v7GmND3+zN7clzXZIpLqxk7ZQmeIB5/b9Ybqqr6gqpmq2p2ampqc741kz5ayZKiPaQlteKpi22c3Rhz6CIjXDw3ZhDJcVF8vaqEvwXx8sD+lPsmIKPe63TfsZAxJ38rr3y31jvOftkgkuKinY5kjAlRRyW22j/+/sdPCvhx7U6HEzXMn3KfD/QQkW4iEg2MBmY0bazAOXCc3ZbxNcYcqVN7ted/fOPvN7+9MCjH3xstd1V1AzcBc4AVwBRVzReRiSKSAyAix4lIEXAx8A8RCYq1MmvcHm56e5GNsxtjAm7cmT3J7uIdfx/3bvCNv/s15q6qs1W1p6oeraqP+I5NUNUZvo/nq2q6qsaraoqqBsUiLU/OWcmSjbttnN0YE3D7xt+T4qL4sqCEF74pdDrSfwjbJ1S/WLmNF79ZS4RLfH8ANs5ujAmsTkmtePoS7/Z8T80pYNGGXQ4n+rewLPfi0irG+cbZbz+rF4O72Di7MaZpnNbbO+Tr9ig3v72I0qpapyMBYVjuHo/y+ymL2VlRw8k92vHbU7o7HckYE+buHN6LfmmJFO3ayx/eX4aq8+PvYVfuL3+7lu/W7KBtfDR/vGQALpeNsxtjmlZMpHd7vrjoCGYs2cz0xc7PFg+rcl+5tZQn5xQA8MSo/rRPiHU4kTGmpejWLp4HzvPOJblver7j2/OFTbnXuD38/p0l1NR5uGxoZ87o28HpSMaYFubi7HTOzuxAebWbO6fmOTo9MmzK/fkv1rB8Symd28Zx7zl9nI5jjGmBRIRHL+hHSnw0cwt38Ma89Y5lCYtyL9haxl+/9K6x/ORF/YmPse3yjDHOSGkdwyMX9APgiY9Xsmn3XkdyhHy5ezzK3dPyqK1TLh/amaHdU5yOZIxp4YZndWREVkcqauqYMH2ZIxlCvtzfW1jEwg27SU2I4a4RvZ2OY4wxADyYk0lCTCSfr9zG5yuKm/39Q7rcK6rdPP6xd3bMPef0pk2s7apkjAkO7dvEctuZPQF4aOZyatzNu5dRSJf7y9+uZXt5NQMzkjh/4GHtH2KMMU3m18O6cHRqPOt2VPLO/A3N+t4hW+6lVbW8+LV3oZ7xI3rbomDGmKATFeHijrN7AfDcv9ZQVVvXbO8dsuX+5rz1lFW7GdY9hePtJqoxJkidndmRvke1oaSsmmkLm+/J1ZAsd3edh9e/984fvfGXRzucxhhjDk5E9vfUK9+tbbZ1Z0Ky3L9eXcLW0iq6psRxSo92TscxxpifNSKrI+1ax7BmWzkL1jfPssAhWe4zl2wB4OLsDBtrN8YEvagIF6MGeyd9zMzb0izvGXLlXudRPl+5DfD+a2iMMaFgRNZRAHySv7VZhmZCrtwLtpaxZ28t6cmt6J7a2uk4xhjjl/5piSTFRbF5TxVFu5p+SQK/yl1EhotIgYisEZHxDXw+RkTe8X3+BxHpGuig+ywp2g3AsZ1tdyVjTOhwuYTBvt5avHF3079fYyeISATwPDAC6AuMEZG+B5x2LbBLVY8BngEeD3TQfdZtrwCgR3u7ajfGhJZjfL21r8eakj9X7kOANapaqKo1wGRg5AHnjARe8308FThdmuhO59bSKsC7Ma0xxoSSfb21r8eakj/lngZsrPe6yHeswXNU1Q3sAZrkyaKKau8TXq1jbVlfY0xoae1bjryypumfVG3WG6oicoOI5IpIbklJyWF9jbjoCJLiokiwNduNMSEmITaS5LgoYiKbvnr9achNQEa91+m+Yw2dUyQikUAisOPAL6SqLwAvAGRnZx/WXKDnxgw6nF9mjDGOOyuzI2dlNs8Ubn/++ZgP9BCRbiISDYwGZhxwzgzgKt/HFwH/0uZ6xtYYY8x/afTKXVXdInITMAeIAF5R1XwRmQjkquoM4GXgDRFZA+zE+w+AMcYYh/g1cK2qs4HZBxybUO/jKuDiwEYzxhhzuELuCVVjjDGNs3I3xpgwZOVujDFhyMrdGGPCkJW7McaEIXFqOrqIlADrDzjcDtjuQJxDEewZgz0fBH/GYM8HljEQgj0fNJyxi6qmNvYLHSv3hohIrqpmO53j5wR7xmDPB8GfMdjzgWUMhGDPB0eW0YZljDEmDFm5G2NMGAq2cn/B6QB+CPaMwZ4Pgj9jsOcDyxgIwZ4PjiBjUI25G2OMCYxgu3I3xhgTAEFX7iIyUETmichi38YeQ5zOdCARuVlEVopIvog84XSegxGRcSKiItLO6SwHEpEnfb+HeSLyvogkOZ0JGt8M3mkikiEiX4jIct/3361OZ2qIiESIyCIRmel0loaISJKITPV9D64QkWFOZ6pPRMb6/nyXicjbIhJ7qF8j6ModeAJ4UFUHAhN8r4OGiJyKd8/YAaqaCTzlcKQGiUgGcBawweksB/EpkKWq/YFVwN0O5/F3M3inuYFxqtoXOB743yDMCHArsMLpED/jWeBjVe0NDCCIsopIGnALkK2qWXiXWj/kZdSDsdwVaOP7OBHY7GCWhvwPMElVqwFUdZvDeQ7mGeBOvL+fQUdVP/HttwswD+8OX07zZzN4R6nqFlVd6Pu4DG8pHbinsaNEJB34FfCS01kaIiKJwCl496FAVWtUdbezqf5LJNDKt7NdHIfRg8FY7rcBT4rIRrxXxY5f0R2gJ3CyiPwgIl+JyHFOBzqQiIwENqnqEqez+Oka4COnQ+DfZvBBQ0S6AoOAH5xN8l/+hPfCwuN0kIPoBpQAr/qGjl4SkXinQ+2jqpvwdt8GYAuwR1U/OdSv48gu0yLyGdDQRoL3AqcDY1X1PRG5BO+/rmcEUb5IoC3eH4mPA6aISPfm3lawkYz34B2ScdTPZVTVD3zn3It3qOGfzZkt1IlIa+A94DZVLXU6zz4ici6wTVUXiMgvnc5zEJHAscDNqvqDiDwLjAfuczaWl4gk4/2JsRuwG3hXRK5Q1TcP5es4Uu6qetCyFpHX8Y7XAbyLAz/aNZLvf4BpvjL/UUQ8eNd/KGmufHDwjCLSD+83xRIRAe9wx0IRGaKqW5sx4s/+PgKIyNXAucDpQbLnrj+bwTtORKLwFvs/VXWa03kOcCKQIyLnALFAGxF5U1WvcDhXfUVAkaru+4lnKt5yDxZnAGtVtQRARKYBJwCHVO7BOCyzGfiF7+PTgNUOZmnIdOBUABHpCUQTRIsPqepSVW2vql1VtSveb+Rjm7vYGyMiw/H+6J6jqpVO5/HxZzN4R4n3X+yXgRWq+rTTeQ6kqnerarrve2808K8gK3Z8fxc2ikgv36HTgeUORjrQBuB4EYnz/XmfzmHc8HXkyr0R1wPP+m4kVAE3OJznQK8Ar4jIMqAGuCpIrjpDzV+AGOBT308Y81T1RicDHWwzeCczNeBE4EpgqYgs9h27x7fPsfHfzcA/ff+IFwK/cTjPfr6hoqnAQrxDlos4jCdV7QlVY4wJQ8E4LGOMMeYIWbkbY0wYsnI3xpgwZOVujDFhyMrdGGPCkJW7McaEISt3Y4wJQ1buxhgThv4f8LDbO0r3VpoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "import matplotlib.pylab as plt\n",
    "%matplotlib inline\n",
    "\n",
    "class LiftingLine:\n",
    "    def __init__(self, N, b, c_root, c_tip,\n",
    "                 cl_alpha=2*np.pi, CL_target=1.0, Cb_target=1.0):\n",
    "        \"\"\"\n",
    "        Instantiate the lifting line theory class.\n",
    "        \n",
    "        This creates the data needed to perform the lifting line theory\n",
    "        analysis for a simple, tapered, unswept wing with a fixed lift\n",
    "        and optionally specified wing root bending moment.\n",
    "        \n",
    "        Args:\n",
    "            N (int): The number of Fourier coefficients\n",
    "            b (float): The wing span\n",
    "            c_root (float): The root chord\n",
    "            c_tip (float): The tip chord\n",
    "            cl_alpha (float): The lift curve slope of the airfoil\n",
    "            CL_target (float): The target coefficient of lift\n",
    "            Cb_target (float): The target bending moment coefficient\n",
    "        \"\"\"\n",
    "        \n",
    "        self.N = N\n",
    "        \n",
    "        self.alpha = 0.0\n",
    "        self.b = b\n",
    "        self.S = 0.5*self.b*(c_root + c_tip)\n",
    "        self.cl_alpha = cl_alpha\n",
    "        self.CL_target = CL_target\n",
    "        self.Cb_target = Cb_target\n",
    "        \n",
    "        # Set the theta locations and compute their corresponding y-locations\n",
    "        delta = 0.5*np.pi/self.N\n",
    "        self.theta = np.linspace(delta, np.pi - delta, self.N)\n",
    "        self.y = 0.5*self.b*np.cos(self.theta)\n",
    "\n",
    "        u = np.absolute(2*self.y/self.b)\n",
    "        self.chord = (1.0 - u)*c_root + u*c_tip \n",
    "        \n",
    "        return    \n",
    "\n",
    "    def solve_lifting_line_model(self, x):\n",
    "        \"\"\"\n",
    "        Solve for the lifting line coefficients based on the values of\n",
    "        the design variable twists.\n",
    "        \n",
    "        This code assembles and solves the governing equations for the\n",
    "        lifting line theory. The resulting coefficients can be used\n",
    "        to compute the lift, drag and bending moment coefficients.\n",
    "        \n",
    "        Args:\n",
    "            x (np.array): An array of the twists of dimension self.N\n",
    "        \n",
    "        Returns:\n",
    "            coef (np.array): An array of the Fourier coefficients\n",
    "        \"\"\"\n",
    "        A = np.zeros((N, N))\n",
    "        b = np.zeros(N)\n",
    "        \n",
    "        for m in range(self.N):\n",
    "            for n in range(self.N):\n",
    "                cval = 0.25*(self.chord[m]/self.b)*self.cl_alpha\n",
    "                A[m,n] = np.sin((n+1)*self.theta[m])*np.sin(self.theta[m]) + cval*np.sin((n+1)*self.theta[m])\n",
    "\n",
    "            b[m] = 0.25*self.chord[m]*np.sin(self.theta[m])*self.cl_alpha*(self.alpha + x[m])\n",
    "            \n",
    "        return np.linalg.solve(A, b)\n",
    "\n",
    "    def compute_drag(self, x):\n",
    "        \"\"\"\n",
    "        Compute the drag coefficient based on the design variable values.\n",
    "        \n",
    "        Args:\n",
    "            x (np.array): An array of the twists of dimension self.N\n",
    "            \n",
    "        Returns:\n",
    "            CD (float): The drag coefficient of the wing at the design point\n",
    "        \"\"\"\n",
    "        AR = self.b**2/self.S\n",
    "        coef = self.solve_lifting_line_model(x)\n",
    "    \n",
    "        s = 0.0\n",
    "        for n in range(self.N):\n",
    "            s += (n+1)*coef[n]**2\n",
    "        return np.pi*AR*s\n",
    "\n",
    "    def compute_lift_constraint(self, x):\n",
    "        \"\"\"\n",
    "        Compute the lift constraint based on the design variable values.\n",
    "        \n",
    "        The lift constraint enforces that the lift must be equal to the\n",
    "        target lift coefficient\n",
    "        \n",
    "        CL = CL_target\n",
    "        \n",
    "        Args:\n",
    "            x (np.array): An array of the twists of dimension self.N\n",
    "            \n",
    "        Returns:\n",
    "            con (float): The value of the lift constraint\n",
    "        \"\"\"\n",
    "        AR = self.b**2/self.S\n",
    "        coef = self.solve_lifting_line_model(x)\n",
    "        return np.pi*AR*coef[0] - self.CL_target\n",
    "\n",
    "    def compute_moment_constraint(self, x):\n",
    "        \"\"\"\n",
    "        Compute the wing root bending moment constraint based on the\n",
    "        design variable values.\n",
    "        \n",
    "        This constraint enforces that the bending moment must be less\n",
    "        than the target bending moment coefficient:\n",
    "        \n",
    "        Cb_target - Cb >= 0.0\n",
    "        \n",
    "        Args:\n",
    "            x (np.array): An array of the twists of dimension self.N\n",
    "            \n",
    "        Returns:\n",
    "            con (float): The bending moment constraint value\n",
    "        \"\"\"\n",
    "        AR = self.b**2/self.S\n",
    "        coef = self.solve_lifting_line_model(x)\n",
    "        \n",
    "        s = 0.0\n",
    "        for n in range(self.N):\n",
    "            if n == 1:\n",
    "                cval = 0.0\n",
    "            else:\n",
    "                cval = (2*np.sin(0.5*np.pi*(n+1)) - np.sin(np.pi*(n+1)))/(4.0 - (n+1)**2)\n",
    "            s += cval*coef[n]\n",
    "            \n",
    "        return self.Cb_target - AR*s # >= 0.0\n",
    "    \n",
    "    def plot_gamma_dist(self, x):\n",
    "        \"\"\"\n",
    "        Plot the circulation distribution as a function of the span.\n",
    "        \n",
    "        Args:\n",
    "            x (np.array): An array of the twists of dimension self.N\n",
    "        \"\"\"\n",
    "        coef = self.solve_lifting_line_model(x)\n",
    "        \n",
    "        size = 100\n",
    "        theta = np.linspace(0, np.pi, size)\n",
    "        y = 0.5*self.b*np.cos(theta)\n",
    "        Gamma = np.zeros(size)\n",
    "        \n",
    "        for n in range(self.N):\n",
    "            Gamma[:] += 2*self.b*coef[n]*np.sin((n+1)*theta)\n",
    "            \n",
    "        plt.figure()\n",
    "        plt.plot(y, Gamma, linewidth=2)\n",
    "        \n",
    "        return\n",
    "\n",
    "# Set the number of coefficients to use as well as the geometric\n",
    "# properties of the wing.\n",
    "N = 5\n",
    "b = 15.0\n",
    "c_root = 1.0\n",
    "c_tip = 0.65\n",
    "\n",
    "# Create the lifting line class.\n",
    "lifting_line = LiftingLine(N, b, c_root, c_tip, Cb_target=0.175)\n",
    "\n",
    "# Run the initial design problem: Minimize induced drag subject to a lift\n",
    "# constraint.\n",
    "x = np.ones(N)\n",
    "res = minimize(lifting_line.compute_drag, x, method='SLSQP',\n",
    "               constraints={'type': 'eq',\n",
    "                            'fun': lifting_line.compute_lift_constraint})\n",
    "lifting_line.plot_gamma_dist(res.x)\n",
    "\n",
    "# Run the modified problem: Minimize induced drag subject to a lift\n",
    "# constraint and a requirement that the root bending moment be bounded\n",
    "# by the target bending moment coefficient.\n",
    "x = np.ones(N)\n",
    "res = minimize(lifting_line.compute_drag, x, method='SLSQP',\n",
    "               constraints=[{'type': 'eq', \n",
    "                             'fun': lifting_line.compute_lift_constraint},\n",
    "                            {'type': 'ineq',\n",
    "                             'fun': lifting_line.compute_moment_constraint}])\n",
    "lifting_line.plot_gamma_dist(res.x)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}