{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Optimization\n", "\n", "\n", "## Overview, Objectives, and Key Terms\n", " \n", "In [Root Finding](#root_finding), the bisection method and Newton's method were used to solve nonlinear equations. Now, those methods (and tools like `fsolve`) are used to solve **optimization problems**. The subject is rich, so we'll touch on only the basics, but you'll have tools at your disposal to tackle such problems in varied applications.\n", " \n", "### Objectives\n", "\n", "By the end of this lesson, you should be able to\n", "\n", "- Find local extrema of a function $f(x)$ using the bisection and Newton methods.\n", "- Use `minimize` to solve nonlinear optimization problems.\n", "\n", "### Key Terms\n", "\n", "- optimization\n", "- uncontrained optimization\n", "- constrained optimization\n", "- extremum\n", "- critical point\n", "- objective function\n", "- linear program\n", "- `scipy.optimize.minimize`\n", "- `scipy.optimize.linprog`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Some Nomenclature\n", "\n", "Optimization problems often defined by minimizing an **objective function**\n", "\n", "$$\n", " \\min_{x} f(x) \\qquad \\text{or} \\qquad \\min_{\\mathbf{x}} f(\\mathbf{x}) \\, ,\n", "$$\n", "\n", "subject to the **inequality constraints**\n", "\n", "$$\n", " g(x) \\ge 0 \\qquad \\text{or} \\qquad \\mathbf{g(x)} \\ge \\mathbf{0} \\, ,\n", "$$\n", "\n", "and/or the **equality constraints**\n", "\n", "$$\n", " h(x) = 0 \\qquad \\text{or} \\qquad \\mathbf{h(x)} = \\mathbf{0} \\, .\n", "$$\n", "\n", "An optimization problem is *unconstrained* if it is subject to no inequality constraints and no equality constraints. Otherwise, the problem is *constrained*. In practice, unconstrained problems are often easier to solve using a variety of techniques for finding *extrema*. We'll start with unconstrained problems and then consider an important class of constrained problems known as *linear programs*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Finding Extrema\n", "\n", "One of the great applications of differentiation is its ability\n", "to identify the *critical points* of a function $f(x)$, which\n", "include its *minima* and *maxima*, i.e., its *extrema*. Suppose\n", "$f(x)$ represents some quantity of interest, perhaps the cost \n", "of materials for some component or the signal-to-noise ratio of\n", "some sensor. If we want to mininimize those costs or maximize those\n", "ratios, we need to use *optimization*. Although optimization problems\n", "and the techniques to solve them are quite diverse, we will stick to\n", "1-D problems in which the *objective function* to be minimized or \n", "maximized is continuous (as opposed to discrete).\n", "\n", "You've had calculus (and we've reviewed some its topics). Hence, you might recall that a function $f(x)$ exhibits an **extremum**, i.e., a **minimum** or **maximum**, at any point $x$ for which $f'(x)$ is zero and $f''(x)$ is nonzero. When $f'(x)$ *and* $f''(x)$ vanish, the point is often called a **saddle point**, though that term is more meaningful in multiple dimensions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "**Exercise**: Given some fencing of *fixed* length `L`, what is the largest area that fencing can enclose if arranged in a rectangle of sides `a` and `b`? \n", "\n", "*Solution*:\n", "\n", "Here's the problem. We want to maximize $A = ab$ where $2a+2b = L$ or $b = L/2 - a$. We know $L$ and $b$ depends on $a$, so there is just one free parameter (i.e., we're in 1-D where we'll stay throughout). Note that $a$ (and $b$) can range from 0 to $L/2$.\n", "\n", "Formally, we see a solution to\n", "\n", "$$\n", " \\max_{a \\in [0,\\, L/2]} a\\left (\\frac{L}{2}-a \\right)\n", "$$\n", "\n", "which reads \"find the value of $a$ in the range $[0, L/2]$ that maximizes the quantity $f(a) = a(L/2-a)$. Here, $f(a)$ is the **objective function**.\n", "\n", "Differentiation of $f(a)$ gives\n", "\n", "$$\n", " f'(a) = L/2 - 2a \n", "$$\n", "\n", "which, when set to zero, requires that $a = L/4$. Consequently, $b = L/2-a = L/4 = a$. That is to say, our fence must form a square to maximize the area enclosed.\n", "\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Connection To Nonlinear Functions\n", "\n", "\n", "As the fence exercise illustrates, if we can write down the derivative and explicitly find the point at which it vanishes, our problem is solved. That is rarely the case. Often, the functions of interest have complicated derivatives that make $f'(x) = 0$ a nonlinear equation. We generally can't solve that problem directly. In many cases, we don't even have $f'(x)$, and at best, we can approximate it numerically using finite differences. \n", "\n", "Further complicating matters is that even if we can solve $f'(x)=0$, the result might not be an extremum (consider $f(x)=x^3$ and the zeros of its derivative $3x^2$), and it might not be unique. \n", "\n", "Consider the following example (for $f(x) = \\sin(x^22)-3)^2$): " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8FPX9+PHXexOgXF/lDHKWS6sWOQKJCFaCYBQU6sGhFW2/asRq7a9aFWurCB6AB4pSUCqirRLgy6lSEDGIIOUsiQQSwIDKIUdBMGCAZN+/P3ZDNyGBZLO7s8f7+XjMIzszn5l5f3aTfWc+85nPiKpijDHGFHM5HYAxxpjwYonBGGNMCZYYjDHGlGCJwRhjTAmWGIwxxpRgicEYY0wJlhhM1BKRliKSLyJxTsdiTCSxxGDCgogsFpFRZSwfKCLfiUh8Zfepqt+oah1VLQpMlMEjIj8VEfWnnsYEmiUGEy6mAcNEREotHwa8p6qFldlZNH7BRmOdTHiyxGDCxTygPnBl8QIRqQdcD7zrne8vIv8WkaMi8q2IjPQpW/wf910i8g3waen/wkWkqYgsEJFDIrJdRO7x2X6aiDzjM99LRHb5zD8mIrtF5AcRyRWRq8uqhIjUEJEXReQbEdknIpNFpKbPPv7lE899IpItIj8Blnt38b23+au7iPxaRFaKyHgROQSM9G73vyKyRUQOe8+0WvkcX0XktyKyzRvraBFpKyKrvO/bTBGpXvz+isiHInLAu68PRaR55T86E20sMZiwoKo/AjOBO3wWDwZyVDXTO3/Mu/58oD9wn4j8stSurgIuBlLLOMx0YBfQFLgFeK68L3hfInIR8ADQTVXreve9s5ziY4ELgU5AO6AZ8KR33QvASeDPItIeeA64XVULgF94y5zvbf5a5Z1PBvKAxsCz3vr+CbgJaAR87q2Xr2uBROBy4FHgTeBXQAvg58Ct3nIu4G2gFdAS+BF4/Vzvh4kBqmqTTWExAT2BI0BN7/xK4A9nKf8KMN77+qeAAm181hcvi8fzpVgE1PVZ/zwwzft6GvCMz7pewC7v63bAfqAPUO0s8Qie5NXWZ1l3YEepmA4BW4DHy4rVZ9mvgW9KHeOfwF0+8y7gONDKO69AD5/164HHfOZfAl4pJ/5OwGGnfw9scn6yMwYTNlR1BXAAGCgibYBuwPvF60UkWUQyvE0fR4DhQMNSu/m2nN03BQ6p6g8+y77G8x/9ueLaDvw/PE05+0UkXUSallG0EVALWC8i34vI98Ai7/Life0EMvAkgonnOjZn1qcV8KrP/g/hSUi+9djn8/rHMubrAIhILRF5Q0S+FpGjeJqzzrdeXMYSgwk37+JpLhoGfKyqvl9q7wMLgBaqeh4wGc+Xoq/yhgveA9QXkbo+y1oCu72vj+H5Ui/WpMROVd9X1Z54vpgVT5NRaQfxfPFeqqrne6fzVLVOcQER6YfnLGIpnqalc8Vdevm3wL0++z9fVWuq6hflbH82DwMXAcmq+j/8tzmr9HtqYowlBhNu3sXTZHMP8E6pdXXx/NdfICJJwG0V3amqfgt8ATwvIj8RkcuAu4D3vEU2Av1EpL6INMFzhgB4rjGISG8RqQEU4PnyP6MLrKq6gSnAeBFp7N22mYikel83BN4C7gbuBG7wJgrwnCm5gTbnqMpk4HERudS7z/NEZFBF34dS6nrr8r2I1Aee8nM/JspYYjBhxdvU8gVQG8/Zga/fAqNE5Ac8F3RnVnL3t+JpwtkDzAWeUtUl3nV/BzLxXFT+GJjhs10NYAyeM4Lv8FwI/lM5x3gM2A78y9s88wme/8rBcxF4vqouVNX/4ElMfxORBqp6HHgWWOltJrq8rJ2r6lw8Zyvp3v1vAq6rzJvg4xWgprde/8LT7GUMomoP6jHGGPNfdsZgjDGmhIAkBhGZKiL7RWSTz7L6IrLEe6PNEu/NSmVte6e3zDYRuTMQ8RhjjPFfoM4YpuG5qcbXCGCpqrbH0wNjROmNfC54JQNJwFPlJRBjjDGhEZDEoKrL8fSn9jWQ//YqeQcofYcqeO4gXaKqh1T1MLCEMxOMMcaYEArmoFwJqroXQFX3FnffK6UZJW/g2UU5NxyJSBqQBlCzZs3EFi1a+BWU2+3G5apcPjx8+DAHDhygbdu2xMWFz70//tQlHEVLPcDqEq4CWZcffviBvXv30qpVK2rUqBGQfVZUVeuxdevWg6ra6JwFA3ULNZ5ugJt85r8vtf6MW+2BR4A/+8z/BXj4XMdKTExUf2VkZFR6m6VLlyqgH3/8sd/HDQZ/6hKOoqUeqlaXcBXIuowaNUoBPX78eMD2WVFVrQewTh0eEmOfiFwA4P25v4wyu/CMYVOsOZ4+5mGlc+fOAGzYsMHhSIwxTtu+fTvNmzenZs2aTocSNMFMDAvw3N2J9+f8MsosBq7xDv9bD7jGuyys1KtXj5/+9KeWGIwxbN++nXbt2jkdRlAFqrvqdGAVcJGI7BKRu/DcKdpXRLYBfb3ziEhXEfkbgKoeAkYDa73TKO+ysNOlSxdLDMaYmEgMAbn4rKq3lrPqjLHuVXUdnrFiiuenAlMDEUcwde3alTlz5nD48GHq1bMetcbEoqNHj7J///6oTwzR0eUgBJKSkgBYu3atw5EYY5zy1VdfAVhiMB5du3ZFRFizZo3ToRhjHLJ9+3bAEoPxOu+88/jZz35micGYGFacGNq2betwJMFliaESkpKSWLNmTfE9F8aYGPPVV1/RpEkT6tSpc+7CEcwSQyUkJSWxb98+vvnmG6dDMcY4YPv27VF/tgCWGCql+AK0NScZE5u2bt0a9dcXwBJDpVx22WVUr17dEoMxMejo0aPs3buXiy666NyFI5wlhkqoXr06nTt3tsRgTAzaunUrAD/72c8cjiT4LDFUUnJyMuvWraOwsNDpUIwxIZSbmwtgZwzmTElJSRw/fpzNmzc7HYoxJoRyc3NxuVx28dmcyS5AGxObcnNzad26dcifweAESwyV1K5dO+rVq8fq1audDsUYE0K5ubkx0YwElhgqTUS4/PLL+eKLL5wOxRgTIm63m61bt1piMOXr2bMnmzdv5tChsBwh3BgTYLt27eLHH3+0xGDK16NHDwA7azAmRsRSjySwxOCXbt26ER8fz8qVK50OxRgTApYYzDnVqlWLxMREVqxY4XQoxpgQyM3NpW7dujRp0sTpUELCEoOfevTowdq1azlx4oTToRhjgqy4R5KIOB1KSAQ1MYjIRSKy0Wc6KiL/r1SZXiJyxKfMk8GMKVB69uzJiRMnWL9+vdOhGGOCLJa6qkKAnvlcHlXNBToBiEgcsBuYW0bRz1X1+mDGEmhXXHEFACtXrjz92hgTfY4dO8Y333wTU4khlE1JVwNfqerXITxm0CQkJNC+fXu7AG1MlNu2bRsQOxeeIbSJYSgwvZx13UUkU0T+KSKXhjCmKunRowcrV660J7oZE8VirUcSBLkpqZiIVAcGAI+XsXoD0EpV80WkHzAPaF/GPtKANPD8t75s2TK/YsnPz/d729IaNmzIwYMH+fvf/07Lli0Dss/KCGRdnBQt9QCrS7iqSl0WLVoEwN69ezl8+HAAo6q8kH0mqhr0CRgIfFzBsjuBhmcrk5iYqP7KyMjwe9vScnJyFNApU6YEbJ+VEci6OCla6qFqdQlXVanLrbfeqq1atQpYLFVR1c8EWKcV+B4OVVPSrZTTjCQiTcTbB0xEkvA0b/0nRHFVyYUXXkiTJk3IyMhwOhRjTJBs3ryZSy+NmBbugAh6YhCRWkBfYI7PsuEiMtw7ewuwSUQygQnAUG9mC3siQq9evcjIyLDrDMZEoaKiInJycrjkkkucDiWkgp4YVPW4qjZQ1SM+yyar6mTv69dV9VJV7aiql6tqRA1AlJKSwt69e08/9s8YEz3y8vI4ceKEnTGYyklJSQGw5iRjolB2djaAnTGYymnXrh3NmjWzxGBMFCp+hO/FF1/scCShZYmhikSElJQUli1bZtcZjIky2dnZtGzZkrp16zodSkhZYgiAlJQU9u/ff/q/C2NMdIjFHklgiSEg7DqDMdEnVnskgSWGgGjdujWtWrWyxGBMFNmxYwcFBQV2xmD817t3bzIyMigqKnI6FGNMAMRqjySwxBAwqampHD58mLVr1zodijEmAIqvGVpiMH7r06cPIsLixYudDsUYEwDZ2dm0aNEi5nokgSWGgGnQoAHdunWzxGBMlIjVHklgiSGgUlNTWb16teND8xpjqqaoqIgtW7bEZDMSWGIIqNTUVNxuN5988onToRhjqmDnzp0x2yMJLDEEVHJyMuedd541JxkT4WK5RxJYYgio+Ph4rr76ahYvXmzDYxgTwb788ksAO2MwgZGamsquXbvYsmWL06EYY/yUlZVFmzZtYrJHElhiCLjrrrsOgA8//NDhSIwx/srMzKRjx45Oh+EYSwwB1qJFCzp37syCBQucDsUY44fjx4+zbds2LrvsMqdDcYwlhiAYMGAAq1at4sCBA06HYoyppM2bN+N2uy0xmMAaMGAAbrebhQsXOh2KMaaSMjMzASwxBJOI7BSRL0Vko4isK2O9iMgEEdkuIlki0iXYMQVb586dadasmTUnGROBsrKyqF27Nm3atHE6FMeE6owhRVU7qWrXMtZdB7T3TmnApBDFFDQiwoABA1i8eDEFBQVOh2OMqYSsrCw6dOiAyxW7DSrhUPOBwLvq8S/gfBG5wOmgqmrAgAEcO3bMntFgTARRVbKysmK6GQkgPgTHUOBjEVHgDVV9s9T6ZsC3PvO7vMv2+hYSkTQ8ZxQkJCSwbNkyv4LJz8/3e9vKcLlc1KxZk0mTJlGzZs2gHCNUdQm2aKkHWF3CVUXrcuDAAQ4dOkTNmjXDsu4h+0xUNagT0NT7szGQCfyi1PqPgJ4+80uBxLPtMzExUf2VkZHh97aVNWjQIG3cuLEWFhYGZf+hrEswRUs9VK0u4aqidfnoo48U0M8//zy4Afmpqp8JsE4r8L0d9KYkVd3j/bkfmAsklSqyC2jhM98c2BPsuEJh0KBB7N+/n+XLlzsdijGmArKysgDo0KGDw5E4K6iJQURqi0jd4tfANcCmUsUWAHd4eyddDhxR1b1EgX79+lGrVi1mzpzpdCjGmArIysqiVatWnHfeeU6H4qhgnzEkACtEJBNYA3ykqotEZLiIDPeWWQjkAduBKcBvgxxTyNSuXZvrr7+eOXPmUFhY6HQ4xphz2LhxY0wPhVEsqIlBVfNUtaN3ulRVn/Uun6yqk72vVVXvV9W2qtpBVc+41yGSWXOSMZHh2LFj5OTk0KVLxN9KVWXh0F01qllzkjGRITMzE1W1xIAlhqCrVauWNScZEwHWr18PYIkBSwwhMWTIEA4cOGCP/DQmjG3YsIGEhASaNm3qdCiOs8QQAv3796devXq8++67TodijCnHhg0b6NKlCyLidCiOs8QQAjVq1GDo0KHMnTuXI0eOOB2OMaaUgoICsrOzrRnJyxJDiNx5550UFBQwa9Ysp0MxxpTy5ZdfUlRUZInByxJDiCQlJXHRRRdZc5IxYWjDhg2AXXguZokhRESEO++8k88//5y8vDynwzHG+NiwYQP16tWjVatWTocSFiwxhNCwYcMQEd555x2nQzHG+Fi/fr1dePZhiSGEmjdvTmpqKm+99Zbd02BMmDh58iRffvmlNSP5sMQQYsOHD2f37t188MEHTodijAGys7M5efKkJQYflhhCrH///jRv3pxJkyL+CabGRIU1a9YA0K1bN4cjCR+WGEIsPj6etLQ0lixZwrZt25wOx5iYt2bNGho2bEibNm2cDiVsWGJwwN133018fDyTJ092OhRjYt7q1atJSkqyC88+LDE44IILLuDGG2/k7bff5tixY06HY0zMOnr0KJs3byY5OdnpUMKKJQaHPPjggxw+fJi3337b6VCMiVnr1q1DVUlKKv3E4dhmicEhPXr0oHv37rz88svWddUYhxRfeLbEUJIlBoeICI888gg7duxgzpw5TodjTExavXo17du3p379+k6HElaClhhEpIWIZIjIFhHJFpHfl1Gml4gcEZGN3unJYMUTjgYMGMCFF17IuHHjUFWnwzEmpqjq6QvPpqRgnjEUAg+r6sXA5cD9InJJGeU+V9VO3mlUEOMJO3FxcTz88MOsX7+ejIwMp8MxJqbs3r2bvXv32oXnMgQtMajqXlXd4H39A7AFaBas40WqO+64gwsuuICRI0faWYMxIbR69WoASwxlkFB8GYnIT4HlwM9V9ajP8l7AbGAXsAf4o6pml7OPNCANICEhITE9Pd2vWPLz86lTp45f2wbL3LlzmTBhAi+88AJdu3at8HbhWBd/REs9wOoSrsqqy+TJk5kzZw4ffvgh1atXdyiyyqnqZ5KSkrJeVc/9JaOqQZ2AOsB64KYy1v0PUMf7uh+wrSL7TExMVH9lZGT4vW2wFBQUaMuWLTUpKUndbneFtwvHuvgjWuqhanUJV2XVpWfPnpqcnBz6YKqgqp8JsE4r8B0b1F5JIlINzxnBe6p6RtcbVT2qqvne1wuBaiLSMJgxhaMaNWrw5JNPsmbNGj766COnwzEm6hUUFLBmzRp69uzpdChhKZi9kgR4C9iiqi+XU6aJtxwikuSN5z/Biimc3XHHHbRt25YnnniCoqIip8MxJqqtX7+ekydPWmIoRzDPGHoAw4DePt1R+4nIcBEZ7i1zC7BJRDKBCcBQ7+lOzKlWrRrPPvssWVlZTJ061elwjIlqK1euBDw3mpozxQdrx6q6AjjrqFSq+jrwerBiiDSDBw9m4sSJPPHEEwwaNIjzzz/f6ZCMiUorVqzgwgsvpFGjRk6HEpbszucwIiK8+uqrHDx4kNGjRzsdTkTKP5nP4FmDyT+Z73QoJky53W5WrlxpzUhnYYkhzHTu3Jm7776bCRMmkJWV5XQ4EWdp3lJmbZ7Fpzs+dToUE6ZycnI4dOiQJYazsMQQhp5//nnq1avHXXfdZQPsVdLcnLmen1vmOhyJCVfF1xcsMZTPEkMYatCgAa+99hrr1q3j1VdfdTqciKGqfLj1QwA+2PqB3UluyrRixQoaN25Mu3btnA4lbFliCFODBw9mwIAB/PnPfyYnJ8fpcCLC5gObKSgsAODHwh/ZcnCLwxGZcLRixQp69OhhT2w7C0sMYUpEmDRpErVr12bo0KEUFBQ4HVLYW7htIYVuT9ObW90s3LbQ4YhMuNmzZw95eXnWjHQOlhjCWNOmTZk2bRqZmZk8+uijTocT9mZmz+RE0QkACgoLmJk90+GITLgpHsW4V69ezgYS5oJ2H4MJjOuvv57f//73vPrqq/ziF7/glltucTokx9w842bm5JT/UKPqcSUHQsvcl4k8XX5zwU0/u4nZQ2YHLD4T/jIyMqhXrx4dO3Z0OpSwZmcMEWDs2LF0796dO+64gw0bNjgdjmPG9BlDpyadqF2tdpnrTxadPOt8sdrVatO5SWfG9BkT8BhNePv000+56qqriIuLczqUsGaJIQLUqFGDuXPn0rBhQwYMGMCePXucDskR7Ru0Z90963i619PUjK+JSyr36+sSFzXjazIqZRTr0tbRvkH7IEVqwtHOnTvZsWMHvXv3djqUsGdNSREiISGBBQsW0LNnT6655ho+++wzp0MKmePHj5OVlcXGjRvJzs4mLy+Ppgeaktc1D+oBFRlK/yRUO1aN5F3J7P1+L7O+mUW3bt1o3bq19U6JEcXXFywxnJslhgjSqVMnFixYQL9+/UhNTWXkyJFV2t+XX8KkSZCRAXl5cPIkVK8ObdpASgrcdx906BCY2CtKVfn6669Zvnw5n332GV988QVbt27F7XYDULduXdq0aUOH1h24/ifXk1s7l6WnlnKKU+XuM554kjWZRt80Iu+rPCZ8MIGTJz3NTI0aNaJv375cc801pKam0qRJk5DU04ReRkYGjRs35pJLynrCsPFliSHC9O7dm9mzZ/PLX/6Shx9+mKSkJBo3blypfeTlwbBhsHEjnDgBvqN8nzwJOTmwbRu88w506gR//7snWQTL8ePHmT17NgsWLCAjI4Nvv/0WgHr16tGjRw+GDBlCp06d6NSpE61atSrxH/78nPl8Me8LTp0oPzHUqlGLR4c+yoCLBnjreJJNmzaxdu1aPv/8c5YsWcL777+PiNCrVy+GDh3KzTffTIMGDYJXaRNSqsqnn35KSkqKnSFWgCWGCNS/f3/mzZvHzTffTI8ePVi0aBFt27at0LYzZ8JvfnNmQiitqAiOH4fVqz1nDW+/DYMHB6gCeB7E/sEHHzB//nyWLl3KqVOnqF+/Pr179+bRRx/lqquu4tJLL8XlOvt1hLk5c/nhxA9nLfPDiR+Yu2Xu6cRQvXp1unTpQpcuXbj33ntxu91kZmYyf/58pk+fzr333svvfvc7Bg8ezP33309ycrJ9mUS43bt3s3v3blJSUpwOJSLYxecI1b9/f1566SUOHTpEt27dWLjw3DdzzZwJv/615wu/os8CKk4Qv/61Z3t/qSpZWVmMHj2abt260bx5c+677z62bdvGL3/5S5YtW8a+ffuYNWsWDzzwAB06dDhnUigeAkP579AXxReYfS9MK3rWITJcLhedO3dm5MiR5OTksGHDBtLS0pg/fz7du3ena9euzJgxwx6gFMGKe/PZ9YWKscQQwS699FLWrFlDy5Yt6d+/PyNGjCj3Dum8PM+Zwo8/+nesH3/0bL9jR8W3OXXqFEuXLuXBBx+kdevWdOzYkaeeeor4+Hiee+45srOz2bZtG7/97W+56qqriI+v3Ans5gOb+bHwvxWqVa0WHRM6Mn/ofDomdCzRrbWiQ2SICJ07d+a1115j9+7d/PWvf+X48eMMHTqUDh06MH36dEsQEWjdunW0aNHCxkeqIEsMEa5t27asWrWKu+++m7Fjx9KpUyeWL19+RrlhwzzNR1Vx4gTcfvvZyxw+fJjp06dz22230ahRI/r06cOUKVO47LLLmDJlCnv27GHVqlU8/vjjXHLJJVVqolm4bSFF7qLTZwmjU0azLm0dfdv2Ze09a0t0ay1yF1V6iIy6dety3333sWnTJtLT03G5XNx222106tSJJUuW+B23Ca1Tp06xYcMGrr32WmsSrCBLDFGgZs2aTJkyhUWLFlFQUMBVV13FgAEDTj/PISvLc6G5qv/oFhXBv//t6c1UTFXJycnhxRdfpFevXjRq1IjbbruNTz75hJtuuol58+Zx8OBBFixYwN133x3QXj8zs2dyyn2KjgkdyRyeyUPdHzrdhBTniuPhKx4mc3gmlyVcxin3Kb+HyIiLi2PIkCFkZWUxc+ZMjh8/zjXXXMMNN9xAbm5uwOpjgmPVqlUcO3aM6667zulQIkbQE4OIXCsiuSKyXURGlLG+hojM8K5fLSI/DXZM0So1NZXs7GyeffZZli9fTseOHbn66qt55JHtnDgRmCGoT56El17K57333uOee+6hXbt2XHzxxTzyyCN8//33jBgxglWrVrF3716mTp3KwIEDqV277DuVq6pJnSa80PeFs96sVnxT3Lg+40ionVCl47lcLgYNGsTmzZsZO3Ysn332GR06dODpp58+3f3VhJ9FixYRFxdn1xcqQ1WDNgFxwFdAGzy3IWUCl5Qq81tgsvf1UGDGufabmJio/srIyPB723Bztrr85z//0eeee05btmypsFlBAzhlK6Dnn3++DhgwQCdNmqRff/11UOoRzr777ju99dZbFdBLL71UV61aFbF1KUu01KVz58562WWXOR1GQFT1MwHWaQW+u4N9xpAEbFfVPFU9CaQDA0uVGQi84339f8DVYg2BVVa/fn0ef/xx8vLyiI+/KKD7jou7iA0bNnDw4EHmz5/P8OHDadmyZUCPEQkSEhJ4//33+eCDDzhy5AhXXHEFkyZN4kRVL+aYgPnuu+/497//TXJystOhRJRg38fQDPjWZ34XUPoTOl1GVQtF5AjQADjoW0hE0oA08PxBLlu2zK+A8vPz/d423FS0LoWFvQJ63KIiF0eOHOHzzz8PyP4i/TOpU6cOkydP5s0332TmzJls2LCBv/zlLxGfLCP9cwFPMxLAz3/+84ivC4TwM6nIaYW/EzAI+JvP/DDgtVJlsoHmPvNfAQ3Otl9rSvKoaF2qVw9kM5Jnf07UIxI888wz2qBBA61Zs6a++eab6na7nQ7Jb9HwuQwZMkSbNGmin376qdOhBES0NCXtAlr4zDcHSg8NerqMiMQD5wGHghxXTAn0cBYVvMk6JvXo0YOsrCx69OhBWloat99+O8ePH3c6rJhUWFjIxx9/TGpqqnVTraRgJ4a1QHsRaS0i1fFcXF5QqswC4E7v61uAT72ZzQRISgoEavj5uDjP/kz5mjZtyuLFi3nmmWeYPn063bt3Jy8vz+mwYs6KFSs4fPgwN9xwg9OhRJygJgZVLQQeABYDW4CZqpotIqNEZIC32FtAAxHZDjwEnNGl1VTN8OFQo0Zg9lW9umd/5uxcLhdPPPEECxcu5NtvvyUxMZF//vOfTocVU+bNm0eNGjVITU11OpSIE/T7GFR1oapeqKptVfVZ77InVXWB93WBqg5S1XaqmqSq9q9VgF12mWeU1KqeNcTFQefOoR+KO5Jde+21rFu3jlatWtG/f39efPHFcsdsMoGjqsyfP5++fftSp04dp8OJOHbnc4z4+9+rftZQowb84x+BiSeWtGnThi+++IJbbrmFRx55hPvuu4/CwkKnw4pqWVlZ7Ny5k4EDS/eONxVhiSFGtGnjGTq7Zk3/tq9Z07N969aBjStW1KpVi/T0dEaMGMEbb7zB9ddfz9GjR50OK2rNmzcPEbHrC36yxBBDBg+GadOgVq2KNyvFxXnKT5sW2OcxxCKXy8Xzzz/PlClT+OSTT+jZs+fphxKZwJo3bx7du3cnIaFqw6DEKksMMWbwYM8geMnJnrOA8hJEXJxnfXIybNpkSSGQ7r77bv75z3/y9ddfc8UVV7Bly7mHAzcVt2PHDjZu3GjNSFVgiSEGtWkDK1d6ns52771w8cWe3kYinp8XX+xZvnq1p5w1HwVe3759Wb58OadOneLKK69k7dq1TocUNWbMmAHAoEFuCk1RAAAQQElEQVSDHI4kctmjPWNYhw4wcaLTUcSujh07snLlSvr27UtKSgrz58/n6quvdjqsiDdjxgySk5Npbf/R+M3OGIxxUNu2bVm5ciVt2rShX79+zJ492+mQIlpubi4bN25k6NChTocS0SwxGOOwCy64gM8++4yuXbsyePBg3n77badDilgzZsxARKwZqYosMRgTBurVq8eSJUvo06cP//u//8ubb77pdEgRR1VJT0/nyiuvpFmzZk6HE9EsMRgTJmrVqsX8+fPp378/9957LxPtAlClfPnll2zZsoUhQ4Y4HUrEs8RgTBj5yU9+wuzZsxk4cCAPPPAA48ePdzqkiDFt2jSqVavGYOtbXWWWGIwJMzVq1GDWrFncfPPNPPTQQ4wbN87pkMLeqVOn+Mc//sENN9xAw4YNnQ4n4lliMCYMVatWjenTpzNkyBAee+wxnnnmGadDCmsfffQRBw4c4De/+Y3ToUQFu4/BmDBVrVo1/vGPf1CtWjX+8pe/UFRUxFNPPeV0WGFp2rRpNGnShGuvvdbpUKKCJQZjwlh8fDzTpk0jPj6ekSNH4na7GTlypD2RzMf+/fv56KOP+MMf/kB8vH2lBYK9i8aEubi4OP72t78hIowaNQpV5emnn7bk4DVt2jQKCwutGSmALDEYEwGKk4PL5WL06NG43W5Gjx4d88mhqKiISZMm0atXLy6++GKnw4kalhiMiRAul4s333wTEeHZZ59FVXnmmWdiOjksXLiQnTt38sILLzgdSlQJSmIQkReAG4CTwFfAb1T1+zLK7QR+AIqAQlXtGox4jIkWLpeLN954A5fLxXPPPYfb7ea5556L2eQwceJEmjVrZkNsB1iwzhiWAI+raqGIjAUeBx4rp2yKqh4MUhzGRB2Xy8WkSZNwuVyMGTMGt9vNmDFjYi45bN26lcWLFzNq1CiqVavmdDhRJSiJQVU/9pn9F3BLMI5jTKxyuVxMnDgREWHcuHG43W7GjRsXU8lhwoQJVKtWjXvuucfpUKKOqGpwDyDyATBDVc94jLyI7AAOAwq8oarljhwmImlAGkBCQkJienq6X/Hk5+dTp04dv7YNN9FSl2ipB4S+LqrKhAkTmDdvHoMHD2b48OEBSw7h/LkcPnyYoUOHcvXVV/Poo4+es3w416UyqlqPlJSU9RVqsldVvybgE2BTGdNAnzJPAHPxJqAy9tHU+7MxkAn8oiLHTkxMVH9lZGT4vW24iZa6REs9VJ2pi9vt1t/97ncK6B/+8Ad1u90B2W84fy5/+tOfVEQ0JyenQuXDuS6VUdV6AOu0At+xfjclqWqfs60XkTuB64GrvQGVtY893p/7RWQukAQs9zcmY2KRiPDqq6/icrkYP348breb8ePHR22z0tGjR5k4cSI33XQTF110kdPhRKVg9Uq6Fs/F5qtU9Xg5ZWoDLlX9wfv6GmBUMOIxJtqJyOlk8Morr6CqvPLKK1GZHCZNmsSRI0cYMWKE06FErWD1SnodqAEs8f5i/ktVh4tIU+BvqtoPSADmetfHA++r6qIgxWNM1BMRXn75ZVwuFy+//DJut5sJEyZEVXI4cuQI48aNIzU1la5drXd7sASrV1K7cpbvAfp5X+cBHYNxfGNilYjw4osvIiK89NJLuN1uXn/99ahJDi+++CKHDh3iueeeczqUqGZ3PhsTZUSEF154AZfLxQsvvICq8vrrr+NyRfYo+/v27ePll19myJAhdOnSxelwopolBmOikIgwduxYXC4XY8eORVWZOHFiRCeH0aNHc+LECUaPHu10KFHPEoMxUUpEeP7553G5XDz//PO43e7Td0xHmszMTCZNmsS9995L+/btnQ4n6lliMCaKFQ+453K5ePbZZ3G73afHWooUbreb+++/n/r169uT7ELEEoMxUU5EGD169OkhuwsLC5kyZUrEPNTm3XffZeXKlUydOpX69es7HU5MiIzfDGNMlYgITz/9NPHx8Tz11FPs37+fGTNmhP0wEd999x1//OMf6d69O3feeafT4cSMyDmfNMZUiYjw5JNP8sYbb7Bo0SJ69erFvn37nA6rXKpKWloa+fn5vPXWWxHV/BXp7J02JsakpaUxf/58tmzZQvfu3cnNzXU6pDJNnTqVDz74gDFjxtjT2ULMEoMxMej6668nIyOD/Px8Lr/8chYvXux0SCVs3LiRBx54gN69e/Pggw86HU7MscRgTIxKSkpi9erVtGzZkn79+p2+3+GcmjQBkcBNTZqU2P2hQ4e46aabaNCgAdOnT7cmJAfYO25MDGvdujVffPEFt9xyCyNGjGDo0KEcO3bs7BsF+rqEz/4KCgq4+eab2b17N7Nnz6Zx48aBPZapEOuVZEyMq127Nunp6XTp0oXHH3+czMxM3n///ZDHUVRUxK9+9SuWLVvGe++9R3JycshjMB52xmCMQUR47LHHWLp06enrDunp6bjd7pAc/+TJk9x+++3MmTOHV155hdtuuy0kxzVls8RgjDktJSWFrKwsBgwYwBtvvMGVV17Jpk2bgn7cG2+8kfT0dMaMGcPvf//7oB/PnJ0lBmNMCfXr12fWrFmMGDGC3NxcOnfuzIgRIzh69GjQjrlo0SImTZrEY489FrRjmIqzxGCMOYOIkJqaSk5ODsOGDWPs2LG0bduW8ePHB+V4ixcvZvjw4UHZt6k8SwzGmHI1bNiQqVOnsmbNGjp16sRDDz0UlOP06XPWR8ibELPEYIw5p27durFkyRIyMjKcDsWEQNASg4iMFJHdIrLRO/Urp9y1IpIrIttFxJ7ubUwY69Wrl9MhmBAI9n0M41X1xfJWikgcMBHoC+wC1orIAlXdHOS4jDHGlMPppqQkYLuq5qnqSSAdGOhwTMYYE9OkQmOj+LNjkZHAr4GjwDrgYVU9XKrMLcC1qnq3d34YkKyqD5SxvzQgDSAhISExPT3dr7jy8/PDfgz6ioqWukRLPSA26tIrJSXgx1oW5GsX0fK5VLUeKSkp61W16zkLqqrfE/AJsKmMaSCQAMThOSt5FphaxvaDgL/5zA8DXjvXcRMTE9VfGRkZfm8bbqKlLtFSD9UYqQsEfnKqLhGmqvUA1mkFvturdI1BVSvUx0xEpgAflrFqF9DCZ745sKcqMRljjKmaYPZKusBn9kY8ZxKlrQXai0hrEakODAUWBCsmY4wx5xbMXknjRKQToMBO4F4AEWmKp/mon6oWisgDwGI8zU5TVTU7iDEZY4w5h6AlBlUdVs7yPUA/n/mFwMJgxWGMMaZynO6uaowxJsxYYjDGGFOCJQZjTOUkJIT3/kyV2aM9jTGV8913TkdggszOGIwxxpRgicEYY0wJlhiMMcaUYInBGGNMCZYYjDHGlGCJwRhjTAmWGIwxxpRgicEYY0wJlhiMMcaUYInBGGNMCZYYjDHGlGCJwRhjTAmWGIwxxpRgicEYY0wJQRl2W0RmABd5Z88HvlfVTmWU2wn8ABQBharaNRjxGGOMqbigJAZVHVL8WkReAo6cpXiKqh4MRhzGGGMqL6gP6hERAQYDvYN5HGOMMYET7GsMVwL7VHVbOesV+FhE1otIWpBjMcYYUwGiqv5tKPIJ0KSMVU+o6nxvmUnAdlV9qZx9NFXVPSLSGFgC/E5Vl5dTNg1IA0hISEhMT0/3K+78/Hzq1Knj17bhJlrqEi31AKtLuIqWulS1HikpKesrdC1XVYMy4Wmm2gc0r2D5kcAfK1I2MTFR/ZWRkeH3tuEmWuoSLfVQtbqEq2ipS1XrAazTCnzHBrMpqQ+Qo6q7ylopIrVFpG7xa+AaYFMQ4zHGGFMBwUwMQ4HpvgtEpKmILPTOJgArRCQTWAN8pKqLghiPMcaYCgharyRV/XUZy/YA/byv84COwTq+McYY/9idz8YYY0qwxGCMMaYESwzGGGNKsMRgjDGmBEsMxhhjSrDEYIwxpgRLDMYYY0qwxGCMMaYESwzGGGNKsMRgjDGmBEsMxhhjSrDEYIwxpgRLDMYYY0qwxGCMMaYESwzGGGNKsMRgjDGmBEsMxhhjSrDEYIwxpgRLDMYYY0qoUmIQkUEiki0ibhHpWmrd4yKyXURyRSS1nO1bi8hqEdkmIjNEpHpV4jHGGFN1VT1j2ATcBCz3XSgilwBDgUuBa4G/ikhcGduPBcaranvgMHBXFeMxxhhTRVVKDKq6RVVzy1g1EEhX1ROqugPYDiT5FhARAXoD/+dd9A7wy6rEY4wxpurig7TfZsC/fOZ3eZf5agB8r6qFZylzmoikAWne2XwRKSshVURD4KCf24abaKlLtNQDrC7hKlrqUtV6tKpIoXMmBhH5BGhSxqonVHV+eZuVsUz9KPPfFapvAm+Wt76iRGSdqnY9d8nwFy11iZZ6gNUlXEVLXUJVj3MmBlXt48d+dwEtfOabA3tKlTkInC8i8d6zhrLKGGOMCbFgdVddAAwVkRoi0hpoD6zxLaCqCmQAt3gX3QmUdwZijDEmRKraXfVGEdkFdAc+EpHFAKqaDcwENgOLgPtVtci7zUIRaerdxWPAQyKyHc81h7eqEk8FVbk5KoxES12ipR5gdQlX0VKXkNRDPP+4G2OMMR5257MxxpgSLDEYY4wpISYSQ3lDd4jIT0XkRxHZ6J0mOxlnRVR1GJJwJSIjRWS3z2fRz+mYKktErvW+99tFZITT8fhLRHaKyJfez2Gd0/FUhohMFZH9IrLJZ1l9EVniHXpniYjUczLGiiqnLiH5O4mJxEA5Q3d4faWqnbzT8BDH5Y+qDkMSzsb7fBYLnQ6mMrzv9UTgOuAS4FbvZxKpUryfQ6T1/Z+G5/ff1whgqXfonaXe+UgwjTPrAiH4O4mJxHCWoTsiTlWGITFBlQRsV9U8VT0JpOP5TEwIqepy4FCpxQPxDLkDETT0Tjl1CYmYSAzn0FpE/i0in4nIlU4HUwXNgG995s86xEiYekBEsryn0BFxuu8jGt7/Ygp8LCLrvUPRRLoEVd0L4P3Z2OF4qirofydRkxhE5BMR2VTGdLb/2vYCLVW1M/AQ8L6I/E9oIi6fn3Wp1BAjTjhHvSYBbYFOeD6XlxwNtvLC/v2vhB6q2gVPs9j9IvILpwMyp4Xk7yRYg+iFnD9Dd6jqCeCE9/V6EfkKuBBw9IJbEIchcVRF6yUiU4APgxxOoIX9+19RqrrH+3O/iMzF00xW1vW5SLFPRC5Q1b0icgGw3+mA/KWq+4pfB/PvJGrOGPwhIo2KL9CKSBs8Q3fkORuV3845DEk48/7BFrsRz0X2SLIWaO99+FR1PB0BFjgcU6WJSG0RqVv8GriGyPssSluAZ8gdiPChd0L1dxI1ZwxnIyI3Aq8BjfAM3bFRVVOBXwCjRKQQKAKGq6ojF3sqqry6qGq2iBQPQ1KIzzAkEWKciHTC0/yyE7jX2XAqR1ULReQBYDEQB0z1Dg0TaRKAuSICnu+H91V1kbMhVZyITAd6AQ29w/U8BYwBZorIXcA3wCDnIqy4curSKxR/JzYkhjHGmBJiuinJGGPMmSwxGGOMKcESgzHGmBIsMRhjjCnBEoMxxpgSLDEYY4wpwRKDMcaYEv4/fAA1ZTY++WQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib\n", "%matplotlib inline\n", "import nonlinear_plots\n", "nonlinear_plots.extrema()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of the three roots, two correspond to minima---but the left one (blue circle) is greater than the right one (red square), so is it really a minimum? *Yes*, if we make the distinction between *local* and *global* minima. For this problem, the right minimum is the *global* minimum, while the central root (green star) represents a local maximum ($f(x)$ is *unbounded*, i.e., there is no value $M < \\infty$ such that $|f(x)| < M$ for all possible values of $x$). Given the choice, a global optimum is usualy preferred, but for many cases, we're at best guaranteed a local minimum. There are techniques to increase our chances of finding a global optimum, but that's outside the present scope!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving Unconstrained Optimization Problems\n", "\n", "\n", "### When the Derivative is Available\n", "\n", "If we have $f(x)$ and can evaluate $f'(x)$, then either Newton's method or the secant method can be applied to $f'(x) = 0$. The same goes for multi-variable problems, where the gradient $\\nabla f(\\mathbf{x})$ replaces $f'(x)$, and the *hessian* matrix $\\mathbf{H}(\\mathbf{x})$ replaces $f''(x)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "**Exercise**: Find the value of $x$ that maximizes the function $f(x) = 2\\sin(x) - x^2/10$ starting with $x_0 = 2.5$. Then, confirm your result by showing it on a graph of $f(x)$.\n", "\n", "*Solution*. We can immediately differentiate $f(x)$ to find $f'(x) = 2\\cos(x) + x/5$ and $f''(x) = -2\\sin(x) + 1/5$. Then, Newton's method gives" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEWCAYAAAByqrw/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4FNXbxvHvSYFA6EV6r6EJSO9IFamG0A1IExBRAUUFBWki/OjYQKogEEoAQYL03kJHuoBIKNIJKaQ97x+z5I0YSNtkdjfnc117JbszO3Nve3b2zJk5SkTQNE3THIeT2QE0TdM069KFXdM0zcHowq5pmuZgdGHXNE1zMLqwa5qmORhd2DVN0xyMLuyAUupzpdRPZuewJqXUKKXUYrNzPKOUSqeU+lUp9UgptSKW6UopNV8p9UApdchyW1ql1BmlVO54LH+QUmpCHPM8UUoVTfyj0DT7oAs7ICLjRaS32TkcXHsgF5BdRLximV4HaALkF5Fqltv6ArtE5FY8lj8b6KaUeuVFM4hIBhG5nJDQSqk8Sql1SqkbSilRShV+ybwFLV8eMS+ilBoSY573lVJXlFKPlVL+Sqk6CckTz8yNlFLnlFLBSqntSqlCMaalVUrNs6z/llJqcBLWc1Up1dg6qV+6ntlKqfNKqSilVI943iebUuqOUmpPjNvKWJ7zB5bLFqVUmWTImyLP/8vowq6llELABRGJeMn0qyISFOO2d4Gf47NwEQkFNgLeSUr5X1GAH+AZjwzXLF8eGUQkA1Decv9VAEqp6sAEjC+5zMBcwFcp5ZzQUJaiWjiW23MAq4EvgGyAP7A8xiyjgBIYz3dD4BOlVPOErj+FnQAGAEcTcJ9vgLPP3XYD47nPBuQA1gHLEhPI5p9/EbHZC3AV+Bg4CQRhfBByYXyAA4EtQNYY868AbgGPgF1AWcvtaYDjwPuW687AXuBLy/VRwGLL/4UBAd4B/gYeAP2AqpYcD4FZMdYZfd/n7u9iub4DGAvsA54AvwLZgSXAY+AwUPgFj98PGPjcbSeAtyz/T7dkfAwcAerGlgtoAFyP5bltbPnfCfgU+BO4B/gA2SzT3IDFltsfWvLmekFeD8vjfQj8AbS23P4VEAaEW56DXs/drxcQCkRapn8FFARCYjyPL30NLbd1Bba/5P0kQHHL/y2AMxjvowBgaBzvRRfL/WN9rV5wn5Ex8wAdgUMxrrtblpnHcj0t8D/gGnAb+AFI95LPxn+yYPzK2ffcOkKA0pbrAUDTGNPHAMtesI5iwDbLa3/X8p7NYpn2M8aXVojlNfsklvsPAw7EeA37W94XbomsB3uAHvGYryawH+MzvOclr+d7QHCM22zq+U/KxaoLs3o448k7gFHM8wH/YHxrV7K8CNuAkTHm7wlktEybBhyPMa0cRpH2AIZblutsmTaK/xb2HzCKWlOMorMGeCVGjvrP3/e5+8cs7JcsH5LMGMXkAtDY8uZaBMx/weP3BvbGuF4Go2imtVzvhvEl4QIMwfhSc4vlMTXg5YX9Q8vzkd/y3P0ILLVMexfjyyg9RjF9DcgUS1ZXy+P8HKMIv45RNEvF9jzFcv8exPgQAm8Cfzw3zwtfQ8v0ysD9l6wjZmG/ieWLEMgKVI7jvZiYwv4nMQoRkAnjC7i65bl8HzgGKMv0aRhbkdkw3se/Al+/5LMRW2GZDnz/3G2nMX5xZLU8hlwxprUHTr1gHcUxmsfSAjkxNpamxfYeesH9nSz3GYWxlfoAqBRj+sOXXD6NZXlxFnbL83rU8j7913vqufVGYHwxjYhxu009/0m5uGD7ZorIbQCl1G7gHxE5ZrnuCzR6NqOIzHv2v1JqFPBAKZVZRB6JyGml1FjAF+OLopqIRL5kvWPE+Hn/u1IqCKPQ/RMjRyVgZzwfw3wR+dNy341AGRHZYrm+AuNbOza+wPdKqUIi8hfGFulqEXlqebwxd45OVkqNAEphbNUnxLsYvwyuWzKNAq4ppd7G2MrOjlEQT2IUptjUADIAE0QkCtimlFoPdMb4YCdUFowvhmjxeA0DMb484yMcKKOUOiEiDzCKjtUopepaMq58Lt8qjAKlMArMGyIiSikF9AEqiMh9yzLGA78AnyVg1RmAO8/d9gijUGWIcf35af8hIpcwvqwB7iilpmD8CokXEYlSSnljFNqOwMRnn13L9CzxXVYCDAIOisgRpVT5F+TKopRyB7oDf4Gx8x4be/6Twh7a2G/H+D8klusZAJRSzkqpCUqpP5VSjzG+UcFoS3tmIcYW9W8ictEa642nRC1LRAKBDUAny02dMH4OA6CUGqKUOmvpafIQo6jl+O+S4lQIo633oWU5ZzGaRXJh/OTeBCyz7ECcqJRyjWUZeYG/LUX9mb8wfuEkxgNif8O/7DXMyL8/NC/jidEc85dSaqdSqmYic75Id2CViDyJcVtvjF+VZTF+1XQD1iul8mJsEacHjsR4Hfwstz/bMfswxrSCwMkYt3WxrOMJxi+DmDJhfKk8iXH9+Wn/oZR6RSm1TCkVYPlMLSaB7y8RuQpsx3jNvk3IfRPK8jwOwvg1F1euIIxf5YssO9xt7vlPCnso7PHVBWiD0cSRGeONBMaW0TPfAeuBZlbsjRCE8YZ4Js6ueQm0FOhsKTzpMD4kz7YIhwEdMPYzZMEoaiqWZfwro2VnXc4Y0//G2HLMEuPiJiIBIhIuIl+JSBmgFtCS2HdQ3gAKKKVivqcKYrQpJsZJoKhS6vlflS97DT2I568VETksIm0wmtfWYOxXsAqlVDrAC+NLKKZXgV9F5IKIRImIH0aTUC2MNuwQjP1Cz16DzGLshEWMHbPRrw9GO3CFGLf9YlnHH5b1PMvijtEM+Ifll8nNmNMt///xgofyNUbTQQURyYTxRRTz/RXnqWGVUi0w2ry3ApOem/Z8D6KYl8/jWnYsqgF5gDNKqVsYzSLVLL1PYttB7YTxuciHbT7/ieZIhT0j8BRjR096YHzMiZZmhWftboOAhUqphGx1v8hxoJ7lGz0zCfvZFh+/YWxRjwaWx9gizojRTngHcFFKfcl/txSeuQC4KaXetGxtj8BoN33mB2Dcs25ZSqmcSqk2lv8bKqXKWz4YjzGaMGJrwjqI8QXyiVLKVSnVAGhFInsdWJqFLmJ8WLFkies1rI+xY/2llFJplFJdLc104ZbH9cJmOaWUG///fKW1XH+ZdhjNLNufu/0w8KZSqqgyNAFKAqctr+scYKplCxKlVD6lVLO4Hs9zfIFySilPS84vgZMics4yfREwQimVVSlVGqP5YcELlpURYyvzoVIqH0ZHhphuAy88LkAZPUTmYvxS6Q60shR6ILr76Ysu42MsJ43lsSjAVSnl9twGxDMbMTboKlouX2Lsw6goIpFKqSZKqUqWX/eZgCkYvwzP2ujzn2iOVNgXYfz0D8DYQXng2QSlVEGMHSPeIvLE8u3qD0xN6kpFZDNGd6Zn7c/rk7rM55b/FKP7VGOM9r5nNmG8kS9gPO5QjC3v2JbxCKO72E8Yz08QcD3GLNMxdhr9rpQKxHjuqlum5cZoJ36M0USzE+Mn+fPrCANaA29gbP18h/F8n3t+3gT4EXgb4n4NLR+iFvx3K/lF3gauWpoY+mFsjb7Is54fAOcs17Gs9wel1A/Pzd8dWCSWvWMxLML4otuB8XzOAN6N8RwNw2jTPmDJtQVjn0m8icgdjGamcRhFqzr/35QHRhv5nxjvmZ3AJMsvh9h8hbFD+hFGk+Dq56Z/jVGkHiqlhsZy/9nAWhH5TUTuYfR++kkplT0hjwn4HeM5r2VZZghQD8DyBf0HGJ8VEbn17GLJHS7/fxxEFoxfwI8sz0FxoLllXxrY3vOfaOq/7z1Nsw1KqbQYW1yNRORmHPO+DxQQkU9SJJym2TBd2DVN0xyMIzXFaJqmaejCrmma5nB0Ydc0TXMwphx5miNHDilcuLAZq9Y0TbNbR44cuSsiOeOaz5TCXrhwYfz9/c1YtaZpmt1SSv0Vn/l0U4ymaZqD0YVd0zTNwejCrmma5mDs4bS9mqbZmfDwcK5fv05oaGjcM2v/4ebmRv78+XF1je1EqnHThV3TNKu7fv06GTNmpHDhwigV2wlHtRcREe7du8f169cpUqRIopaR5KYYpVQBZQzYelYp9YdS6oOkLlOzsokTYfvzJxpMpO3bjeVp2kuEhoaSPXt2XdQTQSlF9uzZk/Rrxxpt7BHAEBHxwBhF5z2VDCN/a0lQtSp06JD04r59u7GcqlWtk0tzaLqoJ15Sn7skF3YRuSkiRy3/B2Kc2jWxo+bEafnp5aw5t4azd84SFhmWXKtxLA0bgo9P0or7s6Lu42MsT9M0m2XVNnalVGGMsUAPxjKtL8YI3hQsWDDR6xi2ZRh/PTL66DsrZ4pkLYJXGS/GNzLOyx8aEYqbS1zjIKRCMYt7QouzLuqaZles1t3RMpLNKuBDEXn8/HQRmS0iVUSkSs6ccR4R+0In+5/kUO9D/NzuZz6r8xmVclciQxpjEJ2wyDBemfQK1eZU49Mtn7Lp0iaCwoISvS6Hk5gtd13UNTs1Y8YMPDw86Nq1K2vWrGH06NEvnX/o0KFs27YthdIlMxFJ8gVwxRjRZ3B85n/ttdckOTwKfSRfbPtC6syrIy6jXYRRiOtoV5l7dG6yrM9ubdsmkiOH8dca82nac86cOWN2BClVqpRcvnxZRERq1qwpd+7ceen8V69elSZNmqREtHiJ7TkE/CUeNTbJTTHKaOWfizFu4JSkLi8pMqXNxOiGxrdyUFgQe//ey6ZLm6ia19jZt/PqTuYdn0e38t14vcjrODvFNr5tKhCfZhm9pa5ZyYcffsjx48etusyKFSsybdq0F07v168fly9fpnXr1nTr1o20adOSI0cOANq0aYOnpyfe3t78+OOP7Nq1iyVLllCoUCHu3bvHrVu3yJ3b2mPSpyxrNMXUxhg/8nWl1HHLpUVcd0pu7mncaVqsKZObTaZ8rvIAXH5wmTXn1tB0cVPyT83PiG0juPXkVhxLclAva5bRRV2zcz/88AN58+Zl+/btvPLKK1SuXDl62uzZsxk9ejS7d+9m8uTJzJw5M3pa5cqV2bt3rxmRrSrJW+wisgdj9HCb906ld+hcvjMbLmxg0clFjN89nkUnFnHlgyupc+s9ti13XdQ1K3vZlnVKuHnzJjH36+XKlYvRo0fTsGFDfH19yZYtW/S0V155hRs3bpgR06pS3ZGnbi5ueJbxxLOMJxfvXeTi/Ys4OzkTGRVJv/X9ePvVt6lbsK5D9sENCgri3Llz3Lp1i8ePH0dfsjdvTocWLdhdrhz1/viDdd268fDcOTLfvEnmzJnJkSMHpUuXJnPmzGY/BE1LsHTp0vHo0aN/3Xbq1CmyZ8/+nyIeGhpKunTpUjJeskh1hT2mEtlLUCJ7CQAu3LvAmvNr+OnYT9QuUJuJTSZSq0AtkxMmTmRkJCdOnODkyZOcOXOGP/74gzNnznD16tVY53dycuKuqyuf+PszVim+mDMn1vny589PuXLloi/ly5fn1Vdfxdk5Ff7a0eyGh4cHixcvjr5+6NAhNm7cyLFjx6hfvz5NmzaNPnT/woULeHl5mRXVapSxozVlValSRWxxoI2Q8BDmH5/P2F1jufnkJp4ensxpNYes6bKaHS1O165dY/Pmzfz+++9s3bqVe/fuAZAmTRpKly5NmTJlKFu2LGXKlCF//vxkzpyZTJkykSlTJtIfPIjq2BH690e+/56nixbxsFIlHj16xKNHj7h9+zZnzpzh9OnTnD59mrNnz/L06VMAsmbNStOmTXnjjTdo1qyZ3e900qzj7NmzeHh4mJrh2YA+6dOnp2rVqpw+fZqwsDCqVavG/PnzqVy5MuvWrWPq1Kls27aNiIgIKlSowKlTp3BxMX+bN7bnUCl1RESqxHnn+HSdsfYlubo7WsuTp09k9I7RUn1OdQmPDBcRif5rK6KiouTAgQMyaNAgKVWqlAACSJ48eaR79+6yePFiOX/+vISHx5H7+S6N8ejiGB4eLufOnZNffvlFevToIblz545ef6VKleTzzz+XEydOWPHRavbGFro7xjRo0CDZvHnzS+dZvXq1jBgxIoUSxS0p3R11YX+JyKhIERF5GPJQik4vKl/v/lqeRjw1NdP169dlwoQJUrp0aQHEzc1NmjdvLlOmTJHTp09LVFRU/Bf2oiKewP7rkZGRcuzYMRk/frzUrVtXnJ2dBZBq1arJTz/9JIGBgQl4hJojsLXCfuvWLVm7du1L5/Hx8ZEHDx6kUKK46cKezK4/ui6tl7YWRiHlvisnB/4+kKLrDw0NlWXLlknz5s3FyclJAKlTp47MnTtXHj16lLiFxlW8k3Bw0t27d2XatGlSpkwZASRjxozy7rvvir+/f+KyanbH1gq7PdKFPYWsO7dO8k/JL2qUkkG/DZKwiLBkXV9gYKBMmTJF8uXLJ4AUKFBARowYIRcvXkzaglPoyNOoqCjZu3evdO/eXdKlSyeA1KxZUzZt2pSwXxaa3dGFPel0YU9Bj0IfycANA6X10tbJVpzu3bsnX331lWTLlk0AqV+/vvz2228SERGR9IUntFhb6bQCDx48kBkzZkiBAgUEkNq1a8vWrVt1gXdQurAnnS7sJoiINIrsn/f/lB5resi94HtJXubNmzdl6NChkiFDBgGkVatWsm/fviQvN1pii7QVzxkTGhoq3377bfSvkPr168uOHTuSvFzNtujCnnS6sJtoyckl4jraVfJPyS87r+5M1DKCg4Nl7Nix4u7uLk5OTtKlSxc5efKkdYMmtThb+YRgISEhMmPGjOgeNU2aNJGzZ89aZdma+XRhTzpd2E3mH+AvxWcUF6evnOSLbV/Eu2tkVFSULFmyJLp5ol27dnLhwgXrB7RWUU6Gsz0GBwfLlClTJEuWLOLq6iqff/65BAUFWW35mjlsobBPnz5dSpcuLXnz5pWRI0dG3z516lRZuHDhS+/bsWPH5PksJoAu7Dbgcehj6e7bXRiFfL376zjn37dvn1SvXj267/f27duTL9w331ivGG/bZizPym7duiVvv/22AFK4cGFZv3691dehpRxbKOzPTts7f/786MIeHh4u5cuXj/P4jh07dkjv3r1TIOWLmXraXs2QMW1GFrRdQJtSbWharCkAj58+JlPaTP+a786dO3z00UcsWbKEPHnyMH/+fN5+++3kPSz/k0+st6yGDZPl5GC5cuVi0aJF9OzZkwEDBtCyZUvatWvH9OnTKVCggNXXp6WsBgsa/Oe2DmU7MKDqAILDg2mx5L8nhO1RsQc9KvbgbvBd2vu0/9e0HT12vHR9z5+2N0MGYzCebdu2UblyZVxcXIiIiKBmzZpMmjSJBg0a8Nlnn+Hk5MS4ceOoW7cuPXr0ICIiwiaOQk0oq42gpBnaebTDPY07weHBVP+pOgN/G0h4ZDgiwtKlSylTpgw+Pj6MGDGCCxcu0KNHD32ulRgaNGjA8ePHGT9+PH5+fnh4eDB37lzj56WmxVPM0/YOGzaMoUOHArB3715ee+01AFxcXFiwYAH9+/dn8+bN+Pn5MXLkSMA4f1Lx4sU5ceKEaY8hSeKzWW/tiyM2xTwvLCJMhm4aKoxCqv9QXZq0bSKAVK9eXU6fPm12PLtw5coVadiwoQDi6ekp9+4lveeRljJsoSmmUKFC/xk1qU+fPrJ06dJ/3TZu3DhJmzatHD169F+3d+nSRdatW5fsOV8kKU0xeos9mbg6uzKxyUR6Z+nNwb8PsqXoFj6c+CF79+6lbNmyZsezC4ULF2bz5s1MmDCBtWvX8uqrr7Jjxw6zY2l2LF26dISGhv7rtlOnTpElSxZu3779r9vt+RS+urAnk+vXr9OkSRN++vAnXjvxGnnz5OVIjiM4OemnPCGcnZ0ZNmwY+/fvJ126dLz++ut89tlnhIeHmx1Ns0MeHh5cunQp+vrq1au5d+8eu3btYtCgQTx8+DB62oULF+x3Iyw+m/XWvjh6U8zatWslW7Zs4u7uLj/88INERkbKnaA7cv3RdRERCXwamOynI3BEgYGB0qtXLwGkSpUq0QMVa7bHVptirl69KnXr1hURkTt37kiJEiXk2rVrImJ0j/T29hYRo5dW1apVUzbwc3R3RxsREhIiAwcOjO7CeP78+f/MExUVJS1/aSmNFjaShyEPTUhp/1auXClZsmSR7Nmzy9atW82Oo8XiP0XJhrrctm3bNs4+6lOmTJGffvop0euwBt3GbgPOnTtHjRo1mDVrFh9++CH79++nZMmS/5lPKUV7j/bs/GsndebX4e9Hf5uQ1r55enpy+PBhXnnlFZo2bcrMmTN1rxlbV7Vq7AOnJ9SzMXmrVk30IiZMmMDNmzdfOk+WLFno3r17otdhuvhUf2tfHGmLPSoqSubOnSvp06eXHDlyxPvAmi1/bpFMX2eSvJPzyrGbx5I5pWN69OiRtG7dWgDp2bOnhIaGmh1Js4i1KcbGTmth63RTjElCQkLknXfeEUAaNmwoAQEBCbr/qdunpMCUAlJ8RnGbG6HJXkRGRsoXX3wRfUrgGzdumB1Jk5e0sdvAiejshemFHZgH/AOcjs/8jlDYr127JlWqVBFAvvjii0SfUjfgcUD0Frs+hW3irVixQtKnTy958+aVw4cPmx0n1XvpzlOTTh1tb2yhjX0B0NxKy7J5u3btokqVKpw/fx5fX19Gjx6d6KNH82bMS8XcFQH4bOtnjNg2QrcXJ0L79u3Zv38/rq6uNGjQgE2bNpkdSXuRhg3Bxyd+be7P2tR9fJLlVBaOyiqFXUR2AfetsSxbJiLMmjWLRo0akSVLFg4ePEjbtm2ttux7wfcYt3scAzYMIEqirLLc1KRChQocOHCA4sWL07JlS5YsWWJ2JO1F4lPcU7ioT5s2jeDg4OjrLVq0+Fe/drsSn836+FyAwrykKQboC/gD/gULFrTub5YUEBISIj169BBAWrZsKQ8fWr+rYlRUlAzbPEwYhXRe2Vn3dU+khw8fRp+K4H//+5/ZcVKlePdjt9KA6tYQW793M5nexi7xKOwxL/bWxn779m2pUaOGAPLll19KZGRksq7v691fC6OQdsva6Xb3RAoNDRUvLy8BZMiQIcn+mmn/lqADlJ4v4lYs6pMnT5ayZctK2bJlZerUqXLlyhUpVaqUeHt7S/ny5cXT01OCgoJk+vTp4urqKuXKlZMGDRqIyP8X+mf36dWrl5QtW1a6dOkimzdvllq1aknx4sXl4MGDIiIycuRImTRpUvS6y5YtK1euXIn3/Z+nC3syOnv2rBQpUkTc3Nxk5cqVKbbe7w59J8tOLUux9TmiiIiI6APGunXrJmFh+hdQSknwkafPivkXX1itqPv7+0u5cuXkyZMnEhgYKGXKlJGjR48KIHv27BERkXfeeSe6GD+/xR6zsDs7O8vJkyclMjJSKleuLO+8845ERUXJmjVrpE2bNiLy8sIen/s/zxZ2njqkHTt2ULNmTYKCgtixYweenp4ptu7+VfvTsVxHAPwu+XE3+G6KrdtRODs7M2PGDMaNG8fixYtp06bNf04ApdmIhg2hf38YM8b4a4U29T179tCuXTvc3d3JkCEDb731Frt376ZAgQLUrl0bgG7durFnz544l1WkSBHKly+Pk5MTZcuWpVGjRiilKF++PFevXk32+yeUVQq7UmopsB8opZS6rpTqZY3lmmnhwoU0bdqUPHnycPDgQapXr25KjgchD+i4siP15tfjZuDLj5bT/kspxeeff86cOXPw8/OjdevWhISEmB1Le9727fD99/DFF8bfpB6hCi/sXaaUeun12KRNmzb6fycnp+jrTk5OREREAMb53aOi/r/TQ8yNiPjc35qs1Sums4jkERFXEckvInOtsVwziAhffvklPXr0oG7duuzbt4/ChQublidruqys67SOa4+u0WBhA24E3jAtiz3r3bs38+bNY8uWLbRq1epfvR80k8Xs/TJ6dPy7QsahXr16rFmzhuDgYIKCgvD19aVu3bpcu3aN/fv3A7B06VLq1KkDQMaMGQkMDEz0+goXLszRo0cBOHr0KFeuXElS/qTQTTExhIWF4e3tzZgxY3jnnXfYuHEjWbJkMTsW9QvXx6+bHzcCb9BgQQMCHgeYHcku9ejRg4ULF7J9+3befPNNgoKCzI6kxdalMSH93F+icuXK9OjRg2rVqlG9enV69+5N1qxZ8fDwYOHChVSoUIH79+/Tv39/APr27csbb7xBw0Q2A3l6enL//n0qVqzI999/H+u5olJMfBrirX2xxZ2nT548kebNmwsgY8aMscneKHuv7ZWM4zPGa7Bs7cWWLFkiTk5OUq9ePQkMDDQ7jkOK187TuHq/JEOXxytXrkjZsmWttrzkZBO9YhJysbXCfvfuXalevbo4OTnJnDlzzI7zUn/e/zP6S8cWv3zsxbJly8TZ2Vlq164tjx8/NjuOw4mzsMe3aFu5uOvCnkoK+7Vr18TDw0PSpk0rvr6+ZseJt3N3zkmtubXk6oOrZkexWz4+PuLs7Cy1atXSW+5Wps8Vk3S6u2MinT17llq1ahEQEMCmTZusdnqAlBAYFsiZO2eov6A+Vx9eNTuOXfLy8mL58uXRp4bQXSGtS2LrlZKY0wRYqc3dnsT63CVAqi3sBw4coE6dOkRERLBr1y7q169vdqQEqZK3Clu9t/Lo6SMaLWqkd6gmkqenJ/PmzWPr1q106tRJj6VqJW5ubty7d+/fBSop535JRcVdRLh37x5ubm6JXoZK6jdDYlSpUkX8/f1TfL3PbN26ldatW5M3b142bdpE0aJFTcuSVIcCDtF4UWPyZMzDrh67yJUhl9mR7NKsWbN4//336datGwsXLtSDjidReHg4169fj/4VlP7gQfINHkzAlCkEJ+GYEGstx9a5ubmRP39+XF1d/3W7UuqIiFSJ6/4uyZbMRv366694eXlRsmRJNm/eTK5c9l0Iq+WrxsauGxm/ZzzpXdObHcduDRw4kEePHjFixAgyZcrErFmz4nXgihY7V1dXihQp8v83/PorrFpFoaQeUerhAYUKUejwYeN/LVapaovdx8eHrl27UqlSJfz8/MiWLVuKZ0huQWFBhEWGkTX9+cqRAAAgAElEQVRdVrOj2B0RYdiwYUyaNInPPvuM8ePHmx1J0/5Fb7E/Z/78+fTu3ZvatWuzfv16MmXKZHYkqxMR2ixrQ2BYIJvf3kymtI73GJOTUopvvvmGR48e8fXXX5M5c2aGDRtmdixNS7BU0ZD47bff0rNnTxo3boyfn59DFnUwCtP71d7n6M2jvPnLmwSF6SMrE0opxXfffUenTp349NNPmTdvntmRNC3BHL6wT5w4kYEDB9KmTRvWrVtH+vSO3Q7dpnQbfnnrF/b9vY92y9vxNOKp2ZHsjrOzM4sWLaJp06b07dsXPz8/syNpWoI4dGH/6quvGDZsGJ07d2bFihX/OsOaI/Mq68Xc1nPZfHkzgzYOMjuOXXJ1dWXlypWUL1+e9u3bR5/cSdPsgUO2sYsII0eOZMyYMXTv3p25c+cmerBpe9WjYg8ioyKpV6ie2VHsVsaMGfntt9+oUaMGLVq0YP/+/f/u6aFpNsrhtthFhOHDhzNmzBh69erFvHnzUl1Rf6ZX5V6UyF4CEeHX878m+Wi21ChPnjz4+fkRFhbGG2+8wb1798yOpGlxcqjC/qy72tdff827777L7Nmz9YEmwMozK2m9rDVjdo0xO4pd8vDwYN26dVy9elUP1KHZBYepeiLCkCFDmDRpEu+99x7ff/+9LuoWnmU86VGxByN3jGT6gelmx7FLderUYfHixezfv59u3boRGRlpdiRNeyGHqHwiwgcffMDUqVP54IMPmDlzpj5qMAYn5cScVnN4y+MtPtz0IQuPLzQ7kl1q3749U6ZMYfXq1XzyySdmx9G0F4vPKSCtfbHmaXujoqLkvffeE0AGDx6sz1H+EqHhodJ4UWNxG+smNwNvmh3Hbr3//vsC2Py5+zXHQzxP22vXvWJEhEGDBvHtt99GN8PoLfUXS+uSltUdVnPs1jFyZ8htdhy7NWXKFC5cuED//v0pVqxYoodS07TkYrdNMWJpfpk1axaDBw/WRT2eMqbNGN0FcuWZlZy8fdLkRPbHxcWF5cuXU7JkSTw9Pbl48aLZkTTtX+yysIsIH330ETNnzuSjjz7if//7ny7qCRQSHsKQ34fQfHFzrjwwbzR1e5U5c2Z+/fVXnJ2dadmyJQ8ePDA7kqZFs0phV0o1V0qdV0pdUkp9ao1lvoiIMHjwYKZPn86HH37I5MmTdVFPhHSu6fity2+ERoTSdHFT/gn6x+xIdqdo0aL4+vpy9epV2rdvrwfp0GxGkgu7UsoZ+BZ4AygDdFZKlUnqcmMjli6N06ZN44MPPmDKlCm6qCdB2VfKsqHLBgIeB/DGkjd4/PSx2ZHsTp06dZgzZw7btm3j/fff1weBaTbBGlvs1YBLInJZRMKAZUAbKyz3P4YPH87UqVMZNGgQU6dO1UXdCmoWqMnKDis5cesEy08vNzuOXfL29ubTTz/lxx9/ZObMmWbH0TSr9IrJB/wd4/p14D9jViml+gJ9AQoWLJioFTVo0ICwsDC9o9TKWpRowcn+JymTM1l+aKUK48aN49y5cwwePJiyZcvSqFEjsyNpqViSR1BSSnkBzUSkt+X620A1EXn/Rfcxe8xT7cVO3DrB6rOrGdVglP7yTKDAwEBq1KjB7du3OXz4sD5hmGZ18R1ByRpNMdeBAjGu5wduWGG5mglWnFnB6F2jmbBngtlR7E7GjBlZu3YtkZGRtG3blqAgPdCJZg5rFPbDQAmlVBGlVBqgE7DOCsvVTDC64Wi6lu/K59s+Z/6x+WbHsTvFixdn2bJlnD59mp49e+qdqZopklzYRSQCGAhsAs4CPiLyR1KXq5nDSTkxr808mhZrSp9f+7D+wnqzI9mdZs2aMWHCBHx8fPjmm2/MjqOlQkluY08M3cZu+wKfBtJwYUMKZC6Ab0dfs+PYHRGha9euLFu2jPXr19OiRQuzI2kOIL5t7Lqway90P+Q+GdJkII1zGrOj2KXg4GBq167NlStXOHToECVLljQ7kmbnUnLnqeagsqXLRhrnNNwJuoPXCi9uBOp94gmRPn161qxZg6urK23btiUwMNDsSFoqoQu7FqeAwAD8LvnRfHFzHoY+NDuOXSlUqBA+Pj6cP3+eXr166Z2pWorQhV2LU8XcFVndYTXn7p6j7bK2hEaEmh3JrjRs2JAJEyawYsUKpk6danYcLRXQhV2LlybFmrCw7UJ2/rWTt33fJjJKDw2XEEOHDuWtt97ik08+YceOHWbH0RycLuxavHUu35kpTadw/NZx7gbfNTuOXVFKMX/+fEqUKEHHjh0JCAgwO5LmwHRh1xLko5ofcezdY+TKkEu3FydQpkyZWL16NcHBwXh5eREWFmZ2JC0F3Q+5T3hkypzaWRd2LcEypMlAeGQ476x9Rx+dmkAeHh7MmzeP/fv3M2TIELPjaCnk8oPL1Jxbk/c3vvAUWlalC7uWKIJw88lNfXRqInh5eTFkyBBmzZrF4sWLzY6jJTP/G/7UnFuTO0F36FahW4qsUxd2LVHSOKdhpddKKuauSIcVHdj/936zI9mVCRMmUL9+ffr27cupU6fMjqMlk98u/kb9BfVJ75qefb32UadgnRRZry7sWqJlTJuR37r+Rr5M+Wi5tCVn75w1O5LdcHFxYdmyZWTOnBlPT08eP9ajVzmah6EP6bq6K6VzlGZ/r/2UzlE6xdatC7uWJK+4v8KmbpvIni47t4Numx3HruTOnZvly5dz+fJlevfurXdGO4hnr2MWtyz4dfVjZ4+d5M6QO0Uz6MKuJVnRrEU5894ZGhRuAEBEVIS5gexIvXr1GD9+PCtWrNDD6jmA8Mhweq7rycyDxmtZPX91MqTJkOI5dGHXrMLFyRhlcebBmby+8HVCwkNMTmQ/hg4dSuvWrRkyZAj79+t9FfYqKCyItsvbsuD4Ah6EPjA1iy7smlXlzpCbPdf20GlVJ73lHk9OTk4sWLCAAgUK0KFDB+7e1Qd/2Zu7wXdptKgRfpf8+OHNH/iy/pem5tGFXbMqr7JezGoxi3Xn19Hn1z663TiesmbNysqVK7lz5w7dunUjKirK7EhaPIVGhFJ3fl2O3zrOqg6reLfKu2ZH0oVds74BVQcwqv4oFhxfwCebPzE7jt2oXLkyM2bMYNOmTYwbN87sOFo8ubm4MbDqQDa/vZm2pduaHQcAF7MDaI7py/pfcj/kPgUyF4h7Zi1anz592L17NyNHjqRmzZo0btzY7EjaC+y8upNIieT1Iq/zXrX3zI7zL3oEJS1FPAx9SBa3LGbHsAtBQUFUq1aNu3fvcvz4cfLkyWN2JO05q8+upsuqLlTMXZH9vfajlEqR9eoRlDSbcTjgMIWnFWbNuTVmR7EL7u7urFixgidPntC5c2ciIvROaFsy+8hsvFZ4USlPJTZ02ZBiRT0hdGHXkl2ZnGUonaM0nVZ2YufVnWbHsQtlypTh+++/Z+fOnYwaNcrsOBrGgUdjdo7h3fXv0qxYM7a8vYXs6bObHStWSSrsSikvpdQfSqkopVScPw+01Mk9jTsbumygaNaitFraCv8buhkuPry9venZsyfjxo3Dz8/P7DgacPH+Rd6u8DZrO63FPY272XFeKElt7EopDyAK+BEYKiLx+sTqNvbUKeBxAHXm1+Hx08cc7nOYolmLmh3J5gUHB1O9enVu3rzJ8ePHyZ8/v9mRUp2nEU+5G3yXfJnyEREVgZNywkmZ09iRIm3sInJWRM4nZRla6pEvUz62vL0F7wreFMxc0Ow4diF9+vSsXLmSp0+f0qlTJ8LDU2agBs0Q+DSQlktb0nBhQ0IjQnFxcjGtqCeE7SfUHEqxbMWY2nwqLk4u3Hpyi4DHeoi4uJQqVYrZs2ezd+9eRowYYXacVOOfoH9osLAB269sZ3jd4bi5uJkdKd7i7MeulNoCxHZqsuEisja+K1JK9QX6AhQsqLfWUrsoiaLFkhaERoSy651d5Eifw+xINq1z587s3LmTiRMnUq9ePd58802zIzm0yw8u02xxMwIeB7Cu8zpalGhhdqQEsUo/dqXUDnQbu5ZAO6/upPmS5pTNWZat3lvJ7JbZ7Eg2LTQ0lJo1a3Lt2jWOHz9OgQL64K/k0nppa/Zc28OGLhuoWaCm2XGi6X7sms2rX7g+qzqs4sTtE7Ra2oqgsCCzI9k0Nzc3fHx8CAsL0+3tyeTZhu7c1nPZ23OvTRX1hEhqd8d2SqnrQE1gg1Jqk3ViaalFixItWPLWEvb+vZePN39sdhybV6JECebMmcO+fft0e7uVrTqzirbL2xIWGUZO95x45PQwO1KiJelcMSLiC/haKYuWSnUo2wE3FzdqFahldhS70KlTJ93ebmXfHf6Ogb8NpHr+6gSHB5PGOY3ZkZJEN8VoNqF1qdbkSJ+DsMgwvtnzDWGRYWZHsmlTp07l1Vdfxdvbm7///tvsOHZLRBi+dTjv/fYeb5Z8k63eWx3inEa6sGs2ZcvlLXy69VM6rexEeKRuQ34RNzc3VqxYodvbk2jYlmGM3zOe3pV649vRl/Su6c2OZBW6sGs2pUWJFsxoPgPfc7508+2mR2F6iZjt7V988YXZcexSl/JdGNtwLLNbzY4e3tEROM4j0RzG+9XfJzwqnCG/D8HVyZWFbRfi7ORsdiyb9Ky9/ZtvvqFevXq0aGFf/a3NcCfoDj5/+PBetfeomLsiFXNXNDuS1ektds0mDa45mK8bfc36C+u5/OCy2XFsmm5vj78/7/9J7Xm1Gbp5qEO/r3Rh12zWp3U+5ex7ZymRvQRgHK2q/dez9nZ9PpmXO3j9IDXn1uReyD22em916JPQ6cKu2bQ8GY3RgybsmUB7n/a6t8wL6P7tL7fm3BoaLmxIxrQZ2ddzn8N3rdWFXbML7q7u+J7zxdPHk6cRT82OY5M6depEv379mDhxIhs2bDA7jk0JjwyPHsauVI5SZsdJdnrMU81u/OD/A/039Kd58eas7rCadK7pzI5kc0JDQ6lRowZ///13qj+fTGRUJEduHqFavmrR1+19J7w+V4zmcPpV6cdPrX5i06VNtFnWRre5x+JZe3t4eDgdO3ZMte3tIeEhdFjZgdrzanPh3gUAuy/qCaG7O2p2pVflXrg4uRAWGWYXAx6Y4Vl7e6dOnRg+fDgTJ040O1KKuhl4kzbL2uB/w5+pzaZSMntJsyOlOF3YNbvTvWL36P+3X9lO6Rylo3eyaoaOHTuyc+dOJk2aRL169WjZsqXZkVLE8VvHabW0FQ9CHuDb0Zc2pduYHckUepNHs1tBYUF0XNmR2vNq8+f9P82OY3OmTJlCxYoV8fb25q+//jI7TopYd34dAHt67km1RR10YdfsmHsad9Z3Wc+jp4+oPa82J26dMDuSTXnW3h4REUHHjh0JC3PMrqIiwvXH1wEYUW8Ex9897pBHkyaELuyaXauWrxp73tmDq7Mr9RfUZ8+1PWZHsinFixdn3rx5HDx4kGHDhpkdx+qeRjyl57qeVPqxEjcDb+KknMiePrvZsUynC7tm9zxyerC3515yZciFzx8+ZsexOe3bt2fQoEFMmzaN1atXmx3Ham49uUWjRY1YcHwB71V9j9wZYhuaOXXS/dg1h3E/5D6Z02bG2cmZe8H39JZbDGFhYdStW5dz585x9OhRihUrZnakJDlw/QCePp48CHnA/Dbz6Viuo9mRUoTux66lOtnSZcPZyZm7wXepPLsygzYOIjIq0uxYNiFNmjQsX74cZ2dnvLy8CA0NNTtSksw4OIO0zmnZ32t/qinqCaELu+ZwsrplxdPDk5mHZtJ2eVuehD0xO5JNKFy4MIsWLeLYsWN89NFHZsdJsKcRT7kZeBOA2a1m49/Xn1dzv2pyKtukC7vmcJydnJnSbArftviW3y7+Rr359Qh4HGB2LJvQsmVLPvnkE3744Qd++eUXs+PE243AGzRY2IA3lrxBRFQEGdJkIFu6bGbHslm6sGsOa0DVAazvvJ6L9y/y8eaPzY5jM8aOHUudOnXo27cvZ8+eNTtOnLZe3krlHytz6vYpRtQb4VAjHSUXvfNUc3inbp8iX6Z8ZEuXjacRT0nrktbsSKYLCAigcuXKZM+enUOHDpEhQwazI/1HRFQEo3eOZuyusZTOUZoVXiso+0pZs2OZKkV2niqlJimlzimlTiqlfJVS9j+8t+ZwyucqT7Z02QiLDKPRokZ86Pdhqh8oO1++fCxdupTz58/Tp08fzNjAi0tEVARrz6+le8XuHO5zONUX9YRIalPMZqCciFQALgCfJT2SpiUPhaJK3ipMPzidRosacevJLbMjmer1119n3LhxLFu2jFmzZpkdJ9rWy1t5/PQxbi5u7H5nN/PbzMc9jbvZsexKkgq7iPwuIs+GkT8A5E96JE1LHq7OrkxrPo0lby3B/4Y/lX+szL6/95kdy1SffPIJrVq1YsiQIezfv9/ULGGRYXy25TMa/9yY8bvHA5ApbSZTM9krq7WxK6V+BZaLyOIXTO8L9AUoWLDga6nlpESabTp5+yRvLX+LNM5pONX/VKo6V/fzHjx4QJUqVQgLC+Po0aPkzJkzxTOc/uc03r7eHLt1jD6V+zC9+XQ9kEos4tvGHmdhV0ptAWI7Vne4iKy1zDMcqAK8JfH4ptA7TzVb8CDkAXeD71IiewmCw4MJCgsip3vKFzVbcOzYMWrWrEndunXx8/PD2TnlvuhWnVlFl9VdyJw2M7NbzaZt6bYptm57Y7WdpyLSWETKxXJ5VtS7Ay2BrvEp6ppmK7Kmy0qJ7CUAGLZ5GOW/L8/GixtNTmWOSpUq8d1337FlyxZGjRqVIut8Vi6q5atGh7IdOD3gtC7qVpLUXjHNgWFAaxEJtk4kTUt5fV/rS073nLT4pQXv//Y+IeEhZkdKcT179qRXr16MHTuW9evXJ9t6RIQf/X/E08cTEaFA5gL83O5nXnF/JdnWmdoktVfMLCAjsFkpdVwp9YMVMmlaiiufqzyH+xzmw+ofMuvwLF6b/Rqn/zltdqwUN3PmTCpXrky3bt24ePGi1Zd/9s5ZGi5sSL8N/QgMCyQwLNDq69CS3iumuIgUEJGKlks/awXTtJTm5uLG1OZT+b3b70RKJG4ubmZHSnHp0qVj9erVuLi40K5dO548sc55dkLCQxi+dTiv/vAqJ2+fZHbL2Wzqtkn3ekkm+pQCmvacJsWacGbAGYpnK46I8N6G99j852azY6WYQoUKsWzZMs6ePUuvXr2scvBSpETy88mf6VK+C+cHnqfPa330YOTJSD+zmhaLZ90f7wbfZfPlzTRd3JQuq7qkmoOaGjduzNdff42Pjw+TJ09O1DKuPbrGBxs/4GnEUzKkycDJ/idZ0HZBqu15lJJ0Yde0l8jpnpOT/U8ysv5IVp1dRelZpfnB/weiJMrsaMnu448/pn379gwbNoytW7fG+373gu8x9PehlJxZkjlH5+B/w+janMVNn3EkpeiTgGlaPJ2/e57+G/pz5s4Z/hz0Z6o4zD0wMJAaNWpw+/Ztjhw5QqFChV44b3hkOJP3T2bCngkEhgXS/dXufNXgKwpkLpCCiR2bHkFJ06ysVI5SbPXeyt6ee3FP4054ZDjdVndj91+7zY6WbDJmzIivry/h4eF4enoSEvLibqAuTi6sPruaeoXqcaLfCea1maeLukl0Yde0BFBKUSybMV7opfuX2HZlG/UW1OPNX97k+K3jJqdLHiVLluTnn3/myJEj9OvXL3pnanB4MN8d/o7XZr/G3eC7KKXY1n0b6zqvo9wr5UxOnbrpM9ZrWiJ55PTg0qBLzDw4kwl7J1Dpx0q0K92OeW3mOVx7cuvWrfnqq68YOXIkxcoXI6pKFN8e/pa7wXeplq8aAY8DyJE+BxnS2N553VMj3cauaVbwIOQBk/dPZudfO9nVYxdKKfxv+FPulXIO0x8+KiqKNp3bsL74ekgDrUu1ZmjNodQpWAellNnxUgWrnQQsOejCrjkqEUEpRVBYEPmm5CONcxr6VelHvyr9yJsxr9nxEuxByAN8/vDhysMrTGg8gaCgIEq+XZLHxx7j7+dPqVKlzI6Yquidp5pmgmdbruld07Oyw0qq56/O2F1jyT8lPw0WNGDvtb0mJ4xbeGQ4v57/Fa8VXuSenJt+G/qx8dJGwiPDcXd3Z9/UfaQLSkebNm14+PCh2XG1WOjCrmnJQClF46KN+bXzr1x4/wIj64/kn6B/ogdiPnLjCN8f/p7bT26bnNRwI/AGQWFBAEw/OJ3Wy1qz8+pO+r3WD/8+/hx/9ziuzq6AcWTqqlWr+PPPP+ncuTORkZFmRtdioZtiNC0FPWuq+WLbF4zdPRaA0jlKU6dAHeoUrEPXCl2ji39yuht8l8MBh9l8eTO///k7f9z5g2Wey+hYriM3Am9w9OZRmhVrFl3MYzN79mzeffddhg4dyqRJk5I9s6bb2DXNpokIp/85zYaLG9j79172XtuLq7Mrt4bcQinFrEOzuBF4gyJZilA0a1GKZC1CgUwFXlponxcRFcGDkAf89egvTv9zmkKZC9GwSENuBN4g35R8AKR1Tku9QvVoWqwpnh6eFMlaJEGPY+DAgXz77bcsXLgQb2/vBN1XS7j4Fnbd3VHTTKCUonyu8pTPVR6AKIniRuCN6Db6Xy/8yrYr24iIioi+T438NdjfyxiX1GuFF7ef3CatS1rSOqclrUtaauavydBaQwEoOLUgfz/++1/r7P5qdxoWaUieDHmY1mwaZV8pS+0CtZM0BN3UqVM5c+YMffr0oUiRItStWzfRy9KsRxd2TbMBTsqJ/Jn+fyz4Td02ERkVyfXH17ny8ApXHlwhvWv66OnpXdPj7ORMUFgQ9yPv8zTiKXky5Ime3qNiD5yUE9nTZSdvxryUz1WeYlmNA6uUUnxQ4wOr5HZ1dWXlypXUrFmTdu3aceDAAYoXL26VZWuJp5tiNE1LskuXLlGjRg2yZ8/O/v37yZYtm9mRHJLu7qhpWoopXrw4a9as4erVq3h6ehIWFmZ2pFRNF3ZN06yiTp06zJs3jx07dtC3b1+rDNChJY5uY9c0zWq6du3KpUuXGDVqFCVLluTzzz83O1KqpAu7pmlW9eWXX3Lx4kWGDx9OsWLF6Nixo9mRUh1d2DVNsyqlFHPnzuWvv/6ie/fu5MmTh3r16pkdK1XRbeyaplld2rRpWbNmDUWKFKF169acOnXK7EipSpIKu1JqjFLqpFLquFLqd6WU/Z2+TtO0ZJE9e3b8/Pxwd3fnjTfe4Nq1a2ZHSjWSusU+SUQqiEhFYD3wpRUyaZrmIAoVKoSfnx9PnjyhefPm3L9/3+xIqUKSCruIPI5x1R3Q/Zs0TfuX8uXLs3btWi5fvkyrVq0IDg42O5LDS3Ibu1JqnFLqb6ArL9liV0r1VUr5K6X879y5k9TVappmR+rXr8+SJUvYv38/nTt3JiIiIu47aYkWZ2FXSm1RSp2O5dIGQESGi0gBYAkw8EXLEZHZIlJFRKrkzJnTeo9A0zS74OnpyaxZs1i3bh0DBgzQBzAlozi7O4pI43gu6xdgAzAySYk0TXNYAwYM4MaNG4wbN46sWbMyYcIEPV5qMkhSP3alVAkRuWi52ho4l/RImqY5sjFjxvDw4UMmTpyIu7s7X36p+1xYW1IPUJqglCoFRAF/Af2SHknTNEemlGLGjBkEBQUxcuRI0qdPz9ChQ82O5VCSVNhFxNNaQTRNSz2cnJz46aefCAkJ4eOPPyZ9+vQMGDDA7FgOQ59SQNM0Uzg7O/Pzzz8THBzMe++9R/r06enRo4fZsRyCPqWApmmmcXV1xcfHhyZNmtCrVy+WL19udiSHoAu7pmmmcnNzw9fXl1q1atGtWzd8fX3NjmT3dGHXNM107u7ubNiwgSpVquDl5YWPj4/ZkeyaLuyaptmETJky8fvvv1OrVi06d+7M4sWLzY5kt3Rh1zTNZmTMmJGNGzfSoEEDvL29mTdvntmR7JIu7Jqm2RR3d3fWr19P06ZN6dWrF99//73ZkeyOLuyaptmcdOnSsWbNGlq1asWAAQOYNm2a2ZHsii7smqbZJDc3N1auXImnpycfffQR48aN0ycOiydd2DVNs1lp0qRh2bJldO3alREjRjBw4EAiIyPNjmXz9JGnmqbZNBcXFxYtWkSePHn43//+x61bt1i8eDHp0qUzO5rN0lvsmqbZPCcnJyZNmsTUqVPx9fWladOmepi9l9CFXdM0u/Hhhx+ybNkyDh06RJ06dfQA2S+gC7umaXalQ4cO+Pn5ERAQQK1atTh16pTZkWyOLuyaptmdhg0bsnv3bkSEWrVqsXbtWrMj2RRd2DVNs0sVKlTg4MGDlC5dmrZt2zJ69GiioqLMjmUTdGHXNM1u5c+fn127dvH2228zcuRIvLy8ePLkidmxTKcLu6Zpdi1dunQsXLiQKVOmsGbNGmrWrMnly5fNjmUqXdg1TbN7Sik++uij6J2qVatWZcuWLWbHMo0u7JqmOYwmTZpw+PBh8uTJQ7Nmzfjyyy+JiIgwO1aK04Vd0zSHUqxYMQ4cOIC3tzdjxoyhfv36XL161exYKUoXdk3THE6GDBmYP38+S5cu5fTp01SsWDFVjadqlcKulBqqlBKlVA5rLE/TNM0aOnXqxPHjx/Hw8KBTp0707NkzVfSaSXJhV0oVAJoA+theTdNsTpEiRdi1axfDhw9nwYIFvPbaa+zdu9fsWMnKGlvsU4FPAH2iZE3TbJKrqytjx45l+/btPH36lDp16tCvXz8ePnxodrRkkaTCrpRqDQSIyIl4zNtXKeWvlPK/c+dOUlaraZqWKPXr1+f06dMMHjyYOXPm4OHhwcqVKx1uAI84C7tSaotS6nQslzbAcODL+KxIRGaLSBURqZIzZ86k5tY0TUuUDBkyMHnyZA4dOkTevHnx8vKidevWDnWmyDgLu4g0FpFyz1+Ay0AR4IRS6iqQHziqlMqdvJE1TdOS7rXXXuPgwYNMnjyZbdu2UaZMGRmMOqkAAAYWSURBVMaPH09wcLDZ0ZIs0U0xInJKRF4RkcIiUhi4DlQWkVtWS6dpmpaMXFxcGDx4MH/88QeNGjVi+PDhFC9enDlz5tj1gU26H7umaale4cKFWbt2Lbt27aJw4cL07duXcuXK4evra5ft71Yr7JYt97vWWp6maVpKq1u3Lnv37sXX1xelFG+99Ra1a9dm69atdlXg9Ra7pmlaDEop2rZty6lTp5gzZw5//fUXjRs3pnLlyvz888+EhYWZHTFOurBrmqbFwsXFhd69e/Pnn3/y008/ERYWhre3N0WKFGHChAk8ePDA7IgvpAu7pmnaS7i5udGrVy9Onz7Nxo0bKVu2LJ999hn58+enf//+7N271+aaaXRh1zRNiwelFM2bN+f333/nxIkTdOjQgYULF1KnTh2KFi3K8OHDOXPmjNkxAVBmfNNUqVJF/P39U3y9mqZp1hQYGMiaNWtYsmQJmzdvJioqildffZXOnTvTrFkzKlSogJOT9baflVJHRKRKnPPpwq5pmpZ0t2/fZvny5SxZsoRDhw4BkCNHDl5//XUaNWpE48aNKVq0aJLWoQu7pmmaSQICAti6dWv0JSAgADD6y8+bN4+GDRsmarnxLewuiVq6pmma9kL58uXD29sbb29vRITz58+zdetWtmzZQr58+ZJ9/XqLXdM0zU7Ed4td94rRNE1zMLqwa5qmORhd2DVN0xyMLuyapmkORhd2TdM0B6MLu6ZpmoPRhV3TNM3B6MKuaZrmYEw5QEkpdQf4K5F3zwGktpGa9GNOHfRjTh2S8pgLiUjOuGYypbAnhVLKPz5HXjkS/ZhTB/2YU4eUeMy6KUbTNM3B6MKuaZrmYOyxsM82O4AJ9GNOHfRjTh2S/THbXRu7pmma9nL2uMWuaZqmvYQu7JqmaQ7Grgq7Uqq5Uuq8UuqSUupTs/MkN6XUPKXUP0qp02ZnSQlKqQJKqe1KqbNKqT+UUh+YnSm5KaXclFKHlFInLI/5K7MzpRSllLNS6phSar3ZWVKCUuqqUuqUUuq4UipZRxqymzZ2pZQzcAFoAlwHDgOdReSMqcGSkVKqHvAEWCQi5czOk9yUUnmAPCJyVCmVETgCtHXw11gB7iLyRCnlCuwBPhCRAyZHS3ZKqcFAFSCTiLQ0O09yU0pdBaqISLIfkGVPW+zVgEsicllEwoBlQBuTMyUrEdkF3Dc7R0oRkZsictTyfyBwFkj+ASJNJIYnlquulot9bG0lgVIqP/Am8JPZWRyRPRX2fMDfMa5fx8E/9KmZUqowUAk4aG6S5GdpkjgO/ANsFhGHf8zANOATIMrsIClIgN+VUkeUUn2Tc0X2VNhVLLc5/JZNaqSUygCsAj4Ukcdm50luIhIpIhWB/EA1pZRDN7sppVoC/4jIEbOzpLDaIlIZeAN4z9LUmizsqbBfBwrEuJ4fuGFSFi2ZWNqZVwFL/q+9O1aNIgrDMPx+IGJYBRFEBAtTiFdgY7ogIhbWFprGNjdgL9h5B3aiYqEQEBQhWoiCaTQgegHBwsrSIvJbzIBtXJwc9+R9YNjt9mv24+fMmTNV9bR1nv1UVT+AN8CVxlGmtgJcG9ecHwOrSR60jTS9qvo2fn4HnjEsL09ikYp9CziXZDnJYeA6sNE4k/6h8UbifeBLVd1rnWc/JDmZ5Pj4fQm4BHxtm2paVXW7qs5U1VmG//FmVd1oHGtSSWbjhgCSzIDLwGS73Ram2KtqF1gHXjLcVHtSVZ/bpppWkkfAe+B8kp0kt1pnmtgKcJNhgvs4Xldbh5rYaeB1km2G4eVVVR2I7X8HzCngbZJPwAfgeVW9mOrHFma7oyRpbxZmYpck7Y3FLkmdsdglqTMWuyR1xmKXpM5Y7JLUGYtdkjpjsUtAkgtJtsfz0Wfj2ehdn9mifvmAkjRKcgc4AiwBO1V1t3EkaS4WuzQazyDaAn4CF6vqV+NI0lxcipH+OAEcBY4xTO7SQnJil0ZJNhiOkV1meEXfeuNI0lwOtQ4g/Q+SrAG7VfVwfL/uuySrVbXZOpv0t5zYJakzrrFLUmcsdknqjMUuSZ2x2CWpMxa7JHXGYpekzljsktSZ3xl4wjwuE+7dAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "f = lambda z: 2*np.sin(z) - z**2/10 # objective function\n", "fp = lambda z: 2*np.cos(z) - z/5 \n", "fpp = lambda z: -2*np.sin(z) - 1/5 \n", "\n", "x = 2.5 # initial guess\n", "while abs(fp(x)) > 1e-12:\n", " x = x - fp(x)/fpp(x)\n", "\n", "x_vals = np.linspace(0, 5)\n", "plt.plot(x_vals, f(x_vals), 'k', label='f(x)')\n", "plt.plot(x_vals, fp(x_vals), 'g--', label=\"f'(x)\")\n", "plt.plot(x, f(x), 'rx', ms='20', label='optimum') # X marks the spot...\n", "plt.xlabel('x')\n", "plt.title('maximum values of f(x) is {:.2e} at x={:.2e}'.format(f(x), x))\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "***\n", "\n", "**Exercise**: Consider the objective function $f(\\mathbf{x}) = ax_0 + bx_1 - (2x_0^2 + x_0 x_1 + 2x_1^2)$. Compute $\\nabla f(\\mathbf{x})$ and $\\mathbf{H}(\\mathbf{x})$ (which is the jacobian of $\\mathbf{\\nabla g}(x)$. Then, for $a = 5$ and $b = 10$, determine the values of $x_0$ and $x_1$ that maximize $f(\\mathbf{x})$ starting with $\\mathbf{x} = [1, 2]^T$.\n", "\n", "*Solution*:\n", "\n", "The gradient is\n", "\n", "$$\n", "\\nabla f(\\mathbf{x}) = \n", "\\left [ \\begin{array}{c}\n", " \\partial f / \\partial x_0 \\\\\n", " \\partial f / \\partial x_1 \n", "\\end{array} \\right ] =\n", "\\left [ \\begin{array}{c}\n", " a - 4 x_0 - x_1 \\\\\n", " b - 4 x_1 - x_0\n", "\\end{array} \\right ] \\, .\n", "$$\n", "\n", "The hessian matrix is the jacobian matrix for $\\mathbf{\\nabla f}(\\mathbf{x})$, or\n", "$$\n", "\\mathbf{H}(\\mathbf{x}) =\n", "\\left [ \\begin{array}{cc}\n", " \\partial \\nabla f_0 / \\partial x_0 & \\partial \\nabla f_0 / \\partial x_1 \\\\\n", " \\partial \\nabla f_1 / \\partial x_0 & \\partial \\nabla f_1 / \\partial x_1 \\\\\n", "\\end{array} \\right ] =\n", "\\left [ \\begin{array}{c}\n", " -4 & -1 \\\\\n", " -1 & -4\n", "\\end{array} \\right ] \\, .\n", "$$\n", "\n", "Now, apply Newton's method:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEaCAYAAAAboUz3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXu0XVV97z9fzjnh5MFDDBoSuIVbqUNFjTWDq4NqEcGLiOC7cK9Fb7WRVkZrW6tSbtWq3Gpp1Tp8puJFbhW1IpVqVKDVqtz6CBBeBRW9WMMBYngkhOSQ8/jdP/baZGVnP9baez3mXPv3GWONrL1ec+59Tr77d77zN39TZobjOI7TPA6ouwOO4zhOObjAO47jNBQXeMdxnIbiAu84jtNQXOAdx3Eaigu84zhOQ3GBdxzHaSgu8E4pSDpc0tWSHpB0saS/lPSmuvsFIOkHkp5Sdz8cp2zkE52cMpD0fmDazH5f0uHAZuAJZra75q4h6VXAb5nZy+vui+OUiUfwTlmcDPxDsv9aYGMI4p5wJfA8SUfU3RHHKRMXeKdQJC2RtB14KvBPkm4GXgj8a+qav5J0Rer1RZL+WdLUEO3lfpaZzQLXAS/I257jxMRk3R1wmoWZ7ZH0bOCbZvZ4AEm/BH6Uuux9wE8lrQWeBZwK/IaZzQ3R5LDPug14+hDtOU40uMA7ZbAWuDH1+lDgofYLM7tP0geBS4FDaAny9vZ5SRcCzwXuBc4xs129Gur1LEmHAFcDTwaeZWa3dNz6EOAWjdNo3KJxyqBT4B8ADuq45gZaNs75ZvaL9kFJxwG/ambPAa4BfidDe92etQt4EfDFHvccBDyY4dmOEy0u8E4ZPJ19Bf4m4NfaLyQ9FfgY8Gn2F/DnAF9L9r8G/Ea/hno9y8zmzOyXfW59UkcfHadxuMA7ZdAp8BuB3wSQtAb4J+Bc4PeBp0o6MXXtY4C2XbMdOCy57xJJl6QbyfCsrkg6EHgmLQvHcRqLC7xTKJJW0RLp21OHLwVOS3zxjcD7zezKxFu/CLgwde0DtLx0kn/vT/aPAq5NtXNwhmf14gzgW2Y2k/f9OU5M+EQnpxIk/S9gq5l9cMB1bS/9v0laDxwIfILWXwRPy5tpk0T9f50eZJX0feB1XQZeHadRuMA7wSHpL2l571tpZdE8PORzNtIa8P058Akzu6SwTjpOBLjAO47jNJRKPHhJd0q6WdJmSZuqaNNxHKdqJH1K0lZJaUvwnZLuSvRvs6TTetx7qqQfSbpD0tsK6U8VEbykO4F1Zrat9MYcx3FqQtJzgZ3ApWZ2XHLsncBOM/vrPvdNAD8GTgG2AD8Ezjazfx+lP55F4ziOUxBm9m32Zn7l4XjgDjP7mZntAT4HnDlqf6oSeAOuknRdkhnhOI4zTpwn6abEwnlMl/NrgF+kXm9Jjo1EVbVoTjCzGUmPA66WdHvyTQdAIvrrAaaWTjzzsGMOLqUT84tx/MEyZxOlPn+hoM9hYVEj3b847P1578t7fQ7XUos1XJfhmkHPGXy+/4eghQHn+92/2KfxhX7nFvY7tMPu32Zmh/ftzABOed5Su+/+bD+gG27acyswmzq0wcw2DLjtY8C7af1mvRv4G/afwd3tl3Rk/7wSgW9PKDGzrUlp1+OBb6fObwA2AKx6ymH26s+eUkW3HuWXj6yotL08bJtdXtqz759dVshzduyeHun+2d1Lct+zsHuIX91d2b84J2azfylM7sp27UTGaviTGa6b7Fl+bS9Tu/rrw9Su/qI2+fD+grrP/Tvn+56f2PlIz3Pa2edNPrSz56nF7Tv2eX3VI5/9ed9OZOC++xf5ztdXZbp2xer/mDWzdXmeb2b3tvcl/R3wlS6XbaE1ma/NkcDIE/FKD2klLZd0UHufVg3uoCaYHH7gzq5bCKycfpiV00OlgQ/ksOldHDadQSkGcPDSWQ5eOjv4wh5ML92T+56JpfNMLO0vMPuxrL9gpVmYNhamswVQ88uyXbewtLUNfF6WazJ8N88t6//FM7es/3//+eX9vxDnVvT/kl1YcWDPc7aiz5s8qHfAdcAh5fx1XyYdC8u8lO7690PgWEnHSFoCnEVrYZqRqCKCfzxwhaR2e581s69X0O7IdBP5uqL9tMgXHdW3RX7UiL4t8sNE9G2RzxvNTyydzxfNt0U+YzS/MG2Zovm2yGeJ5heWDo7m55cOjuTnlw2O5OeWqW8kP7fsgL6R/Pzyib6R/NyKyb6R/MKKA3tG8rZiae9I/qAVPSP5Aw45eL9IPhQkXQacCKyUtAV4B3Bisl6BAXcCb0iuXQ180sxOM7N5SecB3wAmgE+Z2a2j9qd0gTezn9GghRU6Rb8OwW+LfchCP6xtM710z1AiDzltm2ULhYs8tITeRX5fxknkzezsLocv7nHtDHBa6vVGWvWVCiOOUceAqdPWads3RVs4RVg3o9g200v3DG3b5CIAy2bgs5YOtmzmlw22bNyuGU9c4AumLsFvqtDnJbc3v2wht9BnoUiRh2J8+SaKvNMfF/iSqVrsy4jqixL6YQg1ms/C/DLLJPTjJPL9cJEvHl+TtULSIl+Fd1+0V3/Y9K6R/PlRB2ErGYCtyZfP4slDMb58nZ58Pz8ehvfkR2HOjJmFYdZ7Dx+P4Guiysi+yIg+tmh+KMsmI1kjechm2VSZRllnJN/PqnGKxQU+AKoS+5CEvg5vPjM5fPk6B18Ht9f/fKgi71ZNcQRn0UxpgVVLil3s/p49hxb6vDKpwsYp0roZNbVyWNtmmLz5ECyb0NIoQ7Vr6rBqmshYRPCrljzYdwuVsqP6oiP6Uagqmg/BsvFIfi8eyZfLWAj8IEIX/rItnKKEvi7bphLLJiMu8j3ud5GvBRf4PoQo+OMi9HkZdgA2Mzl9+Sy4yO9laJF3+uICn4OQovwyo/oihX5Yqormy8qXd5HPj4t88bjAj0gIgl+m0I9KXdF8Hsry5bNm2BQ5ISp0kS9ttqvTleCyaGInLfJVZ++0Rb7I7JuiMm5GmSQ1TPGy0rNsasqwaUJ2zaDiZFWzhwOYmW+mn+8RfInUFdmXEdEXYduMEs0HadnUNPja9Ejeo/jicIGviDqsnJCFflhc5JNnRSLyfe91kS8dF/iaqFLsyxL6URg1ms9LEzNsYhB5z6ypFxf4AIhV6OuM5oexbPKmUpY5+JqFKkV+EC7y2ZD0REmbU9sOSW/quOZESdtT17y9rP74IGtAVDVAW/Rg7KgDsaOUOxh2ALbuwdfQBl6LWBmqb/sDShr0Y1AFypAwsx8BawEkTQB3AVd0ufQ7ZnZ62f0JWuBXT2WLamfm4qk1k5W22Mcm9KNk2wybaeMiX53I96NJmTUF8Xzgp2b287o6EJzAT2k+s7C3yXt9TF8IVUT1RQp9XdH8MEXLXOS7tRNm+mRMUXyKs4DLepx7tqQbgRngzUUssN2N4AS+Cvp9IYQs/mVH9UULfQzRvIt8t3bGS+TnbIK75h+T9fKVkjalXm8wsw2dF0laApwBnN/lGdcDv2JmOyWdBvwjcGzObmdiLAW+H93EPzTRj0Xoi4jmqxJ5yD4pykV+dJHvR+CR/DYzW5fhuhcC15vZvZ0nzGxHan+jpI9KWmlm24rsKHgWTSZWTz243xYCZadaFpVxM0qmzbDplGWnUo5Ddk2M6ZMBcTY97BlJqyQp2T+elg7fV0YnXOCHJDTBL0voi0qtHDWlcliRHyaVMitNEPnB7fQ/X6bIx4qkZcApwJdSx86VdG7y8hXALYkH/yHgLDMbLj1pAC7wBRGK2Mcg9MNS1cSocRL5KiZC9X12CeUM6sbMdpnZY81se+rYx83s48n+h83sKWb2dDN7lpn937L64gJfAiFE92UK/ai4yA+mSSJfR80ap4ULfAXUKfZlCH0R0fwolk1VvnwsIl8ERcx27UeZJYad3rjAV0xdYl+W0I9KHb58HkorbVCgyDdh0HUQLvLD4QJfI3UJfZGEEM3npUyRhxzRvIv8vs9u6KBrnfjXYgCkRb6KnPsy8ugPP3BnIXnzw+TMD5MvH8yEqALz5GPJke/bdg3lDOZsMri5LkXhEXxgVBnVF23bFBXND8MwvnzeNMpxiuQHt9P/fB2Drs7+uMAHSpVefRlCPwoh+/J5RT4zgYl83YOufe91kc+MC3wEVCn0RdFkXz6G7JoqRN79+PCJ6qtwzeQDma/NUTwoGtoiX6ZfWKQ/X0RNm1B9+TyefFmLeQ8iiyc/8Bk11qwZw/LChROcwE9pIZeQ96LJXwZVCX1Rg7CjDsC6yO9PUcXJYh90dfpTmUUjaULSDZK+UlWbWVkz+UDXLXTKtm6K9ObrHHzNSxB2TQaKmghVd/rkKH68058qP9k/BG6rsL2RiUX0qxD6Ihg1y2asRD4wP75s3I8vh0r+xpF0JPAi4ELgj6tosyw6RT4ke6dM66Zob35Yy2bYGvNR2jWB5ciH6sePypxNlLo0Zp1UFcF/EHgL0PWnJ2m9pE2SNm2/r5wfYlmEGN2XGdHHbNl4JB9+Zk3fZ3skn5vSBV7S6cBWM7uu1zVmtsHM1pnZukMeG+8PMTQ7pyyhL8qbr8OycZEf3bd3Pz4eqvg0TwDOkHQn8DngJEl/X0G7tROK2Dc5mneRLx7345tD6QJvZueb2ZFmdjStVcb/xcxeXXa7oVG30JcZzReBi3wGMoi8WzX1IulUST+SdIekt3U5f6Ckzyfnvy/p6DL7438PVUzdUX0ZQl+kZTMsLvJ7aYLIx2jVSJoAPkJrwe0nA2dLenLHZa8DHjCzJwAfAN5XZp8q/RTN7FtmdnqVbYZM3UJfNC7y+1Na7ZoBxOLH9yNCq+Z44A4z+5mZ7aFlSZ/Zcc2ZwKeT/S8Cz28vwF0G8X1NNpC6ovqyovlRGVeRD3HQtWw/vmH1atYAv0i93pIc63qNmc0D24HHltUhn+sbGG2RrzK/fvXUg4XmzheRMz9KHZthShsMkyefh6x58kXnyBdB3fnxZTO/eECe37OVkjalXm8wsw3Jfrdvq843luWawvAIPlCqjuibFs1XEcnnXRkqK+Pox/cjMD9+WzulO9k2pM5tAY5KvT4SmOm4/9FrJE0ChwD3l9XZoD45Z3/qEPoicZHfS+yDrgOfMYLIl7mea4X8EDhW0jGSltDKGryy45orgdck+6+glVXoEfy4U6XQFx3Nu8jvpa7iZE3w40Mn8dTPA75Bq+7WF8zsVknvknRGctnFwGMl3UGrbMt+qZRF4h58ZFTp0RfpzRdRfnjYGjZVePJl1a3JxJj48f1q1YSCmW0ENnYce3tqfxZ4ZVX98Qg+UqqK6IuO5EeN5sctkh83q2YQkVg1weCfVuRUIfShWTYu8j2oSOTdqomHaC2a1ZOjVSMEmJkffpWh0Fgz+UDptk1Ilk3Idk0Z5FoRagBZV4Pq+wy3aqIgaIEvQsSHfX6M4l+FP+8iP5gm+PEh1I93Ric4gV/CYunCnoVufYhF9MuO5sdV5PNQhsgXOQmqqEVCyqSqKH7OJkr9XagT9+BzsHpy5z5byJTtzRfpy8fiyY+jHz/wGT4BKmj8ExqBGMS+igHYInCRz07Ri3f3o4pZrv3wAdfRcIEviJDFvopovgjqEvm8DFOcLAuFV54MKHWy//P7n3eRHx4X+BIIVezHQeSHoez0ySZYNaOmTpZp1Ti9cYEvmdDEvsxoPgSRrzJHPg+1inxFuFUTHi7wFRKa0JfBOIl8CH58JtyqGVtc4GsgFKF3kd+fskU+K6FaNQOf4VZNUASXBz9OtEW+zvz6siZHFZUrP0qe/LA58nmpexJUkbNcBxFybvywLCweUPtM5bLwCD4AQojoy4jmY43kGzvo6gOuY4cLfEDULfQu8nsJadC1UCKxavrhXnx2XOADpE6hb6rID0Nj/fgCqHvA1cmGC3zA1CX0IYv8sIQ46FrWSlADaYBVE2MUL+nzkjYn252SNve47k5JNyfXbep2TVZc4CPARX4vdWTWlEnRVk1oWTVlEpvIm9lvmdlaM1sLXA58qc/lz0uuXTdKmy7wkVBHNF/GpKgYRd6tmt74gGt+JAl4FXBZ2W25wEdGE6J5F/l98QHXkZuIjecA95rZT3qcN+AqSddJWj9KQ54HHyF15M8XXWO+iDz5IhbyDoWs+fHjmBtf9sIgC4vKPDcBWNnhi28wsw3tF5KuAVZ1ue8CM/tysn82/aP3E8xsRtLjgKsl3W5m387awTQu8BGzenLn2Iv8sAwzCSrvIiFlLfVXqMgXtDjIIEZd/SkgtvXzxc3s5H43S5oEXgY8s88zZpJ/t0q6AjgeGErg3aKJnKq9+dDsGrdqCiCDVTMIT5vMzMnA7Wa2pdtJScslHdTeB14A3DJsYy7wDcFFvjrKnASVVeSrHnANIW2yIZxFhz0jabWkjcnLxwPflXQj8APgq2b29WEbc4umQVRp2ZS97mtehvXjq6hXU7dVk4mKrJpRaMIi3Wb22i7HZoDTkv2fAU8vqj2P4BtG3TVthqXOzJqmWzVNSpt08uEC30Cq8uVDs2pGIcRJUFkIcYZrmTTIi68EF/gGM44iX6Uf71F8dzyKD4fSBV7StKQfSLpR0q2S/qLsNp29xCjyozJuVk2IJYXLxKP47FQxyPoIcJKZ7ZQ0RWuE+Gtm9r0K2naoZvC1yEHX2PLjnf0ZNPkppLz4xUUVN1gdGKVH8NaiHRpNJVvYVYwaSGyRvFs1LTyK745H8dmoxIOXNJGUxtwKXG1m36+i3TSrJ6Yyb03FRT4bIS4QEiPuxddPJXnwZrYArJV0KHCFpOPM7NHZWUlBnfUAq9f0z8XNwqgi3ev+mYW5kZ4bAlWXN3C6U8Y6rjGWMBgFj+IHU2kWjZk9CHwLOLXj+AYzW2dm6w47bLguVRGBNyXKLzuS9yi+PgpNmywAj+LrpYosmsOTyB1JS0lqMRT1/LoEN3axd5EfTNm58cGnTUbgxTv9qSKCPwL4pqSbgB/S8uC/MupDQxLXkPqSh5hEPhY8it8fj+Lro4osmpvM7Blm9jQzO87M3jXqM0MV01iFPgY8ivco3slPVDNZYxHQWPoJcUXxdS/cnZUQlvjzKD4Hi2Jh92SmLTaiEfhYBDNNLEIfk8iPQqhRfB48infyEIXAxyCS/Yih/7FUoWxqFO8MjuKd/AQv8DGIYxZiiObLFHmP4vszDrNbR131yQdb8xO8wDeN0EW+TIoS+VjKGHgU79RN0ALfVDEM+X3FYtXUQZOj+ExUsHZrk1MmJb0yqai7KGld6vjRknZL2pxsH+9x/2GSrpb0k+TfgdX9ghX4kEWwCEK2bGKwajyKL5Yq124dY24BXgZ8u8u5n5rZ2mQ7t8f9bwP+2cyOBf45ed2XYAV+XAhV5J3ujH0UXwDjGsWb2W1m9qMRHnEm8Olk/9PASwbdEKTAj5vohfh+xyGKr4pGRfFNTJlcVKuwWpYNVkralNrWF9SLYyTdIOlfJT2nxzWPN7O7AZJ/HzfoofFl7jeU1RNTwVWrbHrlyVVLHuSePfkXFil7UZA8lSazkrXSZCgMWhCkRraZ2bpeJyVdA6zqcuoCM/tyj9vuBv6Tmd0n6ZnAP0p6ipntGLWzQUbw40qIkXxZeBRfDKFOfBpEU1MmzezkpCRL59ZL3DGzR8zsvmT/OuCnwK91ufReSUcAJP9uHdQfF3inL03PqmlCXnyIRGfT1EhScXci2f/PwLHAz7pceiXwmmT/NUDPL402wQn8lOpbQCAEPIrPTyxRfFmM62BrbEh6qaQtwLOBr0r6RnLqucBNkm4Evgica2b3J/d8MpVS+V7gFEk/AU5JXvdlaA9e0lvN7H3D3u/0JjQ/3r34Ylg5/TDbZpeX3s6oxLLiU6g2TS/M7Argii7HLwcu73HP61P79wHPz9Nm5ghe0hdS2z8Arx94U8BMfWQHB1xbjH95wLWzTH1k5PGQfRiXSD7mKD5Gmya0KN4plzwWzQ4ze1WyvRK4pqxOVcHC2iVMv2HbyCJ/wLWzTL9hGwtr48lQGAb34oshlsHWTAQws9XpTx6Bv7Dj9QVFdqRqFk+YZvYTK0cS+ba4z35iJYsnDF4UOS8exeejiVF8XfjM1mYwUOAlfVCSzOz/pY+3BwFiZhSRL1vc24Qk8k2P4kNkHGa21o7BxKwybbGRJYLfCVwpaTmApBdIurbcblXHMCJflbg71RK7TVM0oeTEu00zPAMF3sz+J3AZ8C1J3wX+hAxFbmIij8jXIe4hRfFl4TZNbzwn3hmWLBbN84HfBR4GDgf+wMy+U3bHqiaLyHvk7jZN6LhN46TJYtFcAPy5mZ0IvAL4vKSTyurQnNX3bd9P5F3cxwO3afbFbZq4yWLRnGRm3032bwZeCLyn7I7VRTeRD0Hc3abJjts04eE2TT3kLlWQlKnMNZsqNtIiP/VXD9Yu7qHhNk3YuE3jtBmqFo2ZhVnIs0AWT5hm7pwVLPngDubOWeHiPka4TTMEbtMESXDFxoAg6rAccO0sU5fuZM+bDmbq0p2FlTUYBbdpstNEm6Yumj7pSYswuUuZttgIUuDrJu25z73l0JFnvDYRt2mqx2vTOHkJVuDriuK7DagWUdbAcbrRqNo0BeA2TbEEK/B10C9bxkV+vKjKh4+FqtIlnWIJWuCrjOKzpEK6yFeD+/DNJVQfvqkELfBVkSfP3UXeqRP34eNF0kWSbpd0k6QrJB2aHD9F0nWSbk7+7TqRVNI7Jd0laXOynTaozeAFvuwofphJTC7yLXygtRga5cN7umQ/rgaOM7OnAT8Gzk+ObwNebGZPpbXW6v/p84wPmNnaZNs4qMHgBb5MRpmh6iLffNyH35ei0iXHFTO7yszaH+L3gCOT4zeY2Uxy/FZgWtKBRbQZhcCXEcUXUX7ARd5x8hO5D79S0qbUtn7I5/wO8LUux18O3GBmj/S477zE4vmUpMcMaiQKgYfiRX5i855Cyg+0RX5is/uYIdLEgdbY69KEhhZhYne2DdhmZutS24Z9niVdI+mWLtuZqWsuAOaBz3Tc+xTgfcAbenT1Y8CvAmuBu4G/GfTeBiydPjqSjgIuBVYBi8AGM/vbYZ41szBX2GzOuTceXMhzoCXyXsqgWNZMPsBd8wMDFGcABy+dZcfuwb+b00v3MLs7nHWF55dZ35mjC0sfFdygMLOT+52X9BrgdOD5ZntL50o6ErgCOMfMftrj2femrv874CuD+lNFBD8P/ImZPQl4FvBGSU+uoF3HiYZY6tJ4PvzwSDoVeCtwhpntSh0/FPgqcL6Z9VwtT9IRqZcvBW4Z1GbpAm9md5vZ9cn+Q8BtwJphnxdCnZq6GOf37jSLyH34YfkwcBBwdZLm+PHk+HnAE4A/T6VAPg5A0iclrUuu+6sklfIm4HnAHw1qsHSLJo2ko4FnAN/vOL4eWA+wes3EwOcUadU4Tj9WLXmQe/YcWnc3nAZgZk/ocfw99Fhjw8xen9r/7bxtVjbIKmkFcDnwJjPbkT5nZhvagxaHHZatS+MWzY7b+3X64xOenCxUIvCSpmiJ+2fM7EtFPddFzwmVkEoWeOGx8aV0gZck4GLgNjN7f9HPn1mYa7zQN/39OQ3EB1qDoAoP/gTgt4GbJW1Ojv1Zlmm2eXBf3nGqYWLpPAu7R5eOhWljYrb+RTS0CJMBplwWQekCnyzYXclPsYki79H7+LBy+mG2zS6vuxtOg4hmJmtWmmTZhP4+ZuZX1N0Fx3H60DiBbxO6OA4i9v47jlM/jRV4iDeaj7HPTvWMQ6rkoEwapz+NFvg2MQl9LP10nIF4Jk3tjIXAtwlZ6EPumxM/oebCj2nJgsoYK4FvE5qYhtSXrPgA63jji3/EQaW1aEIjLaxVp1fGKOqO00gWYbKhJfbHWuDTdApuWYLvwp4NrwXvOKPjAt+DbkKcV/SbKuZuzzhOHLjA56Cpgu04ITNodSenN2M5yOo4jjMOuMA7uXB7xnGGQ9K7Jd2UrNh0laTVyXFJ+pCkO5Lzv97j/mcmKzrdkVw/8M8aF3gnOHyA1WkoF5nZ08xsLa0Fs9+eHH8hcGyyrQc+1uP+jyXn29eeOqhBF3gnMx69O7nx2ayP0rGS3XKgPcvrTOBSa/E94NCOBbbbC24fbGb/ZmYGXAq8ZFCbPsjqNJqZOV9PNXRCqQufkZWSNqVebzCzDVlvlnQhcA6wndbC2QBrgF+kLtuSHLs7dWxNcrzzmr64wDuZ8Oh9fJheuofZ3Uvq7kZlaBGmdmUumbDNzNb1fJZ0DbCqy6kLzOzLZnYBcIGk84HzgHfQfb2Mzg5luWY/XOCdgVQp7u6/OzFjZidnvPSzwFdpCfwW4KjUuSOBmY7rtyTH+12zH+7BO04f7tnjFo9TDJKOTb08A7g92b8SOCfJpnkWsN3M0vYMyeuHJD0ryZ45B/jyoDY9gnf64tG74xTGeyU9EVgEfg6cmxzfCJwG3AHsAv5H+wZJm5OsG4DfAy4BlgJfS7a+uMA7PYndd/cBVickzOzlPY4b8MYe59am9jcBx+Vp0y0apyuxi7szPOM0wNp0XOCd/ahD3Jtmz/zyEf+CdOonOIHfwwEePdZIUz77IuwZH2B1YidYD35mfgWrJ3fW3Q2nApoWvQ/LttnldXdhLGnlwS/W3Y1SCFbgYW806UJfPk2J3MeJ+2eX1d0FJ3CCs2i6MTO/wgWoROr8bMuI3j17xnFaRCHwbVzoi8c/z+64/+40gagEvo0LfTHU/Rk21Xv3DBonFIL24AfhHv1w1C3sUJ64uz3TnR27p+vuglMDUUbwnXhEnx3/nMLEM2icMog6gu/EI/rehCTsoUfv7r8XyK6Junsw1jRK4NukxWzcxT4kYYfm+u5tqvLfPUWyOLRoTD7czJWnGinwacZV7EMTdihX3D16d5z9abzAp2m62Ico6m2aHrk3hToKjUW0XF90lC7wkj4FnA5sNbNcpS7LpEliH7KwV0EomTPD2jM+wOqURRUR/CXAh2mtAh4knQIZuuDHJuixRO9NtWfKSJFc2D1Wf/xHS+k/JTP7tqSjy26nSLoJaJ2iH5ugpylb3EOJ3qvEB1jjRNK7gTNprei0FXitmc1I+lPgvyeXTQJPAg43s/s77r8E+E1ge3LotWa2uV+b/jWckX574SGCAAALDElEQVQiW5T4xyzk3YhJ3EeN3t2ecTJwkZn9OYCkPwDeDpxrZhcBFyXHXwz8Uae4p/hTM/ti1gaDEHhJ64H1AI9fHUSXctE0YS6CWGwZx6kKM9uRerkcsC6XnQ1cVlSbQcxkNbMNZrbOzNYd8lifGBEzd80/phJxDyl6D5k8/rsv1ZeJlZI2pbb1eW6WdKGkX9CyZN7ecW4ZcCpweZ9HXCjpJkkfkHTgoPbiC5edYKkqag/Nd6/SnonKf49kFqsWjKmd81kv32Zm63o+S7oGWNXl1AVm9mUzuwC4QNL5wHnAO1LXvBi4to89cz5wD7AE2AC8FXhXv85WkSZ5GXAirW++LcA7zOzistt1qiVWS6bJ0XtZeAZNb8zs5IyXfhb4KvsK/Fn0sWfM7O5k9xFJ/xt486BGqsiiOTvP9XM2wV3zj2HN5ANldckpmCrFvSnRu9MiyySnyV3NmAgl6Vgz+0ny8gzg9tS5Q2hlyLy6z/1HmNndkgS8BLhlUJvBfhW3RcOFPlyqjtqLFvc6o/cq7BkvERwc75X0RFppkj8Hzk2deylwlZk9nL5B0kbg9WY2A3xG0uGAgM0d93clWIFv40IfJi7uzYrefYC1fMzs5X3OXUJrUmjn8dNS+yflbTN4gW/jQh8GdXjtodkyjhML0Qh8m7TAuNhXR12DqGWIe93Re8zZM5kHWCPJoGk6QeTBD0tVOdfjjot7fLj/7kCEEXw3PKovhzq/PEO2ZaqO3sumLv89lAwaLRoTOx8pvZ06aITAp3GvfjRC+IuoLHGPNXoPxZ5x4qNxAt/Go/p8hCDsEL64xxC9l2XP+ASn+BiLn1ineLngtwhF1Ns0WdyHpezovRR7xgdYg2EsBL6TcY/uQxN2CNtzL4IQvXen+YylwKcZh+g+REFvU7awxxy9D0Ms2TO+Dms1jL3Ad9JNDGMT/ZAFPc24iPuw0XtI9kwd+e9NqUFTJy7wGeglmCEIfyxi3kks4h4TsUTvTnW4wI9AVnEd5osgVuEeRBVee5Hi3tTo3UmxuIh27q67F6UQnMDPWXBdGpmminVeXNzLI2/0Xoo9k5Gi/PeJZmpyoQSppm0hWD31YM09cYqgqgyZkMR9FKKN3t1/D44gBb5NWhhc7OOjytTH0Dz3kNMivTTw+BC0wKfxqD4eqs5pL1rc67RmhoneyxxcDdWecbIRXTXJmblDH92c8Bh3cR+FaK0ZqHz2asz+u6Q3SzJJK5PXJ0raLmlzsr29x33HSPq+pJ9I+rykgX+KRSfwaVzow6Hqn8U9ew4NUtyrtmZiGlzNSpP9d0lHAacA/9Fx6jtmtjbZ3tXj9vcBHzCzY4EHgNcNai9qgW/jUX091PW5l+G31y3uUUfvTh4+ALwFsDw3JQttnwR8MTn0aVoLb/clGg8+Kz4wWz51fpGGKu51UGb0nouM9kxD/PeVkjalXm8wsw1ZbpR0BnCXmd3Y0ut9eLakG4EZ4M1mdmvH+ccCD5rZfPJ6C7BmUJuNE/g0nULkgj88df91VFaWTFHi3qSBVQjbnincf19YhId2Zr16m5mt63VS0jXAqi6nLgD+DHhBl3PXA79iZjslnQb8I3Bs56O73Dfwr4BGC3wnHt3no25Rb+PiPjoevVeDmZ3c7bikpwLHAO3o/UjgeknHm9k9qfs3SvqopJVmti31iG3AoZImkyj+SFrRfl/GSuDTeHTfnVBEvU2TxX1Ymhq9Nxkzuxl4XPu1pDuBdWa2TdIq4F4zM0nH0xobva/jfpP0TeAVwOeA1wBfHtSu/yQTxlXwQxP0NmVOXApF3KOP3gumFnsmDF4B/J6keWA3cJaZGYCkjcDrzWwGeCvwOUnvAW4ALh70YBf4HvQSvtiFP1RBb1P2jNTYxT2o6N3tmaExs6NT+x8GPtzjutNS+z8Djs/TTnACP2dhL/cVi/CHLuTdiEXcR6VKcY8lenfKITiBh33/o69aEpZw9iKvoA7zhRCjaGehijoyRYp7yHVmRqHO6H2M7ZlSCVLg07T/88ci9FlpqljnITZhh7isGY/eneAFvk2MUb3TnaoqP7q4Z6eM6D0rtZcmWFhgcfuOevtQEtEIfBoX+zipsqRvU8Q9ZoocXHV7ZjiiFPg0TbVwmkLVddrLGEitU9ybHr075RK9wLfpFBIX/HqpYwEOF/dwfPciB1ed4WmMwHfiNk711LWqUojCPipl57u3iSF6d3tmeCoReEmnAn8LTACfNLP3VtFuGxf78qh7qbyQxb1q371UayYHHr2HQ+kCL2kC+AitIvdbgB9KutLM/r3strvhVs7o1C3qUN6kpRDEPUhrxqP3KKkigj8euCOZZoukzwFnArUIfCcu+NkIQdSh3NmosYr7MHj0Ph5UIfBrgF+kXm8B/kv6AknrgfXJy0f+Zu0XbqmgX1WzklbJz6bRxPfVxPcEzXxfTxz1ATvs/m9c9chnV2a8PKrPrwqBH1ioPlkRZQOApE39CurHir+veGjie4Jmvq+O1ZWGwsxOLaIvIVLFmqxbgKNSrzMVqnccx3FGowqB/yFwrKRjJC0BzgKurKBdx3GcsaZ0i8bM5iWdB3yDVprkp7osKJsm0wK2EeLvKx6a+J6gme+rie+pMJQsHOI4juM0jCosGsdxHKcGXOAdx3EaSlACL+lUST+SdIekt9XdnyKQ9ClJWyU1Jrdf0lGSvinpNkm3SvrDuvtUBJKmJf1A0o3J+/qLuvtUFJImJN0g6St196UoJN0p6WZJm4tIl2wiwXjwSUmDH5MqaQCcXVdJg6KQ9FxgJ3CpmR1Xd3+KQNIRwBFmdr2kg4DrgJc04GclYLmZ7ZQ0BXwX+EMz+17NXRsZSX8MrAMONrPT6+5PEUi6E1hnZlFNPqqSkCL4R0samNkeoF3SIGrM7NvA/XX3o0jM7G4zuz7Zfwi4jdaM5aixFjuTl1PJFkYENAKSjgReBHyy7r441RKSwHcraRC9aDQdSUcDzwC+X29PiiGxMjYDW4GrzawJ7+uDwFuAxbo7UjAGXCXpuqTcidNBSAI/sKSBExaSVgCXA28ys0YsamlmC2a2ltaM6+MlRW2rSTod2Gpm19XdlxI4wcx+HXgh8MbEDnVShCTwXtIgIhKP+nLgM2b2pbr7UzRm9iDwLSD2OiUnAGckfvXngJMk/X29XSoGM5tJ/t0KXEHL5nVShCTwXtIgEpLByIuB28zs/XX3pygkHS7p0GR/KXAycHu9vRoNMzvfzI40s6Np/Z/6FzN7dc3dGhlJy5MBfiQtB14ANCZTrSiCEXgzmwfaJQ1uA74woKRBFEi6DPg34ImStkh6Xd19KoATgN+mFQ1uTrbT6u5UARwBfFPSTbQCjqvNrDFphQ3j8cB3Jd0I/AD4qpl9veY+BUcwaZKO4zhOsQQTwTuO4zjF4gLvOI7TUFzgHcdxGooLvOM4TkNxgXccx2koLvCO4zgNxQXeCZKkHPEpyf57JH2o7j45TmyUviar4wzJO4B3SXocrWJmZ9TcH8eJDp/o5ASLpH8FVgAnmtlDyZT0jwJ7gG+Z2Wdq7aDjBI5bNE6QSHoqrdIBjyQ15wFeBnzRzH4Xj+gdZyAu8E5wJCtGfYbWgi8PS/qvyakj2btmwEIdfXOcmHCBd4JC0jLgS8CfmNltwLuBdyant9ASefDfXccZiHvwTjQkHvyHgVngu+7BO05/XOAdx3Eaiv+Z6ziO01Bc4B3HcRqKC7zjOE5DcYF3HMdpKC7wjuM4DcUF3nEcp6G4wDuO4zQUF3jHcZyG4gLvOI7TUP4/1YTzLuXA+88AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "f = lambda z: 5*z[0] + 10*z[1] - (2*z[0]**2 + z[0]*z[1] + 2*z[1]**2)\n", "gf = lambda z: np.array([5-4*x[0]-x[1], \n", " 10-4*x[1]-x[0]]) # gradient of f\n", "H = np.array([[-4.0, -1.0],\n", " [-1.0, -4.0]]) \n", "\n", "# Apply Newton's method\n", "x = np.array([1.0, 2.0])\n", "while sum(abs(gf(x))) > 1e-8:\n", " x = x - np.linalg.solve(H, gf(x))\n", " \n", "# Plot the results and show where the maximum exists\n", "x0, x1 = np.linspace(0, 5), np.linspace(0, 5)\n", "x0, x1 = np.meshgrid(x0, x1)\n", "plt.contourf(x0, x1, f([x0, x1]), 25)\n", "plt.colorbar()\n", "plt.plot(x[0], x[1], 'rx', ms=20)\n", "plt.xlabel('$x_0$')\n", "plt.ylabel('$x_1$')\n", "plt.title('$f(x_0, x_1)$');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### When the Derivative is Unavailable\n", "\n", "If we have $f'(x)$, we can apply Newton's method to $f'(x) = 0$. Assuming we don't, how can we proceed? One option is to adapt the bisection method. Now, the search is not based on the change of the sign of $f(x)$ but rather it's magnitude relative to the left and right points. Again, we have to have isolated our (single) target value in some finite range, which may be no easy task (plotting helps). Alternatively, Newton-like methods can be used by approximating both $f'(x)$ and $f''(x)$ (or their multi-dimensional analogs) with finite differences." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using `scipy.optimize.minimize`\n", "\n", "An alternative to Newton's method is the `minimize` function in `scipy.optimize`. From `help(minimize)`, we find\n", "\n", "```\n", "Help on function minimize in module scipy.optimize._minimize:\n", "\n", "minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), \n", " tol=None, callback=None, options=None)\n", " Minimization of scalar function of one or more variables.\n", "```\n", "\n", "At a minimum, one must define (1) the function `fun` to minimize, and (2) an initial guess `x0` for the solution. The function minimize returns an *object* with several attributes, one of which is the solution `x`. Therefore, if one calls `sol = minimize(fun, x0)`, the solution is `sol.x`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "**Exercise**: Use `minimize` to find the *maximum* of $f(x) = 2\\sin(x) - x^2/10$ starting with $x_0 = 2.5$. \n", "\n", "**Solution**: Note that any maximization problem can be turned into a minimization problem negating the objective function, i.e., $f(x)$ becomes $-f(x)$. Then:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.42755179]\n" ] } ], "source": [ "import numpy as np\n", "from scipy.optimize import minimize\n", "f = lambda x: -2*np.sin(x) + x**2 / 10\n", "sol = minimize(f, x0=2.5)\n", "print(sol.x) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example**: Use `minimize` to find $\\mathbf{x} = [x_0, x_1]^T$ that minimizes $f(x) =x_0^2 - x_0 - 2 x_1 - x_0 x_1 + x_1^2$. Use an initial guess of $x_0 = x_1 = 0$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.33333353 1.66666666]\n" ] } ], "source": [ "f = lambda x: x[0]**2 - x[0] -2*x[1] - x[0]*x[1] + x[1]**2\n", "x_initial = [0.0, 0.0]\n", "sol = minimize(f, x0=x_initial)\n", "print(sol.x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving Constrained Optimization Problems\n", "\n", "Optimization with constraints is a much more challenging problem. Although `minimize` can be used for such problems, we limit our focus to a specific class of constrained optimization problems known as [linear programs](https://en.wikipedia.org/wiki/Linear_programming). The standard form of a linear program is\n", "\n", "$$\n", " \\min_{\\mathbf{x}} \\mathbf{c}^T \\mathbf{x}\\, ,\n", "$$\n", "\n", "subject to \n", "\n", "$$\n", " \\mathbf{A}_{ub}\\mathbf{x} \\le \\mathbf{b}_{ub} \\, ,\n", "$$\n", "\n", "and \n", "\n", "$$\n", " \\mathbf{A}_{eq}\\mathbf{x} = \\mathbf{b}_{eq} \\, .\n", "$$\n", "\n", "Note, the objective function is *linear*, and $\\mathbf{c}$ represents a vector of costs associated with each element of $\\mathbf{x}$. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The ideal tool for solving such problems is `scipy.optimize.linprog`, for which help provides\n", "\n", "```\n", "linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='simplex', callback=None, options=None)\n", " Minimize a linear objective function subject to linear\n", " equality and inequality constraints.\n", " \n", " Linear Programming is intended to solve the following problem form::\n", " \n", " Minimize: c^T * x\n", " \n", " Subject to: A_ub * x <= b_ub\n", " A_eq * x == b_eq\n", "```\n", "\n", "Here, `A_ub` and `A_eq` are two-dimensional arrays, while `c`, `b_ub`, and `b_eq` are one-dimensional arrays. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "**Exercise**. Use `linprog` to find $\\mathbf{x} = [x_0, x_1]^T$ that *maximizes* $x_0 + x_1$ subject to the constraints\n", "\n", "$$\n", "x_0 \\ge 0 \\\\\n", "x_1 \\ge 0 \\\\\n", "x_0 + 2 x_1 \\leq 4 \\\\\n", "4 x_0 + 2 x_1 \\leq 12 \\\\\n", "-x_0 + x_1 \\leq 1 \\\\\n", "$$\n", "\n", "*Solution*: Note first that the cost vector $\\mathbf{c}$ should be $[-1, -1]^T$ so that the sum of $x_0$ and $x_1$ is *maximized*. All of the inequality constraints can be cast in $\\mathbf{A}_{ub} \\mathbf{x} \\le \\mathbf{b}_{ub}$ form by noting that if $x_0 \\ge 0$ then $-x_0 \\le 0$. Therefore, these constraints can be written as\n", "\n", "$$\n", "\\left [ \\begin{array}{cc}\n", "-1 & 0 \\\\\n", " 0 & -1 \\\\\n", " 1 & 2 \\\\\n", " 4 & 2 \\\\\n", "-1 & 1\n", "\\end{array} \\right] \n", "\\left [ \\begin{array}{cc}\n", "x_0 \\\\\n", "x_1 \n", "\\end{array} \\right] \\leq\n", "\\left [ \\begin{array}{c}\n", " 0 \\\\\n", " 0 \\\\\n", " 4 \\\\\n", "12 \\\\\n", "1 \n", "\\end{array} \\right] \n", "$$\n", "\n", "Finally, there are no equality constraints. From here, `linprog` can be used via " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[2.66666667 0.66666667]\n" ] } ], "source": [ "from scipy.optimize import linprog\n", "import numpy as np\n", "c = np.array([-1, -1])\n", "A_ub = np.array([[-1.0, 0.0],\n", " [0.0, -1.0],\n", " [1.0, 2.0],\n", " [4.0, 2.0],\n", " [-1.0, 1.0]])\n", "b_ub = np.array([0, 0, 4, 12, 1])\n", "sol=linprog(c, A_ub, b_ub)\n", "print(sol.x) # like minimize, linprog returns an object \n", " # with attribute x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "***\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "**Exercise**. Suppose the Acme Concrete Company is supplying concrete for three large-scale projects. The company has production facilities in four different locations. Projects 0, 1, and 2 require 250, 400, and 500 cubic meters of concrete each week. Plants A, B, C, and D produce 350, 200, 400, and 450 cubic meters each week. Although the plants, in total, produce more than enough to supply all three projects, the cost to provide the concrete depends on the distance traveled between the plant and the project. The cost per cubic meter to ship from any plant to any project is summarized n the table below:\n", "\n", "| | 0 | 1 | 2 |\n", "|----|---|---|\n", "| **A** | 10 | 20 | 30 |\n", "| **B** | 15 | 5 | 20 |\n", "| **C** | 25 | 20 |15 |\n", "| **D** | 30 | 10 | 20 |\n", "\n", "How much should each plant supply to each project? Convert this problem into a linear program, and then use `linprog` to determine the solution.\n", "\n", "***\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Further Reading\n", "\n", "Checkout out the [SciPy documentation](https://docs.scipy.org/doc/scipy/reference/optimize.html) on optimization and root finding." ] } ], "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.8" } }, "nbformat": 4, "nbformat_minor": 2 }