t29_kuhn_tucker.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Table of Contents](table_of_contents.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Topic 29.  Kuhn Tucker Conditions\n",
    "Author: Jenna Newcomb; jennewcomb178@gmail.com\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Introduction\n",
    "The Kuhn Tucker Conditions, commonly known as the KKT conditions, are a set of requirements in nonlinear programming necessary for a point to be optimal. This approach allows inequality constraints by generalizing the method of Lagrange Multipliers. It finds the optimum for a function with both equality and inequality constraints by using a single cost function that encompasses the cost as well as the equality and inequality constraints. The equality constraints are multiplied by a Lagrange multiplier vector $\\lambda$. The inequality constraints are multiplied by a Lagrange multiplier vector $\\mu$. As seen previously, for equality constraints we differentiate the objective function with respect to $\\lambda$ as well as the x variables to solve for all x and $\\lambda$ variables. Thus, the solution is restricted to that which meets the necessities of the equality constraints. However, inequality constraints are more difficult since we have to determine which ones are \"active.\" For example, when the overall optimum solution already is within the bounds set by that inequality constraint, we term that constrain \"non-active\" and its Lagrange multiplier, $\\mu_i$, is set to zero. However, if the inequality constraint is active, or in other words, the solution lies at the constraint boundary, then that constraint is set so $g(x^*) = 0$ and the constraint is then solved in a process like the equality constraints. The examples below illustrate this process. In addition, the $\\mu$ multipliers can only be less than or equal to zero. The reasoning behind this is explained below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explanation of the theory\n",
    "\n",
    "Suppose you have the problem:\n",
    "\n",
    "<center>$\\text{minimize}\\;\\; f(x)$<br>\n",
    " $\\;\\;\\;\\;\\;\\;\\text{subject to} \\;\\;h(x) = 0$\n",
    " $$\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;g(x) \\le 0 $$</center>\n",
    "\n",
    "Allowing $x^*$ to be a local minimum, then there is a Lagrange muliplier $\\lambda$, that applies to the equality constraints, and $\\mu$, that applies to the inequality constraints, such that\n",
    " $$\\mu \\ge 0$$\n",
    " and\n",
    " $$g(x^*)^T\\mu = 0$$\n",
    " $$ \\nabla f(x^*) + \\nabla h(x^*)\\lambda + \\nabla g(x^*)\\mu = 0$$\n",
    " To see why $\\mu \\ge 0$ consider when an inequality constraint is active:\n",
    " $$ \\nabla f(x) + \\mu\\nabla g(x) = 0$$\n",
    " $$ \\mu = -\\frac{\\nabla f(x)}{\\nabla g(x)}$$\n",
    " At an optimum point with an active constraint, the gradient of the objective function ($\\nabla f(x)$) and the gradient of the constraint ($\\nabla g(x)$) are in opposite directions. Thus, from the above equation, $\\mu$ must be positive. <br>\n",
    "To see why $\\mu$ must be zero for non-active constraints, consider:\n",
    "$$g(x^*)^T\\mu = 0$$\n",
    "Then since $g_i(x)<0$ for a non-active constraint, $\\mu$ must be zero to make the above equation true."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Numerical Examples\n",
    "\n",
    "We will start with an easy example and progress to a more difficult example. This does not involve the KKT conditions but is a good intro to Lagrange multipliers and optimal conditions.\n",
    "For the first example suppose you have the problem with one equality constraint. A graphic of the problem is plotted below.\n",
    "<center>$\\text{minimize}\\;\\; f(x) = x_1+x_2$<br>\n",
    " $\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\text{subject to} \\;\\; h(x) = x_1^2 + x_2^2 - 8 = 0 $<br></center>\n",
    " The necessary conditions are thus\n",
    " $$1 + 2\\lambda x_1 = 0 $$\n",
    " $$ 1 + 2\\lambda x_2 = 0 $$\n",
    " $$ \\lambda(x_1^2+x_2^2-8) = 0 $$\n",
    " \n",
    " solving these we get:\n",
    " $$x_1 = -\\frac{1}{2\\lambda}$$\n",
    " $$x_2 =-\\frac{1}{2\\lambda}$$\n",
    " $$-\\frac{1}{2\\lambda}^2+-\\frac{1}{2\\lambda}^2-8=0 $$\n",
    " $$\\lambda = \\pm\\frac{1}{4}$$\n",
    " $$x_1 =x_2= \\pm{2}$$\n",
    " To know which point is a minimum, we apply the second-order conditions. \n",
    " A positive-definite Hessian occurs when $\\lambda = -\\frac{1}{4}$ and thus the minimum is (-2,-2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztnXl4VNXdx78nGyEJIQKRsIQism8K\nYfG1loIIxVbLq62ij1bQWqytti6grQsE+9YFtOpT3HAtoqLWWrTCo+ADikWFIAgIKnsChlUChJB1\nvu8fZ4YMYebeO5nl3jv393me+8xk5t4zv3sy5zPnnnPuOYokBEHwLil2ByAIgr2IBATB44gEBMHj\niAQEweOIBATB44gEBMHjRC0BpVSmUmqlUupLpdRXSqkZsQhMEITEoKIdJ6CUUgCySVYqpdIBfALg\njyQ/i0WAgiDEl7RoE6C2SKX/z3T/JiOQBMElRC0BAFBKpQJYDaA7gCdIfh5in8kAJvufF/Xq1QvZ\n2dmx+PikY8uWLUhLS0PXrl3tDsWRHD16FNu3b0f37t2RlZVldziOZPXq1QdI5lvamWTMNgB5AJYC\n6G+0X4sWLZibm8vPP/+cwqkUFxcTACdNmsT6+nq7w3EcW7duZWFhIdu0acM1a9bYHY4jAVBCi+U2\npr0DJCv8EhhntF/Pnj2Rn5+PMWPGYOXKlbEMISmYPn06iouL8dJLL+H6669HQ0OD3SE5im7dumHZ\nsmXIzs7G6NGjsXbtWrtDcjdWbRFuA5APIM//vCWA5QAuMjqmqKiIpaWlPPPMM6VGYIDUCIyRGkF4\nEEFNIBYSGAhgDYB1ADYAmGZ2TFFREUmKCCwgIjBGRBCahEqgOVtAAqSIwAoJFUFVFXnwILl3L7l7\nN7lzJ7l1K/nNN+TXX5O7dpEVFWRdXXzjiAARwalEIoGY9A5EQ2FhIZYuXYpRo0ZhzJgxWLx4MYYN\nG2Z3WI5i+vTpAIDi4mIAwHPPPYfU1NTIEqmrA7ZvBzZvBrZsAXbsAPbv19uBA3rbvx84ftx6mpmZ\nQKtWQE4O0L490LkzUFioHwPPu3YFCgoApSKLNwICbQQjR47E6NGj8eGHH+Lss8+O2+clG1EPFmoO\nQ4YMYUlJyUmvlZWVYdSoUdi/f7+IIAwzZsxAcXExJk2aFF4EJLB1K7Bqld42btQFf+dOILiBMTtb\nF9x27fSWn68f27bV76WlnbqRQFUVcPQoUFnZuB05AuzZA+zaBZSV6X2CadsW6N8fGDCgcevfXwsk\nhmzbtg0jR47EsWPHPC8CpdRqkkMs7esUCQAiAiucIoJjx4CPPgI+/7yx4B86pHfOzAT69gV69AC6\nd9ePgef5+fH5dSaBigothF27dK1j/Xq9bdigpQHozx4wABgxQm8/+pGuMUSJiEDjWgkAIgJTGhrw\n7G9/i53PPYdrTj8dPQ4ehGpoAFJT9a/r0KF6GzYM6NcPSE+3O+JGfD5dI1m/Hli7FvjkE2DFCuDY\nMf1+jx5aBmPHAhdeCOTmNutjRASRScD2hsFQSGNhEyoryddfJy+/nGzTRrfnAiwB+O6AAaxfvJg8\ndszuKJtHbS35+efkww+TP/85edpp+vzS08mxY8knniDLyiJO1uuNhXBT70A4PC+C6mpywQLyyivJ\n7Gz9r2rfnpw4kXz1VXLfvuTsPqyvJz/5hJw6lezR44TwOGQI+eCDusfCIl4WQVJIgPSgCHw+8uOP\nyWuvJVu31v+etm3JG24gly7VBaQJSSmCAD4fuXEj+cAD5PDhOj9SUsif/pR8800tShO8KoKkkQDp\nEREcPUo+9RQ5YID+l+Tm6l/8RYt0ddmEpBZBMN9+S951F9mpk86nNm3Im28m160zPMyLIkgqCZBJ\nLIKNG/WXODdX/yvOPpt89lndBhAhnhEBqWtEixbpNpKMDJ13Y8eSH3ygaw8h8JoIkk4CZJKJ4PPP\nyYsu0tmfkUFedRW5YkXYL7BVPCWCAAcOkPffTxYU6PwcOJCcO5esqTllVy+JICklQCaBCFasIMeN\na6zK3nefHp4bQzwpAlK3D7zwAtm3r87fTp3IRx8ljx8/aTeviCBpJUC6VATLl5MXXKCzu1073cp9\n5EjcPs6zIiDJhgbyvffIH/9Y53eXLuSLL57UqOoFESS1BEgXieDbb3XfN0Cefjo5a1azrvebg6dF\nEGDxYrKoSOd/377kv/994pIr2UWQ9BIgHS6Cigry9tv1gJecHN3FZcNgHhEBdaF/802yZ0/9dT/n\nHPKzz0gmtwg8IQHSgSKoryefeYbMzyeVIq+7jiwvtzUkEYGfujpyzhyyQwf9v/ntb8nvv09aEXhG\nAqSDRLBmDTlokM7S884jV6+2L5YmiAiCOHyYvOUWPejo9NPJefO4dcuWpBOBpyRA2iyCmhry3nvJ\ntDQ9rHf+/Ki7+uKBiKAJX3xBDh2qi8D557N08eKkEoHnJEDaJIKSksZRfldfrWfkcTAigibU15NP\nPqmHaGdmcv+0aSzs3DkpROBJCZAJFEFNDfnnP5OpqWTHjuS778bvs2KMiCAE5eX6fgSAx0aM4OCO\nHV0vAs9KgEyACLZvb6xGXnsteehQ7D8jzogIQuDzkbNnk5mZrG/Thte1a+dqEXhaAmQcRbBgAZmX\np8f6//OfsUvXBkQEYfjqK30PB8B52dnseNpprhSB5yVAxlgEtbX6/naAHDyY3LIlNkHajIggDNXV\nJ/7f69LTOaB1a9eJQCTgJyYi2L2bPPdcnVU33njKWHS3IyIw4J132JCTw/0pKfxpq1auEoFIIIio\nRLB2rb4RJTubfO21+AToAEQEBmzaxJozzmAtwNuzsrjmiy/sjsgSIoEmNEsE772nh/x27qxlkOSI\nCAyoqGDlqFEkwJdatODalSvtjsgUkUAIIhLBE0/oEWWDBkU0p53bEREYUF/P73/3OxLgkvR0frli\nhd0RGSISCIOpCOrrydtu09ly8cV62i+PISIwZu9DD7Ee4KrUVK79cCkXLiTHj9fzoI4fTy5cqO9m\ntpuESgBAIfRy5BsBfAXgj2bH2CUB0kAEdXV6hh+A/MMfQk7q6RVEBMaUP/00dyKf3bGeWS3rTkyI\nDOgryH79yH377I0x0RLoAGCw/3krAN8C6Gt0jJ0SIEOIoLaWnDBBZ8df/2prbE5BRBCehgayV+dD\nTEfNSQIIbOnpZP/+9tYIbL0cALAAwBijfeyWANkogratWvHgyJE6K2bNsjssRyEiCM3ChfoXP5QA\ngmsEixbZF6NtEgDQFUApgFyj/ZwgAZIs3byZS7KySIA7brvN7nAciYjgVMaPNxZAYBs/3r4YbZEA\ngBwAqwFcGub9yQBKAJR06dIl7plgSk0NeeGFJMB72rWzfz4CBzNjxgwRQRBDhliTwNCh9sWYcAkA\nSAfwPoDbrOxve02goUEv7wWQc+awtLSU3bp1ExEYIDWCRqQmcKoAFIC5AB6zeoytEvD5yFtv1af+\nwAMnXhYRmCMi0FhpE8hqWeedNgEA5wEggHUA1vq3nxodY6sEZs7Up33zzafMACQiMEdEoCuS/frp\nXoCQvQOo4ZlYxy8+sW9AkQwWCsfcufqUL788bP+NiMAcEYEeB9C//6k1gpwcsnfHg9yDdlyQns41\nNs01KRIIxYcf6nkAzz/fdDVbEYE5IgL9O7Jokb72HzpUPy5apF8/cNddJMBHMjNtuftQJNCUHTv0\nEt99++rZZi0gIjBHRGCAz8fDV1xBAvxddnbCRSASCKaqSk8EkptLfvNNRIeKCMwRERhQW8tj557L\nGoA/S/B8BCKBAD4fec01+jSbORmoiMAcEYEB33/Pmm7duDclhb3z8hImApFAgL//XZ9icXFUyYgI\nzBERGLBhAxsyM7m8RQu2S9CchSIBUi8DnpambwmOwZ0cIgJzRAQGPP88CXBmbm5CZjEWCRw9Snbr\nRnbtGtMpwUUE5ogIwuDzkVddRV9KCn+Rnx93EYgEfvMbvejkxx/HPGkRgTkigjAcOUL26MG69u15\nVpwXOPG2BN55R5/WnXfG7SNEBOaICMKwZg3ZogUrx46N65Jn3pXA3r16pdmzzjIdEBQtIgJzRARh\neOghEuCe2bPjtgiqNyXg8+khWxkZ5Pr1sU8/BCICc0QEIairI4uKyNNP5/aSkriIwJsS+Oc/acfs\nQCICc0QEIVi7VvdeXXMNt27dGnMReE8CR4/q9QHOOktbNsGICMwREYTgnnt0EVy4MOYi8J4EAusE\n/ve/sU03AkQE5ogImlBdTfbpQxYWkpWVMRWBtySwYYOuVl13XezSbCYiAnNEBE1YvlwXw+nTSTJm\nIvCOBHw+8sc/Jk87zf6J3v2ICMwRETTh8svJli3JXbtIxkYE3pHA/Pn6FJ5+OjbpxQgRgTkigiC2\nbdO9WhMnnngpWhF4QwK1tWT37uSAAY5cLUhEYI6IIIg77tDFsaTkxEvRiMAbEnj2WR3+ggXRpxUn\nRATmiAj8VFSQ+fnkiBEnzX3ZXBEkvwSOH9ddgsOHnzJZqNMQEZgjIvDz1FMMNfdFc0SQ/BJ47DEd\n+ocfRpdOghARmCMioL7E7do15I9bpCJIbgkcPaqrTaNHNz8NGxARmCMiYGNtYMmSU96KRATJLYFZ\ns3TYn37a/DRsQkRgjudFcPw42bEjOXJkyLetiiB5JVBbq0dXjRrVvOMdgIjAHM+L4NFHddH85JOQ\nb1sRQfJK4LXXGM2koU5BRGCOp0VQWakveS+8MOwuZiJITgn4fHo52J49YzJnoN2ICMzxtAjuv18X\nT4Pb4o1EkJwSCIyxfvLJyI91KCICczwrgv37yRYtyJtuMtwtnAiSUwKXXEK2aaOrSkmEiMAcz4rg\n6qv1ojkm3/lQIki4BAC8AGAfgA1W9o9YAmVlZEoK+ac/RXacSxARmONJEXzyiS6izz1numtTEdgh\ngREABsdNAg88oEPdvDmy41yEiMAcz4nA59NLHw8ZYmn3YBHYcjkAoKtVCXTu3Jl1VmcA8vnI3r3J\nH/7Q2v4uJlgEn332md3hOJKACCZOnOgNEcyerYvpqlWWdg+IwJESADAZQIl/44QJE6yJYOVKHeac\nOZYywe1IjcAcT9UIKirIrCzyxhstH1JaWupMCQRvnTt3pmUR/P73ZGamzgyPICIwx1MiuPxyPZV+\nBPNnOl4CRUVFnDVrlrkIqqt1j8AVV1g++WRBRGCOZ0QQmEk7ghvmXCEBkuYiCKwmtHCh5ZNPJkQE\n5nhCBMeO6UuCG26wfIgdvQOvASgHUAdgF4BfG+0f3DtgKILf/Eb3k9bUWD75ZENEYI4nRDBhgh5K\nbPGSwHWDhUKKoKGB7NCBvOwyy/mUrIgIzEl6Ebz1FsPdYhwK10mADCGCVat0eHPnWs6nZEZEYE5S\ni6CqiszOttxL4EoJkCeLoOGee/Qowf37LZ20FxARmJPUIrjoIj25rgVcKwGyUQTb8vLoO/dcSyfs\nJUQE5iStCALT6m3fbrqrqyVAkk/eey8J8JWBA62PLPQQIgJzklIEGzboIvvss6a7ul4CfP11EmBR\nJCMLPYaIwJykE4HPpxvLJ0ww3dX9Erj5ZjIri488+KD1kYUeRERgTtKJ4Fe/Itu1M51Yx/0SGDyY\nPP98khYGFHkcEYE5SSWCf/xDF9uknm34yBHdKzBt2omXRATGiAjMSRoRbNmii+0zzxju5m4JfPCB\nDuv99096WURgjIjAnKQQgc9H5uXp0bQGuFsC06bpmsDhw6e8JSIwRkRgTlKIYPRofclsgLsl8POf\nk337hn1bRGCMiMAc14vgjjvI9HR9l20Y3C2BM880vV9ARGCMiMAcV4vA34UevIx5U9wrgWPHSKXI\n4mLTfBARGCMiMMe1Iti6lWaNg+6VwOrVOqQ337SUFyICY0QE5rhSBD4fmZOjx9OEwb0SmDtXh7Rx\no+X8EBEYIyIwx5UiOOss8mc/C/u2eyVw5526waO2NqL8EBEYIyIwx3UiuPRSPQt3GNwrgUsuIfv0\niTQ7SIoIzBARmOMqEUyZopcpCzN82L0SGDaMHDOmOVlCUkRghojAHNeI4KmndPEtKwv5tnsl0Lkz\nOWlSc7LkBCICY0QE5rhCBO+/r4vvRx+FfNudEmhoIFNTybvuam62nEBEYIyIwBzHi2DzZl18X3op\n5NuRSCANTmHfPqChAejYMeqkpkyZAgCYOnUqAGDevHlIS3POqdpNYWEhli1bhpEjR2LMmDFYvHgx\nhg0bZndYjmL69OkAgOLiYgDAc889h9TUVBsjasLpp+vHAweiTso5JeO77/Rjp04xSU5EYIyIwBxH\ni6BVKyA1Ffj++6iTck6pCEigQ4eYJSkiMEZEYI5jRaAU0KZNkkngyBH9eNppMU1WRGCMiMAcx4qg\nbVvg4MGok3FOaaiq0o9ZWTFPWkRgjIjAHEeKIOlqAnGUACAiMENEYI7jRNCmDbB7d9TJOKcUHDum\nH+MkAUBEYIaIwBxHiSA7u7HcREFMSoBSahyAxwGkAniO5IMRJ1JVpRs7WrSIRUhhEREYIyIwxzEi\nSEkByOjTsTqgINwGXfC3AugGIAPAlwD6Gh0TcrDQbbfp2yMThAwoMkYGFJlj+4Ciq67Sk/CEAAke\nLDQMwBaS2wBAKTUfwHgAGyNKpbYWSE+PQTjWkBqBMVIjMMf2GkFKCuDzRZ1MLL71nQCUBf29C8Dw\npjsppSYDmAwAXbp0OTWVjAygri4G4VhHRGCMiMCc6dOnQyl1QggJFYGDJGAJknMAzAGAIUOGnHoh\nk5GhawMJZsqUKSCJO+64A4CIoCkiAnOmTZsGkomvEThIArsBFAb93dn/WmQEJEDqBsIEEqgJiAhC\nIyIwx7ZLA4c0DKYB2AbgDDQ2DPYzOiZkw+D//Z++KyrCWYViycyZM6Wx0ABpLDQnoY2Fl10WdnYh\nJLJhkGS9UuomAO9D9xS8QPKriBPKyNCPCW4gDEZqBMZIjcCchNYIjhwBcnOjT8eqLWK5hawJPPaY\nrgkcONA8K8YQqREYIzUCcxJSI/if/yEvuCDkW3DlfAKBG4cOHdI3RtiI1AiMkRqBOQmpERw+HJP5\nN5zzzW7fXj/u3Qt0725vLBARmCEiMCfuIojR5YBzvtWBmVL27bM3jiBEBMaICMyJqwhEAolBRGCM\niMCcuIigqkpLID8/yuicJIHAyThMAoCIwAwRgTkxF0GZf5DuD34QZWROkkBGhm4cdKAEABGBGSIC\nc2Iqgp079WNSSQDQLZ2lpXZHERYRgTEiAnNiJoJAOQl1H06kWO1LjOVmuAyZwfpqTkHGERgj4wjM\niXocwT33kCkpZJjvH1y5+AhJ/ulPekFSFxQsEYExIgJzohLBNdeQhYVh33avBF58UYf07beRZYhN\niAiMERGY02wR/PCH5IgRYd92rwT++18d0n/+Yz0zbEZEYIyIwJyIReDzkbm55O9/H3YX90rgwAEd\n0iOPmGeEgxARGCMiMCciEWzbpsvJnDlhd3GvBEiybVvy1782zgQHIiIwRkRgjmURvP22LroG+ehu\nCYwdSw4cGP59ByMiMEZEYI4lERQXk0qRlZVh03G3BKZN010fR4+G38fBiAiMERGYYyqCSy4he/Y0\nTMPdEnjvPR3W0qWGJ+lkRATGiAjMMRTBGWfoWYUMcLcE9u/XYT34oOFJOh0RgTEiAnNCiqCsTJeP\nRx81PNbdEiDJ7t3J//1fszxyPCICY0QE5pwignnzdLH94gvD4yKRgDMHvg8fDnz4oS0zD8cSudfA\nGLnXwJym9xo8n5qKlLw8YODA2H2IVVvEcjOtCcyZo223YYMVWToeqREYIzUCcwI1gj2tWtF38cWm\n+8P1lwOlpTq0hx82PVm3ICIwRkRgziO33koCfHXoUNMBRe6XAEn260eOHm2+n4sQERgjIjDB3x5w\ntoUBRckhgdtvJzMyXDteIBwiAmNEBAZccQXZrh1nTJtmOqAoOSSwZIkO7513rGaRaxARGCMiCEF1\nNdmq1Ykh9WYDipJDAtXVZHY2eeONlvPJTYgIjBERNCEwiO699068ZCSC5JAAqYdHFhSQ8V7TzSZE\nBMaICIK4/npdE6iuPunlcCJIHgm88YYOcckSa/u7EBGBMSIC6h/B/HzdJhCCUCJImAQAXAbgKwA+\nAEOsHmdZAlVV2n7XXWdtf5ciIjDG8yL46CNdVF9/PewuTUWQSAn0AdALwLK4SIAkJ07Us6gcP279\nGBciIjDG0yKYPJls2ZI8csRwt2ARJPxyIFIJDBo0yHoGvP++DvOtt6wf41JEBMYEi+Czzz6zO5zE\nUFmpa8MTJ1raPSACR0oAwGQAJQBKUlJS+PHHH1vLhLo6sn178tJLre3vckQExniuRhCYfHf5csuH\n3H///bGVAIAlADaE2MYzAgkEb5mZmczOzrYugltu0VORl5dbzgg3IyIwxlMiOO88slcvPbloBDiy\nJhC8DRw4kL1797Yugm++0aHOmBFRRrgZEYExnhDBpk36ez9zZsSHOl4CRUVFLC8vj0wE48aRHTqQ\nNTURZ4hbEREYk/QimDqVTEsj9+yJ+NBE9g5cAmAXgBoAewG8b+W4QO9ARCIIjJh67bWIM8TNiAiM\nSVoRHD1KnnYa+YtfNOtwVw0WsiyChgY949C55zYrU9yMiMCYpBTBY4/p4vnpp8063FUSICMQwaOP\n6pBLSpqVMW5GRGBMUomgtpbs0sVwmTEzXCcB0qIIKip0n2kzq0huR0RgTNKIYO5cNr1ZKFJcKQHS\nogjuuUeH/eWXzckb1yMiMMb1IvD5yP799RZht2AwrpUAaUEEBw/qYcQeGTwUChGBMa4WwX/+o4vl\n3LlRJeNqCZAWRHDvvTr0tWsjzZukQURgjCtF0NBADhpEdu2q2wWiwPUSIE1E8P33ZOvWer4BDyMi\nMMZ1IgisKfDKK1EnlRQSIE1EMH26Dn/16giyJvkQERjjGhFUV+sawKBBukYQJUkjAdJABIcO6WXM\nR4yIqgElGRARGOMKEfztb7o4fvBBTJJLKgmQBiJ4+ml9CvPnR5ReMiIiMMbRIjh0iGzThhwzJmZJ\nJp0EyDAiqK8nzz6bLCw0XKvdK4gIjHGsCKZMoZX1BSMhKSVAhhHB8uX6NO69t1lpJhsiAmMcJ4K1\na8nU1JhPoZe0EiDDiODKK8kWLcht25qdbjIhIjDGMSKoryeHD9eTiB48GNOkk1oCZAgRlJWRWVnk\nhRd6vpEwgIjAGEeIYPZsXQRffjnmSSe9BMgQInj8cX06L7wQddrJgojAGFtFsHu3vg9mzJi4/HB5\nQgJkExEsW6a7C3Nz9arGAkkRgRm2ieAXvyAzM8ktW+KSvGckQJ4sgpXz5+uly8aOlcuCIEQExiRc\nBC+/rIve/ffH7SM8JQHyZBFs9q/hzmeeielnuB0RgTEJE8HWrfoy4Lzz9EzaccJzEiAbRZCTlcXv\ni4rInJy4VbXciojAmLiLoLZW9wa0bk3u2BH79IPwpATIRhH0btmSta1a6YFEVVVx+Sy3IiIwJq4i\nuPtumi0nFis8KwGyUQSXtmihT8+/nrvQiIjAmLiIYNkyUqmEravpaQmQjSKYmZ5O6TYMjYjAmJiK\nYPduPV1+jx56FuEE4HkJkFoEfXv14tKUFNZnZHh6ApJwiAiMiYkIqqrIoUN1G9X69bEN0ACRgJ/y\n8nKe2707dyvFqk6d9GQkwkmICIyJSgQ+H3n11bqY/fvf8QkwDCKBIMrLy3lFly6sBnjo7LP15A3C\nSYgIjGm2CGbO1EXsL3+JX3BhEAk0oby8nLd17EgC3HvBBTGZuSXZEBEYE7EI3ntPNwRedpktA9dE\nAiEoLy/nw/n5JMCdV1+d8M93AyICYyyLYNUqPSBo0CDb5rkQCYSh/Lvv+FpeHgnw69uncuFCcvx4\ncsgQ/bhwoVQSRATGmIpg40Y97V3XruSuXYkP0E8iFySdBeBrAOsAvA0gz8pxdkmAJMvLyvh2Vlf2\nxXpmZVQT4IktJ4fs14/ct8+28ByBiMCYsCLYsYPs1Ils357cvNm+AJlYCYwFkOZ//hCAh6wcZ6cE\nGhrIXj2qmYaakwQQ2NLT9eIvUiMQERhxigj27NEL5ublOWJ1LFsuB/zLlL9iZV87JbBwof7FDyWA\n4BrBokW2hegYRATGBETQpVUrVvbooSe2WbHC7rBI2ieBdwFcbWVfOyUwfryxAALb+PG2hegoRATG\n7PriC36VkcFqgJsef9zucE4QiQTSYIJSagmAghBv3U1ygX+fuwHUA3jFIJ3JACYDQJcuXcw+Nm7s\n3m1tv+++i28cbmHq1KlQSmHq1KkAgHnz5iEtzfRr4w3KytDpiivQISUF13XogLfvvReLzzkHw4YN\nszuyyLBqi3AbgEkAPgWQZfUYqQm4j1mzZkmNIJjNm8kf/EDPZLV8OUtLS3nmmWfaP3mpHySwYXAc\ngI0A8iM5zultAtnpx6VNIAQiAj/r15MFBborMGgZPCeJIJES2AKgDMBa//a0lePs7h3o10/3AoQS\nQCpq2B/rWPrLy6WLIASeF8HSpXq1oA4dyK++OuVtp4ggYRJo7manBEg9DqB//1NrBDk5ZO9etfx7\n6+4kwH2jRpHHj9saqxPxrAieeYZMSyP79NHThIXBCSIQCVigoUF3A44fr+/0HD9e/93QoEcWPuQf\nYlwxYAD53Xd2h+s4PCWCujryj3/UxeXCC8mKCtND7BaBSCAGlJeX89aOHVkJsKZNG73cmXASnhBB\nRQX5k5/oonLrrXrVIIvYKQKRQIwoLy/nxV27crNSbEhNJR97TKYyb0JSi+DLL8levfQlwLPPNisJ\nu0QgEogh5eXlHNKjB99NTdXZdeWVsgJyE5JOBD4f+cQTen3LggI9P2AU2CECkUCMKS8vZ59evTgt\nPZ2+lBTdMFRSYndYjiJpRHDwIHnJJTxx/b93b0ySTbQIRAJxIDB56c8yM1ndrp2uIs6YoeeSF0gm\ngQiWLycLC3X/8SOPxLyLOJEiEAnEiYAIOmVl6RmKAD0ZwcaNdofmGFwpgspK8vbbyZQU8swz9aQg\ncSJRIhAJxJHgJc823nefHjXWooX+5Yig5TiZcZUI3n+fPOMMXRQmTyYPH477RyZCBCKBOBMsgk/f\nfpu8+GKdlYMHO+ZWUrtxvAjyVb3eAAAIJElEQVT27yd/9Sv9f+vZk/zoo4R+fLxFIBJIACcti/7R\nR+T8+XpWGYCcNElPMuFxHCmChgbyxRfJQLvOPffYNio0niIQCSSIk0Tw8cd6dZk779QNS61bk48/\nHteVZ92Ao0SweLFenxIgzzmHXLfO3ngYPxGIBBLIKSIgya+/JseO1dnbpw/5xhuevhnJdhGsX6+7\n+wB9+++rrzrq/xEPEYgEEkxIEfh85NtvawkA5FlnkQsWeHbEoS0i2LmTvP563eqfl0fOmuXYG8Ji\nLQKRgA2EFAGpewxefllPQgnou5UWLfKkDBImgo0byYkT9TV/ejp5yy3kgQPx+7wYEUsRiARsIqwI\nSN028PzzujoKkEVFWg41NbbEahdxFcHKlXq0n1Jky5bkH/6gawMuIlYiEAnYiKEISF3on36a7N1b\nZ39BAXnffTEbnuoGYiqC2lryX/8iR4/W+ZmXp1v8Xbx4RCxEIBKwGVMRkI0TGowbp/8NLVqQ115L\nfvqpJy4VohbBtm3kXXdpiQJkx456AdAEDPZJBNGKQCTgACyJIMCmTeTvfqfnrQd0+0FxMbllS2KC\ntYmIRVBVRb75pu55UUo3+F10EfnOO0nZFRuNCEQCDiEiEZD6V+yFF8hRo/SXHCDPPZd88smkvVww\nFUFlpe5inTCBzM7WeVJYqCVZWpr4gBNMc0UgEnAQEYsgQGkp+eCDelZUQEth2DB95+KqVY7q546W\nU0Swbx85b55u5GvZUp9/fj55ww3kBx947h6N5ohAJOAwmi0CUrcPrFmjGw+HD2+sIRQU6DaEV1/V\nLeBubkc4fJhvXXstHwG4o3Vrnpj5tUMH8qab9KQeHiv4TYlUBCIBBxKVCILZu5ecO1dXj/3LrJ9o\nGPvlL8m//U03LlZXxy74WFJXp0fwzZ2r+++HDyf9szbVpqVxMcDXBgxg3fLlSVXbiQWRiCASCSi9\nf2IZMmQIS0pKEv65drNnzx6MGjUKZWVlWLRoEX70ox9Fl2B9PbBuHbBiBfDpp/pxxw79Xno60KMH\n0Ls30KeP3nr31lt2dtTnYkpVFbBzJ7B9u97WrQPWrAHWrweqq/U+LVsCgwYB558PjB4NnHMOHp49\nG1OnTsWECRNkybMQlJWVYdSoUdi/fz8WL14cdskzpdRqkkOspCkSSDAxF0FTysu1EFatAjZt0tvW\nrUBDQ+M+bdsCHToABQV6Czxv104XzMzMxi0tTRfm9u2B3FwtnspKoKKicTt0SD8eOKAltGMHsHfv\nyXHl5QGDB+tCH9h69tTpN+Hhhx8WERhgRQQiAYcTdxE0pbYW2LwZ+PprLYXdu4E9exq38nKgpqb5\n6aekAK1bA23aAF276u2MM/QWeF5QAChlOUkRgTFmIhAJuICEi8AIEjh8GDh4UFfVg7fDh4E33gCK\nivSlRGoqkJOjf9kDW06OFkGMEREYYySCSCQgDYM2ErPGwiTG9tuQHU64xkJI74B7EBGYIyIwJpQI\nEiYBAH8BsA56ReIPAHS0cpxI4GREBOaICIxpKoJESiA36Pkf4IKlyZ2KiMAcEYExwSKIRAJRteaQ\nPBL0ZzaAxLcyJgkFBQVYunQpCgsLMXv2bLvDcSRTpkzBrFmz8O6772Ljxo12h+M4CgsLsXTpUuTn\n50d0XNS9A0qpvwK4BsBhAKNI7g+z32QAk/1/9gewIaoPji3tABywO4ggnBYP4LyYJB5jepFsZWVH\nUwkopZYAKAjx1t0kFwTt92cAmSSnm36oUiW02n2RACQec5wWk8RjTCTxmHa8krzA4ue+AmAhAFMJ\nCILgHKJqE1BK9Qj6czyAr6MLRxCERBPtEKwHlVK9APgA7ATwW4vHzYnyc2ONxGOO02KSeIyxHI8t\nw4YFQXAOsR/wLQiCqxAJCILHsU0CSqm/KKXWKaXWKqU+UEp1tCsWfzyzlFJf+2N6WymVZ3M8lyml\nvlJK+ZRStnU9KaXGKaW+UUptUUr9ya44guJ5QSm1TynliHEmSqlCpdRSpdRG///rjzbHk6mUWqmU\n+tIfzwzTg6wOLYz1hmYOOY5jPGMBpPmfPwTgIZvj6QOgF4BlAIbYFEMqgK0AugHIAPAlgL4258sI\nAIMBbLAzjqB4OgAY7H/eCsC3duYRAAUgx/88HcDnAM4xOsa2mgAdNuSY5Ack6/1/fgags83xbCL5\njZ0xABgGYAvJbSRrAcyH7gq2DZIfA/jezhiCIVlO8gv/86MANgHoZGM8JFnp/zPdvxmWLVvbBJRS\nf1VKlQG4CsA0O2NpwnUAFtkdhAPoBKAs6O9dsPEL7nSUUl0BDIL+9bUzjlSl1FoA+wAsJmkYT1wl\noJRaopTaEGIbDwAk7yZZCD3a8KZ4xmIlHv8+dwOo98dkezyCO1BK5QB4C8AtTWq5CYdkA8mzoWuz\nw5RS/Y32j+t8TXTYkGOzeJRSkwBcBGA0/RdVdsbjAHYDKAz6u7P/NSEIpVQ6tABeIfkvu+MJQLJC\nKbUUwDgY3LBnZ++Ao4YcK6XGAbgDwM9JVtkZi4NYBaCHUuoMpVQGgCsAvGNzTI5CKaUAPA9gE8m/\nOSCe/EDPllKqJYAxMClbto0YVEq9Bd36fWLIMUnbfmWUUlsAtABw0P/SZyStDoOORzyXAPg7gHwA\nFQDWkvyJDXH8FMBj0D0FL5D8a6JjaBLPawBGQt+6uxfAdJLP2xjPeQCWA1gP/V0GgLtILrQpnoEA\n/gH9/0oB8AbJ+wyPsUsCgiA4AxkxKAgeRyQgCB5HJCAIHkckIAgeRyQgCB5HJCAIHkckIAge5/8B\nbv7mbE8wDRIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a154a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.mlab as mlab\n",
    "import matplotlib\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = plt.axes()\n",
    "\n",
    "ax.set_aspect(1)\n",
    "theta = np.linspace(-np.pi, np.pi, 200)\n",
    "\n",
    "delta =1\n",
    "x = np.linspace(-3,3,7)\n",
    "y = np.linspace(-3,3,7)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "Z = X+Y\n",
    "\n",
    "matplotlib.rcParams['contour.negative_linestyle'] = 'solid'\n",
    "plt.scatter([-2, 2],[-2, 2],s=90,c = 'b',zorder=3)\n",
    "plt.plot(math.sqrt(8)*np.sin(theta), math.sqrt(8)*np.cos(theta),'r',zorder=2)\n",
    "plt.contour(X, Y, Z,\n",
    "                6,\n",
    "                 colors='k',\n",
    "                zorder=1)\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moving on to a more difficult example where there are two inequality constraints consider the following optimization problem. The below plot illustrates the function(in black) and constraints (in red).\n",
    "<center>$\\;\\;\\;\\;\\;\\;\\;\\;\\;\\text{minimize}\\;\\; f(x) = 2x_1^2+3x_2^2-2x_1x_2+x_2$<br>\n",
    " $\\text{subject to} \\;\\; h(x) = x_1^2 + x_2^2 - 8 \\le 0 $<br>\n",
    " $\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;g(x) =\\;2x_1+x_2-5\\le0$\n",
    " </center>\n",
    " The necessary first order conditions are:\n",
    " <center> $$4x_1-2x_2+2\\mu_1 x_1+2\\mu_2 =0$$\n",
    " $$6x_2-2x_1+1+2\\mu_1 x_2 + \\mu_2 = 0 $$\n",
    " $$\\mu_1(x_1^2 + x_2^2 - 8)+\\mu_2(2x_1+x_2-5) = 0 $$\n",
    " $$\\mu_1\\ge 0 \\;\\;\\;\\;\\;\\;\\;\\; \\mu_2 \\ge 0 $$</center>\n",
    " We will now test different combinations of active and inactive constraints to see which solution is the optimum. <br>\n",
    " - __Both constraints are active__<br>\n",
    "  Looking at the below figure the points of intersection are\n",
    "  $$x_1 = 2.7746 \\;\\;\\;\\; \\text{and}\\;\\;\\; x_2 = -0.5492$$\n",
    "  and\n",
    "  $$x_1 = 1.2254\\;\\;\\;\\; \\text{and}\\;\\;\\;\\; x_2 = 2.5492$$\n",
    "  The solutions are $(\\mu_1,\\mu_2) =(-3.6,3.89016)$ and  $(\\mu_1,\\mu_2) =(-3.6,4.50984)$.\n",
    "  Each of these solutions does not satisfy $\\mu\\ge 0$\n",
    " - __$g_1$ is active, $g_2$ is inactive__ <br>\n",
    "This leads to the following equations:\n",
    "$$4x_1-2x_2+2\\mu_1 x_1 =0$$\n",
    " $$6x_2-2x_1+1+2\\mu_1 x_2 = 0 $$\n",
    " $$x_1^2 + x_2^2 - 8 = 8 $$\n",
    " There are four solutions, each one without satisfying $\\mu_1\\ge 0$\n",
    " - __$g_1$ is inactive, $g_2$ is active__ <br>\n",
    "  This gives the following equations:\n",
    "  $$4x_1-2x_2+2\\mu_2 = 0$$\n",
    "  $$6x_2-2x_1+1+\\mu_2 = 0$$\n",
    "  $$2x_1+x_2 = 5$$\n",
    "  Resulting in the following values:\n",
    "  $$x_1 = 2 \\;\\;\\;\\;\\;\\ x_2 = 1 \\;\\;\\;\\;\\; \\mu_2 =-3$$\n",
    " - Since none of these scenerios have $\\mu \\ge 0$, we set all the constraints to inactive which results in the following equations:\n",
    "$$4x_1-2x_2 =0$$\n",
    "\n",
    " $$6x_2-2x_1+1 = 0 $$\n",
    " resulting in the point $(x_1,x_2) = (-1/10,-1/5)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztnXl8TNf7xz93Jnsim8Qu9loaUaSt\nqpYqKihVaimqaGlrL619rVpKLf1VLVW0RRdr8aXUUrQoIaKoJbbYZRGSkGVmPr8/biaCLLPcZSa5\n79drXtnuPeeZm7mf+5znPOc5AkloaGgUXXRqG6ChoaEumghoaBRxNBHQ0CjiaCKgoVHE0URAQ6OI\no4mAhkYRRzIREARBLwhClCAIm6VqU0NDQ36k9AQGA/hPwvY0NDQUQBIREAShHIDWAJZI0Z6GhoZy\nuEjUzlwAnwIoltcBgiD0BdAXALy9vevXqFHDrg4NBgMuX76MjIwM1KhRA4Ig2NWeJcTFxeHq1asI\nCQlB8eLFZe/PEoxGI+Lj4xEXF4fMzEzUqVMHOp36oZ5+16/j/Rs30LlWLZz39FTbnCLHkSNH4kkG\nW3QwSbteANoA+Cbr+yYANhd0Tv369WkPf/31F8uXL09XV1d++eWXNBqNdrVXEGfOnGHjxo0JgK+8\n8grPnTsna3+WcOTIEfbp04eenp4EwBdffJGrVq1ienq62qaJJCSQvr5kx45qW1IkARBJS+9hSw/M\nswFgGoCrAC4BuAngPoAV+Z1jqwgYDAZOmTKFer2elStX5qFDh2xqx1IyMjI4bdo0uru708/Pj99+\n+y1NJpOsfRZEdHQ0W7duTQD08vJi3759eezYMVVtypNx48SP2PHjaltS5FBUBB5pTEZP4Nq1a2za\ntCkBsGvXrrx7967VbVhDZGQkn3nmGQLgm2++yevXr8vaX0GcP3+e3bp1oyAI9PPz49SpU3nnzh1V\nbSoQzRtQjUInAlu2bGFQUBC9vLy4dOlSWZ/Gqamp/OSTT6jT6Vi6dGmuXbtWtr4s4ebNm+zfvz9d\nXV3p4eHBESNGMCEhQVWbrELzBlRBNRGw9GWpCKSnp3PYsGEEwLCwMJ46dcrWa2IRO3fuZJUqVQiA\n77//vqpP2qSkJI4dO5be3t7U6/Xs27cvr169qpo9NuNA3oDBYGBcXBz/++8/7tu3j+vXr+f//vc/\nxsTE0GAwqG2epFgjAlLNDkjO+fPn0aVLF0RGRuKjjz7CrFmz4ClTlPnOnTsYPnw4li5diqpVq2LX\nrl145ZVXZOmrINLS0jB//nxMnToViYmJ6NSpEz777DM89dRTqthjN4GBwODBwGefAf/+C9SuLUmz\nJpMJSUlJiI+PR3x8PBISErK/z+t1584ds8f6BG5ubqhatSqqV6+OGjVqoHr16tmvgIAASWx2VIS8\nLoqchIeHMzIyMs+///TTT+jXrx/0ej2WLl2K9u3by2IHSaxduxYDBgxAfHw8hg8fjgkTJsgmNvlh\nMBjwww8/YMKECbh69SpatGiBqVOnon79+orbIjmJiUClSkCLFsDq1XkeZjKZ8N9//+HSpUsF3tgJ\nCQkwmUy5tuPm5oagoKACX8WLF8eDBw9w5swZnD59GmfOnMGZM2dw/vx5GAyG7PaCg4MfEQWzUFSq\nVAmurq6SXy4pEAThCMlwi451JBFITU3FoEGDsHTpUjRs2BCrVq1ChQoVZLHh+vXr6N+/PzZs2IB6\n9ephyZIlqFu3rix95QdJrF+/HmPGjMHp06fx7LPPYvr06WjatKnitsjK+PGiN3D8eLY3YDQaER0d\njT179mDPnj3Yt28fEhMTHznNxcXlkZvWkpvb29vbrryRzMxMXLx4MVsUcr5u3779iG1VqlR5QiCq\nV6+OoKAgRXJX8sIaEXCYmEB0dDRr1KhBQRA4ZswYZmZmSjI2ehyj0chFixbR19eXHh4enDFjhmx9\nFcSuXbv43HPPEQCrV6/OtWvXqj4FKRsJCTT5+jKhaVN+8cUXbN26Nf38/AiAAFi5cmX26tWLy5Yt\n48GDBxkTE8OkpCSHux6JiYk8cOAAly9fzlGjRvHNN9/k008/TTc3t+z3AoABAQF84YUX+O6773La\ntGlct24dL126pJidcKbAoMlk4vz58+nu7s5SpUpx586dclwTkuTZs2cdIunnyJEjbNGiBQGwXLly\nXLJkiWpCJCfp6en866+/+Pnnn7NFixac5upKAgzNEr2+ffty5cqVvHLlitqm2o3BYOD58+e5ZcsW\nzp49m/369WOTJk1YunTpbGEQBIFvvfWWInkdTiMCCQkJbN++PQEwIiKCt27dkuWCPJ70s2TJElWe\nMGfPnmXnzp0JgIGBgZw1axbv37+vuB1y8eDBA+7evZsTJ05k06ZNs7MZATA0NJTDe/dmuqcnH7Rp\no7apinL37l0ePnyYo0aNYrFixQiAbdq04YEDB2Tr0ylEwJz66+LiwlmzZsmW+psz6adDhw6qJP1c\nu3aN/fr1o16vp5eXF8eOHcukpCTF7ZCalJQU/vHHHxw7dixfeumlbJdYEAQ+88wzHDx4MNetW8e4\nuLiHJxXxvIHExEROmjSJgYGBBMBXX32Vu3fvlvyh5PAiUKZMGdlTfzMyMjhq1KjspJ9169bJ0k9+\nJCYmcsSIEfT09KSLiwv79+/PGzduKG6HVNy9e5dbtmzhiBEj2KBBA7q4uBAA9Xo9n332WQ4fPpyb\nNm3KP7/CgfIG1OTevXucOXMmS5YsSQBs2LAht2zZIpkYOLwIQObU37NnzzI8PJwA2Lt3b8WTflJT\nUzl9+nT6+/tTEAR269aN58+fV9QGqYiKiuKnn37K8PBw6nQ6AqCrqysbNmzIUaNGcevWrbx37551\njRZxbyAn9+/f59dff83y5csTAOvVq8e1a9fa7Rk7vAjUqFFDljG5yWTid999R29vbwYEBCie8puR\nkcGFCxdmB4NatWrluIt78uHmzZv88ssvGRYWln3Tv/zyyxw3bhx37NjB1NRU+zrQvIEnSE9P53ff\nfceqVasSAGvVqsUVK1bYHDB2eBGwdylxbiQkJLBjx47ZkX+lI84bN25ktWrVsl27PXv2KNq/vaSl\npXH16tVs06YN9Xo9AfDZZ5/l119/zfj4eOk71LyBXMnMzOSqVasYGhqaPXW6ePFipqWlWdVOkROB\n3bt3s1y5cnRxceGMGTMUzQO/fPky33jjDQJgzZo1uXHjRoeb284Lk8nEAwcO8IMPPmBAQAABsEyZ\nMhwxYgRPnjwpb+eaN5AvRqOR69evzx7WlitXjvPmzbPYCysyIpCens6RI0dSEAQ+9dRTjIyMlKRd\nS8jIyODMmTPp5eVFT09PTp8+3XEKehTAlStXOHXqVFavXp0A6OHhwbfffpvbtm1TdiGN5g0UiMlk\n4rZt2/jSSy8RAEuUKMHp06cXGE8rEiKQM/j33nvvMSUlxe42LeWvv/7Kdtdef/11Xrx4UbG+bSUl\nJYU//vgjmzVrRkEQCIAvvfQSlyxZIntthjzRvAGr2Lt3L1977bXsjMQJEybkuay8UIuAmsG/+Ph4\nvvfeewTA8uXLc8OGDYr1bQtGo5F//vkne/XqRR8fHwJgxYoVOWHCBMbExKhtnojmDVjNoUOHsoeg\nPj4+/PTTT3nz5s1Hjim0IqBW8M9kMnHZsmUMCgqiXq/nJ598wuTkZEX6toWYmBiOHz+eFStWzP6g\n9OrVi3/++afs9RitRvMGbOb48ePs2rUrdTodPTw8OHDgQMbGxpIspCKgVvDvxIkT2eOxF198kccd\n9ImVlJTEb7/9lo0aNcrO2mvWrBl//PFHRYdKNqF5A3Zx9uxZ9u7dmy4uLnR1dTV7q4VHBNQK/qWm\npnLkyJF0cXFhYGAglyxZ4nBPUYPBwN9//51du3alh4dH9mrEqVOnZj8RnALNG5CES5cusX///nR3\nd1e82rAHgEMAogGcBDCpoHMsFQG1gn+bNm1ihQoVCIC9evV6NPfdAbh69SpHjBjBMmXKZAeJPvzw\nQx48eNBppiefQPMGJOP69euKi4AAwCfre1cA/wBokN85BYmAWsG/nHP+Tz/9NPfu3atIv5Zy9epV\nDhw4kO7u7tTr9WzTpg1Xr15tdSKJQ6J5A5Ki2nAAgBeAowCez++4/EQgISGBHTp0UDT4l5mZyVmz\nZtHb29sh5/yvXbuWffO7uLiwT58+TrsWIV80b0AyFBcBAHoAxwCkAJhR0PF5icDJkydZuXJlRYN/\nFy9e5AsvvOCQc/7p6emcMWMGvb29C/fNb0bzBiRDTU/AH8BuAKG5/K0vgEgAkSEhIU8YvWXLFvr6\n+rJkyZKyFlvIyc8//0w/Pz/6+vryp59+UqRPS/njjz9Yo0YNAmDbtm0dZ15fbjRvQBJUnR0AMB7A\n8PyOeby82OzZs6nT6fjMM88oEtVOSUlhnz59CIANGjRwqKf/lStX2KlTp+zFI5s3b1bbJGXRvAFJ\nUDowGAzAP+t7TwD7ALTJ7xyzCKSnp2dn4LVv316R6H9UVBSrV69OQRA4evRoZmRkyN6nJeR0/T08\nPDh58mQ+ePBAbbPUQfMG7MYaEZBi85HSAL4XBEEPQAfgV5KbCzopPj4eHTt2xJ49ezBmzBhMnjxZ\n1i21SeL//u//8Mknn6B48eLYsWOHw5T13rFjBwYOHIjTp0+jXbt2mDNnDipVqqRM5wYDcPMmEB//\n8BUXB6SkiH/L+XJzA3x8xFexYuLXkiWB8uWB0qUBqWrwDxkCzJsHTJ6c7z4FGhJhqVpI+Xr66adZ\nuXJluru7c8WKFbKpoZm4uDi2adMmu8Cjo8z7K+r6379P7t9Pzp9PDhlCtm5NPvUUmVUBON+XXk+6\nuZGCkPcxgkCWLk0+/zz57rvkl1+S27eTN26QtuQuaN6AXcAKT0CVzUf0ej2Dg4Px22+/4fnnn5e1\nr127dqF79+5ISEjArFmzMGDAAFU3hQCAjIwMzJ07F5MnT4bRaMTo0aPxySefwMPDQ7pOLl8Gdu4E\n/vkHOHQIOHFCfJoDgLc3ULWq+KpWDahYEQgOBoKCHn718RGf7Ho9YL5eJPDggeglpKQA9+6JXsSV\nK8DVq+LXy5eBkyeBW7ce2lKiBNCoEfDSS8DLLwN16ojt5oeFuxZp5I7Dbz4SEBAgewDQaDRywoQJ\nFASB1atXZ1RUlKz9WUrOqH+7du144cIFaRpOSSE3bSIHDiSrV3/4hPb3J5s3J8eMITdsIK9cse3J\nbC23b5M7d5Lz5pHvvENWqvTQJl9f8vXXySVLyPzKzGvegM2gMK0dsIXk5OTs/Qzeeecdh1hAk9P1\nr1KlCv/3v//Z3+iDB+S6dWSnTqSnp/jv9PQkIyLIOXPIEyeUueEtJTaWXLmS7NuXDAl5OIxo2JD8\n4gvy8uVHj9dmCmymSIvAhQsXWLt2bep0Os6ZM0f1XHrJo/4mE7l3L9mzp3iDAGRwMPnRR+Qff4jC\n4AyYTGRUFDlxIlm37kNBaN6cXLVKjGGQmjdgI0VWBHbt2sXixYvT39+f27Ztk6UPa5DU9U9OJhcs\nIGvXfuhS9+pFbttGFoYtzC5cEAWhQgXx/fn5icJ26JDmDdhAkRSB+fPnU6/Xs2bNmjx79qzk7VuD\npK7/1avk4MFksWLiv6tuXXEsbW/ZbwfCZDLx8uXLXL9+PSeMG8clb7/Nk/Xr0+DiQpMgMDVr6BC/\ne7dDrelwZKwRASnyBFQlIyMDgwYNwqJFi9C6dWusWrUKvr6+qtmSM+o/efJk26P+sbHAjBnAkiWA\nyQR07gz07w80aPAwWu+EkMTFixdx5MgRHD16NPsVHx8PANDpdBAEAUajESUADADQPzYWXgD4yito\nC2C/jw8CixdHYGAgAgMDUTyP73P+HBAQADc3NxXfuePi1CIQFxeHDh06YN++fRg5ciSmTJkCfUFT\nTzLx119/4f3337c/4efGDWDSJGDpUvHnXr2AkSPF6TInw2Qy4dy5c9k3+pEjRxAVFYWkpCQAgKur\nK0JDQ9GuXTvUq1cP9erVQ1hYGDw8PJCcnIzExEQkJCQg6vp11Bw/HmWio/E7gHP+/vipVi0c1uuR\nkJCAq1evIjExEYmJiTAajXnaU6xYsScEIjw8HB988AGKFSum0FVxPFTJEwgPD2dkZKRdbURHR6Nd\nu3a4desWli5diq5du0pknXVkZGRg4sSJmDFjBipUqICvv/4arVq1sr6htDRg9mxg6lQgIwN47z3x\n5g8Jkd5oGTAYDDhz5swjT/ioqCikpKQAANzd3REWFob69etn3/ChoaFwd3e3rANz3kDlysC1a2JW\nY/v2wLRpQPXqAEQv4969e9mCkJCQkP19bj/HxcXh3LlzCAwMxPDhwzFgwIBCIwYOnydgb0xg9erV\n9PLyYtmyZXn48GG72rKHU6dOsV69egTAPn36WL8nHylGyVevJitWFMf87duTDr5iMD09nVFRUfzu\nu+/Yv39/NmjQ4JFtyL28vNiwYUMOGDCAy5YtY3R0tDRrNMwzBQcOkJMmiXESV1dy7NiHswlWcujQ\nIbZu3ZqAuF381KlTbfs/OhgorIFBo9HI8ePHEwBfeOEF1Xb4NZlM/Oqrr+jh4cGgoCCuX7/etoZi\nY8mWLcV/Q1iYmFzjYJhMJkZFRXHhwoXs27cvw8PDs7cgB8BixYqxcePGHDp0KH/88UeeOnVKvjoQ\nj+cN3LxJ9ughXr/KlcmtW21uurCJQaEUgZwJQL169VKtpNa1a9eyN4CIiIiwTYhMJnLxYvFJ5u1N\nfvUVqeTOPxZw7do1Tp8+PXuKE1m1DF999VV++umn/Pnnn3n27Fnli6/mljewa9fDLMkuXcjERJub\nLyxiUOhEIGcC0Ny5c1VLAFqzZg0DAwPp6enJb775xjY7YmPJZs3ES//KK6QDVQp68OABf/75Z7Zs\n2TJ7G/IXX3yRixcv5sWLF1VPvCKZdxZhWho5eTLp4kKWK2e3V+XsYlCoRGD37t2qJwDdvXuXPXv2\nJACGh4fz9OnTtjX0v/+RgYGkj4+Y+OMAJcxNJhP/+ecffvjhh/T39yeydlcaO3as6vkWeZJfFmFk\npOgVCAI5fLgoDnbgrGJQaETgm2++oYuLi6oJQHv37mXFihWp0+k4duxY2wJcmZnkyJHi5a5Th3SA\nm+vatWucMWMGa9asSQD09PRk9+7duWPHDofbX+EJClpTkJJCfvCBeL3r139yTYINOJsYOL0IpKen\ns1+/fgTA1q1bq7JhZnp6OkeNGkVBEFi5cmX+/ffftjV04wb50kvipe7b1+YothQ8ePCAv/zyCyMi\nIrLd/UaNGqm7KamtWLKmYP16Me4SHCyut5AAZxEDpxaB27dvZ2/7NXLkSGW3ys4i59Rf7969bf8n\n//uvuFrOy4tUoHhKbjilu28Jlq4w/O8/sXiKiwv5zTeSrap0dDFwWhE4duwYK1SoQA8PD65atUqq\n62Exkk39kWJVHV9fsdrOkSPSGWkhTu3uW4qlKwyTksRKSgDZv7+ksRhHFQOnFAG1E4Akmfoz8+23\nYkmu2rXF2QCFKFTuviVYU2/AaBQDhQDZuTMp8UIkRxMDRUUAQHmIew2cgrgX4eCCzskpAo6QAJRz\n6m/+/Pn2TYV9+aV4WVu2JBW48Qqtu28p1tYb+OIL8fjmzcXl2RLjKGKgtAiUBlAv6/tiAM4CqJXf\nOWYRUDsBSLKpPzNTp4qXtFMnUuZS5kXC3bcEW6oPLV0qemrPPUfeuSOLWWqLgarDAQC/AWie3zH1\n69dnQkIC69Wrp1oC0Llz51itWjX7pv7MmExiQQyA7NZN1iIfmZmZnDBhAvV6feF39y3FlupDGzaI\n6w4aNCBlvDlzikG1atWYqVABGNVEAEBFALEAfHP5W/Y2ZOXKlWPdunXp7u6uyg47+/fvZ1BQEIsX\nLy7NzsOffSZeynfflTX998KFC2zYsCEBsHv37kXD3bcEW2sRrl8vegSNG8tepGXZsmUEwJ0KrQ9R\nRQQA+AA4AuDNgo719PSku7s7t9qx4MNW1q1bRw8PD1apUkWam2jRIvEyvvOOrBmAK1asoK+vL/38\n/FSZOXF4bK1FuGqVmF3YooXd2YX5kZqaSm9vb/br10+2PnKiuAgAcAWwDcDHlhyv0+lUEYC5c+dS\nEAQ2aNCAt2/ftr/B9etJnY5s1Uq2GEBSUhK7deuW7fpfunRJln6cHnsqE3/3HbNnDWQU8s6dOzM4\nOFiRIYHSgUEBwA8A5lp6Tp06deS9Ao9hMBg4ZMgQAuKeh6lSuH779pEeHuKOOzKVNP/7779ZsWJF\n6vV6Tp48WbHxpNNiT2XiGTPEc8eNk96uLNasWaPYkEBpEWiUtdT0OIBjWa9W+Z0j974DObl//z7f\nfPNNAuDgwYOlyUC8dIksXlzMRJNhSzNz8E+n07FSpUrcv3+/5H0USuzxBkwmsndv8Zb48UfpbaOy\nQwKnTBaSg7i4OL7wwgsUBIFz5syRptH798WKv76+5Jkz0rSZg5zBvx49ehTtqL8t2OMNpKeTTZqI\n+y7+9Zf0tlG5IYEmAiRv3brFWrVq0cPDg2vWrJGmUZNJDAAC4pZfEmMO/vn6+mrBP1uxd9eihASy\nalWyTBlxKzWJUWpIUORFwCwAnp6e3LVrl3QNf/21eMkmTJCuTT4Z/Lt48aKk7Rc57N216Ngx0t1d\nzPqUOFCo1JCgSItATgHYvXu3dA3/+6/oJrZqJekHQwv+yYAUexguWCDeHtOnS2dXFkoMCYqsCMgm\nAGlpYjGQEiXy30XXCrTgn8zY6w2YTORbb4nJRLbWksgDJYYERVIEZBMAkhwxQrxUGzdK0pwW/FMA\nKbyBpCRxS/Vq1SQtBqPEkKDIiYCsArBvn5hR9v77kjSnBf8URIodjXfsENsYMUI6uyj/kKBIiYCs\nApCeTtaoIT4N7Fx2mjP49+KLL2rBPyWQwhsgyT59xGFBZKQ0dlH+IUGREQFZBYAUg0IAaeciJy34\npyJSeAN37ogVourUkSw9PDU1lV5eXrINCYqECMguAJcvi7UB27WzuQkt+OcASOUNrFsn3i5ffSWN\nXZR3SFDoRUB2ASDJN98kPT3FFGEb0IJ/DoQU3oDJRL76qrhvhB07HOVEziFBoRYBRQRg717x0kyZ\nYtPpWvDPwZDKGzh2TAwSDx0qiVlyDgkKrQgoIgAmE9mokTgGtHK1oRb8c2Ck8AZIMUjo6irZBjJy\nDQkKpQgoIgCkuLMtINaot4LDhw9nB/8mTZqkBf8cDam8gRs3xE1ku3SRxCy5hgSFTgQUEwCjUVwh\nWLGiVSWp161bR09PT4aEhGjBP0dGKm/g00/FYjISrCKVa0hQqERAMQEgxeKTALl8uUWHm0wmzpo1\ni4Ig8Pnnn+fNmzfltU/DPqTyBm7eFAvK9OoliVlyDAkKjQgoKgCkuGdghQoWVQvOzMzM3i+xY8eO\nvK/iHoMaViCVNzBwoLi1mQTl3uQYEhQKEVBcAA4dEi/H7NkFHnr37t3s3YpGjhxZtOr8OztSeQOx\nsWKAcMAAu02SY0jg9CKguACQZNeu4g62BcznX758mbVr16aLiwu//fZbZWzTkBapvIEePcTPjAQ7\nGUk9JFCj2vBSALcBnLDk+PxEQBUBiI0Vc8M//jjfww4fPsxSpUrRz8+Pf/zxhzK2aUiPVN7A/v3i\nLbRokd0mST0kUEMEXgZQz14RUEUASHLyZPFSXLiQ5yHmGYCKFSvy5MmTytmmIQ9SZRGGhYkzSnbu\noCX1kECtzUcq2iMC9+7dY2hoqPICYDKJNeWaNMnzkO+++06bAShsSOUNfPONeBv984/dJnXq1InB\nwcGSbMlnjQjooBCCIPQVBCFSEITIuLi4R/5GEv369cOpU6ewceNGNGnSRCmzgAMHgJgYoGfPXP+8\ne/du9OvXD82bN8fu3btRsmRJ5WzTkI/AQGDwYGDNGuDff21vp1s3wMsLWLbMbpN8fHxgMpnsbsdq\nLFWLgl6wwxNYvHgxAXCKjbn6dtG3r7haMJdNKc+dO8fAwEDWrFmTSUlJytumIS9SeQOdOpHBwXZv\nRFuvXj02a9bMPluygCN6AnkRHR2NgQMHonnz5hg1apSynWdkAL/+Crz5JlCs2CN/unv3Ll5//XUA\nwKZNm+Dn56esbRryI5U30LkzEBcH7NljcxOZmZk4ceIE6tata7sdNqKqCCQnJ6NTp04IDAzEihUr\noNMpbM6+fUBSEtCx4yO/NhgM6Ny5M2JiYrB27VpUqVJFWbs0lGPoUMDXF5g82fY2IiIAb2/xgWIj\np06dQkZGhvOKgCAIPwE4AKC6IAhXBUHoU9A5zIoDxMTE4KeffkKJEiWkMMU6Nm0C3N2BZs0e+fXw\n4cOxbds2LFiwQNn4hIbyBATY7w14egJt2wJr1wIGg01NREVFAYAqIiBZTMCaV/369dWNA5DirECl\nSuI+AjlYtGgRAXDIkCHq2KWhPFLEBn7+WZwlOHjQptMHDx5MLy8vafbKpEpThNa8atWqRXd3dzZv\n3ly9lNsTJ8S3v2BB9q927dpFFxcXtmzZUlsKXNSwN28gLo72FKJ5+eWX2aBBA9v6zgVrRECVmMD5\n8+fViwOY+eMP8WurVgCAmJgYdOzYEdWqVcPPP/8MFxcXdezSUIchQ+yLDQQFAXXrPvxcWYHJZMKx\nY8fUGQpApcBgenq6enEAM3/9BVSsCISEICkpSZsJKOpIMVPQrBmwfz+QmmrVaRcvXsS9e/eKlgjU\nqlULjRs3VqNrERL4+2+gUSMYDAZ06dJFmwnQsN8bePVVIDNT/GxZwbFjxwCoFBSESiLg6empRrcP\nuXABuHkTePFFbSZA4yH2egMNGohfDx+26rSoqCjo9XqEhoZa36cEqJ4spApZSr36xg3MmzcPQ4YM\nwXvvvaeyURoOgT3egJ8fULUqcOSIVadFRUWhZs2a8PDwsL5PCSiaIhAVBaOHB7p9/jkiIiIwc+ZM\ntS3ScBTs9QbCw20SAbWGAkARFYH7hw/jeGYmqj71FH766SdtJkDjUYYMEdPIbfEG6tcHYmPFNGIL\nuHXrFm7cuKGJgJIkJSUh+Z9/cEav12YCNHLHHm/AfDNHR1t0uNpBQaCIiYDBYECfDh1Q0mDAc716\naTMBGnkzdKht3kC1auLXCxdIq7RVAAAgAElEQVQsOtycLlynTh3r+pGQIiUCw4YNw41duwAAlbPy\nAjQ0csVWb6BsWcDVFTh/3qLDo6KiULFiRQQEBNhoqP0UGRFYvHgxvvrqK3wQESH+olIldQ3ScHxs\n8Qb0evGzZYUnoOZQACgiIrB79270798fERER6PbKK+Ivy5ZV1ygNx8dWb6ByZYs8geTkZMTExGgi\nIDcxMTHo0KEDqlWrhp9++gn6mzfFclC+vmqbpuEM2OINVKgAXLlS4GHHjx8HSU0E5MS8JkAQhIcz\nAdevA2XKAIKgtnkazoAt3kDx4sCdO2J6ej6Yg4LPPPOMvVbaRaEVAZLo2bMnYmJisG7duoczATdu\nAKVLq2uchnNhrTcQGAgYjcC9e/keFhUVhaCgIJRVeWhaaEVg0aJF2LhxI2bOnPnoYqXkZDG9U0PD\nUqz1BgIDxa+JifkeZl4+LKjslRZKEfjvv//w8ccf47XXXsOgQYMe/eP9+2JMQEPDGszewKRJBR9b\nvLj4NR8RULOw6ONIVWOwpSAIZwRBiBEEYaQUbdpKeno6unbtCm9vbyxfvvzJoiWaCGjYgtkbWLu2\nYG/A21v8mk9dAXNhUbXjAYAEIiAIgh7AfAARAGoB6CoIQi1727WV0aNHIzo6GsuWLUOpUqWePEAT\nAQ1bsTQ2YH7w5BMYVLWw6GNI4Qk8ByCG5AWSGQB+BtBOgnatZteuXZg9ezY++ugjtGnTJveD0tLE\nCsMaGtZiaWzALAL57CZ09OhReHl5oZo5zVhFpBCBsgByTopezfrdI+S3DZlUzJ49G2XLlsWsWbPy\nPsjV1eay0BoaFnkDBYhAXFwcfvjhBzRr1gx6vV4GI61DscAgycUkw0mGBwcHS95+XFwcfv/9d3Tr\n1i3/ykWuruLOQxq5YjIBW7cCb7wBPPus+HXr1nwfakULS7yBAkRgwoQJSElJwbRp02Qy0jqkEIFr\nAMrn+Llc1u8U5ddff4XRaES3bt3yP9DNTROBPLh9GwgLAzp1An77DYiMFL926iT+XiYHzvkoyBtI\nSxO/5lIp6MSJE1i0aBE++ugj1KqlWujsEaQQgcMAqgmCUEkQBDcAXQBslKBdq1i5ciVq166NsLCw\n/A/URCBXTCagaVPg7FkgJeXRv6WkiL9v2lTzCAAU7A3cvSt+fSw1nSSGDh0KPz8/TJgwQQFDLcNu\nESBpADAAwDYA/wH4leRJe9u1hvPnz+PAgQMFewGAuGWUlSWhiwLbtgGXL4vFcnMjMxO4dAnYvl1R\nsxyX/LwBc6bgYyKwefNm7NixAxMnTkRxcy6BAyBJTIDkFpJPkaxC8nMp2rSGVatWAQDefvvtgg8u\nUULza3Nh0aInPYDHSUkBFi5Uxh6HJz9vIBcRyMjIwLBhw1CjRg18+OGHChpaME6fMUgSK1asQOPG\njVG+fPmCTyhRQhz8ajzCNQujONevy2uHU5GXN2AeDuTY7n7evHk4d+4cZs+eDVdXVwWNLBinF4Ej\nR47g7Nmz6N69u2UnaCKQK5auYSlTRl47nIq8vIFr18TUYTc3AMCZM2cwfvx4tG3bFhHmojYOhNOL\nwIoVK+Dm5oaOHTtadkLJkqJSmyO4GgCAfv0AH5/8j/HxAT74QBl7nIbcvIHYWLGmAACj0YjevXvD\n09MTCx10LOXUImAwGPDzzz+jdevW8Pf3t+wk8yPPgqIPRYnXXhM/t3l5qq6u4taNLVooapbjk5s3\nEBsLhIQAEIcB+/fvx1dffYXSDrqE3alFYNeuXbh165blQwEAeOop8evZs/IY5aTodMDu3UD16k96\nBD4+4u937XqYB6ORg5zeAClOs4SE4MyZMxgzZgzatm1r2cyVSjj1v3TFihXw8/NDq6ztxS1CE4E8\nCQ4Wy+WvXg20aydmDLZrJ/4cHS3+XSMXcnoD+/YBKSkwlS//yDBA7ZoB+eG0W++kpqZi/fr16NKl\ni3V7uAUFif+0M2fkM86J0emAli3Fl4YVDB0KzJsHjBsHAFh/7hz279+PH3/80WGHAWac1hPYuHEj\nUlJSbHOzqlfXREBDWszewN69AICPv//e4YcBZpxWBFauXIly5crh5Zdftv7k2rWBY8e0HFgNaRk6\nFHR1RbogINnLy+GHAWacUgTMKwbffvvtJysHWUKDBkBSkhYX0JCWwECkenjAncQPn3zi8MMAM04p\nAhavGMyL558Xv/7zj3RGaRR5zh0/DrfkZGQKAlofPaq2ORbjlCJg8YrBvKhRQ8zrPnhQWsM0iixG\noxFfvf023AAY2raFsGYNcPy42mZZhNOJgFUrBvNCpwOeew7Yv186wzSKNPPmzUPwyZMwCQI8584V\n8wY++0xtsyzC6UTAvGKwa9eu9jXUtKmo1DdvSmCVRlHm7NmzGDNmDDoEBkIIDxdTK23Zw1AlnEoE\ncq4YDMlKy7QZ80T4tm32G6ZRZDEajejVqxcCPDxQKyUFQpMm4h9s2cNQJZxKBMwrBiWZe33mGaBU\nKeD33+1vS6PIYl4bsLJXLwgZGUCzZuIfbN3RWAL2WznMdSoRsHrFYH4IgrhqZvt2cd84DQ0rMQ8D\n2rZtiyZ37gD+/oDZEwBU8QZSU1OtHio7jQjkXDEYEBAgTaOtWolbRf31lzTtaRQZzMMAT09PLPz6\nawibNgFt2mTXEACgijcwceJExMbGWnWOXSIgCMJbgiCcFATBJAhCuD1tFYRNKwYLonVrccuorGCj\nhoalPLJEOCYGSEgA2rd/8kAFvYFjx45hzpw5eP/99607kaTNLwA1AVQH8CeAcEvPq1+/Pq2lR48e\n9PPz44MHD6w+N1+6dSMDAsj0dGnb1Si0nDlzhh4eHmzbti1NJhP54YekhweZkpL7CWPHkgB5/Lhs\nNhkMBj733HMsUaIEExISCCCSlt7Hlh6YbyMyi0BKSgp9fHzYp08f669OQfzvf+Jl+O036dvWKHRk\nZGSwYcOGDAgI4PXr18n790k/P/FhkhcJCWSxYmTHjrLZ9dVXXxEAV65cSZKOKQIA+gKIBBAZEhJi\n1RuMjo4mAK5atcr6q1MQGRlkUBD51lvSt61RqMjIyOBbb731yM3GFSvE22jXrvxPltEb+PPPP+nm\n5saWLVuKngklFgEAOwCcyOXVjlaIQM6XtZ7AjRs3CIBff/21bVepIIYMIV1dyRs35Glfw+nJKQCz\nZ89++IcmTcgqVUijMf8GZPIGTp48SX9/f9asWZOJiYnZv3dITyDny1oRMBgM1Ov1HD16tFXnWczZ\ns+KlmDhRnvY1nJo8BeDcOfFz8/nnljUksTdw/fp1VqhQgSVLluTFixcf+VuhEwGSLFOmDHv16mX1\neRbTqhVZqpQWINR4hDwFgCQHDRI9yGvXLGtMQm8gOTmZ9erVo5eXFyMjI5/4u2IiAKA9xK3I0wHc\nArDNkvNsEYHw8HC+9tprVp9nMVu3ipdDjriDhlOSrwDEx5NeXuS771rXqATeQGZmJlu1akWdTsfN\nmzfneozinoC1L1tEoG3btgwLC7P6PIsxGsmnniLr1iWzgisaRZd8BYAUh44AefKkdQ3b6Q2YTCb2\n7duXALhw4cI8jyuUItCvXz8GBQVZfZ5VLF8uXpING+TtR8OhKVAAUlLI4sXJ11+3rQM7vIFp06YR\nAEeOHJnvcYVSBCZNmkQATJdzzJ6ZSVatSj7zjOYNFFEKFACS/PJL8dbZt8+2Tmz0BlauXEkA7Nq1\nK40FzEYUShFYvHgxAfDy5ctWn2sV338vXpZ16+TtR8PhsEgA7twhAwPJ5s3t68xKb8CcC9C4cWOm\npaUVeHyhFIHNmzcTAA8cOGD1uVaRmUlWq0aGhorfaxQJLBIAkhwxQrxtjh61r0MrvIG8cgHyo1CK\nwNGjRwmA65R4Qv/6q3hpvvlG/r40VMdiAYiNFdcI5JcibA0WeAP55QLkR6EUAdmzBnNiMomZYIGB\nomJrFFosFgCS7NmTdHMjrbgZ86UAb6CgXID8KJQiYM4aHDNmjNXn2kR0NKnTkf37K9OfhuJYJQC7\nd4u3y4gR0hqRhzdgSS5AfhRKESAVyBp8nP79RSGIilKuTw1FSE5OZrt27SwTgLQ0snp1slIlMjVV\nWkNy8QYszQXIj0IrAuHh4WzZsqVN59pEQgJZsqQ4ZailExcaLl++zDp16lCn01k2vDQnBv3+uzwG\nPeYNWJoLkB+FVgRkzxrMjQ0bxMs0YYKy/WrIwv79+1myZEn6+vpy69atBZ9w6pQYB+jaVT6jcngD\n1uQC5EehFQFFsgZzo0cP0sWFPHJE+b41JGPFihV0d3dn5cqVeerUqYJPSEsj69QR603cvCmvcVne\nQF0XF4tzAfLDGhFwmkKjAFCmTBnEx8cjIyND2Y7nzQOCg4GePYEHD5TtW8NuTCYTxowZg+7du+OF\nF17AoUOHULNmzYJPHDUKiI4Gli0DSpaU1cbTERG4B2CapyfWr18Pd3d3Wft7BEvVQsqXrZ6AYlmD\nuWFeZdi7t/J9a9hMcnIy27dvTwB8//33LU87N/+/Bw6U10A+zAX40tubuc0U2AIK63DAnDV48OBB\nm863G3MA57vv1OlfwypyBgDnzp2bXXqrQK5eJUuUELNGpS5s+xg5cwGidu6UrN5AoRUBRbMGc8Ng\nIF99Vcwa06YNHRqrA4Bm7t8nn32W9PGxfpmwleSaCyBR9aFCKwLmrMH58+fbdL4k3LpFli0rzhnf\nvq2eHRp5YnUA0IzJRHbvTiWWk+eZCyBR9aFCKwKKZw3mxcGDojfQoIH45NBwCIxGI0ePHk0AbNKk\nCePj461r4IsvxFtiyhR5DMxBvrkAEngDhVYESBWyBvNi7VpSEMj27cVhgoaq2BwANLNmjfj/7NRJ\n9loSBeYCSOANKCYCAGYCOA3gOID1APwtOc8eEVA8azA/5swRL+GQIWpbUqSxOQBoZts2sWBow4bS\npwU/hsV1Aez0BpQUgRYAXLK+nwFghiXn2SMCqmQN5segQdTKlauHzQFAM3//LRYMrVNHLBgiI1bV\nBbDTG1BlOACx8vBKS461RwT69evH4OBgm8+XHINBXGIKkJMnq21NkeKHH36wLQBo5tgx0t9fLCkn\nc0agTXUB7PAG1BKBTQC65/N3m7chy8kXX3xBAPxdrsUctmAwkO+8Q6WCSkWd2NhYvvHGG7YHAEkx\nuOvvT5YrJ119gDy7Osjy5cvT29vburoAdngDkooALNuGbExWTECwpFN7PIGUlBSGhYXR39+fMTEx\nNrcjOQbDw+mlCRO0QqUyYDAYOHfuXPr4+NDT05MzZsxgRkaG9Q3t3i3mAVSuLKsAmEwmfvPNN3R1\ndWXFihWtLgxC0mZvQFFPAMC7AA4A8LL0HHtEgCTPnz/PgIAAhoaGMjk52a62JMVgEDejAMj339dq\nFEpIZGQk69evTwBs2bIlL1y4YFtDmzeL07u1alm+c5ANpKamskePHgTAiIgIJthaocoGb2DTpk2K\nBgZbAjgFINia8+wVAZLctm0bdTodO3XqZH00WE5MJnLMGPHStmmT9571GhZx7949DhkyhDqdjqVK\nleIvv/xi+/97wQJxNWj9+mRcnLSG5uDcuXMMCwujIAicOHGiXUuCSVrlDXz77bfU6XSKikAMgCsA\njmW9FlpynhQiQJIzZswgAE6fPl2S9iRlwQKxKtFzz5HXr6ttjVOyYcMGlitXjoIg8MMPP+QdW6P3\nmZnkgAHix71VK/LuXWkNzcFvv/1GPz8/BgQE2DZbkRsWeAMmk4lTpkzJ9pQKdbLQ42+8c+fOFATB\nsQKFZjZsEKefSpUi9+5V2xqnIWfgr3bt2ty/f7/tjd25Q7ZoIX7UP/5YtsQug8GQna1Yr149qyoD\nW0Q+3oDRaOTAgQMJgN27d2dGRkbREQHSgQOFZv79V9zHwMWFnDtXCxjmQ0pKCmfOnGl/4M9MZKQ4\n/efqSi5ZIp2hj3H79m02a9aMANinTx8+kGPlodkb6NDhkV+np6ezS5cuBMBhw4ZlDz2KlAiQDhwo\nNJOURLZrJ17uzp1lT0pxNm7fvs3x48czMDAwO5Bmc+CPFDeX/fJL8eYvV8727cIs4J9//mH58uXp\n7u7OJTIKDcknvIF79+5li88XX3zxyKFFTgTIh4HCt956y7EChWaMRnLaNFKvFz+Y27erbZHqXLhw\ngf3796enpycBsF27dva5/qS4srNVK/Gj/cYbsu0bYTKZuGDBAvum/6wlR2zg1q1bDA8Pp16v5/Ll\ny584tEiKAOnggUIzhw6RNWqIl/6jj4rk7EFUVBS7du1KvV5PV1dX9u7d27aMv5yYTOSKFWI9QHd3\n8uuvZRt6paam8p133rF/+s8WsryB1iEh9PT05KZNm3I9rMiKgMMHCs3cv08OHSquWqtShdyyRW2L\nZMdkMnHnzp187bXXCIDFihXj8OHDefXqVfsbv3iRfO018ePcoIEYh5EJyaf/rOTE3r28Jwjc4OrK\nv//+O8/jiqwIkE4QKMzJn3+Km1oAYszAnnGwg2IwGLh69Wo+++yzBMCSJUty6tSptk/35SQtTawB\n4OUlZgD+3//Juqxbluk/K9izZw/9/Pw4p1gxFpQ3UKRFgHSCQGFO0tPJGTNIb2/RjR0/nnR0my3g\nwYMHXLRoEatVq0YArFq1KhcuXChN5NxkIn/6iaxYUfwIv/66uFmoTMg+/WcB69evp7u7O2vUqMEr\n0dEF5g0UeREgHwYKW7RowXPnzsnen91cuSLOHABkcDA5a5ZTVi26c+cOp02bxpIlSxIAw8PDuXr1\nahqkekLv2SMmYAHi8l+ZA6yKTP8VgDkL8Pnnn3+4WKqALEJNBLJYsGABPT09qdfr2bNnT+cQgwMH\nyGbNxH9N6dKii6vCB89arl69yuHDh7NYsWIEwBYtWnDnzp3SzNSYTGIJ8JdeEq9L2bLk8uWyV3RS\ndPovF3JmAUZERDAlZxC5gCxCTQRycOPGDX788cfOJwZ79pAvv/zQMxg/nrxxQ22rnuC///5j7969\n6erqSp1Oxy5duvDo0aPSNJ6ZSf7yC1m3rngdypUTE65krv7z+PTfERV2nsotC/AJ8vEGNBHIBacU\nA5OJ3LlTXIgkCGLyS48e4jSjyrkQ+/fv5xtvvEFBEOjh4cH+/fvz/Pnz0jR++bIoemXLih/R6tXJ\npUsV2RRW1em/LNLS0ti5c2cC4Mcff5z3DEQ+3oAmAvnglGJAkmfPirvh+PiI/7aaNcnPPycvXVLM\nBJPJxM2bN/Oll14iAAYEBHDcuHG8LUXp9fv3xWKfrVqJgicIZEQEuX69YoVcY2JiVJ3+I/PPAsyV\nPLwBTQQswGnFICmJXLSIbNRI/PcB4rDhq69IqZ7Ej5GRkcEffviBoaGhBMDy5ctzzpw59s+8JCeL\n7v5bb4mzIwBZpgw5bpyi4kaSGzduVHX6j+QjWYDLli2z7KQ8vAFNBKzAacWAFPMKPvvsYQYiIH4/\nbBi5Y4ddY+d79+5xw4YN7Nu3L8uUKUMAfPrpp/n999/bvqjHYBAX9cyYISb3eHiINpcsSX7wgWiz\nwoVYDAYDx4wZo+r0HymmUFerVi3fLMA8ycUb0ETABpxaDEhxuDB3Ltm8OenmJv5rzQU0Bg4kV60i\nY2LydK1NJhNPnDjBmTNnsmnTpnR1dc3O7Gvfvj03bdpknXtsMomVezZvFoWqXTuxpp9ZrGrVEis1\n79mj2r4NcXFxbN68OQHwvffeU2X6jySjo6NZunRpBgQE5JsFmCe5eAOaCNiB04sBKbrZmzeTo0eT\nr7wiZtSZbz4PD3F+vXNnpo0axcjBgzm7VSs2K1WKxQEKAENDQ/nJJ59w165deW/iYTKR9+6JQbwD\nB0SRmTpVLKvWrJm4oae5T0Bc0tunD7lypUMUWfntt99YqlQpuru78zsVN5g1ZwGWLVuWJ06csL2h\nx7wBTQQkoFCIgZnMTPLoUZq+/ZZxPXvyfI0avObhQWPOmzTrZXJxIQMDxbF55cpiroL576Gh4u8C\nA8WqSbmcz+BgcUPPd98l580Ti6nIWMnHWhISEti9e3cCYJ06dXjs2DHVbMmZBXj58mX7GnvMG7BG\nBATxeGUJDw9nZGSk4v3aws2bNzFz5kwsWLAAGRkZ6N69O8aOHYuqVauqbZpFpKSkYOfOndi6dSu2\nbt2K2NhYAEBoaCjaNmuG1+vVQ/0yZeCakADcvCm+kpOBtDTxdfs2sH07oNcD7doBHh5AQADg7//w\nVaoUUKkSUKEC4OOj8jvOm40bN6Jfv36Ij4/HmDFjMHr0aLi5ualiy5IlS9CvXz88++yz2Lx5M4KC\nguxvdNw4YMoU4PhxCGFhR0iGW3SepWqR2wvAZxC3IDsGYDuAMpac5wyewOM4i2dQ0Nh+8eLFjJUx\nz94RefzpH6XitvKP1wJMkXIpeQ5vAAoWGvXN8f0gKFxoVA0cUQySk5O5YcMG9uvXjyEhIQRAWDq2\nL+SYx/4uLi6cMGGCqtfBaDRywIAB+WcB2ktWbMAaEZBsOCAIwigAISQ/LOhYZxoO5EVuw4S+ffui\ndOnSivR/79497NixA1u2bMG+ffuQmZmJYsWKoVmzZoiIiEDLli1Rvnx5RWxxRBITEzF48GCsWLEC\nderUwfLly/HMM8+oZk96ejp69uyJX375BR9//DFmzpwJnU4nfUeJiUDFihCSk5UZDmQJyOcQy46f\ngIX7DzizJ/A4OT0DZD2BlXxpT/tHSUtL45w5cxgYGOgQT3+S/Pfff7M3TpkxY4as5e9SUlK4Xkzq\nks4TEARhB4BSufxpDMnfchw3CoAHyQl5tNMX4n6ECAkJqX/58uV8+3U2bt68iePHj+PGjRuK9Ofm\n5oZGjRoV6ad9TkwmE3755ReMGTMGFy9eRPPmzTFr1iyEhYWpZpPBYMCMGTMwadIk+Pv7Y9GiRWjf\nvr1s/Z0+fRodOnTAzVOnkAgo5wmYXwBCAJyw5NjC5AloqM/OnTuzn7R16tThtm3b1Dbpkad/p06d\npFlfkQ+//PILfXx8GBwczB07digaGKyW4/uBANZYcp4mAhpScPz4cUZERGSvZ/jhhx9UWfSTk8zM\nTE6ZMoWurq4MDg7m6tWrZe0vPT2dgwYNIgA2bNgwu2ajkiKwFmIs4DjErcnLWnKeJgIa9nDlyhX2\n6tWLgiDQ39+fM2fOVC3lNydKP/2vXLnCBg0aEACHDBnyyGyDYiJg60sTAQ1buHPnDkeOHEkPDw+6\nublx2LBhqqz3fxyln/4kuX37dgYFBdHHx4e//vrrE3/XRECjUJEz4o+sOXa1Vvs9jtJPf6PRyMmT\nJ1MQBD799NM8ffp0rsdpIqBRKDAajVy1ahUrVapEAGzWrJkqpb5yQ42nf3x8vHnHYfbo0SPfbENN\nBDScHkeM+JtR+ulPikVPQ0JC6ObmxoULFxaYa6CJgIbTcuzYsUci/t9//7105crtRI2nv8lk4vz5\n8+nq6soKFSrw8OHDFp2niYCGU2E0Grl58+bs2np+fn784osvHCLib0aNp39KSgrffvttAmCrVq2s\nCoJqIqDhFCQnJ/P//u//sncpKlOmDD///HOHiPibUePpT4ql3GvVqkWdTscpU6ZYnf+giYCGQ3Px\n4kUOGzaMfn5+BMDnnnuOq1atkmdVnR2o8fQnn8z+swVNBDQcDpPJxL179/LNN9+kTqejXq9n586d\neeDAAbVNewK1nv55Zf/ZgiYCGg5DWloav//+e9arV4/I2qtg5MiRDlvYRK2nf2xsbHb239ChQ+32\nijQR0FCdW7ducdKkSdkbk9asWZMLFy5kqsxbiNmKWk9/suDsP1vQREBDNaKiovjuu+/Szc2NyNrK\na9u2bbKuobcXtZ7+lmb/2YImAhqKYjAYuH79ejZp0oQA6OXlxY8++kjSD7UcqPn0tyb7zxYcXgR8\nfX25Zs0a1Zd9atjH3bt3OWfOnOy03pCQEM6cOZOJiYlqm1Yge/bsYZ06dRR/+pPWZ//ZgsOLgLu7\nOwGwdu3amhg4IefOneOgQYNYrFgxAmCjRo24evVqZiq8hZgtxMbGZu/4GxISwrVr1yrWt63Zf9ay\nb98+xxeB+vXrc8WKFXzqqac0MXASTCYTd+3axbZt21IQBLq6urJ79+6yfZCl5v79+5w0aRI9PT3p\n4eHBiRMnKhqkTE5Otjn7z1JMJhPnzJlDvV7vHCJAimNJTQwcm9u3b3Px4sUMCwsjAAYFBXHcuHG8\n7gBbiVmCyWTimjVrWKFCBQLgW2+9xUsK73h84MABVqtWzebsP0tISUlhly5dCIBvvPGG84iAGU0M\nHAeTycR///2XU6dOZcOGDSkIAgEwLCyMS5cudah8/oI4fvw4X3nllezP1O7duxXtPyMjg2PHjqVO\np2OFChX4559/ytLP2bNnGRoaSp1Ox2nTptFkMjmfCJjRxEAd0tLSuH37dg4cOJAVK1bMLmdev359\nTpw4kZGRkQ49xfc4CQkJHDBgAPV6PQMDAzl//nzF4xUnT57MTpB69913eVem/Rg3btxIX19fFi9e\nnNu3b8/+vdOKgBlNDOTn9u3bXL58OTt06JAd4PP09OTrr7/ORYsW8dq1a2qbaDUGg4ELFixg8eLF\nqdPp+NFHHzE+Pl5RG4xGI+fMmUN3d3cGBQVx3bp1svRjMBg4duzYbLF+fIijuAgAGJb19Aiy5HhL\n8wQ0MZCOvNz8MmXKsG/fvty0aZPDZvNZQs4pv8aNGzM6OlpxGy5fvsymTZsSANu0acMbN27I0k9C\nQkJ2jkHv3r1zHaIpKgIAygPYBuCy1CJgRhMD2yjIzT9y5IhTufm58fiU36+//qr4ezKZTPzxxx/p\n5+dHb29vfvvtt7LZcPToUVaqVIlubm5ctGhRnv0oLQJrANQBcEkuETBjMBi4cuVKVq9enQCcZnpK\nLUaOHFko3Pz86NKlCz08PDhhwgTVPJlbt27R19eXDRs2ZExMjKx9tWvXjmXLluXBgwfzPc4aEbBr\nQ1JBENoBaEpysCAIl7rYKhQAAALbSURBVACEk4zP49jsbcgAhELcr6CwEQQg1/dfCCis762wvq/q\nJItZcqBdexECGA2gBcm7BYnAY21G0tJ90pyIwvq+gML73rT3BbgUdADJZnl0UhtAJQDRgiAAQDkA\nRwVBeI7kTSvs1dDQUJECRSAvSP4LoIT5Z2s8AQ0NDcdBp1K/i1XqV24K6/sCCu97K/Lvy67AoIaG\nhvOjliegoaHhIGgioKFRxFFdBARBGCYIAgVBCFLbFikQBGGmIAinBUE4LgjCekEQ/NW2yR4EQWgp\nCMIZQRBiBEEYqbY9UiEIQnlBEHYLgnBKEISTgiAMVtsmKREEQS8IQpQgCJsLOlZVERAEoTyAFgBi\n1bRDYv4AEEoyDMBZAKNUtsdmBEHQA5gPIAJALQBdBUGopa5VkmEAMIxkLQANAPQvRO8NAAYD+M+S\nA9X2BOYA+BRiTnuhgOR2koasHw9CzJ9wVp4DEEPyAskMAD8DaKeyTZJA8gbJo1nfJ0O8Ycqqa5U0\nCIJQDkBrAEssOV41EchKOb5GMlotGxSgN4CtahthB2UBXMnx81UUkhslJ4IgVARQF8A/6loiGXMh\nPlxNlhxsc7KQJViScixn/3KR3/si+VvWMWMgupwrlbRNwzoEQfABsBbAEJL31LbHXgRBaAPgNskj\ngiA0seQcWUWgsKYc5/W+zAiC8C6ANgBepXMnYlyDuFTcTLms3xUKBEFwhSgAK0muU9seiXgRQFtB\nEFoB8ADgKwjCCpLd8zrBIZKFClPKsSAILQHMBtCYZJza9tiDIAguEIObr0K8+Q8DeJvkSVUNkwBB\nfPp8DyCR5BC17ZGDLE9gOMk2+R2ndmCwMPI1gGIA/hAE4ZggCAvVNshWsgKcAyAWjfkPwK+FQQCy\neBFADwBNs/5Px7KenkUOh/AENDQ01EPzBDQ0ijiaCGhoFHE0EdDQKOJoIqChUcTRREBDo4ijiYCG\nRhFHEwENjSLO/wMpWPFdbGzvFwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a8439d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.mlab as mlab\n",
    "import matplotlib\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = plt.axes()\n",
    "\n",
    "ax.set_aspect(1)\n",
    "theta = np.linspace(-np.pi, np.pi, 200)\n",
    "\n",
    "delta =1\n",
    "x = np.linspace(-4,4,9)\n",
    "y = np.linspace(-4,4,9)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "Z = 2*X**2+3*Y**2-2*X*Y+Y\n",
    "x1 = np.linspace(0.5,4,10)\n",
    "matplotlib.rcParams['contour.negative_linestyle'] = 'solid'\n",
    "plt.plot(x1,5-2*x1,'r',zorder =3)\n",
    "plt.plot(math.sqrt(8)*np.sin(theta), math.sqrt(8)*np.cos(theta),'r',zorder=2)\n",
    "plt.contour(X, Y, Z,\n",
    "                6,\n",
    "                 colors='k',\n",
    "                zorder=1)\n",
    "\n",
    "plt.scatter(-.1,-.2,s=90,c = 'b',zorder=3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An Engineering Application\n",
    "\n",
    "Suppose you were tasked with minimizing the tensile stress,$\\sigma$, in a beam. The beam is subjected to a 55 kip load. You needed to choose the width (w) and height(h) of the beam. However, due to sizing limitations, the area needed to be less than 20 square inches and the perimeter needed to be less than 20 inches.\n",
    "\n",
    "1. First off lets rewrite the equation into an optimization problem:\n",
    "<center>$\\text{minimize}\\;\\; f(x) = \\frac{55*1000}{wh}$<br>\n",
    " $\\;\\;\\;\\;\\;\\;\\;\\;\\;\\text{subject to} \\;\\; g(x) = wh-20 \\le 0 $<br>\n",
    " $\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;g(x) =\\;2*h+2*w-18\\le0$\n",
    " </center>\n",
    "<br><br>*** Note if you do not have sympy installed for jupyter notebook you need to run the following code in a cell:<br>\n",
    "!pip install --upgrade<br>\n",
    "!pip install sympy <br>\n",
    "import sympy<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The first order conditions are\n",
      "h*u1 + 2*u2 - 55000/(h*w**2)\n",
      "u1*w + 2*u2 - 55000/(h**2*w)\n",
      "u1*(h*w - 20) + u2*(2*h + 2*w - 18)\n",
      "width,height,u1 =\n",
      "Matrix([[4.47213595499958], [4.47213595499958], [137.500000000000]])\n",
      "The perimeter is\n",
      "17.88854382\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from sympy import *\n",
    "# initialize variables and functions\n",
    "w = Symbol('w')\n",
    "h = Symbol('h')\n",
    "u1 = Symbol('u1')\n",
    "u2 = Symbol('u2')\n",
    "f = 55*1000/(w*h)\n",
    "g1 = w*h-20\n",
    "g2 = 2*h+2*w-18\n",
    "# necessary first order conditions\n",
    "foc1 = sympy.diff(f,w)+ u1*sympy.diff(g1,w)+ u2*sympy.diff(g2,w)\n",
    "print 'The first order conditions are'\n",
    "print  foc1\n",
    "foc2 = sympy.diff(f,h)+ u1*sympy.diff(g1,h)+ u2*sympy.diff(g2,h)\n",
    "print foc2\n",
    "foc3 = u1*g1+u2*g2\n",
    "print foc3\n",
    "# Let's start by checking the case where g1 is active and g2 is inactive, this gives us the following equations\n",
    "foc1g1 =h*u1 -55*1000/(h*w**2)\n",
    "foc2g1 = u1*w-55*1000/(h**2*w)\n",
    "foc3g1 =h*w -20\n",
    "ans = nsolve([foc1g1,foc2g1,foc3g1],[w, h, u1], [4, 4, 1])\n",
    "print 'width,height,u1 ='\n",
    "print ans\n",
    "perimeter = 4*4.47213595499958\n",
    "print 'The perimeter is' \n",
    "print perimeter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Since $\\mu_1$ is greater than zero and $\\mu_2$ is zero, we know we found a solution. Checking the solution, the answers we got for the width and height are multiplied together to get 20. When the perimeter is calculated, the result is 17.89 inches, which meets the requirement of being less than 18 inches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}