{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Unconstrained optimization\n", "\n", "In this chapter, we'll cover how to apply `scipy.optimize.minimize` to unconstrained optimization problems. As a reminder, unconstrained optimization considers:\n", "\n", "$$\\mathop {\\min }\\limits_x f\\left( x \\right)$$\n", "with $x$ the design variable of length $n$ and $f$ the objective function.\n", "\n", "## Method\n", "In this course, we're making use of the function `scipy.optimize.minimize`. The documentation of this function is available here: \n", "https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html. In this course we'll cover only the relevant parts.\n", "\n", "For unconstrained optimization we need to run at least `scipy.optimize.minimize(fun, x0, ...)` with:\n", "- `fun`, the objective function $f(x)$ to be minimized. `fun` is a callable. The `scipy.optmize.minimize` function takes care of defining and inputing our design variable $x$.\n", "- `x0`, the initial guess for our design variable $x$. It needs to be a `ndarray` with length $n$\n", "\n", "The function `scipy.optimize.minimize` outputs an object `scipy.optimize.OptimizeResult`. with:\n", "- `scipy.optimize.OptimizeResult.x` the optimized solution of the design variable `x`. It is a `ndarray` with length $n$\n", "- `scipy.optimize.OptimizeResult.success`, a indication whether or not the optimizer was executed succesfully. It is a `bool`, indicating `True` or `False`\n", "- `scipy.optimize.OptimizeResult.message`, a message describing the cause of termination of the optimizatino algorithm. It is a `str`.\n", "- `scipy.optimize.OptimizeResult.fun`, the values of the optimized objective function $f$. It is a `int` or `float`\n", "- `scipy.optimize.OptimizeResult.nit`, the number of iteration performed by the optimizer. It is a `int`\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 1\n", "\n", "### Problem\n", "It is desired to determine the number of bathymetry maps $n$ of a local area that should be produced to maximize the profit of a company. The total cost of production and distribution is €$75$ per unit $n$. The revenues are proportional to the number of units multiplied by its price: $Revenues = n \\cdot Price $\n", "\n", "The demand depends on the price ($Price = 150 - 0.01n^2$), as shown in the graph:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAG4CAYAAABVU6MoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABYxklEQVR4nO3deVhUhf4G8PfMzr7LooCouIKIoOYOqai5lZlrLqWlmQuZpmX9whYxK7OrpWWWXpe0cms1sQgXTBEEBVcUZRHEBVkEhgHO7w+vcy+hCDhwZuD9PA+PzZkzh3e+UbyeOYsgiqIIIiIiIrovmdQBiIiIiIwZyxIRERFRFViWiIiIiKrAskRERERUBZYlIiIioiqwLBERERFVgWWJiIiIqAosS0RERERVYFkiIiIiqgLLEhEREVEVjKIsHThwAMOGDYObmxsEQcDu3bsrrXPmzBkMHz4cNjY2sLKywmOPPYbU1FT981qtFrNnz4ajoyMsLCwwfPhwpKen1+O7ICIioobIKMrSnTt34Ofnh9WrV9/3+YsXL6JXr15o27Yt/vrrLyQkJOCtt96CRqPRrxMaGopdu3Zh27ZtOHToEAoKCjB06FCUlZXV19sgIiKiBkgwthvpCoKAXbt24cknn9QvGzt2LJRKJTZt2nTf1+Tm5sLJyQmbNm3CmDFjAABXr16Fu7s7fv31VwwcOLA+ohMREVEDZBR7lqpSXl6OX375Ba1bt8bAgQPRpEkTdOvWrcJHdbGxsdDpdAgJCdEvc3Nzg4+PD6Kjox+4ba1Wi7y8PP1Xbm4url+/DiPrj0RERCQhhdQBHiY7OxsFBQVYtmwZ3nvvPXzwwQfYu3cvRo4cicjISPTt2xdZWVlQqVSws7Or8FpnZ2dkZWU9cNvh4eFYsmRJpeVbt26Fubm5wd8LERER1Y0RI0bU2baNviyVl5cDuDuEV155BQDQqVMnREdHY+3atejbt+8DXyuKIgRBeODzr7/+OubNm6d/nJeXB3d3dwQHB8PBwcFA76Bx0ul0iIiIwIABA6BUKqWOY7I4R8PhLA2HszQMztFwdDpdnW7f6MuSo6MjFAoF2rdvX2F5u3btcOjQIQCAi4sLSkpKkJOTU2HvUnZ2Nnr06PHAbavVaqjV6krLlUolf3ANhLM0DM7RcDhLw+EsDYNzNH5Gf8ySSqVCly5dcO7cuQrLz58/D09PTwBAQEAAlEolIiIi9M9nZmYiMTGxyrJERERE9DBGsWepoKAAycnJ+scpKSmIj4+Hvb09PDw8sGDBAowZMwZ9+vRBcHAw9u7di59++gl//fUXAMDGxgZTp07Fq6++CgcHB9jb22P+/Pnw9fVF//79JXpXRERE1BAYRVk6fvw4goOD9Y/vHUc0efJkbNiwAU899RTWrl2L8PBwzJkzB23atMGOHTvQq1cv/Ws++eQTKBQKjB49GkVFRejXrx82bNgAuVxe7++HiIiIGg6jKEtBQUEPPV3/+eefx/PPP//A5zUaDVatWoVVq1YZOh4RERE1YkZ/zBIRERGRlFiWiIiIiKrAskRERERUBZYlIiIioiqwLBERERFVgWWJiIiIqAosS0RERERVYFkiIiIiqgLLEhEREVEVWJaIiIiIqsCyRERERFQFliUiIiKiKrAsEREREVWBZYmIiIioCixLRERERFVQSB3AGI1Zdwxqc0vIBAEyAZDJBChlMijkApRyGZRyASqFDBqFHGqlDGqFHBqlHOaqu18WasXdP1UKWGoUsNYoYaVRwNpMCUu1AnKZIPVbJCIiompiWbqP89cKIFOX19n2rTQK2JmrYGeuhK25CvYWKtiZq+BopYKjhRqOVio4WKjRxFoNR0s1lHLuACQiIpIKy9J9fDbeD1bWthABlJeLKCsXUVouQldWDl2ZiNKycpSUlaNYV4Zi3X//LNKVobCkFHe09/4sRYG2FHnFpcgr0kFbereA5ReXIr+4FKm3Hp5FEAAHCzWcrdVwttbA2VqDprYauNqYwdVWg6a2ZnCx0UCtkNftUIiIiBoplqX76NHCAQ4ODgbfrra0DPnFpbhdqMPtwhLkFOqQU1iC24UluHmnBDcLSnCjQHv3K//uP5eWi/plSVfz7rtdQQCcrTRwtzeDu505mtmbw8PeHF6O5vB0sICDhQqCwI/+iIiIaoNlqR6pFXKoLeVwtFRXa/3ychG3CktwLa8Y2XlaZOUVIyu3GFdvFyHzP39ezS1Csa787nN5xYi5nFNpO1ZqBTwdzeHlaImWThZo6WSJlk6WaOFkAY2Se6SIiIiqwrJkxGQyAY6Wd49b6uB2/3VEUcTNOyVIu1WItJwipN0qRHpOIa7cvPt1NbcI+dpSJGbkITGj4p4pQQCa2ZmhdRMreDtboY2LJbybWKFVE0uWKCIiov9gWTJxgvDfQuXvYVfp+WJdGdJuFSLlxh2k3LiDi9cLcPH63T9vF+qQdqsIabeK8MfZbP1r5DIBXo4WaOdqjXauVmjnao0OrtZoYq2pz7dGRERkFFiWGjiNUg5v57t7jv7pZoEWF7ILcOFaPs5dy8f5awU4fy0ftwt1SM4uQHJ2AX5K+O/6TlZq+LhZw6epDTq42cC3mQ3cbDQ8HoqIiBo0lqVGzMFSDQdLNR5r8d+D2UVRRHa+Fqcz83AmMw9nMvNxJjMPl64X4Hq+FpHnriPy3HX9+o6WanRyt4FfM1v4udvCr5ktbMyVUrwdIiKiOsGyRBUIgqC/REFwmyb65UUlZTiTlYfEjFwkZuTiVEYezl/Lx40CLfafycb+M//9GK9VE0t0amYDxW0B3tcK0NbNFjJeiJOIiEwUyxJVi5lKjs4eduj8P8dFFevKkHQ1FwlpuUhIv434tNu4crNQ/xEeIMe21dGw1igQ2Nwegc3t0LW5PXyb2fC6UEREZDJYlqjWNEo5AjztEeBpr192s0CLE6m3cfzyTeyPv4iMIgXyikvx59ls/Pmfg8hVChk6udvisRYO6N7CAf4etjz7joiIjBbLEhmUg6Ua/ds7o6+3PdrpLiBk4OO4cKMIx1Ju4fjlHMRcvoWbd0pwLOUWjqXcwr/+uACVQgb//5SnXt6O6ORuy1u8EBGR0WBZojqlkMvQsZktOjazxbTedw8gv3TjDo6l3MKRizfx96WbyM7X4mjKLRxNuYVP/7gAC5UcXb3s0bOVI3p5O6KNsxXPuCMiIsmwLFG9EgRBfwXxcV09IIoiUm7cwZFLNxF98SaOXLyJW3dKKpx152ytRm9vJ/Rp7YTerRxhZ6GS+F0QEVFjwrJEkhIEAS2cLNHCyRITunmivFzEmaw8HE6+gUPJN3Es5Sau5WnxQ2w6fohNhyAAHZvZIriNEx5v2wQ+bjY8046IiOoUyxIZFZlMQAe3uxe9fLFPSxTryhBz+RYOnL+OA+dv4Ny1fCSk3UZC2m2s3H8BjpZqBLVxQr+2TdC7tRMs1fyRJiIiw+JvFjJqGqUcvb2d0NvbCYuHAFm5xYg6f/fMukMXbuBGwX/3OqnkMjzW0gED2jVBv3bOcLM1kzo+ERE1ACxLZFJcbDQY08UDY7p4oKS0HDGXb+HPs9n448w1XL5Z+J89UNfx1p4kdHCzxqAOLhjo4wLvJpY8SJyIiGqFZYlMlkohQ89WjujZyhFvDmmHi9fvYP+Za9h/+hpiU3OQdDUPSVfz8HHEebRwtEBIBxcM8nGBXzMbFiciIqo2o7iYzYEDBzBs2DC4ublBEATs3r37getOnz4dgiBg5cqVFZZrtVrMnj0bjo6OsLCwwPDhw5Genl63wcloCIKAVk0sMaNvS/zwUg/ELO6PD572xeNtm0All+HSjTtYG3URT352GL0+iMT7v5zGidQciKIodXQiIjJyRlGW7ty5Az8/P6xevbrK9Xbv3o2jR4/Czc2t0nOhoaHYtWsXtm3bhkOHDqGgoABDhw5FWVlZXcUmI+ZoqcaYLh74ekoXxL7VH6vG+WNIR1eYq+TIuF2EdQdT8NTn0fridDL9NosTERHdl1F8DDd48GAMHjy4ynUyMjIwa9Ys/P777xgyZEiF53Jzc7F+/Xps2rQJ/fv3BwBs3rwZ7u7u2L9/PwYOHFhn2cn4WWmUGObnhmF+bijWleGvc9n45VQW/jhzTV+c1h1MQXMHc/16rZ2tpI5NRERGwijK0sOUl5dj4sSJWLBgATp06FDp+djYWOh0OoSEhOiXubm5wcfHB9HR0Q8sS1qtFlqtVv84Ly8PAKDT6aDT6Qz8LhqXe/MztjnKAfRr44h+bRxRrGuHAxdu4NdT1/DHuWxcvlmIVX8mY9WfyWjjbIlhHV0x3M8VrjYayfIa6xxNEWdpOJylYXCOhqPT6aBUKuts+yZRlj744AMoFArMmTPnvs9nZWVBpVLBzs6uwnJnZ2dkZWU9cLvh4eFYsmRJpeWRkZEwNzd/tNAEAIiIiJA6wkOFWAF9/YHEHAFxNwScuS3g3LUCnIu4gI8jzqOltYhARxF+DiLMJfovxhTmaCo4S8PhLA2DczSMESNG1Nm2jb4sxcbG4tNPP0VcXFyNz2ASRbHK17z++uuYN2+e/nFeXh7c3d0RHBwMBweHWmemuy0/IiICAwYMqNO2b0hP/efP3CIdfk+6hj0JmTh2OQfJeQKS84CdqTIMaNsET/m7omdLByjq4Wa/pjhHY8VZGg5naRico+HU9d45oy9LBw8eRHZ2Njw8PPTLysrK8Oqrr2LlypW4fPkyXFxcUFJSgpycnAp7l7Kzs9GjR48HblutVkOtVldarlQq+YNrIKY4S0elEhO6e2FCdy9k3C7CnvgM7D6RgfPXCvBLYhZ+ScyCk5UaI/2b4umAZvVyfJMpztFYcZaGw1kaBudo/IzibLiqTJw4ESdPnkR8fLz+y83NDQsWLMDvv/8OAAgICIBSqaywKzMzMxOJiYlVliWih2lqa4aZQa3we2gf/Dy7F6b0aA47cyWu52vxxYFLCPnkAJ787DC+PZaK/GIed0BE1BAZxZ6lgoICJCcn6x+npKQgPj4e9vb28PDwqPSRmFKphIuLC9q0aQMAsLGxwdSpU/Hqq6/CwcEB9vb2mD9/Pnx9ffVnxxE9CkEQ4NPUBj5NbfDGE+0QeS4bP8SmI/JsNuLTbiM+7Tbe+ek0nvB1xdiu7gj0tOOFL4mIGgijKEvHjx9HcHCw/vG944gmT56MDRs2VGsbn3zyCRQKBUaPHo2ioiL069cPGzZsgFwur4vI1IipFDIM7OCCgR1ccD1fi51x6dh+PA2Xrt/Bjrh07IhLR6smlhjX1QNPd24KW3OV1JGJiOgRGEVZCgoKqtEFAS9fvlxpmUajwapVq7Bq1SoDJiOqmpOVGtP7tsSLfVogLjUH22PS8FNCJpKzC/Duz6fxwd6zeMLHBeO6eqCrlz33NhERmSCjKEtEpk4QBAR42iPA0x5vDW2P3fFXsfVoKs5k5mF3/FXsjr+K1s6WmPiYJ570bworDQ/mJCIyFSxLRAZmpVFi4mOeeLabB06m52Lr0VT8mHAV568V4K09SVj221k81bkpnn3ME21drKWOS0RED8GyRFRHBEGAn7st/Nxt8caQdtgZl45Nf1/Bpet3sPnvVGz+OxXdvOzxXM/m6N/OuV6u20RERDXHskRUD2zMlHiupxem9GiOI5duYvPfV/B70jUcTbmFoym34GajwbPdPTG2iwfsLXhAOBGRMWFZIqpHgiCgR0tH9GjpiMzcImz++wq+PZaGq7nFWL73HD7dfwFP+TfF8728eDNfIiIjwbJEJBFXGzMsGNgWsx/3xs8nM7EhOgWJGXnYFpOGbTFp6O3tiCndPVBe/RNFiYioDrAsEUlMo5RjVEAzPN25KWIu5+DrQynYdzoLBy/cwMELN+BsJkehczqeDvSARsnrhhER1TeWJSIjIQgCunrZo6uXPdJuFWJD9GVsi0nFtaIyLN5zGp/8cRFTenhiQjdP2PG4JiKiesPTb4iMkLu9Od4a2h4H5/fFk55lcLXR4EaBFh/tO48ey/7E23sSkXarUOqYRESNAssSkRGz0igQ7Cbij1d64dOxndDe1RpFujJsPHIFQR/9hdBtJ3AmM0/qmEREDRo/hiMyAUq5DCM6NcVwPzccTr6JLw5cxMELN/RXBw9u44QZfVuiWwuHh2+MiIhqhGWJyIQIgoBe3o7o5e2IxIxcrIm6iN9OZSLy3HVEnruOLs3tMDO4FYJaO/E+dEREBsKP4YhMlE9TG3w2vjP+fDUI47t5QCWXIeZyDp77JgbDVh/C3sRMlPO6A0REj4xlicjENXe0wNKnfHFwYTCm9fKCmVKOxIw8zNgch4ErD2BPfAbKWJqIiGqNZYmogXC21uDNoe1xeNHjmBXcClZqBS5kF2DutngM+CQKu06ko7SsXOqYREQmh2WJqIGxt1Bh/sA2OPz645g3oDVszJS4dP0OXtmegAGfHMCOWJYmIqKaYFkiaqCsNUrM6eeNQwuDsWBgG9iaK5Fy4w5e/T4BIZ/w4zkioupiWSJq4Kw0Srwc3AqHFj6OhYPaws5ciUs37mDutngMXHkAP5+8ygPBiYiqwLJE1EhYqhV4KaglDi58HAsGtoGNmRLJ2QWYtfUEBn96EL8nZUEUWZqIiP6JZYmokbFUK/BycCscXBiMeQNaw0qjwLlr+Zi+KRZPfnYYB85fZ2kiIvofLEtEjZT+mKbX7p49Z66SIyE9F5O+PoaxX/6N45dvSR2RiMgosCwRNXI25krMH9gGUQuC8VzP5lDJZTiacguj1h7B1A0xOJeVL3VEIiJJsSwREQDAyUqNt4d1QOSCIIzt4g65TMAfZ7Mx6NMDmPddPNJzCqWOSEQkCZYlIqqgqa0Zlj3dEfte6YPBPi4QRWBnXAYe/ygK7/x0Gjl3SqSOSERUr1iWiOi+WjpZYs2zAdj9ck90b+GAkrJyfH04BX0+jMTaqIso1pVJHZGIqF6wLBFRlTq522LrC92w8fmuaOdqjfziUiz77Swe/+gv7IxL5zWaiKjBY1kioocSBAF9Wzvh59m98PEzfnCz0eBqbjHmfZeAoasOITr5htQRiYjqDMsSEVWbXCbg6YBm+HN+EBYOagsrtQKnM/Mw/qujmLYxBhevF0gdkYjI4FiWiKjGNEo5XgpqiajXgjGlR3PIZQL2n8nGwE8OIOzHJB4ETkQNCssSEdWavYUKYcM74PfQPujfrglKy0VsiL6Mvh9GYv2hFOjKyqWOSET0yFiWiOiRtWpiia8md8GWad3QztUaecWlePfn0xi48gAiz2VLHY+I6JGwLBGRwfRs5YifZ/dC+EhfOFiocOn6HTz3TQymfHMMydk8nomITBPLEhEZlFwmYFxXD0QuCMKLfVpAKRfw17nrGLTyAN79+TTyinVSRyQiqhGWJSKqE9YaJd54oh32vdJXfzzT+kMpePyjKPwQy+szEZHpYFkiojrl5WiBryZ3wYbnuqCFowVuFGgx//sEPL02GqfSc6WOR0T0UEZRlg4cOIBhw4bBzc0NgiBg9+7d+ud0Oh0WLlwIX19fWFhYwM3NDZMmTcLVq1crbEOr1WL27NlwdHSEhYUFhg8fjvT09Hp+J0T0IEFtmmBvaB8sGtwWFio5TqTexvDPDuGNXad4qQEiMmpGUZbu3LkDPz8/rF69utJzhYWFiIuLw1tvvYW4uDjs3LkT58+fx/DhwyusFxoail27dmHbtm04dOgQCgoKMHToUJSV8f5VRMZCpZBhRt+W+HN+EJ7s5AZRBLYeTcXjH/+F7TGp/GiOiIySQuoAADB48GAMHjz4vs/Z2NggIiKiwrJVq1aha9euSE1NhYeHB3Jzc7F+/Xps2rQJ/fv3BwBs3rwZ7u7u2L9/PwYOHFjn74GIqs/ZWoOVY/0xrqsH3tqTiPPXCrBwxyl8eywN7z3pA5+mNlJHJCLSM4o9SzWVm5sLQRBga2sLAIiNjYVOp0NISIh+HTc3N/j4+CA6OlqilET0MN1aOOCXOb3x5pB2sFDJEZ92G8NXH8LbexJ51hwRGQ2j2LNUE8XFxVi0aBHGjx8Pa2trAEBWVhZUKhXs7OwqrOvs7IysrKwHbkur1UKr1eof5+XlAbh7nJROx/9RP4p78+McH01jmePkx9wxqL0Twveexy+nsrDxyBX8eioTi59oiyd8nCEIwiN/j8Yyy/rAWRoG52g4Op0OSqWyzrZvUmVJp9Nh7NixKC8vx+eff/7Q9UVRrPJ/suHh4ViyZEml5ZGRkTA3N3+krHTXPz9CpdppLHMMsQQ82wn4PkWG6wUlCP3uJNb+Xo5RXuVwMjPM92gss6wPnKVhcI6GMWLEiDrbtsmUJZ1Oh9GjRyMlJQV//vmnfq8SALi4uKCkpAQ5OTkV9i5lZ2ejR48eD9zm66+/jnnz5ukf5+Xlwd3dHcHBwXBwcKibN9JI6HQ6REREYMCAAXXa9hu6xjjHJwDM1JXhy4OXsfZgCs7mAssTFZjZtwVe6NUcKkXtjh5ojLOsK5ylYXCOhlPXe+dMoizdK0oXLlxAZGRkpSITEBAApVKJiIgIjB49GgCQmZmJxMRELF++/IHbVavVUKvVlZYrlUr+4BoIZ2kYjW2OSqUS8wa2xVMB7vi/PYk4eOEGVv6RjF9OZSF8pC8Cm9s/0rYb0yzrEmdpGJyj8TOKslRQUIDk5GT945SUFMTHx8Pe3h5ubm4YNWoU4uLi8PPPP6OsrEx/HJK9vT1UKhVsbGwwdepUvPrqq3BwcIC9vT3mz58PX19f/dlxRGR6vBwt8O/nu2JP/FW8+/NpXMguwKi1RzChmwdeG9QWNmb8BUNEdc8oytLx48cRHBysf3zvo7HJkycjLCwMP/74IwCgU6dOFV4XGRmJoKAgAMAnn3wChUKB0aNHo6ioCP369cOGDRsgl8vr5T0QUd0QBAFP+jdF39ZOCP/tDL47no4tR1Ox7/Q1vDO8Awb7ukodkYgaOKMoS0FBQRDFB1+Mrqrn7tFoNFi1ahVWrVplyGhEZCTsLFRYPsoPT/k3w+Jdp3Dpxh28tCUOAzs4450RPnC21kgdkYgaKJO8zhIRNV7dWzrg17m9MSu4FRQyAb8nXUP/FVH49lhqtf5iRURUUyxLRGRyNEo55g9sg59m94JfMxvkF5fi9Z2nMG7d37h8447U8YiogWFZIiKT1c7VGjtn9sSbQ9pBo5Th70u3MOjTA/jq4CWU8T5zRGQgLEtEZNLkMgHTerfAvtC+6NHSAcW6crz3yxk8szYaydkFUscjogaAZYmIGgQPB3NsmdYN4SN9YalWIC71Np7410F8/lcySsvKpY5HRCaMZYmIGgxBEDCuqwf2vdIHQW2cUFJajuV7z+HpNdG4wL1MRFRLLEtE1OC42Zrhmyld8NEzfrDWKJCQnosRnx/B/gyBe5mIqMZYloioQRIEAaMCmiFiXl8Et3GCrkzET6lyjP0qhscyEVGNsCwRUYPmbK3B11O6YNlTHaCRi0hIz8UT/zqIdQd4xhwRVQ/LEhE1eIIg4OnOTfG6Xxn6eDugpLQc7/96BuPW/Y3Um4VSxyMiI8eyRESNhq0a+GpiZyx9yhfmKjmOpdy9LtPWo7z6NxE9GMsSETUqgiBgfDcP7J3bB12b26OwpAxv7DqF5zbEIDuvWOp4RGSEWJaIqFHycDDHty8+hsVPtINKIcNf565j4MoD2JuYKXU0IjIyLEtE1GjJZQJe6NMCP8/uhfau1sgp1GHG5ji8+l0C8ot1UscjIiPBskREjV5rZyvsfrknZga1hCAAO+LSMfjTgziWckvqaERkBFiWiIgAqBQyvDaoLb6b3h3N7MyQnlOEMV8ewYe/n4WOF7IkatRYloiI/keX5vb4bW5vjApoBlEEPou8iKfXRCPlxh2poxGRRFiWiIj+wUqjxEfP+OHzCZ1hY6bEyfRcPPHpQWyP4SUGiBojliUiogd4wtcVe0N7o3sLBxTpyrBwxynM2ByLnDslUkcjonrEskREVAVXGzNsmdYNrw9uC6VcwO9J1zD404M4cvGm1NGIqJ6wLBERPYRMJmB635bYNbMnWjhaICuvGOO/+hsf7zuHUh78TdTgsSwREVWTT1Mb/DS7F0YH3j34e9WfyRj9xRGk3eL95YgaMpYlIqIasFArsHyUH1aN84eVWoG41Nt44l8H8ctJXvmbqKFiWSIiqoVhfm74dW5v+HvYIr+4FC9vjcMbu06hWFcmdTQiMjCWJSKiWnK3N8d307vrr/y99WgqnvzsMJKz86WORkQGxLJERPQIlPK7V/7+9/Nd4WipwtmsfAxbdRjfHU/jNZmIGgiWJSIiA+jt7YRf5/ZGr1aOKNKV4bUfTuLV7xJQWFIqdTQiekQsS0REBtLESoN/P98VCwa2gVwmYOeJDAxffRjnr/FjOSJTxrJERGRAMpmAl4Nb4dsXHoOztRrJ2QUYvvoQvj+eJnU0IqolliUiojrQ1csev8zpjd7ejijWlWMBP5YjMlksS0REdcTRUo2Nz939WE4mADvi0vHUZ9G4eL1A6mhEVAMsS0REdejex3JbX3gMTlZqnLuWj+GrDuHnk1eljkZE1cSyRERUDx5r4YBf5vTCYy3scaekDLO2nkDYj0koKeW95YiMHcsSEVE9aWKlweap3TAzqCUAYEP0ZYz+4ggybhdJnIyIqsKyRERUjxT/uYjl+smBsDFTIj7tNob+6yAOXbghdTQiegCWJSIiCfRr54yfZ/eCT1Nr5BTqMOnro/gsMhnl5bzqN5GxMYqydODAAQwbNgxubm4QBAG7d++u8LwoiggLC4ObmxvMzMwQFBSEpKSkCutotVrMnj0bjo6OsLCwwPDhw5Genl6P74KIqGbc7c3xw4weGBPojnIR+PD3c3hx03HkFumkjkZE/8MoytKdO3fg5+eH1atX3/f55cuXY8WKFVi9ejViYmLg4uKCAQMGID//v1fFDQ0Nxa5du7Bt2zYcOnQIBQUFGDp0KMrKeAdwIjJeGqUcH4zqiGUjfaFSyLD/TDaGrz6EM5l5Ukcjov8wirI0ePBgvPfeexg5cmSl50RRxMqVK7F48WKMHDkSPj4+2LhxIwoLC7F161YAQG5uLtavX4+PP/4Y/fv3h7+/PzZv3oxTp05h//799f12iIhqbGxXD+yY0QNNbc1w5WYhRn4ejR8TeHkBImOgkDrAw6SkpCArKwshISH6ZWq1Gn379kV0dDSmT5+O2NhY6HS6Cuu4ubnBx8cH0dHRGDhw4H23rdVqodVq9Y/z8u7+TU6n00Gn427wR3Fvfpzjo+EcDccUZtnW2Ry7XuqGV747hcMXb2LOtycQf+UWFoR4QyE3ir/bAjCNWZoCztFwdDodlEplnW3f6MtSVlYWAMDZ2bnCcmdnZ1y5ckW/jkqlgp2dXaV17r3+fsLDw7FkyZJKyyMjI2Fubv6o0QlARESE1BEaBM7RcExhlqOcALMiGfZfleHr6Cs4mJiCKa3LYVl3vwtqxRRmaQo4R8MYMWJEnW3b6MvSPYIgVHgsimKlZf/0sHVef/11zJs3T/84Ly8P7u7uCA4OhoODw6MFbuR0Oh0iIiIwYMCAOm37DR3naDimNsuhAPYmXcPCnYm4kAesvmCOz8b5wbepjdTRTG6WxopzNJy63jtn9GXJxcUFwN29R66urvrl2dnZ+r1NLi4uKCkpQU5OToW9S9nZ2ejRo8cDt61Wq6FWqystVyqV/ME1EM7SMDhHwzGlWQ7r1AxtXW3w4qZYpNy4g3FfxWDZ0754yr+Z1NEAmNYsjRnnaPyM50PwB/Dy8oKLi0uF3ZQlJSWIiorSF6GAgAAolcoK62RmZiIxMbHKskREZOy8na2wZ1ZP9GvbBNrScryyPQHv/XwapWW8TQpRfTGKPUsFBQVITk7WP05JSUF8fDzs7e3h4eGB0NBQLF26FN7e3vD29sbSpUthbm6O8ePHAwBsbGwwdepUvPrqq3BwcIC9vT3mz58PX19f9O/fX6q3RURkENYaJdZNCsQn+89j1Z/J+OpQCs5k5WH1uM6ws1BJHY+owTOKsnT8+HEEBwfrH987jmjy5MnYsGEDXnvtNRQVFWHmzJnIyclBt27dsG/fPlhZWelf88knn0ChUGD06NEoKipCv379sGHDBsjl8np/P0REhiaTCXg1pA3au1rj1e8TcDj5JoZ/dgjrJgWirYu11PGIGjSjKEtBQUEQxQdf4l8QBISFhSEsLOyB62g0GqxatQqrVq2qg4RERMZhsK8rWjhZ4oV/H0fqrbvXY1oxuhMG+bhIHY2owTL6Y5aIiKiiNi5W+HFWT/Rq5YjCkjLM2ByLT/df4H3liOoIyxIRkQmyNVdhw3Nd8FzP5gCAT/afx6xv41BYUiptMKIGiGWJiMhEKeQyvD2sA5Y/3RFKuYBfT2Xh6TVHkJ5TKHU0ogaFZYmIyMSN7uKOb194DI6WKpzJzMOTnx1G7JVbUsciajBYloiIGoDA5vbYM6sX2rta40ZBCcZ9eRQ7YtOljkXUILAsERE1EE1tzfDDS90xsIMzSsrK8er3CVj221ke+E30iFiWiIgaEHOVAmsmBGBWcCsAwNqoi3hxUyzuaHngN1FtsSwRETUwMpmA+QPbYOWYTlApZNh/5hpGrT2Cq7eLpI5GZJJYloiIGqgn/Zti24v/PfB7xGeHkZB2W+pYRCaHZYmIqAHr7GGH3S/3RBtnK1zP12LMl0fw26lMqWMRmRSWJSKiBq6ZnTl+eKk7gto4oVhXjpe2xOGzyOQqbzNFRP/FskRE1AhYaZT4alIgpvRoDgD48PdzWPDDSZSUlksbjMgEPFJZ0ul0SEtLw7lz53DrFi+ARkRkzBRyGcKGd8C7IzpAJgA/xKZjyjfHkFukkzoakVGrcVkqKCjAF198gaCgINjY2KB58+Zo3749nJyc4OnpiRdeeAExMTF1kZWIiAxgYvfmWD+lCyxUckRfvImn10Qj7RZvkUL0IDUqS5988gmaN2+OdevW4fHHH8fOnTsRHx+Pc+fO4ciRI3j77bdRWlqKAQMGYNCgQbhw4UJd5SYiokcQ3KYJvp/RAy7WGiRnF+DJzw4jLjVH6lhERklRk5Wjo6MRGRkJX1/f+z7ftWtXPP/881i7di3Wr1+PqKgoeHt7GyQoEREZVns3a+x+uSemboxB0tU8jPvyb6wc0wmDfV2ljkZkVGq0Z+n7779/YFH6X2q1GjNnzsS0adNqHYyIiOqei40G303vjn5tm0BbWo6ZW+Pw1cFLUsciMio8G46IqJGzUCvw5aRATOruCVEE3vvlDMJ+TEIZ7ylHBKAWZSkqKgpjx46973Pp6eno2rUrDh8+/MjBiIio/shlApYM74A3nmgLANgQfRkvbY5FUUmZxMmIpFfjshQWFobx48cDAJYtW4bi4mL9c82aNcOsWbPw1ltvGS4hERHVC0EQ8GKfllg93h8qhQz7Tl/DuHV/40aBVupoRJKqcVmKjY3VH7e0dOlSXL16FR999JG+NPXp0wfHjx83bEoiIqo3Qzu6Ycu0brA1VyI+7TaeXhONyzfuSB2LSDI1Lkv29vbIysrSPxZFEe+88w6uXr0KADh9+jRsbGwMl5CIiOpdl+b22PFSD7jbm+HKzUI8vSYa8bwJLzVSNS5LgwYNwjvvvAOd7u4VXwVBgCiKEAQBiYmJmDt3rv5jOiIiMl0tnSyx46Ue8GlqjZt3SjDuy7/xx5lrUsciqne1OmYpJSUFzZs3R3FxMd59911otVo88cQT6NSpE7p27Yp33323LrISEVE9a2KlwbYXu6NPaycU6crwwr+PY/vxdKljEdWrGpclFxcXHD16FAsXLgQAnDlzBqIoIiQkBAcPHsSWLVugUqkMHpSIiKRhqVZg/eRAjApohnIReHPPafyaJoMo8tIC1DjU6jpLNjY2mDNnDnx8fLBlyxb4+vpi7ty56N69u6HzERGREVDKZfhwVEfMebwVAOD3dBne+vE0SsvKJU5GVPdqVJZSU1MrPI6Li0PLli0RFxeHFi1aVFo/IyPj0dIREZHREAQB80LaYMmwdhAgYvvxDMzYHIdiHa/FRA1bjcpSly5d8MILL+DYsWMPXCc3Nxfr1q2Dj48Pdu7c+cgBiYjIuIzv6o7nWpdDpZBh/5lrmPDVUdwuLJE6FlGdqdGNdM+cOYOlS5di0KBBUCqVCAwMhJubGzQaDXJycnD69GkkJSUhMDAQH374IQYPHlxXuYmISEJ+DiL69w7AjC0nEHslB8+sPYKNz3eFm62Z1NGIDK5Ge5bs7e3x0Ucf4erVq1izZg1at26NGzdu4MKFCwCACRMmIDY2FocPH2ZRIiJq4Lo0t8P3M3rAxVqDC9kFeHpNNJKz86WORWRwNdqzdI9Go8HIkSMxcuRIQ+chIiIT0sbFCjtm9sDE9Udx6fodjFp7BN9M6QJ/DzupoxEZTK3OhiMiIrqnqa0ZfpjRA37NbHC7UIfx644i6vx1qWMRGQzLEhERPTJ7CxW2vvAYens7okhXhmkbY7AnnmdEU8PAskRERAZhoVZg/eQuGNrRFboyEXO3xWNj9GWpYxE9MpYlIiIyGJVChn+N9cek7p4AgLd/TMKn+y/wat9k0kyiLJWWluLNN9+El5cXzMzM0KJFC7zzzjsoL//vlWNFUURYWBjc3NxgZmaGoKAgJCUlSZiaiKhxkskELBneAXP7eQMAPtl/Hkt+Oo3ychYmMk21Lku5ubl48cUX0apVK7Rr1w6ZmZmGzFXBBx98gLVr12L16tU4c+YMli9fjg8//BCrVq3Sr7N8+XKsWLECq1evRkxMDFxcXDBgwADk5/M0ViKi+iYIAl4Z0BpvD2sPANgQfRnzv0+AjrdHIRNU67I0c+ZMnDp1CsuXL8eVK1dQVFQEAAgNDcWnn35qsIAAcOTIEYwYMQJDhgxB8+bNMWrUKISEhOD48eMA7u5VWrlyJRYvXoyRI0fCx8cHGzduRGFhIbZu3WrQLEREVH3P9fTCitF+kMsE7DyRgZc2x/L2KGRyanWdJQD47bff8Mcff8Df3x9yuVy/fNCgQXjzzTcxd+5cgwQEgF69emHt2rU4f/48WrdujYSEBBw6dAgrV64EAKSkpCArKwshISH616jVavTt2xfR0dGYPn36fber1Wqh1Wr1j/Py8gAAOp0OOp3OYPkbo3vz4xwfDedoOJyl4dR0lsN8nWGu8MOc705i/5lsTP76KNZO8Ieluta/ghoE/kwajk6ng1KprLPtP9JPqqWlZaVl3t7eSE5OfpTNVrJw4ULk5uaibdu2kMvlKCsrw/vvv49x48YBALKysgAAzs7OFV7n7OyMK1euPHC74eHhWLJkSaXlkZGRMDc3N+A7aLwiIiKkjtAgcI6Gw1kaTk1n+WJrAevOyXA0JQcjPvkDM9qVwaLufr+ZDP5MGsaIESPqbNu1LktPPPEEtm7dirfffrvC8oKCAgiC8MjB/tf27duxefNmbN26FR06dEB8fDxCQ0Ph5uaGyZMn69f75/cVRbHKLK+//jrmzZunf5yXlwd3d3cEBwfDwcHBoO+hsdHpdIiIiMCAAQPqtO03dJyj4XCWhvMoswxOz8XUf8ch9Y4OG1Jt8c2UADSxUtdRUuPGn0nDqeu9c7UuS+Hh4QgMDATw31JSVFSEd955B507dzZYQABYsGABFi1ahLFjxwIAfH19ceXKFYSHh2Py5MlwcXEBcHcPk6urq/512dnZlfY2/S+1Wg21uvJ/pEqlkj+4BsJZGgbnaDicpeHUZpYBXo74bkZ3PPvVUZzPLsD49THYPLUb3O0b7958/kwav1of4O3u7o7Dhw/j0KFDKCwsRNeuXWFnZ4eDBw/igw8+MGRGFBYWQiarGFUul+svHeDl5QUXF5cKuzJLSkoQFRWFHj16GDQLERE9mtbOVvhhRg+425vhys1CjP7iCJKzC6SORfRAj3TMUqtWrRAREYHU1FQkJCRAqVSiW7dusLMz7A0Uhw0bhvfffx8eHh7o0KEDTpw4gRUrVuD5558HcPfjt9DQUCxduhTe3t7w9vbG0qVLYW5ujvHjxxs0CxERPToPB3N8P70Hnl1/FMnZBRjzxRFsntYN7VytpY5GVEmty9KhQ4dgZWUFPz8/eHh4wMPDw5C5Kli1ahXeeustzJw5E9nZ2XBzc8P06dPxf//3f/p1XnvtNRQVFWHmzJnIyclBt27dsG/fPlhZWdVZLiIiqj0XGw2+m94dk74+isSMPIz98m9sfL4rOrnbSh2NqIJafww3e/ZsxMbGVlp+/vx55ObmPlKof7KyssLKlSv113O6ePEi3nvvPahUKv06giAgLCwMmZmZKC4uRlRUFHx8fAyag4iIDMveQoUt0x5DZw9b5BbpMGHd3zh66abUsYgqqHVZOnfuHPr27VtpeWRkpP6UfiIiooexMVNi09Ru6N7CAXdKyjD5m2OIOn9d6lhEerUuS9bW1rh161al5b1798axY8ceKRQRETUuFmoFvnmuC4LbOKFYV44XNh7HvqQsqWMRAXiEsjR8+HB89NFHlTcok6GkpOSRQhERUeOjUcrxxcRADPZxQUlZOWZuicPPJ69KHYuo9mVp6dKlOHjwIPr374+TJ08CAIqLi/HBBx+gY8eOBgtIRESNh0ohw6px/niykxtKy0XM+fYEdp1IlzoWNXK1PhvO0dERR44cwUsvvYROnTpBrVajtLQUNjY2+OmnnwyZkYiIGhGFXIaPR3eCSiHDd8fTMe+7BJSUlmNMl7o765qoKjUqS88++yzWrVsHMzMzpKWlwdPTE7/++ivS0tJw4sQJ/XWW7O3t6yovERE1AnKZgGUjO0KlkGHz36lYuOMUSkrLMbF7c6mjUSNUo7JkaWkJrVYLMzMzeHp6ws7ODn5+fvDz80OnTp3g5+fH6xoREZFByGQC3h3hA7VCjvWHUvDWniRoS8sxrXcLqaNRI1OjsrR27Vr9P6ekpCA+Ph4JCQmIj4/Hjz/+iMuXL0OhUKBt27ZISEgweFgiImpcBEHAm0PaQa2Q4fO/LuK9X86gtFzEjL4tpY5GjUitjlkqKyvD8ePHERISghEjRuiX5+fnIz4+Xn/ANxER0aMSBAELBraBSiHDyv0XsOy3sygtK8esx72ljkaNRK3Kklwux7PPPoukpKQKH7tZWVmhd+/e6N27t8ECEhERCYKA0P6tIRcEfBxxHh/tO4/SchGh/VtLHY0agVpfOqBr165ISUkxZBYiIqIqze7njYWD2gIAVu6/gI/3nYMoihKnooau1mVpzpw5eOONN5CWlmbIPERERFV6Kagl3hzSDgCw6s9kLP+dhYnqVq2vs/TMM88AADp06IDhw4cjKCgI/v7+8PX1rXCDWyIiIkOb1rsF5DIBS346jTV/XYQoAgsHtYEgCFJHowao1mXp3tlw986ICw8Px+XLlyGXy9G2bVse5E1ERHXquZ5ekAkC3v4xCWujLkKEiEWD2rIwkcHVqixduXIFJ0+ehKurK8+GIyIiyUzu0RyCAPzfniR8EXUJEIFFg1mYyLBqXJa+/fZbTJkyBTqdDoIgwN/fH7/99hucnJx4NhwREdW7Sd2bQwDw1p4kfHHgEkQAr7MwkQHV+ADvJUuWYOLEiUhOTsaff/4JmUyGRYsW1UU2IiKiapnYvTnefdIHAPDlgUtY+usZHvRNBlPjPUuXLl3C77//Dk9PT7Ro0QKbNm1C586dsX79+rrIR0REVC0TH/OETAAW70rEuoMpkAkCP5Ijg6jxnqXS0lKYmZnpH7dp0wbl5eXIysoyaDAiIqKamtDNU7+H6YsDl3hZATKIWl1naePGjYiOjkZBQQEAQKFQoLCw0KDBiIiIamPiY554Z0QHAMCavy7iI164kh5RjT+G69WrF9577z3k5+dDJpPBy8sLxcXFWL9+Pfr374+AgABYW1vXRVYiIqJqmdS9OcrLRYT9dBqfRV6EXBAwL6SN1LHIRNW4LB04cAAAcOHCBcTGxiIuLg6xsbFYs2YNwsPDIZPJ4O3tjTNnzhg8LBERUXVN6emFMhF49+fT+NefyZDJBN5Ljmql1hel9Pb2hre3N8aOHatflpKSguPHj+PEiRMGCUdERPQopvbygiiKeO+XM1i5/wKUchleDm4ldSwyMbUuS/fj5eUFLy8v/a1QiIiIpDatdwuUlotY9ttZfPj7OSjlAl7s01LqWGRCan0jXSIiIlMxo29LvDrg7kdwS389i68PpUiciEwJyxIRETUKs/t5Y87jdz+Ce+fn09h05LK0gchksCwREVGj8cqA1ngp6O5HcG/tScK3x1IlTkSmgGWJiIgaDUEQ8NrANpjWywsA8MauU9h1Il3iVGTsWJaIiKhREQQBi4e0w6TunhBF4NXvEvDLyUypY5ERY1kiIqJGRxAEhA3rgDGB7igXgbnbTmBfEm/bRffHskRERI2STCZg6UhfPNnJDaXlImZtPYGo89eljkVGiGWJiIgaLblMwEfP+OEJXxeUlJXjxX8fR/TFG1LHIiPDskRERI2aQi7DyjH+6N+uCbSl5Zi28Thir+RIHYuMCMsSERE1eiqFDKvHd0Zvb0cUlpRhyjfHkJiRK3UsMhIsS0RERAA0Sjm+mBiAQE875BeXYtLXx5CcnS91LDICJlOWMjIy8Oyzz8LBwQHm5ubo1KkTYmNj9c+LooiwsDC4ubnBzMwMQUFBSEpKkjAxERGZGnOVAl8/1wW+TW1w604Jxq87iis370gdiyRmEmUpJycHPXv2hFKpxG+//YbTp0/j448/hq2trX6d5cuXY8WKFVi9ejViYmLg4uKCAQMGID+ffysgIqLqs9YosfH5rmjtbInsfC3GrzuKq7eLpI5FEjKJsvTBBx/A3d0d33zzDbp27YrmzZujX79+aNny7iXrRVHEypUrsXjxYowcORI+Pj7YuHEjCgsLsXXrVonTExGRqbG3UGHz1G5o7mCOjNtFeHb9Udwo0EodiySikDpAdfz4448YOHAgnnnmGURFRaFp06aYOXMmXnjhBQBASkoKsrKyEBISon+NWq1G3759ER0djenTp993u1qtFlrtf3/48/LyAAA6nQ46na4O31HDd29+nOOj4RwNh7M0nMYySzszOTZOCcDYr2Jw6fodTPzqKDY/HwhrM6VBtt9Y5lgfdDodlErD/Hu5H0EURbHOtm4gGo0GADBv3jw888wzOHbsGEJDQ/HFF19g0qRJiI6ORs+ePZGRkQE3Nzf961588UVcuXIFv//++323GxYWhiVLllRavnXrVpibm9fNmyEiIpOSXQR8miRHgU6Al5WIl9qVQS2XOhX904gRI+ps2yZRllQqFQIDAxEdHa1fNmfOHMTExODIkSP6snT16lW4urrq13nhhReQlpaGvXv33ne799uz5O7ujszMTDg4ONTdG2oEdDodIiIiMGDAgDpt+w0d52g4nKXhNMZZnsnMx4SvY5BfXIperRywdoI/1IpHO5KlMc6xruh0ujrdyWESH8O5urqiffv2FZa1a9cOO3bsAAC4uLgAALKysiqUpezsbDg7Oz9wu2q1Gmq1utJypVLJH1wD4SwNg3M0HM7ScBrTLDt62GPDc13w7FfHcCj5JhbsSMSqcf5QyB/90N/GNEdTZRIHePfs2RPnzp2rsOz8+fPw9PQEAHh5ecHFxQURERH650tKShAVFYUePXrUa1YiImqYAjztsW5SIFRyGX5LzMLrO0/BBD6cIQMwibL0yiuv4O+//8bSpUuRnJyMrVu34ssvv8TLL78M4O7do0NDQ7F06VLs2rULiYmJmDJlCszNzTF+/HiJ0xMRUUPRy9sRq8b7Qy4T8H1sOt7/5QwLUyNgEmWpS5cu2LVrF7799lv4+Pjg3XffxcqVKzFhwgT9Oq+99hpCQ0Mxc+ZMBAYGIiMjA/v27YOVlZWEyYmIqKEZ2MEFHzzdEQDw1aEUfP7XRYkTUV0ziWOWAGDo0KEYOnToA58XBAFhYWEICwurv1BERNQojQpohtwiHd79+TQ+/P0crM2UmPiYp9SxqI6YxJ4lIiIiYzO1lxfmPN4KAPB/exKxJz5D4kRUV1iWiIiIaumVAa0xqbsnRBF49bsERJ7NljoS1QGWJSIioloSBAFhwzpgRCc3lJaLeGlLLGKv3JI6FhkYyxIREdEjkMkEfPSMH4LaOKFYV47nvonBuSzexL0hYVkiIiJ6REq5DJ9P6IzOHrbIKy7FpK+PIu1WodSxyEBYloiIiAzAXKXA11O6oLWzJa7laTFx/VHcKNA+/IVk9FiWiIiIDMTWXIV/P98NTW3NcPlmISZ/fQz5xTqpY9EjYlkiIiIyIBcbDTZN7Qp7CxWSrubhxX/HQltaJnUsegQsS0RERAbWwskSG5/rCguVHEcu3cQr2+NRVs7bopgqliUiIqI64NvMBl9MDIRSLuDXU1lY8lMS7yNnoliWiIiI6kgvb0esGN0JggD8+8gVfBaZLHUkqgWWJSIiojo0zM8Nbw9tDwD4aN95bDuWKnEiqimWJSIiojo2pacXXg5uCQB4Y9cp7EvKkjgR1QTLEhERUT2YH9IGowOboVwEZn97AnGpt6WORNXEskRERFQPBEHA0qd80a9tE2hLyzF98wlcK5I6FVUHyxIREVE9UchlWDXeH53cbXG7SIc1p+W4llcsdSx6CJYlIiKienTvtiheDubIKREw7d9xyONVvo0ayxIREVE9s7dQYf3kzrBWijh7rQAv/vs4r/JtxFiWiIiIJOBuZ47p7cpgoZbj70u38Op3CSjnVb6NEssSERGRRJpZAJ+N6wSlXMDPJzOxbO9ZqSPRfbAsERERSahnSwcsH9URAPDlgUv45nCKxInon1iWiIiIJPaUfzO8NqgNAOCdn0/jt1OZEiei/8WyREREZARe6tsSzz7mAVEE5m6Px/HLt6SORP/BskRERGQEBEHAkuE+6N/OGSWl5Zi68TiSswukjkVgWSIiIjIacpmAVePuXrQyt0iHKd8cQ3Y+L1opNZYlIiIiI2KmkmP95EA0dzBHek4Rpm08jsKSUqljNWosS0REREbGwVKNDc91hZ25EifTczF76wmUlpVLHavRYlkiIiIyQs0dLfDV5C5QK2T442w2lvx0GqLIi1ZKgWWJiIjISAV42mHlmE4QBGDT31ew7uAlqSM1SixLRERERmywryveHNIeALD017P4+eRViRM1PixLRERERm5qLy9M6dEcADDvuwReg6mesSwRERGZgLeGtseA9nevwfTCv4/j8o07UkdqNFiWiIiITIBcJuDTsZ3QsZkNcgp1eG5DDHLulEgdq1FgWSIiIjIR5ioFvpociKa2Zki5cQfTN8VCW1omdawGj2WJiIjIhDSx0uCb57rASq3Ascu38NoPJ3lJgTpmkmUpPDwcgiAgNDRUv0wURYSFhcHNzQ1mZmYICgpCUlKSdCGJiIjqSGtnK6x5NgAKmYA98VexIuK81JEaNJMrSzExMfjyyy/RsWPHCsuXL1+OFStWYPXq1YiJiYGLiwsGDBiA/Px8iZISERHVnV7ejlj6lC8AYNWfydgRmy5xoobLpMpSQUEBJkyYgHXr1sHOzk6/XBRFrFy5EosXL8bIkSPh4+ODjRs3orCwEFu3bpUwMRERUd0Z3cUdLwW1BAAs2nkSRy/dlDhRw6SQOkBNvPzyyxgyZAj69++P9957T788JSUFWVlZCAkJ0S9Tq9Xo27cvoqOjMX369PtuT6vVQqvV6h/n5eUBAHQ6HXQ6XR29i8bh3vw4x0fDORoOZ2k4nKVhGGqOocEtkHK9AHuTrmH6plh8P70rmjtYGCKiydDpdFAqlXW2fZMpS9u2bUNcXBxiYmIqPZeVlQUAcHZ2rrDc2dkZV65ceeA2w8PDsWTJkkrLIyMjYW5u/oiJCQAiIiKkjtAgcI6Gw1kaDmdpGIaYYz8L4IylHFcKdBi/9hBe8SmDRd11B6M0YsSIOtu2SZSltLQ0zJ07F/v27YNGo3ngeoIgVHgsimKlZf/r9ddfx7x58/SP8/Ly4O7ujuDgYDg4ODx68EZMp9MhIiICAwYMqNO239BxjobDWRoOZ2kYhp5jr2AtRn1xFBm3i7HrhhM2TA6ASmFSR9vUWl3v5TSJshQbG4vs7GwEBATol5WVleHAgQNYvXo1zp07B+DuHiZXV1f9OtnZ2ZX2Nv0vtVoNtVpdablSqeT/AAyEszQMztFwOEvD4SwNw1BzdLVT4uspXTFqTTRiLufg/346i4+e6VjlTgOqHpOonP369cOpU6cQHx+v/woMDMSECRMQHx+PFi1awMXFpcKuzJKSEkRFRaFHjx4SJiciIqo/bVys8NmEzpDLBOyIS8eaqItSR2oQTGLPkpWVFXx8fCoss7CwgIODg355aGgoli5dCm9vb3h7e2Pp0qUwNzfH+PHjpYhMREQkiT6tnRA2vAPe2p2I5XvPoYWjBQb5uD78hfRAJlGWquO1115DUVERZs6ciZycHHTr1g379u2DlZWV1NGIiIjq1cTHPHExuwAboi8jdHs8vrc1h28zG6ljmSyTLUt//fVXhceCICAsLAxhYWGS5CEiIjImbw5ph8s37+Cvc9cxdWMM9szqCVcbM6ljmSSTOGaJiIiIakYhl2HVOH+0drZEdr4WUzccxx1tqdSxTBLLEhERUQNlpVFi/eQucLBQ4XRmHl7ZHo/yct50t6ZYloiIiBowd3tzfDkpACq5DPtOX8NH+85JHcnksCwRERE1cAGe9vhg1N2b7n7+10XsjONNd2uCZYmIiKgReMq/GWbeu+nujlOIvZIjcSLTwbJERETUSMwPaYOQ9s4oKSvH9E3HkZ5TKHUkk8CyRERE1EjIZAI+GdMJ7VytcaOgBNM28gy56mBZIiIiakQs1Ap8NTkQjpZqnM3KRyjPkHsoliUiIqJGpqmtGb6YGACVQoaI09fwcQTPkKsKyxIREVEjFOBph2Uj754h91nkReyJz5A4kfFiWSIiImqkRnZuhul9WgAAXvvhJBLSbksbyEixLBERETVirw1qi8fbNoG2tBwvbjqOa3nFUkcyOixLREREjZhcJuDTsZ3g3cQS1/K0eHFTLIp1ZVLHMiosS0RERI2clUaJryYHwtZciYS021i04yREkWfI3cOyRERERPB0sMDnEzpDLhOwO/4qvjxwSepIRoNliYiIiAAAPVo64u1h7QEAy/aeReS5bIkTGQeWJSIiItKb+JgnxnZxhygCc749gYvXC6SOJDmWJSIiItITBAHvjPBBoKcd8otL8cLG48gt0kkdS1IsS0RERFSBSiHDmmcD4GajwaUbdzDn2xMoa8S3RGFZIiIiokqcrNT4clIg1AoZos5fx/K9Z6WOJBmWJSIiIrovn6Y2+PAZPwDAFwcuNdpborAsERER0QMN93PDjL4tAQALd5xEYkauxInqH8sSERERVWnBwDbo29oJxbpyTN8Ui5sFWqkj1SuWJSIiIqqSXCbgX2P90dzBHBm3i/Dy1jjoysqljlVvWJaIiIjooWzMlVg3KRAWKjn+vnQL7/9yRupI9YZliYiIiKrF29kKn4zpBADYEH0Z3x1PkzZQPWFZIiIiomoL6eCC0P7eAIA3dyUiPu22tIHqAcsSERER1cicx70xoL0zSsrKMWNTLK7nN+wDvlmWiIiIqEZkMgErRvuhhZMFsvKKG/wB3yxLREREVGNWGiW+nBgIS7UCx1Ia9gHfLEtERERUK62aWFY44HtHbLq0geoIyxIRERHV2oD2zpjb7+4B32/sOoVT6Q3vCt8sS0RERPRI5vbzRv92TaAtLcf0Tccb3BW+WZaIiIjokchkAlaM6YQWjha4mluMOdtOoLQBHfDNskRERESPzFqjxBcTA2CukuNw8k18uO+c1JEMxiTKUnh4OLp06QIrKys0adIETz75JM6dq/gvQRRFhIWFwc3NDWZmZggKCkJSUpJEiYmIiBofb2crfDjKDwDwRdQl/HIyU+JEhmESZSkqKgovv/wy/v77b0RERKC0tBQhISG4c+eOfp3ly5djxYoVWL16NWJiYuDi4oIBAwYgPz9fwuRERESNy5COrpjepwUAYMEPCTh/zfR/D5tEWdq7dy+mTJmCDh06wM/PD9988w1SU1MRGxsL4O5epZUrV2Lx4sUYOXIkfHx8sHHjRhQWFmLr1q0SpyciImpcFgxsgx4tHVBYUobpm2KRV6yTOtIjMYmy9E+5uXdPS7S3twcApKSkICsrCyEhIfp11Go1+vbti+joaEkyEhERNVYKuQyrxvnDzUaDlBt3MG97AsrLRalj1ZpC6gA1JYoi5s2bh169esHHxwcAkJWVBQBwdnausK6zszOuXLnywG1ptVpotf89vTEvLw8AoNPpoNOZdguW2r35cY6PhnM0HM7ScDhLw2joc7RWy7B6nB/GfhWD/WeuYfWf5/FS3xZ18r10Oh2USmWdbBswwbI0a9YsnDx5EocOHar0nCAIFR6Lolhp2f8KDw/HkiVLKi2PjIyEubn5o4clRERESB2hQeAcDYezNBzO0jAa+hxHegjYdkmOT/ZfQFHGObS1rZs9TCNGjKiT7QKAIIqiyewXmz17Nnbv3o0DBw7Ay8tLv/zSpUto2bIl4uLi4O/vr18+YsQI2NraYuPGjffd3v32LLm7uyMzMxMODg5190YaAZ1Oh4iICAwYMKBO235DxzkaDmdpOJylYTSmOS7enYTvYjNgZ67ErpceQ1NbM4NuX6fT1elODpPYsySKImbPno1du3bhr7/+qlCUAMDLywsuLi6IiIjQl6WSkhJERUXhgw8+eOB21Wo11Gp1peVKpbLB/+DWF87SMDhHw+EsDYezNIzGMMd3nvTFmawCnMrIxZztJ/Hd9O7QKOVSx6o2kzjA++WXX8bmzZuxdetWWFlZISsrC1lZWSgqKgJw9+O30NBQLF26FLt27UJiYiKmTJkCc3NzjB8/XuL0REREjZtGKceaZzvD1lyJk+m5WPLTaakj1YhJlKU1a9YgNzcXQUFBcHV11X9t375dv85rr72G0NBQzJw5E4GBgcjIyMC+fftgZWUlYXIiIiICgGZ25vjXWH8IAvDtsVR8F5MmdaRqM5mP4R5GEASEhYUhLCys7gMRERFRjfVp7YRXB7TGR/vO4809iWjvZg2fpjZSx3ook9izRERERA3DzKBW6Ne2CUpKy/HSlljkFhr/pRNYloiIiKjeyGQCVozuBHd7M6TdKsK87+KN/oKVLEtERERUr2zMlVgzIQAqhQx/nM3GmqiLUkeqEssSERER1TufpjZ4d0QHAMDH+87h0IUbEid6MJYlIiIiksSYLh4YHdgM5SIwZ9sJZOYWSR3pvliWiIiISDLvjPBBe1dr3LpTgplb4lBSWi51pEpYloiIiEgyGqUca58NgJVGgROptxH+2xmpI1XCskRERESS8nAwx4rRnQAA3xy+jF9PZUob6B9YloiIiEhyA9o7Y3rfFgCA1344iUvXCyRO9F8sS0RERGQUFoS0QVcvexRoSzFzSxyKSsqkjgSAZYmIiIiMhEIuw+px/nC0VONsVj7e2pModSQALEtERERkRJpYa/CvcZ0gE4AfYtOxPSZV6kgsS0RERGRcerR0xKshbQAA/7cnCaev5kmah2WJiIiIjM5LfVsiuI0TtKXleHlrHPKLpbvhLssSERERGZ17N9x1s9Eg5cYdLNp5CqIozQ13WZaIiIjIKNlZqLB6QmcoZAJ+OZmJzX9fkSQHyxIREREZrc4edlg0uC0A4N2fz+Bk+u16z8CyREREREZtai8vhLR3RklZOWZuiUNuYf0ev8SyREREREZNEAR8+Iwf3O3NkJ5ThPk/JNTr8UssS0RERGT0bMyU+Hx8AFRyGSJOX8PXhy/X2/dmWSIiIiKT4NvMBm8ObQcAWPbbGcSn3a6X78uyRERERCZj4mOeeMLXBboyEbO21s/xSyxLREREZDIEQcCypzvCw94c6TlFWFAPxy+xLBEREZFJsdYo8dn4zlDJZdh3+ho2/l23949jWSIiIiKT49vMBouH3D1+afnv5+v0e7EsERERkUma1N0Tg33uHr9Ul1iWiIiIyCQJgoAPRnWEu51ZnX4fliUiIiIyWdYaJb6d1qVOvwfLEhEREZk0Z2tNnW6fZYmIiIioCixLRERERFVgWSIiIiKqAssSERERURVYloiIiIiqwLJEREREVAWWJSIiIqIqNLiy9Pnnn8PLywsajQYBAQE4ePCg1JGIiIjIhDWosrR9+3aEhoZi8eLFOHHiBHr37o3BgwcjNbVu70ZMREREDVeDKksrVqzA1KlTMW3aNLRr1w4rV66Eu7s71qxZI3U0IiIiMlEKqQMYSklJCWJjY7Fo0aIKy0NCQhAdHX3f12i1Wmi1Wv3j3NxcAMCtW7fqLmgjodPpUFhYiJs3b0KpVEodx2RxjobDWRoOZ2kYnKPh6HQ6mJubw8rKCoIgGHz7DaYs3bhxA2VlZXB2dq6w3NnZGVlZWfd9TXh4OJYsWVJpeevWreskIxEREdWd7OxsODk5GXy7DaYs3fPPRimK4gNb5uuvv4558+bpH9++fRuenp5ITU2FjY1NneZs6PLy8uDu7o60tDRYW1tLHcdkcY6Gw1kaDmdpGJyj4dybpUqlqpPtN5iy5OjoCLlcXmkvUnZ2dqW9Tfeo1Wqo1epKy21sbPiDayDW1tacpQFwjobDWRoOZ2kYnKPh1MVHcEADOsBbpVIhICAAERERFZZHRESgR48eEqUiIiIiU9dg9iwBwLx58zBx4kQEBgaie/fu+PLLL5GamooZM2ZIHY2IiIhMVIMqS2PGjMHNmzfxzjvvIDMzEz4+Pvj111/h6elZrder1Wq8/fbb9/1ojmqGszQMztFwOEvD4SwNg3M0nLqepSCKolgnWyYiIiJqABrMMUtEREREdYFliYiIiKgKLEtEREREVWBZIiIiIqoCy9J/fP755/Dy8oJGo0FAQAAOHjwodSSjFx4eji5dusDKygpNmjTBk08+iXPnzlVYRxRFhIWFwc3NDWZmZggKCkJSUpJEiU1DeHg4BEFAaGiofhnnWH0ZGRl49tln4eDgAHNzc3Tq1AmxsbH65znL6iktLcWbb74JLy8vmJmZoUWLFnjnnXdQXl6uX4ezrOzAgQMYNmwY3NzcIAgCdu/eXeH56sxMq9Vi9uzZcHR0hIWFBYYPH4709PR6fBfGoapZ6nQ6LFy4EL6+vrCwsICbmxsmTZqEq1evVtiGwWYpkrht2zZRqVSK69atE0+fPi3OnTtXtLCwEK9cuSJ1NKM2cOBA8ZtvvhETExPF+Ph4cciQIaKHh4dYUFCgX2fZsmWilZWVuGPHDvHUqVPimDFjRFdXVzEvL0/C5Mbr2LFjYvPmzcWOHTuKc+fO1S/nHKvn1q1boqenpzhlyhTx6NGjYkpKirh//34xOTlZvw5nWT3vvfee6ODgIP78889iSkqK+P3334uWlpbiypUr9etwlpX9+uuv4uLFi8UdO3aIAMRdu3ZVeL46M5sxY4bYtGlTMSIiQoyLixODg4NFPz8/sbS0tJ7fjbSqmuXt27fF/v37i9u3bxfPnj0rHjlyROzWrZsYEBBQYRuGmiXLkiiKXbt2FWfMmFFhWdu2bcVFixZJlMg0ZWdniwDEqKgoURRFsby8XHRxcRGXLVumX6e4uFi0sbER165dK1VMo5Wfny96e3uLERERYt++ffVliXOsvoULF4q9evV64POcZfUNGTJEfP755yssGzlypPjss8+KoshZVsc/f8FXZ2a3b98WlUqluG3bNv06GRkZokwmE/fu3Vtv2Y3N/YrnPx07dkwEoN/RYchZNvqP4UpKShAbG4uQkJAKy0NCQhAdHS1RKtOUm5sLALC3twcApKSkICsrq8Js1Wo1+vbty9nex8svv4whQ4agf//+FZZzjtX3448/IjAwEM888wyaNGkCf39/rFu3Tv88Z1l9vXr1wh9//IHz588DABISEnDo0CE88cQTADjL2qjOzGJjY6HT6Sqs4+bmBh8fH871IXJzcyEIAmxtbQEYdpYN6gretXHjxg2UlZVVutmus7NzpZvy0oOJooh58+ahV69e8PHxAQD9/O432ytXrtR7RmO2bds2xMXFISYmptJznGP1Xbp0CWvWrMG8efPwxhtv4NixY5gzZw7UajUmTZrEWdbAwoULkZubi7Zt20Iul6OsrAzvv/8+xo0bB4A/l7VRnZllZWVBpVLBzs6u0jr8nfRgxcXFWLRoEcaPH6+/KbEhZ9noy9I9/7xTsSiKdXb34oZo1qxZOHnyJA4dOlTpOc62amlpaZg7dy727dsHjUbzwPU4x4crLy9HYGAgli5dCgDw9/dHUlIS1qxZg0mTJunX4ywfbvv27di8eTO2bt2KDh06ID4+HqGhoXBzc8PkyZP163GWNVebmXGuD6bT6TB27FiUl5fj888/f+j6tZllo/8YztHREXK5vFLLzM7OrtT+6f5mz56NH3/8EZGRkWjWrJl+uYuLCwBwtg8RGxuL7OxsBAQEQKFQQKFQICoqCv/617+gUCj0s+IcH87V1RXt27evsKxdu3ZITU0FwJ/JmliwYAEWLVqEsWPHwtfXFxMnTsQrr7yC8PBwAJxlbVRnZi4uLigpKUFOTs4D16H/0ul0GD16NFJSUhAREaHfqwQYdpaNviypVCoEBAQgIiKiwvKIiAj06NFDolSmQRRFzJo1Czt37sSff/4JLy+vCs97eXnBxcWlwmxLSkoQFRXF2f6Pfv364dSpU4iPj9d/BQYGYsKECYiPj0eLFi04x2rq2bNnpctXnD9/Xn8zbf5MVl9hYSFksoq/IuRyuf7SAZxlzVVnZgEBAVAqlRXWyczMRGJiIuf6D/eK0oULF7B//344ODhUeN6gs6zR4eAN1L1LB6xfv148ffq0GBoaKlpYWIiXL1+WOppRe+mll0QbGxvxr7/+EjMzM/VfhYWF+nWWLVsm2tjYiDt37hRPnToljhs3rtGfWlwd/3s2nChyjtV17NgxUaFQiO+//7544cIFccuWLaK5ubm4efNm/TqcZfVMnjxZbNq0qf7SATt37hQdHR3F1157Tb8OZ1lZfn6+eOLECfHEiRMiAHHFihXiiRMn9GdoVWdmM2bMEJs1aybu379fjIuLEx9//PFGeemAqmap0+nE4cOHi82aNRPj4+Mr/A7SarX6bRhqlixL//HZZ5+Jnp6eokqlEjt37qw//Z0eDMB9v7755hv9OuXl5eLbb78turi4iGq1WuzTp4946tQp6UKbiH+WJc6x+n766SfRx8dHVKvVYtu2bcUvv/yywvOcZfXk5eWJc+fOFT08PESNRiO2aNFCXLx4cYVfRJxlZZGRkff9/+LkyZNFUazezIqKisRZs2aJ9vb2opmZmTh06FAxNTVVgncjrapmmZKS8sDfQZGRkfptGGqWgiiKYs32RRERERE1Ho3+mCUiIiKiqrAsEREREVWBZYmIiIioCixLRERERFVgWSIiIiKqAssSERERURVYloiIiIiqwLJEREREVAWWJSIiIqIqsCwRUaNy8+ZNNGnSBJcvX37gOqNGjcKKFSvqLxQRGTXe7oSIGpX58+cjJycH69evf+A6J0+eRHBwMFJSUmBtbV2P6YjIGHHPEhE1GkVFRVi/fj2mTZtW5XodO3ZE8+bNsWXLlnpKRkTGjGWJiExScnIyBEHAL7/8gn79+sHc3Bxt2rTB0aNHH/ia3377DQqFAt27d3/o9ocPH45vv/3WkJGJyESxLBGRSUpISIAgCPj444/x5ptvIiEhAR4eHli0aNEDX3PgwAEEBgZWa/tdu3bFsWPHoNVqDRWZiEwUyxIRmaSEhATY2Nhg+/btCA4Ohre3N5588klcv379ga+5fPky3NzcqrX9pk2bQqvVIisry1CRichEsSwRkUlKSEjAsGHD4OTkpF926dIltGrV6oGvKSoqgkaj0T/esmULLC0t9V8HDx7UP2dmZgYAKCwsrIP0RGRKFFIHICKqjYSEBCxcuLDCshMnTqBPnz4PfI2joyNycnL0j4cPH45u3brpHzdt2lT/z7du3QKACmWMiBon7lkiIpOTm5uLK1euwN/fv8Ly+Ph4dOrU6YGv8/f3x+nTp/WPrays0KpVK/3Xvb1JAJCYmIhmzZrB0dHR4PmJyLSwLBGRyUlISIBcLoefn59+2ZUrV5CTk1NlWRo4cCCSkpIq7F16kIMHDyIkJMQQcYnIxLEsEZHJSUhIQNu2bSvsCTpx4gRsbW3RvHnzB77O19cXgYGB+O6776rcfnFxMXbt2oUXXnjBUJGJyITxCt5E1Kj8+uuvmD9/PhITEyGT3f/vi5999hn27NmDffv21XM6IjJGPMCbiBqVJ554AhcuXEBGRgbc3d3vu45SqcSqVavqORkRGSvuWSIiIiKqAo9ZIiIiIqoCyxIRERFRFViWiIiIiKrAskRERERUBZYlIiIioiqwLBERERFVgWWJiIiIqAosS0RERERVYFkiIiIiqgLLEhEREVEV/h+FAsGQMZkQhAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pylab as plt\n", "n = np.linspace(0,120,100)\n", "price = 150 - 0.01 * n**2\n", "plt.plot(n,price)\n", "plt.xlabel('$n$ (-)')\n", "plt.ylabel('$Price$ (€)')\n", "plt.ylim(0,160)\n", "plt.xlim(0,120)\n", "ax = plt.gca()\n", "ax.grid(True, which='both')\n", "ax.spines['left'].set_position('zero')\n", "ax.spines['right'].set_color('none')\n", "ax.yaxis.tick_left()\n", "ax.spines['bottom'].set_position('zero')\n", "ax.spines['top'].set_color('none')\n", "ax.xaxis.tick_bottom()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The profit can be estimated as the revenues minus the total costs." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Model\n", "The function for the profit can be found by combining the relations in the problem statement. However, this is the profit which should be maximized. To turn this into a minimization problem, the profit can be multiplied with $-1$. The final model of this problem results in:\n", "\n", "$$ \\mathop {\\min }\\limits_n \\left(75n - \\left( 150 - 0.01n^2 \\right) n \\right) $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Method\n", "\n", "This model is described using `scipy.optimize.minimize` according to the standard structure in this course.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Importing libraries\n", "For this problem, we'll use all three packages `scipy`, `numpy`, `matplotlib`." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import scipy as sp \n", "import numpy as np\n", "import matplotlib.pylab as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Defining the variables\n", "There are very few variables in this problem. In fact, the only variable we have to specify is the initial guess for the optimization algorithm. The objective function will be treated later.\n", "The length of $n$ doesn't have to be specified." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "n0 = 20" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Defining the objective function\n", "In the objective function, the formula given in the model description can be inserted. Or, each individual step can be calculated on a seperate line. Again, note that the profit is multiplied with $-1$ to maximize the profit in the minimization formulation.\n", "This results in:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def negprofit(n):\n", " price = 150 - 0.01 * n**2\n", " revenues = price * n\n", " totalcost = 75 * n\n", " profit = revenues - totalcost\n", " return -profit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Solving the problem\n", "Now, the problem can be solved. The result is stored in the variables `result` which is printed." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " message: Optimization terminated successfully.\n", " success: True\n", " status: 0\n", " fun: -2499.9999999998727\n", " x: [ 5.000e+01]\n", " nit: 8\n", " jac: [ 0.000e+00]\n", " hess_inv: [[ 3.503e-01]]\n", " nfev: 22\n", " njev: 11\n" ] } ], "source": [ "result = sp.optimize.minimize(negprofit,n0)\n", "print(result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Postprocessing\n", "As this case is only one-dimensional and the potential range of values is limited, we can easily check this answer by an exhaustive search, evaluating all possible values for $n$. The plot below shows the $\\text{negative profit}$ for $0" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n_range = np.linspace(0,100,100)\n", "negprofit_result = negprofit(n_range)\n", "plt.plot(n_range,negprofit_result)\n", "plt.plot(result.x,result.fun,'o');\n", "plt.xlabel('$n$')\n", "plt.ylabel('Negative profit');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 2\n", "\n", "### Problem\n", "8 data points are collected where $X_i$ is the explanatory variable and $Y_i$ the response variable. It is believed that the data follows the trend $Y = aX^b+c$ where $a$, $b$ and $c$ are the parameters to be determined to fit the model.\n", "\n", "| $X_i$ | $Y_i$ |\n", "| :----- | :----- |\n", "| 0.97 | 0.97 |\n", "| 0 | 0.06 |\n", "| 0.5 | 0.7 |\n", "| 0.85 | 0.74 |\n", "| 0.7 | 0.2 |\n", "| 0.19 | 0.34 |\n", "| 0.41 | 0.29 |\n", "| 0.78 | 0.94 |" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAGzCAYAAAArAc0KAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAjA0lEQVR4nO3de3BU5f3H8c/uhiSI7DpcjEAiogZFMmJJhkso9aeROMhombaCgyOXQseM1QipWlKqiOM0FWvQyEVUkLGizUDF0SkiW6wQQNuSho4SxkRBAxJkgrobvBDZfX5/YLbGbDQJ2duz79fM/rEnZ5NvzlD33XOe7HEYY4wAAAAs5oz1AAAAAJFG8AAAAOsRPAAAwHoEDwAAsB7BAwAArEfwAAAA6xE8AADAegQPAACwXkqsB4ilYDCoI0eOqG/fvnI4HLEeBwAAdIIxRs3NzRo8eLCczs6du0nq4Dly5IiysrJiPQYAAOiGQ4cOKTMzs1P7xk3w7NixQw8//LCqq6vV2NioTZs2aerUqd/7mu3bt6ukpET79u3T4MGDdc8996ioqKjTP7Nv376STh8wt9t9JuMDAIAo8fv9ysrKCr2Pd0bcBM/nn3+uUaNGac6cOfr5z3/+g/sfPHhQ1113nX71q1/pueee065du3Tbbbdp4MCBnXq9pNBlLLfbTfAAAJBgurIcJW6CZ/LkyZo8eXKn93/iiSd0/vnn69FHH5UkjRgxQnv27NGf/vSnTgcPAABIDgn7V1pvvvmmCgsL22y79tprtWfPHn399ddhX3Py5En5/f42DwAAYL+EDZ6jR48qIyOjzbaMjAydOnVKTU1NYV9TVlYmj8cTerBgGQCA5JCwwSO1v3ZnjAm7vVVpaal8Pl/ocejQoYjPCAAAYi9u1vB01XnnnaejR4+22Xbs2DGlpKSof//+YV+TlpamtLS0aIwHAADiSMKe4Rk/fry8Xm+bbVu3blVeXp569eoVo6kAAEA8ipvgOXHihPbu3au9e/dKOv1n53v37lVDQ4Ok05ejZs6cGdq/qKhIH374oUpKSrR//36tXbtWa9as0V133RWL8QEAQByLm0tae/bs0VVXXRV6XlJSIkmaNWuW1q1bp8bGxlD8SNKwYcO0efNmLViwQCtWrNDgwYNVUVHBn6QDAIB2HKZ1pW8S8vv98ng88vl8fPAgAAAJojvv33FzSQsAACBSCB4AANAjlnnrVLGtPuzXKrbVa5m3LsoT/Q/BAwAAeoTL6VB5mOip2Favcm+dXM7O3/uqp8XNomUAAJDYiguyJUnl35zJKS7IDsVOyaThoa/HAsEDAAB6zLejZ/nr76klEIx57Ehc0gIAAD2suCBbqS6nWgJBpbqcMY8dieABAAA9rGJbfSh2WgLBDhcyRxOXtAAAQI/57pqd1ueSWMMDAAASX7gFyuEWMscCwQMAAHpEIGjCLlBufR4Ixu7mDtxagltLAACQULi1BAAAQBgEDwAgacTzrQ8QWQQPACBpxPOtDxBZLFoGACSNeL71ASKL4AEAJJV4vfUBIotLWgCApBOPtz5AZBE8AICkE4+3PkBkcUkLAJBU4vXWB4gsggcAkDTi+dYHiCyCBwCQNOL51geILG4twa0lAABIKNxaAgAAIAyCBwAAWI/gAQAA1iN4AACA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgAAYD2CBwAAWI/gAQAA1iN4AACA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgAAYD2CBwAAWI/gAQAA1iN4AACA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgAAYD2CBwAAWI/gAQAA1iN4AACA9QgeAABgvbgKnpUrV2rYsGFKT09Xbm6uqqqqvnf/9evXa9SoUTrrrLM0aNAgzZkzR8ePH4/StAAAIFHETfBUVlZq/vz5WrRokWpqajRx4kRNnjxZDQ0NYfffuXOnZs6cqblz52rfvn3asGGD/v3vf2vevHlRnhwAAMS7uAme8vJyzZ07V/PmzdOIESP06KOPKisrS6tWrQq7/1tvvaULLrhAxcXFGjZsmH784x/r1ltv1Z49e6I8OQAAiHdxETwtLS2qrq5WYWFhm+2FhYXavXt32Nfk5+fr8OHD2rx5s4wx+vjjj7Vx40ZNmTKlw59z8uRJ+f3+Ng8AAGC/uAiepqYmBQIBZWRktNmekZGho0ePhn1Nfn6+1q9fr+nTpys1NVXnnXeezjnnHD3++OMd/pyysjJ5PJ7QIysrq0d/DwAAEJ/iInhaORyONs+NMe22taqtrVVxcbHuu+8+VVdXa8uWLTp48KCKioo6/P6lpaXy+Xyhx6FDh3p0fgAAEJ9SYj2AJA0YMEAul6vd2Zxjx461O+vTqqysTBMmTNDdd98tSbr88svVp08fTZw4UQ8++KAGDRrU7jVpaWlKS0vr+V8AAADEtbg4w5Oamqrc3Fx5vd42271er/Lz88O+5osvvpDT2XZ8l8sl6fSZIQAAgFZxETySVFJSoqefflpr167V/v37tWDBAjU0NIQuUZWWlmrmzJmh/a+//nq9+OKLWrVqlQ4cOKBdu3apuLhYY8aM0eDBg2P1awAAgDgUF5e0JGn69Ok6fvy4HnjgATU2NionJ0ebN2/W0KFDJUmNjY1tPpNn9uzZam5u1vLly/Wb3/xG55xzjq6++mo99NBDsfoVAABAnHKYJL7+4/f75fF45PP55Ha7Yz0OAADohO68f8fNJS0AAIBIIXgAAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgAAYD2CBwAAWI/gAQAA1iN4AACA9QgeAABgPYIHAIAoWeatU8W2+rBfq9hWr2XeuihPlDwIHgAAosTldKg8TPRUbKtXubdOLqcjRpPZLyXWAwAAkCyKC7IlSeXfnMkpLsgOxU7JpOGhr6PnETwAAETRt6Nn+evvqSUQJHaigEtaAABEWXFBtlJdTrUEgkp1OYmdKCB4AACIsopt9aHYaQkEO1zIjJ7DJS0AAKLou2t2Wp9L4kxPBBE8AABESbgFyuEWMqPnETwAAERJIGjCLlBufR4ImliMlRQcxpikPbp+v18ej0c+n09utzvW4wAAgE7ozvs3i5YBAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgBJYVmYO1S3qthWr2XffAYKADsRPACSgsvpUHmY6Gn9IDiX0xGjyQBEAx88CCAphPs023CfegvATgQPgKTx7ehZ/vp7agkEiR0gSXBJC0BSKS7IDt2hOtXlJHaAJEHwAEgqFdvqQ7HTEgh2uJAZgF24pAUgaXx3zU7rc4k7VAO2I3gAJIVwC5TDLWQGYCeCB0BSCARN2AXKrc8DQROLsQBEicMYk7T/K+/O7eUBAEBsdef9m0XLAADAegQPAACwHsEDAACsR/AAAADrETwAAMB6BA8AALAewQMAAKxH8AAAAOsRPAAAwHoEDwAAsB7BAwAArEfwAAAA6xE8AADAegQPAACwHsEDAACsR/AAAADrETwAAMB6BA8AALAewQMAAKxH8AAAAOsRPAAAwHoEDwAAsB7BAwAArEfwAAAA6xE8AADAegQPAACwHsEDAACsF1fBs3LlSg0bNkzp6enKzc1VVVXV9+5/8uRJLVq0SEOHDlVaWpouuugirV27NkrTAgCARJES6wFaVVZWav78+Vq5cqUmTJig1atXa/LkyaqtrdX5558f9jXTpk3Txx9/rDVr1ujiiy/WsWPHdOrUqShPDgAA4p3DGGNiPYQkjR07VqNHj9aqVatC20aMGKGpU6eqrKys3f5btmzRTTfdpAMHDqhfv37d+pl+v18ej0c+n09ut7vbswMAgOjpzvt3XFzSamlpUXV1tQoLC9tsLyws1O7du8O+5uWXX1ZeXp6WLl2qIUOGaPjw4brrrrv05ZdfdvhzTp48Kb/f3+YBAADsFxeXtJqamhQIBJSRkdFme0ZGho4ePRr2NQcOHNDOnTuVnp6uTZs2qampSbfddps++eSTDtfxlJWVacmSJT0+PwAAiG9xcYanlcPhaPPcGNNuW6tgMCiHw6H169drzJgxuu6661ReXq5169Z1eJantLRUPp8v9Dh06FCP/w4AACD+xMUZngEDBsjlcrU7m3Ps2LF2Z31aDRo0SEOGDJHH4wltGzFihIwxOnz4sLKzs9u9Ji0tTWlpaT07PAAAiHtxcYYnNTVVubm58nq9bbZ7vV7l5+eHfc2ECRN05MgRnThxIrStrq5OTqdTmZmZEZ0XAAAklrgIHkkqKSnR008/rbVr12r//v1asGCBGhoaVFRUJOn05aiZM2eG9p8xY4b69++vOXPmqLa2Vjt27NDdd9+tX/7yl+rdu3esfg0AABCH4uKSliRNnz5dx48f1wMPPKDGxkbl5ORo8+bNGjp0qCSpsbFRDQ0Nof3PPvtseb1e3XHHHcrLy1P//v01bdo0Pfjgg7H6FQAAQJyKm8/hiQU+hwcAgMSTsJ/DAwAAEEkEDwAAsB7BAwAArEfwAAAA6xE8AADAegQPAACwHsEDAACsR/AAAADrETwAAMB6BA8AALAewQMAAKxH8AAAAOsRPAAAwHoEDwAAsB7BAwAArEfwAAAA6xE8AADAegQPAACwHsEDAACsR/AAAADrdTt4HnvsMUnSu+++q2Aw2GMDAQAA9LSUzu7Y1NSkAQMGhJ7n5ORIkhYsWKD33ntPZ599tkaOHKmcnBzl5ORoypQpPT8tAABANziMMaYzO2ZnZ2vLli266KKLwn7d7/frnXfe0TvvvKPa2lo9+uijPTlnRPj9fnk8Hvl8Prnd7liPAwAAOqE779+dPsNTUFCg8ePH65VXXtHYsWPbfd3tdis/P1/5+fmdnxgAACAKOr2G54knntCCBQtUUFCgl156KYIjAQAA9KwuLVouLS3V6tWrNWPGDD3++OORmgkAAKBHdfqSVqubb75Z/fr10w033KCNGzdq7Nixys3N1ejRo5WdnR2JGQEAAM5Il4Ln008/1YoVK7R8+XINHDhQmZmZevXVV7Vs2TIFAgH17dtXPp8vUrMCAAB0S6eDZ/78+VqzZo369eune++9V/PmzVNaWpok6auvvlJNTY1qamoiNigAAEB3dfrP0i+++GItXLhQs2bNUq9evSI9V1TwZ+kAACSeiP5Z+rvvviuXy9Xt4QAAAGKl03+lRewAAIBExc1DAQCA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgAAYD2CBwAAWI/gAQAA1iN4AACA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgAAYD2CBwAAWI/gAQAA1iN4AACA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFgvroJn5cqVGjZsmNLT05Wbm6uqqqpOvW7Xrl1KSUnRFVdcEdkBAQBAQoqb4KmsrNT8+fO1aNEi1dTUaOLEiZo8ebIaGhq+93U+n08zZ85UQUFBlCYFAACJxmGMMbEeQpLGjh2r0aNHa9WqVaFtI0aM0NSpU1VWVtbh62666SZlZ2fL5XLppZde0t69ezv9M/1+vzwej3w+n9xu95mMDwAAoqQ7799xcYanpaVF1dXVKiwsbLO9sLBQu3fv7vB1zzzzjN5//30tXry4Uz/n5MmT8vv9bR4AAMB+cRE8TU1NCgQCysjIaLM9IyNDR48eDfua+vp6LVy4UOvXr1dKSkqnfk5ZWZk8Hk/okZWVdcazAwCA+BcXwdPK4XC0eW6MabdNkgKBgGbMmKElS5Zo+PDhnf7+paWl8vl8ocehQ4fOeGYAABD/OndqJMIGDBggl8vV7mzOsWPH2p31kaTm5mbt2bNHNTU1uv322yVJwWBQxhilpKRo69atuvrqq9u9Li0tTWlpaZH5JQAAQNyKizM8qampys3NldfrbbPd6/UqPz+/3f5ut1tvv/229u7dG3oUFRXpkksu0d69ezV27NhojQ4AABJAXJzhkaSSkhLdcsstysvL0/jx4/Xkk0+qoaFBRUVFkk5fjvroo4/07LPPyul0Kicnp83rzz33XKWnp7fbDgAAEDfBM336dB0/flwPPPCAGhsblZOTo82bN2vo0KGSpMbGxh/8TB4AAIBw4uZzeGKBz+EBACDxJOzn8AAAAEQSwQMAAKxH8AAAAOsRPAAAwHoEDwAAsB7BAwAArEfwAAAA6xE8AADAegQPAACwHsEDAACsR/AAAADrETxIGMu8darYVh/2axXb6rXMWxfliQAAiYLgQcJwOR0qDxM9FdvqVe6tk8vpiNFkAIB4lxLrAYDOKi7IliSVf3Mmp7ggOxQ7JZOGh74OAMB3ETxIKN+OnuWvv6eWQJDYAQD8IC5pIeEUF2Qr1eVUSyCoVJeT2AEA/CCCBwmnYlt9KHZaAsEOFzIDANCKS1pIKN9ds9P6XBJnegAAHSJ4kDDCLVAOt5AZAIDvIniQMAJBE3aBcuvzQNDEYiwAQAJwGGOS9l3C7/fL4/HI5/PJ7XbHehwAANAJ3Xn/ZtEyAACwHsEDAACsR/AAAADrETwAAMB6BA8AALAewQMAAKxH8AAAAOsRPAAAwHoEDwAAsB7BA+CMLPPWdXjH+opt9Vr2zX3OACCWCB4AZ8TldKg8TPS03uzV5XTEaDIA+B9uHgrgjIS7Y324O9sDQCwRPADO2LejZ/nr76klECR2AMQVLmkB6BHFBdlKdTnVEggq1eUkdgDEFYIHQI+o2FYfip2WQLDDhcwAEAtc0gJwxr67Zqf1uSTO9ACICwQPgDMSboFyuIXMABBLBA+AMxIImrALlFufB4ImFmMBQBsOY0zS/tfI7/fL4/HI5/PJ7XbHehwAANAJ3Xn/ZtEyAACwHsEDAACsR/AAAADrETwAAMB6BA8AALAewQMAAKxH8AAAAOsRPAAAwHoEDwAAsB7BAwAArEfwAAAA6xE8AADAegQPAACwHsEDAACsR/AAAADrETwAAMB6BA8AALAewQMAAKxH8AAAAOsRPAAAwHoEDwAAsB7BAwAArEfwAAAA68VV8KxcuVLDhg1Tenq6cnNzVVVV1eG+L774oiZNmqSBAwfK7XZr/Pjxeu2116I4LQDg+yzz1qliW33Yr1Vsq9cyb12UJ0Iyi5vgqays1Pz587Vo0SLV1NRo4sSJmjx5shoaGsLuv2PHDk2aNEmbN29WdXW1rrrqKl1//fWqqamJ8uQAgHBcTofKw0RPxbZ6lXvr5HI6YjQZkpHDGGNiPYQkjR07VqNHj9aqVatC20aMGKGpU6eqrKysU99j5MiRmj59uu67775O7e/3++XxeOTz+eR2u7s1NwCgY61xUzJpuIoLsts9B7qjO+/fKRGeqVNaWlpUXV2thQsXttleWFio3bt3d+p7BINBNTc3q1+/fh3uc/LkSZ08eTL03O/3d29gAECntEZNubdOy19/Ty2BILGDmIiLS1pNTU0KBALKyMhosz0jI0NHjx7t1Pd45JFH9Pnnn2vatGkd7lNWViaPxxN6ZGVlndHcAIAfVlyQrVSXUy2BoFJdTmIHMREXwdPK4Wh7PdcY025bOC+88ILuv/9+VVZW6txzz+1wv9LSUvl8vtDj0KFDZzwzAOD7VWyrD8VOSyDY4UJmIJLi4pLWgAED5HK52p3NOXbsWLuzPt9VWVmpuXPnasOGDbrmmmu+d9+0tDSlpaWd8bwAgM7paA2PJM70IKri4gxPamqqcnNz5fV622z3er3Kz8/v8HUvvPCCZs+ereeff15TpkyJ9JgAgC4It0C5uCBbJZOGh/3rLSCS4uIMjySVlJTolltuUV5ensaPH68nn3xSDQ0NKioqknT6ctRHH32kZ599VtLp2Jk5c6Yee+wxjRs3LnR2qHfv3vJ4PDH7PQAApwWCJuwC5dbngWBc/JEwkkTc/Fm6dPqDB5cuXarGxkbl5ORo2bJl+slPfiJJmj17tj744AO98cYbkqT/+7//0/bt29t9j1mzZmndunWd+nn8WToAAImnO+/fcRU80UbwAACQeLrz/h0Xa3gAAAAiieABAADWI3gAAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgAAYD2CBwAAWI/gAQAA1iN4AACA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgAAYD2CBwAAWI/gAQAA1iN4AACA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgAAYD2CBwAAWI/gAQAA1iN4AACA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFiP4Okhy7x1qthWH/ZrFdvqtcxbF+WJAABAK4Knh7icDpWHiZ6KbfUq99bJ5XTEaDIAAJAS6wFsUVyQLUkq/+ZMTnFBdih2SiYND30dAABEH8HTg74dPctff08tgSCxAwBAHOCSVg8rLshWqsuplkBQqS4nsQMAQBwgeHpYxbb6UOy0BIIdLmQGAADRwyWtHvTdNTutzyVxpgcAgBgieHpIuAXK4RYyAwCA6CN4ekggaMIuUG59HgiaWIwFAAAkOYwxSftO7Pf75fF45PP55Ha7Yz0OAADohO68f7NoGQAAWI/gAQAA1iN4AACA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFgvqW8t0foh036/P8aTAACAzmp93+7KzSKSOniam5slSVlZWTGeBAAAdFVzc7M8Hk+n9k3qe2kFg0EdOXJEffv2lcPh6LHv6/f7lZWVpUOHDnGPrijhmEcfxzz6OObRxzGPvs4cc2OMmpubNXjwYDmdnVudk9RneJxOpzIzMyP2/d1uN/8DiTKOefRxzKOPYx59HPPo+6Fj3tkzO61YtAwAAKxH8AAAAOsRPBGQlpamxYsXKy0tLdajJA2OefRxzKOPYx59HPPoi9QxT+pFywAAIDlwhgcAAFiP4AEAANYjeAAAgPUIHgAAYD2Cp5tWrlypYcOGKT09Xbm5uaqqqvre/bdv367c3Fylp6frwgsv1BNPPBGlSe3RlWP+4osvatKkSRo4cKDcbrfGjx+v1157LYrT2qGr/85b7dq1SykpKbriiisiO6CFunrMT548qUWLFmno0KFKS0vTRRddpLVr10ZpWjt09ZivX79eo0aN0llnnaVBgwZpzpw5On78eJSmTXw7duzQ9ddfr8GDB8vhcOill176wdf0yHuoQZf95S9/Mb169TJPPfWUqa2tNXfeeafp06eP+fDDD8Puf+DAAXPWWWeZO++809TW1pqnnnrK9OrVy2zcuDHKkyeurh7zO++80zz00EPmX//6l6mrqzOlpaWmV69e5j//+U+UJ09cXT3mrT777DNz4YUXmsLCQjNq1KjoDGuJ7hzzG264wYwdO9Z4vV5z8OBB889//tPs2rUrilMntq4e86qqKuN0Os1jjz1mDhw4YKqqqszIkSPN1KlTozx54tq8ebNZtGiR+etf/2okmU2bNn3v/j31HkrwdMOYMWNMUVFRm22XXnqpWbhwYdj977nnHnPppZe22XbrrbeacePGRWxG23T1mIdz2WWXmSVLlvT0aNbq7jGfPn26+f3vf28WL15M8HRRV4/5q6++ajwejzl+/Hg0xrNSV4/5ww8/bC688MI22yoqKkxmZmbEZrRZZ4Knp95DuaTVRS0tLaqurlZhYWGb7YWFhdq9e3fY17z55pvt9r/22mu1Z88eff311xGb1RbdOebfFQwG1dzcrH79+kViROt095g/88wzev/997V48eJIj2id7hzzl19+WXl5eVq6dKmGDBmi4cOH66677tKXX34ZjZETXneOeX5+vg4fPqzNmzfLGKOPP/5YGzdu1JQpU6IxclLqqffQpL55aHc0NTUpEAgoIyOjzfaMjAwdPXo07GuOHj0adv9Tp06pqalJgwYNiti8NujOMf+uRx55RJ9//rmmTZsWiRGt051jXl9fr4ULF6qqqkopKfynpau6c8wPHDignTt3Kj09XZs2bVJTU5Nuu+02ffLJJ6zj6YTuHPP8/HytX79e06dP11dffaVTp07phhtu0OOPPx6NkZNST72HcoanmxwOR5vnxph2235o/3Db0bGuHvNWL7zwgu6//35VVlbq3HPPjdR4VursMQ8EApoxY4aWLFmi4cOHR2s8K3Xl33kwGJTD4dD69es1ZswYXXfddSovL9e6des4y9MFXTnmtbW1Ki4u1n333afq6mpt2bJFBw8eVFFRUTRGTVo98R7K/w3rogEDBsjlcrWr/2PHjrUr0FbnnXde2P1TUlLUv3//iM1qi+4c81aVlZWaO3euNmzYoGuuuSaSY1qlq8e8ublZe/bsUU1NjW6//XZJp9+MjTFKSUnR1q1bdfXVV0dl9kTVnX/ngwYN0pAhQ+TxeELbRowYIWOMDh8+rOzs7IjOnOi6c8zLyso0YcIE3X333ZKkyy+/XH369NHEiRP14IMPcsY+AnrqPZQzPF2Umpqq3Nxceb3eNtu9Xq/y8/PDvmb8+PHt9t+6davy8vLUq1eviM1qi+4cc+n0mZ3Zs2fr+eef5/p6F3X1mLvdbr399tvau3dv6FFUVKRLLrlEe/fu1dixY6M1esLqzr/zCRMm6MiRIzpx4kRoW11dnZxOpzIzMyM6rw26c8y/+OILOZ1t3zpdLpek/511QM/qsffQLi1xhjHmf3/GuGbNGlNbW2vmz59v+vTpYz744ANjjDELFy40t9xyS2j/1j+pW7BggamtrTVr1qzhz9K7qKvH/PnnnzcpKSlmxYoVprGxMfT47LPPYvUrJJyuHvPv4q+0uq6rx7y5udlkZmaaX/ziF2bfvn1m+/btJjs728ybNy9Wv0LC6eoxf+aZZ0xKSopZuXKlef/9983OnTtNXl6eGTNmTKx+hYTT3NxsampqTE1NjZFkysvLTU1NTeijACL1HkrwdNOKFSvM0KFDTWpqqhk9erTZvn176GuzZs0yV155ZZv933jjDfOjH/3IpKammgsuuMCsWrUqyhMnvq4c8yuvvNJIaveYNWtW9AdPYF39d/5tBE/3dPWY79+/31xzzTWmd+/eJjMz05SUlJgvvvgiylMntq4e84qKCnPZZZeZ3r17m0GDBpmbb77ZHD58OMpTJ65//OMf3/vf50i9hzqM4RwcAACwG2t4AACA9QgeAABgPYIHAABYj+ABAADWI3gAAID1CB4AAGA9ggcAAFiP4AEAANYjeAAAgPUIHgAJLxgM6tJLL9Vvf/vbNttfe+01paamasOGDTGaDEC8IHgAJDyn06nS0lKtWrVKn376qSTpv//9r2688Ub94Q9/0I033hjjCQHEGvfSAmCFU6dOafjw4Zo1a5bmzp2rcePG6ac//alWrFgR69EAxAGCB4A1Vq9erd/97ncaMmSILrjgAm3atEkulyvWYwGIAwQPAGucOHFCAwcO1MUXX6y33npLffr0afP1v//973r77be1YMGCGE0IIFYIHgDWmD17tiorK3XOOefo4MGDSk9Pj/VIAOIEi5YBWOHee+/V3/72N7311ls6deqU1qxZ026fyZMna//+/TGYDkCsETwAEt7TTz+tRx55RK+88opGjRql+fPna+nSpfr666/b7FdfX6/s7OwYTQkglggeAAnt1Vdf1a9//Ws999xzGjdunCTpjjvukN/v15///OfQfj6fT2effbZSUlJiNSqAGCJ4ACSs6upq3XjjjVq6dKl+9rOfhba73W7dcccd+uMf/6hAICBJeueddzRy5MhYjQogxli0DCAprF69Wp988olKS0tjPQqAGOAMD4CksG/fPuXk5MR6DAAxwhkeAABgPc7wAAAA6xE8AADAegQPAACwHsEDAACsR/AAAADrETwAAMB6BA8AALAewQMAAKxH8AAAAOsRPAAAwHoEDwAAsB7BAwAArPf/PM1lxqCxnMcAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "xi = np.array([0.97,0,0.5,0.85,0.7,0.19,0.41,0.78])\n", "yi = np.array([0.97,0.06,0.7,0.74,0.2,0.34,0.29,0.94])\n", "plt.xlabel('$X_i$')\n", "plt.ylabel('$Y_i$')\n", "plt.plot(xi,yi,'x');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Examples of the source of this data can be:\n", " - Stress of a pre-stressed concrete specimen based on the deformation when applied an axial force\n", " - Permeability-coefficient for compacted sand based on the grain size\n", " - Job performance based on the number of working hours\n", " - etc.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Model\n", "The problem is modelled as a minimization of the squared error, where error is the difference between $a X_i^b+C$ and $Y_i$\n", "$$\\mathop {\\min }\\limits_{a,b,c} \\sum\\limits_{i = 1}^8 {{{\\left( {a{X_i}^b + c - {Y_i}} \\right)}^2}} $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The plot below shows the data points and the total error squared." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [] }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "029c3407d1ac46e7a9765bca5443248d", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(FloatSlider(value=2.0, description='a', max=5.0, step=0.2), FloatSlider(value=2.0, descr…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from ipywidgets import widgets, interact\n", "x_range = np.linspace(0.0001,1,100)\n", "def func(a,b,c):\n", " fig, ax = plt.subplots(1, 1)\n", " ax.plot(xi,yi,'x')\n", " ax.plot(x_range,a * x_range**b + c);\n", " ax.set_xlim([-0,1])\n", " ax.set_ylim([-0,1])\n", " y_est = a * xi ** b + c\n", " error = np.sum((yi-y_est)**2)\n", " title = 'Total error squared $=$'+str(round(error,4))\n", " ax.set_title(title)\n", " for i in range(8):\n", " ax.annotate(\"\",xy=[xi[i],yi[i]],xytext=[xi[i],max(min(a * xi[i] ** b + c,1),0)],arrowprops=dict(arrowstyle= '<->'))\n", " return\n", " plt.show()\n", "\n", "interact(func,\n", " a = widgets.FloatSlider(min=0, max=5, value=2, step=0.2, description=\"a\"),\n", " b = widgets.FloatSlider(min=0, max=5, value=2, step=0.2, description=\"b\"),\n", " c = widgets.FloatSlider(min=0, max=1, value=0.5, step=0.05, description=\"c\"));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Method\n", "\n", "Again, this model is described using `scipy.optimize.minimize` according to the [standard structure in this course](./Python_start_3.ipynb)\n", "\n", "#### Importing libraries\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "import scipy as sp\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Defining the variables\n", "In this problem, there are 3 design variables: $a$, $b$ and $c$. In `scipy.optimize`, these variables are part of one array. You don't have to specify this variable beforehand, but the optimization algorithm will define the number of variables on your other input. Furthermore, input data is given which is required to specify explicitly. Additionally, an initial guess is required for the otimization algorithm." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "xi = np.array([0.97, 0 , 0.5, 0.85, 0.7, 0.19, 0.41, 0.78])\n", "yi = np.array([0.97, 0.06, 0.7, 0.74, 0.2, 0.34, 0.29, 0.94])\n", "abc0 = np.array([4., 2., 0.5]) #initial guess" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Objective function\n", "The objective function was alerady defined. Because `a`, `b`, `c` and `xi` are all numpy arrays, we can simply use `*`, `**` and `+` to deal with the multiplication, exponent and summation of those arrays. `np.sum` is used to sum up all components of the array which is the result of `(yi-y_est)**2`" ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [], "source": [ "def squarederror(abc):\n", " a = abc[0]\n", " b = abc[1]\n", " c = abc[2]\n", " y_est = a * xi ** b + c\n", " error = np.sum((yi-y_est)**2)\n", " return error" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Solve the problem\n", "Now we can solve the actual problem" ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " message: Optimization terminated successfully.\n", " success: True\n", " status: 0\n", " fun: 0.3342203980898386\n", " x: [ 8.028e-01 1.231e+00 1.233e-01]\n", " nit: 54\n", " jac: [-2.332e-06 2.570e-06 -3.166e-06]\n", " hess_inv: [[ 6.233e-01 -3.674e-01 -3.660e-01]\n", " [-3.674e-01 1.485e+01 2.251e+00]\n", " [-3.660e-01 2.251e+00 5.564e-01]]\n", " nfev: 276\n", " njev: 69\n" ] } ], "source": [ "result = sp.optimize.minimize(squarederror,abc0)\n", "print(result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Postprocessing\n", "As this problem involves three design variables, a plot of all possible design variables with the objective function is not possible because it would require four dimensions. However, we can plot the data which the obtained solution to see how it looks like" ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAGzCAYAAADNKAZOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8r0lEQVR4nO3daXgV5cHG8X8WEtYEEUR2UUFR6gbKJm4IFjfQqvjaKlq1UhcEqlbqirVSbQVFFneta6lUXHGJqIiiVihYBRVEFBAQAU3CFuCceT+MiSJBk5Bkzjn5/64rH57JOeHOFDl3n3lmnrQgCAIkSZJSTHrUASRJkqqCJUeSJKUkS44kSUpJlhxJkpSSLDmSJCklWXIkSVJKsuRIkqSUZMmRJEkpKTPqAFGKx+MsW7aMBg0akJaWFnUcSZJUBkEQUFhYSPPmzUlP3/58TY0uOcuWLaNVq1ZRx5AkSRWwZMkSWrZsud3vJ0zJeeONN/jb3/7GrFmzWL58OZMnT6Z///4/+Z5p06YxbNgw5s6dS/PmzbniiisYNGhQmf/MBg0aAOFJysnJ2ZH4kiSpmhQUFNCqVauSz/HtSZiSs27dOvbff3/OOeccfvWrX/3s6xctWsSxxx7L+eefzyOPPMJbb73FhRdeSJMmTcr0fqDkElVOTo4lR5KkJPNzS00SpuT07duXvn37lvn1d955J61bt+a2224DoEOHDsycOZO///3vZS45kiQpdSXt3VVvv/02ffr02erYMcccw8yZM9m8eXOp7ykqKqKgoGCrL0mSlJqStuSsWLGCpk2bbnWsadOmbNmyhVWrVpX6npEjR5Kbm1vy5aJjSZJSV9KWHNj2WlwQBKUeLzZ8+HDy8/NLvpYsWVLlGSVJUjQSZk1Oee26666sWLFiq2MrV64kMzOTnXfeudT3ZGdnk52dXR3xJElSxJJ2Jqdbt27k5eVtdezll1+mc+fO1KpVK6JUkiQpUSRMyVm7di1z5sxhzpw5QHiL+Jw5c1i8eDEQXmo666yzSl4/aNAgvvjiC4YNG8ZHH33E/fffz3333cdll10WRXxJkpRgEuZy1cyZMznyyCNLxsOGDQNg4MCBPPjggyxfvryk8AC0bduWKVOmMHToUMaNG0fz5s0ZM2aMt49LkiQA0oLi1bo1UEFBAbm5ueTn5/swQEmSkkRZP78T5nKVJElSZbLkSJKkbYzOm8+YqQtK/d6YqQsYnTe/mhOVnyVHkiRtIyM9jVGlFJ0xUxcwKm8+Gek/vW9UIkiYhceSJClxDO7VDoBR383YDO7VrqTgDOvdvuT7icySI0mSSvXDojP21U/ZFIsnTcEBL1dJkqSfMLhXO7Iy0tkUi5OVkZ40BQcsOZIk6SeMmbqgpOBsisW3uxg5EXm5SpIklerHa3CKx0BSzOhYciRJ0jZKW2Rc2mLkRGbJkSRJ24jFg1IXGRePY/HE3zDBbR3c1kGSpKTitg6SJKlGs+RIkpJaKmw/oKphyZEkJbVU2H5AVcOFx5KkpJYK2w+oalhyJElJL9m3H1DV8HKVJCklJPP2A6oalhxJUkpI5u0HVDW8XCVJSnrJvv2AqoYlR5KU1FJh+4GUVLAcPngCul8CadHc4WbJkSQltVTYfiDlLHkPJv4G1q6A7AbQ+ZxIYlhyJElJbWjv9tv9njM4EfjvQ/D8HyC2CZp0gN0PjyyKJUeSJO242GZ4cTi8d0843vt4OOnOcCYnIpYcSZK0Y9atgn8NhC/eDMdHXgU9L4P0aG/ituRIkqSKWzYnXH+TvwSyGsDJd8Pex0adCrDkSJKkivpgEjx9MWzZAI32gNMfg132jjpVCUuOJEkqn3gMXrkeZowJx3v2hl/dC3UaRplqG5YcSZJUduvXwL/PhYWvhuNDh8FRV0N6RrS5SmHJkSRJZfPVXPjnGfDN51CrLvQbBx1PjjrVdllyJEnSz5v7FDx1IWxeBw3bhOtvdu0YdaqfZMmRJEnbF4/Ba3+B6beG492PgFMegLqNIo1VFpYcSZJUug3fwr/Pg0/zwnG3i+HoEZCRHPUhOVJKkqTqtfLjcP3NmoWQWRtOvAP2Oy3qVOViyZEkSVv76FmYPAg2rYXcVjDgEWh+QNSpys2SI0mSQvE4vD4S3rglHO/WE059EOo1jjRWRVlyJElSuP7myd/BgpfCcdcLofefk2b9TWmSN7kkSaocP15/c/xtcMD/RZ1qh1lyJEmqyT56DiZfEK6/yWkJpz8CzQ+MOlWlsORIklQTpdj6m9JYciRJqml+vP6myyDocyNk1Io0VmWz5EiSVJP8eP3NCbfD/qdHnapKWHIkSaoptnn+zcMps/6mNJYcSZJSXTwGr90E0/8ejlNw/U1pLDmSJKWyDd/Av8//fv+pFHj+TVml/m8oSVJN9dW8cP3NN4sgsw6cOCbp9p/aEZYcSZJS0dzJ8NRFsHkd5LYOn3/TbP+oU1UrS44kSakktgWmjoAZY8Jx28PhlAeg3s7R5oqAJUeSpFSxbjVMOgcWTQvH3QdDr+tqxPqb0tTM31qSpFSz/H34528gfzHUqgv9xkHHk6NOFSlLjiRJye79f8Kzl8KWjbBTWzj9UWi6b9SpImfJkSQpWcU2w0tXwX/uCsd79oZf3QN1doo2V4Kw5EiSlIwKv4InzobFM8LxYVfAEcMhPT3SWInEkiNJUrJZOhMm/gYKl0N2Dpx0F+x9bNSpEo4lR5KkZBEEMOsBmHIFxDdD473C9TeN20WdLCFZciRJSgabN8KUy2D2w+G4wwnQfwJkN4g2VwKz5EiSlOi+XQL/OhOWzYa0dOh1LfQYAmlpUSdLaJYcSZIS2WfTwgf8rV8d3jV1yv2wx1FRp0oKlhxJkhJREMCMO+CV6yCIw677wYBHYKc2USdLGpYcSZISTdFaePoimPdUON7/DDh+FNSqE2msZGPJkSQpkaz6FCb+Gr7+GNIz4Zd/hYPPc/1NBVhyJElKFB8/D5MHQVEB1N8VTnsIWneJOlXSsuRIkhS1eAxeuwmm/z0ct+4Gp/4DGjSNNleSs+RIkhSl9Wvg3+fBwqnhuMvvoc+fIaNWtLlSQEJtcDF+/Hjatm1L7dq16dSpE9OnT//J1z/66KPsv//+1K1bl2bNmnHOOeewevXqakorSdIOWv4+3H14WHAy68DJ90Lfv1pwKknClJyJEycyZMgQrrrqKmbPnk3Pnj3p27cvixcvLvX1b775JmeddRbnnnsuc+fO5YknnuC9997jvPPOq+bkkiRVwJzH4L4+8O1i2KktnPcK7Hdq1KlSSsKUnFGjRnHuuedy3nnn0aFDB2677TZatWrFhAkTSn39O++8w2677cbgwYNp27Ythx56KBdccAEzZ86s5uSSJJXDliJ4big89XvYshHaHQO/ew127Rh1spSTECVn06ZNzJo1iz59+mx1vE+fPsyYMaPU93Tv3p2lS5cyZcoUgiDgq6++YtKkSRx33HHb/XOKioooKCjY6kuSpGqT/yU8cCzMvB9IgyP+BP/3z/BJxqp0CVFyVq1aRSwWo2nTrVeRN23alBUrVpT6nu7du/Poo48yYMAAsrKy2HXXXWnYsCF33HHHdv+ckSNHkpubW/LVqlWrSv09JEnarkVvhOtvvpwJtXPh10/AEX+E9IT4KE5JCXVm0370oKMgCLY5VmzevHkMHjyYa6+9llmzZvHiiy+yaNEiBg0atN2fP3z4cPLz80u+lixZUqn5JUnaRhDAW7fDQ/1g3dfQ9Bfwu2nQrnfUyVJeQtxC3rhxYzIyMraZtVm5cuU2szvFRo4cSY8ePbj88ssB2G+//ahXrx49e/bkxhtvpFmzZtu8Jzs7m+zs7Mr/BSRJKs3GAnj6Qvjo2XC8///BcaMgq260uWqIhJjJycrKolOnTuTl5W11PC8vj+7du5f6nvXr15P+oym+jIwMIJwBkiQpUl9/Avf2CgtOeq2w3PSfYMGpRgkxkwMwbNgwzjzzTDp37ky3bt24++67Wbx4ccnlp+HDh/Pll1/y0EMPAXDCCSdw/vnnM2HCBI455hiWL1/OkCFDOOSQQ2jevHmUv4okqab78El4+mLYvA4aNIcBD0PLzlGnqnESpuQMGDCA1atXc8MNN7B8+XI6duzIlClTaNMm3FJ++fLlWz0z5+yzz6awsJCxY8fyhz/8gYYNG3LUUUdx8803R/UrSJJquthmeOV6eHtsON6tJ5zyANRvEmmsmiotqMHXdgoKCsjNzSU/P5+cnJyo40iSklnhVzDpHPjirXDcYwgcdQ1kJMx8Qsoo6+e3Z16SpB21+B3410BYuwKyGsBJE6DDCVGnqvEsOZIkVVQQwLt3wstXQ3wLNNkbBjwKjfeMOpmw5EiSVDFFa+GZS2Duk+G44ylwwu2QXT/aXCphyZEkqby+ng//OhO+/hjSM6HPX6DLBbCdB9gqGpYcSZLKY+5T8PRFsGkt1N8VTvsHtO4adSqVwpIjSVJZ/Pj28DaHwin3Q4PSn8yv6FlyJEn6OYUr4IlzYPGMcNzjUjjqWm8PT3D+ryNJ0k/5/C144mxYt9Lbw5OMJUeSpNIEAcy4I7xEFcRgl33gtIe9PTyJWHIkSfqxjfnh4uLi3cP3GwDHj4asetHmUrlYciRJ+qGv5sLEM2HNQsjIgl/+FTr/1tvDk5AlR5KkYu//E54dAls2QE5LOO0haNnpJ98yOm8+GelpDO7VbpvvjZm6gFg8YGjv9lUUWD8lPeoAkiRFbksRPDcUJl8QFpw9joIL3vjZggOQkZ7GqLz5jJm6YKvjY6YuYNR3BUjRcCZHklSzffMFPDEQls0G0uDwP8LhV0B6RpneXjyDMypvfsm4uOAM692+1BkeVQ9LjiSp5lqQB0+eDxu+gTo7wcn3Qrujy/1jflh0xr76KZticQtOAvBylSSp5onH4NW/wKOnhAWn+YHh5akKFJxig3u1IysjnU2xOFkZ6RacBGDJkSTVLOtWwSO/gjduCcedz4XfvgQNW+/Qjx0zdUFJwdkUi2+zRkfVz8tVkqSaY8l/wqcXF3wJterCCbfDfqft8I/98Rqc4jHgjE6ELDmSpNQXBPDunfDy1RDfAju3C28Pb7rPDv/o0hYZl7YYWdXPkiNJSm0bC+CZi2He0+F4n/7QbyxkN6iUHx+LB6UuMi4ex+JBpfw5Kr+0IAhq7NkvKCggNzeX/Px8cnJyoo4jSapsX82Ff50Fqz+F9FpwzF/gkN/59OIkV9bPb2dyJEmpac7j4QP+ip9efOqD0OrgqFOpGllyJEmpZfMGeOEK+O9D4XiPo8Ln39TbOdpcqnaWHElS6ljzWXh5asUHQBocMRwOu6zMTy9WarHkSJJSw0fPwVMXQlE+1N0ZfnVvOIujGsuSI0lKbrHN8Mr18PbYcNyqC5zyAOS2iDSWomfJkZQyRn+343NpzyQZM3UBsXjA0N7tI0imKpP/JUz6LSx5Jxx3vQh6j4CMWtHmUkJwWwdJKSMjPY1RefO3eZx+8cPaMtK9bTilLHwV7uoZFpzsHDjtYfjlTRYclXAmR1LKKO0ps6U9jVZJLh6DabfAtJuBAHb9BZz6D9h5j6iTKcFYciSllB8WnbGvfsqmWNyCk0rWfg1PngefvR6ODxoIfW+GWnUijaXE5OUqSSlncK92JTtBZ2WkW3BSxRczwstTn70ebq550l1w4hgLjrbLkiMp5YyZuqCk4GyKxbdZo6MkE4/Dm7fBg8dD4XJovBec/yrsf3rUyZTgvFwlKaX8eA1O8RjcCToprV8DT/0e5r8Yjn9xGhw/GrLrR5tLScGSIylllLbIuLTFyEoSS2fCE2dD/hLIyA7X3nQ62801VWaWHEkpIxYPSl1kXDyOxYMoYqm8ggDevQtevhrim2GntnDaP6DZ/lEnU5JJC4Kgxv5XX9at2iVJ1WRjPjx9MXz0TDjucCL0Gwu1c6PNpYRS1s9vZ3IkSYlh2Zzw8tQ3iyC9FvT5M3QZ5OUpVZglR5IUrSCAmffDi8MhVgS5reHUB6Flp6iTKclZciRJ0SkqhGeHwIeTwnH7vtB/PNRtFGkspQZLjiQpGis+hCcGwupPIS0j3Fiz28VenlKlseRIkqpXEMB/H4IXroAtGyGnBZzyALTuEnUypRhLjiSp+hStheeHwf8mhuN2fcLtGbw8pSpgyZEkVY+v5oZ3T62aH16e6nUNdL8U0t1hSFXDkiNJqlpBALMfhimXh5enGjSDU+6HNt2jTqYUZ8mRJFWdH1+e2vPo8PJUvcbR5lKNYMmRJFWNH1+eOupq6DHEy1OqNpYcSVLl8vKUEoQlR5JUeYrWwnND4YN/heM9esHJd3t5SpGw5EiSKseKD8PLU6sXeHlKCcGSI0naMUEAsx6EF6/87vJU8+8uT3WLOplqOEuOJKniNhbAc0Pgw3+H4z17f3f31M6RxpLAkiNJqqjl74eXp9Z8Fl6eOvo66HaJl6eUMCw5kqTyCQJ471546U8Q2wS5rcLLU60OiTqZtBVLjiSp7DZ8C89cAh89E473Ohb6jXPvKSUkS44kqWyWzoJJZ8O3iyG9FvS+Abr+HtLSok4mlcqSI0n6aUEAb4+DV66D+BZo2AZOfQBadIo6mfSTLDmSpO1btxqe+j0seCkc79MPTrwDaudGm0sqA0uOJKl0X8yASedC4TLIyIZj/gIHn+flKSUNS44kaWvxGLw5Cl67CYI47LwnnPIANNsv6mRSuVhyJEnfK1wBT/4OFk0Lx/udDsfdCtn1o80lVYAlR5IU+nQqTL4A1n0NterCsX+HA38ddSqpwiw5klTTxTbDa3+BN0eH4132De+earJXtLmkHWTJkaSa7NvF4eLipf8Jx53PDRcY16oTbS6pEiTUBiPjx4+nbdu21K5dm06dOjF9+vSffH1RURFXXXUVbdq0ITs7mz322IP777+/mtJKUpKb9wzceWhYcLJz4dR/wPGjLDhKGQkzkzNx4kSGDBnC+PHj6dGjB3fddRd9+/Zl3rx5tG7dutT3nHbaaXz11Vfcd9997LnnnqxcuZItW7ZUc3JJSjKbN8LLV4X7TwG06Ayn3Ac77RZpLKmypQVBEEQdAqBLly4cdNBBTJgwoeRYhw4d6N+/PyNHjtzm9S+++CKnn346n332GY0aVWzPlIKCAnJzc8nPzycnJ6fC2SUpaXw9HyadA199GI57DIGjroaMWpHGksqjrJ/fCXG5atOmTcyaNYs+ffpsdbxPnz7MmDGj1Pc888wzdO7cmVtuuYUWLVrQvn17LrvsMjZs2LDdP6eoqIiCgoKtviSpRggCmP0I3H14WHDqNobf/Bt6j7DgKGUlxOWqVatWEYvFaNq06VbHmzZtyooVK0p9z2effcabb75J7dq1mTx5MqtWreLCCy9kzZo1212XM3LkSEaMGFHp+SUpoW0sgOeGwoeTwnHbw+Hku6HBrtHmkqpYQszkFEv70aPCgyDY5lixeDxOWloajz76KIcccgjHHnsso0aN4sEHH9zubM7w4cPJz88v+VqyZEml/w6SlFC+nAV3HRYWnLQM6HUtnPmUBUc1QkLM5DRu3JiMjIxtZm1Wrly5zexOsWbNmtGiRQtyc7/fJK5Dhw4EQcDSpUtp167dNu/Jzs4mOzu7csNLUiKKx+GdcfDKCIhvhtzW8Kt7oXWXqJNJ1SYhZnKysrLo1KkTeXl5Wx3Py8uje/fupb6nR48eLFu2jLVr15Ycmz9/Punp6bRs2bJK80pSQlu7Eh47FV6+Oiw4HU6EQW9YcFTjJETJARg2bBj33nsv999/Px999BFDhw5l8eLFDBo0CAgvNZ111lklrz/jjDPYeeedOeecc5g3bx5vvPEGl19+Ob/97W+pU8dnPEiqoRa+ChN6wKevQGZtOH40nPYQ1Nkp6mRStUuIy1UAAwYMYPXq1dxwww0sX76cjh07MmXKFNq0aQPA8uXLWbx4ccnr69evT15eHpdccgmdO3dm55135rTTTuPGG2+M6leQpOjENsOrN8Jbt4XjJh3CrRl26RBpLClKCfOcnCj4nBxJKWHNIvj3ueEiY4DOv4VjbvLJxUpZZf38TpiZHElSBXwwCZ4dApsKoXYunHgH7NMv6lRSQrDkSFIyKloLL1wBcx4Nx627wcn3QMNW0eaSEoglR5KSzbI5MOm3sGYhpKXDYVfAYZdDhv+kSz/kfxGSlCzicXhnPLxyfXhreE6L8MnFux0adTIpIVlyJCkZrF0JT/0+vDUcYO/jw/U3dSu2QbFUE1hyJCnRffoKTB4E674On33zy5HQ6RzYzrY3kkKWHElKVFuKYOoN8PbYcLzLPnDK/T77RiojS44kJaJVC8LFxSv+F44PPh/6/Nln30jlYMmRpEQSBDD7YXjhj7B5PdRpBP3Gwd7HRp1MSjqWHElKFBu+gWcvhXlPh+O2h8FJd0FO82hzSUnKkiNJieDzt+DJ30HBUkjPhKOuge6DIT1h9lGWko4lR5KiFNsC0/4K02+FIA6Ndodf3QstOkWdTEp6lhxJisqaRfDk+bD0vXB8wK+h782Q3SDaXFKKsORIUhTenwjP/yHcWDM7F44fBb84JepUUkqx5EhSddqYD89fBh/8Kxy37hZuzdCwdbS5pBRkyZGk6rL4XXjyPPh2MaRlwBFXwqHD3FhTqiL+lyVJVS22Bd74G7xxS7i4uGEbOPkeaN0l6mRSSrPkSFJV+uZz+Pf5sPQ/4Xi/0+HYv0HtnEhjSTWBJUeSqkIQwP/+9YPFxTlw/GgXF0vVyJIjSZVtw7dhuflwUjhu1TVcXLxTm0hjSTVNhR+lefvttwPwySefEI/HKy2QJCW1L2bAnYeGBSctA468Gs5+3oIjRaDMMzmrVq2icePGJeOOHTsCMHToUD799FPq16/PvvvuS8eOHenYsSPHHXdc5aeVpEQV2wyvj4Q3R4eLi3dqGz65uGXnqJNJNVZaEARBWV7Yrl07XnzxRfbYY49Sv19QUMCHH37Ihx9+yLx587jtttsqM2eVKCgoIDc3l/z8fHJyXAQoqYJWfRreGr5sdjj2ycVSlSrr53eZZ3J69epFt27dePbZZ+nSZdvbHnNycujevTvdu3evWGJJSjZBAP/9B7w4HDavh9oN4YTbYN+Tok4miXKUnDvvvJM2bdrQq1cvHnnkEfr371+FsSQpwa1bDc8Oho+fC8dtD4P+d0Jui2hzSSpRroXHw4cP56677uKMM87gjjvuqKpMkpTYFrwCE7qFBSe9FvT+M5z5tAVHSjDlvoX817/+NY0aNeLEE09k0qRJdOnShU6dOnHQQQfRrl27qsgoSYlh8wbIuw7+c1c4brJ3+OTiZvtFm0tSqcpVcr755hvGjRvH2LFjadKkCS1btuSFF15g9OjRxGIxGjRoQH5+flVllaToLP8fPHk+fP1xOD7kAug9AmrViTaXpO0qc8kZMmQI9913H40aNeKaa67hvPPOIzs7G4CNGzcye/ZsZs+eXWVBJSkS8RjMuANevRHim6F+U+g3HtodHXUyST+jzLeQ77nnnlx55ZUMHDiQWrVqVXWuauEt5JJ+0reLYfLv4Ys3w/Fex8GJY6Be459+n6QqVem3kH/yySdkZGRUSjhJSmhBAB88EW7NUFQAtepB37/CgWdCWlrU6SSVUZlLjgVHUo2w4Rt4bhjMfTIctzw43Heq0e7R5pJUbm7QKUnFPns9vDxVuCzcd+rwP0LPP0CG/1RKycj/ciVp8waYegO8Mz4cN9ojvDW8Zadoc0naIZYcSTXb8vfhyQvg64/Ccedzoc+fIatetLkk7TBLjqSaKR6Dt26H124Kbw2vtwv0Gwvtj4k6maRKYsmRVPOsWQSTB8GSd8Lx3sfDCbd7a7iUYiw5kmqOIIDZj8CLV8KmtZDVAPreDAec4a3hUgqy5EiqGdZ+He4a/smUcNy6G5x0J+y0W6SxJFUdS46k1PfxFHjmEli/Ktw1/KirofslkO7zv6RUZsmRlLo2FsBLf4LZD4fjXfaFk++CXX8RbS5J1cKSIyk1ff4WPDUo3H+KtHDm5qirITM76mSSqoklR1Jq2bwRXrsRZowFAmjYGvrfCbv1iDqZpGpmyZGUOpb/DyZfACvnheMDz4RjboLa29+lWFLqsuRISn6xLfDWaHj9rxDfAvWawAljYO9jo04mKUKWHEnJbdWn4ezNlzPD8d7Hw/G3Qf0mkcaSFD1LjqTkFI/De/dC3rWwZQNk58Cxf4P9BvhgP0mAJUdSMspfCk9fBJ+9Ho53PwL6jYPcllGmkpRgLDmSkkcQwPuPwwt/hKICyKwDvW+Ag8+D9PSo00lKMJYcSclh7Up4dgh88nw4bnlweGt44z0jjSUpcVlyJCW+uU/Bc0Nhw5pwW4Yjh0P3SyHDf8IkbZ//QkhKXOvXwJTL4cNJ4bhpx3BTTbdlkFQGlhxJiWn+y+GmmmtXQFoG9BwGh10BmVlRJ5OUJCw5khLLxvzvNtV8JBw3bh+uvWnZKdpckpKOJUdS4lj4Gjx9MRQsBdKg64XQ6xqoVSfqZJKSkCVHUvSK1sIr14UP9wPYaTfoPwHadI80lqTkZsmRFK3P34KnL4RvPg/HB58HR4+A7PqRxpKU/Cw5kqKxaT1MvQHenRCOc1pCv7Gwx5HR5pKUMiw5kqrf4nfgqQthzcJwfNBZ0OcvUDsn2lySUoolR1L12bwBXr0R3h4HBNCgOZx4B7Q7OupkklKQJUdS9VjyHjz1e1i9IBzvfwb8ciTUaRhpLEmpy5IjqWpt3giv/QXeHgtBHOrvCieOgfbHRJ1MUoqz5EiqOktnhrM3q+aH4/1Oh75/hTo7RZtLUo2QHnWAHxo/fjxt27aldu3adOrUienTp5fpfW+99RaZmZkccMABVRtQUtls3ggvXwP39Q4LTv2mcPrjcPJdFhxJ1SZhSs7EiRMZMmQIV111FbNnz6Znz5707duXxYsX/+T78vPzOeuss+jVq1c1JZX0k5a8B3f1hBljwstT+w2AC9+BvY+NOpmkGiYtCIIg6hAAXbp04aCDDmLChAklxzp06ED//v0ZOXLkdt93+umn065dOzIyMnjqqaeYM2dOmf/MgoICcnNzyc/PJyfHW1elHbJ5w3drb8Z9t/amKRx/m+VGUqUr6+d3QszkbNq0iVmzZtGnT5+tjvfp04cZM2Zs930PPPAACxcu5LrrrivTn1NUVERBQcFWX5IqweJ34M5DYcYd383enO7sjaTIJcTC41WrVhGLxWjatOlWx5s2bcqKFStKfc+CBQu48sormT59OpmZZfs1Ro4cyYgRI3Y4r6TvbFoPr/4Z3plA+NybZuHszV6/jDqZJCXGTE6xtLS0rcZBEGxzDCAWi3HGGWcwYsQI2rdvX+afP3z4cPLz80u+lixZssOZpRrr8zdhQnd4ZzwQwAG/CWdvLDiSEkRCzOQ0btyYjIyMbWZtVq5cuc3sDkBhYSEzZ85k9uzZXHzxxQDE43GCICAzM5OXX36Zo446apv3ZWdnk52dXTW/hFRTFBXCK9d/v2N4g+bhc2/a9Y40liT9WEKUnKysLDp16kReXh4nnXRSyfG8vDz69eu3zetzcnL44IMPtjo2fvx4Xn31VSZNmkTbtm2rPLNUIy18FZ65FPK/u+vxoIHQ589QOzfaXJJUioQoOQDDhg3jzDPPpHPnznTr1o27776bxYsXM2jQICC81PTll1/y0EMPkZ6eTseOHbd6/y677ELt2rW3OS6pEmz4Fl6+GmY/HI4btoYTxrhjuKSEljAlZ8CAAaxevZobbriB5cuX07FjR6ZMmUKbNm0AWL58+c8+M0dSFfjkBXhuKBQuD8cHnw9HXw/Z9SONJUk/J2GekxMFn5Mj/YR1q+HFP8IHT4TjRntAv7HQpnu0uSTVeGX9/E6YmRxJCSIIYO5kmHI5rF8FaenQ7WI48k9Qq07U6SSpzCw5kr5XsBye/wN88nw4btIB+o+DFp2izSVJFWDJkRTO3sx+GF66GoryIT0Tel4GPf8AmVlRp5OkCrHkSDXdmkXw7KWwaFo4bn4g9BsHTfeNNpck7SBLjlRTxWPw7p3w6o2weT1k1oYjr4KuF0KG/zRISn7+SybVRF/Ng2cuhi9nhePdesIJt8POe0SbS5IqkSVHqkm2FMH0W2H6KIhvhuyc8InFBw2EUvaJk6RkZsmRaorF78Izl8CqT8LxXsfBcX+HnObR5pKkKmLJkVJdUSG8MuK7DTUDqLcLHPs32KefszeSUpolR0pln7wIzw+Dgi/D8YG/gd5/hrqNos0lSdXAkiOlorUr4YUrwicXA+y0Gxx/mxtqSqpRLDlKaKPz5pORnsbgXu22+d6YqQuIxQOG9m4fQbIEFQQw51F46SrY+C2kZUC3i+CI4ZBVN+p0klStLDlKaBnpaYzKmw+wVdEZM3UBo/LmM8yC873VC8OH+n0+PRzvuh+ceAc0PyDSWJIUFUuOElpxsflh0flhwSlthqfG2bIJZoyBabdArAgy68CRw6HrRT7UT1KN5r+ASng/LDpjX/2UTbG4BafYkvfg2cGwcl443uMoOH50uAZHkmq49KgDSGUxuFc7sjLS2RSLk5WRbsHZmB/uFn5f77Dg1N0ZTr4HfvOkBUeSvmPJUVIYM3VBScHZFIszZuqCqCNF56NnYVyX7597s/8ZcNF7sN9pPvdGkn7Ay1VKeD9eg1M8BmrWjE7+UphyBXzyfDhutHt4W/juh0caS5ISlSVHCa20RcalLUZOafEY/OfucLfwTWshPRN6DIHDLoNadaJOJ0kJy5KjhBaLB6UuMi4ex+JBFLGqz7LZ8OwQWD4nHLfqEs7eNN0nwlCSlBzSgiBI8U+J7SsoKCA3N5f8/HxycnKijiN9r6gQXrsJ3r0TgjjUzoWjR4S7hae7lE5SzVbWz29ncqREEgTw8XPwwh+/32+q4ylwzE3QoGm02SQpyVhypETx7eJwYfH8F8LxTrvBcbfCnkdHGkuSkpUlR4pabDO8MwFeHwmb10N6LehxqQuLJWkHWXKkKC1+B54b+v0Ti1t3D59YvMve0eaSpBRgyZGisH4NvHId/PehcFynEfT5c/hgPxcWS1KlsORI1Skeh/cfg7xrYf3q8NiBZ0LvG6Buo2izSVKKseRI1eWrefD8MFj8djhu0iG8NNWmW7S5JClFWXKkqla0Fqb9Fd4eD0EMatWFI66ErhdCRq2o00lSyrLkSFUlCOCjZ+DF4d8/82bv4+GXf4WGraLNJkk1gCVHqgqrF8ILV8Cnr4Tjhq2h799gr19Gm2sHjM6bT0Z6Wql7hY2ZuoBYPGBo7/YRJJOk0nkbh1SZNm+A10bC+G5hwcnIgsMuhwvfTeqCA5CRnsaovPmMmbpgq+PFm6hmpKdFlEySSudMjlRZ5r8Uzt5883k43uOocPam8Z6Rxqospe3+Xtou8ZKUKCw50o765vNw3c0nU8Jxg+bwy5tgn/6QllqzGz8sOmNf/ZRNsbgFR1LC8nKVVFGbN8LrN8O4LmHBSc8Mt2O4+D3Y96SUKzjFBvdqR1ZGOpticbIy0i04khKWMzlSRXzyIrz4x+8vTbU9DI79OzTZK9JY1WHM1AUlBWdTLM6YqQssOpISkiVHKo81i+DFK2H+i+G4QXM45kbY9+SUnbn5oR+vwSkeAxYdSQnHkiOVxab18OZoeOt2iBWFO4V3uyi8cyq7ftTpqkVpi4xLW4wsSYnCkiP9lCCAeU/Dy1dD/pLw2O5HQt9boEnNeiZMLB6Uusi4eByLB1HEkqTtSguCoMb+y1RQUEBubi75+fnk5OREHUeJZuXH4S3hi6aF49zWcMxfoMMJNeLSlCQlqrJ+fjuTI/3Yhm9h2s3w7l3hXlMZ2XDoEOgxBLLqRhxOklRWlhypWDwOcx6BV0bA+lXhsb2OC595s9NukUaTJJWfJUcCWPKf8NLUstnhuHH7cCPNPXtFm0uSVGGWHNVsBcvglevhfxPDcVYDOOJK6HIBZNSKNJokacdYclQzbd4Ib4+F6aNg8zogDQ78NRx1LTRoGnU6SVIlsOSoZgkC+OgZePka+PaL8FirLuGlqRYHRZtNklSpLDmqOZb/L9xI84s3w3GDZtD7BvjFqd4SLkkpyJKj1Lf2a3jtRpj1DyCAzNrQfXB4W3hWvajTSZKqiCVHqWtLEbx7J7zxdygqCI/tezL0HgENW0ebTZJU5Sw5Sj1BAB8/F27FULxLeLP9w3U3bbpHGk2SVH0sOUoty+aE5ebz6eG4/q7Q61rY//8gPT3SaJKk6mXJUWooWAav3ghzHuP7dTeXhFsx1JBdwiVJW7PkKLltWgcz7oC3bofN68Njvzg1nL1x3Y0k1WiWHCWneCyctXn1Rli7IjzWqgsccxO07BxtNklSQrDkKPksfDV8mN9XH4bjhq3D593s09/n3UiSSlhylDy+mgt518Knr4Tj2rlw2OVwyO8gMzvabJKkhGPJUeIrWAav/SW8PBXEIb0WHHJ+WHDqNoo6nSQpQVlylLg2FoQLit8eB1s2hMf26Qe9roOd94g2myQp4VlylHi2bIJZD8C0m2H96vBYq67Q50ZodXC02SRJScOSo8QRj8PcJ+HVP3//pOKd94Sjr4e9j3dRsSSpXCw5SgyfvQ6vXA/LZofjervAEVfCQWdBRq0ok0mSkpQlR9FaNicsN5+9Fo6z6kOPS6HrhT6pWJK0Qyw5isbqheGD/OY+GY7Ta8HB50LPy6B+k2izSZJSQkLtWDh+/Hjatm1L7dq16dSpE9OnT9/ua5988kl69+5NkyZNyMnJoVu3brz00kvVmFYVUrAMnr0Uxh78XcFJg/0GwCUzoe/NFhwpCYzOm8+YqQtK/d6YqQsYnTe/mhNJpUuYkjNx4kSGDBnCVVddxezZs+nZsyd9+/Zl8eLFpb7+jTfeoHfv3kyZMoVZs2Zx5JFHcsIJJzB79uxqTq4yWb8m3B18zIEw60EIYtCuDwyaDiffDTvtFnVCSWWUkZ7GqFKKzpipCxiVN5+MdG8SUGJIC4IgiDoEQJcuXTjooIOYMGFCybEOHTrQv39/Ro4cWaafse+++zJgwACuvfbaMr2+oKCA3Nxc8vPzycnJqVBu/YyNBfDOBHh7LBQVhMdadws30GzTPdpskiqsuNAM692ewb3abTOWqlJZP78TYk3Opk2bmDVrFldeeeVWx/v06cOMGTPK9DPi8TiFhYU0arT9J+AWFRVRVFRUMi4oKKhYYP28zRvgP/fAm6Nhw5rwWNNfhOWmXW9vB5eSXHGRGZU3n7GvfsqmWNyCo4STEJerVq1aRSwWo2nTplsdb9q0KStWrCjTz7j11ltZt24dp5122nZfM3LkSHJzc0u+WrVqtUO5VYotRWG5uf0AyLsmLDg77wmn3A8XvAHt+1hwpBQxuFc7sjLS2RSLk5WRbsFRwkmIklMs7UcffkEQbHOsNI8//jjXX389EydOZJdddtnu64YPH05+fn7J15IlS3Y4s74T2xyutbmjE0y5DNaugNzW0G8cXPgudPwVpCfUXzdJO2jM1AUlBWdTLL7dxchSVBLiclXjxo3JyMjYZtZm5cqV28zu/NjEiRM599xzeeKJJzj66KN/8rXZ2dlkZ7tbdaWKbYEP/hVuwVD8lOL6u8Jhl4UP8nN3cCklbW9NDuCMjhJGQpScrKwsOnXqRF5eHieddFLJ8by8PPr167fd9z3++OP89re/5fHHH+e4446rjqgqFtsCH06CabfAmoXhsXpN4NBh0PkcqFUn2nySqkxpi4x/uEbnh2MpSglRcgCGDRvGmWeeSefOnenWrRt33303ixcvZtCgQUB4qenLL7/koYceAsKCc9ZZZ3H77bfTtWvXklmgOnXqkJubG9nvkfLiMfhgErxxC6z+NDxWp1H4lOJDzoesetHmk1TlYvGg1EXGxeNYPCFu2pUS5xZyCB8GeMstt7B8+XI6duzI6NGjOeywwwA4++yz+fzzz3n99dcBOOKII5g2bdo2P2PgwIE8+OCDZfrzvIW8HIpnbt74O6z+7rp7nUbQ/RI45HduwSBJqjZl/fxOqJJT3Sw5ZRDbDP+bGJabbxaFx+rs9INy0yDafJKkGiepnpOjBLSlCOY8Fj7n5tsvwmN1dw7LzcHnWW4kSQnPkqOtbVof3go+4w4oXBYeq9ckXHPT+beuuZEkJQ1LjkIb8+G9e+Ht8bB+VXisQXPoMRgOGghZdaPNJ0lSOVlyarrCr+DdCfDefd/vLdWwDRw6FA44w+fcSJKSliWnplqzKLwkNfsRiH23n1eTvaHHEPjFqZDhXw1JUnLzk6ym+fK/MGMMzHsagnh4rEVn6DkM2vd16wVJUsqw5NQEQQCfToUZt8OiN74/vkev8LLUboe6aaYkKeVYclLZ5o3hvlJvj4OvPw6PpWeGm2V2vwR2/UW0+SRJqkKWnFS0blW4kPi9e2Dd1+GxrPrhhpldL4SGraLNJ0lSNbDkpJLl/4N37wz3lipeTJzTEroOCgtObff0kiTVHJacZBfbAp88D+/cCYtnfH+8RSfodhF0OBEyakWXT5KkiFhyklXhCpj1j/DpxMVPJk7PhH36Q5dB0OrgKNNJkhQ5S04yCQL4fDrMvB8+ehbiW8LjdRtDp7Ph4HMhp3mkESVJShSWnGSw9mt4/7Fw5mbNwu+Pt+oabpa5z4k+mViSpB+x5CSqeAwWvgazH4aPn4f45vB4Vn34xSlhufEWcEmStsuSk2hWfhzO2rw/Edau+P54i07hRpkdfwXZ9aPLJ0lSkrDkJIL8L2HuZPhwEiyb/f3xOo3CWZsDz4Rm+0WXT5KkJGTJiUrhV/Dxc/Dhv+GLGUAQHk/PhHZ9wh3A2x0DmVmRxpQkKVlZcqpLEIRbK3wyBT6eAl/O3Pr7rbuFl6L26Q/1m0QSUZKkVGLJqSrxGHw1F5a8G34tfhfyF2/9muYHwb4nhV9utSBJUqWy5FSFpy6EeU/DprVbH8/Iht0Ph736Qvu+kNMsmnySJNUAlpyqENsUFpysBuGTh1t1gVaHQMtDvDNKkqRqYsmpCocOgx5DYJcOkJ4RdRpJkmokS05VaLpP1AkkSarx0qMOIEmSVBUsOZIkKSVZciRJUkqy5EiSpJRkyZEkSSnJkiNJklKSJUeSJKUkS44kSUpJlhxJkpSSLDmSJCklWXIkSVJKsuRIkqSUZMmRJEkpyZIjSZJSkiVHkiSlJEuOJElKSZYcSZKUkiw5lWR03nzGTF1Q6vfGTF3A6Lz51ZxIkqSazZJTSTLS0xhVStEZM3UBo/Lmk5GeFlEySZJqpsyoA6SKwb3aATDquxmbwb3alRScYb3bl3xfkiRVD0tOJfph0Rn76qdsisUtOJIkRcTLVZVscK92ZGWksykWJysj3YIjSVJELDmVbMzUBSUFZ1Msvt3FyJIkqWp5uaoS/XgNTvEYcEZHkqRqZsmpJKUtMi5tMbIkSaoelpxKEosHpS4yLh7H4kEUsSRJqrHSgiCosZ++BQUF5Obmkp+fT05OTtRxJElSGZT189uFx5IkKSVZciRJUkqy5EiSpJRkyZEkSSnJkiNJklKSJUeSJKUkS44kSUpJlhxJkpSSLDmSJCkl1ehtHYof9lxQUBBxEkmSVFbFn9s/t2lDjS45hYWFALRq1SriJJIkqbwKCwvJzc3d7vdr9N5V8XicZcuW0aBBA9LS0irt5xYUFNCqVSuWLFninlhVzHNdPTzP1cPzXD08z9WjKs9zEAQUFhbSvHlz0tO3v/KmRs/kpKen07Jlyyr7+Tk5Of4HVE0819XD81w9PM/Vw/NcParqPP/UDE4xFx5LkqSUZMmRJEkpyZJTBbKzs7nuuuvIzs6OOkrK81xXD89z9fA8Vw/Pc/VIhPNcoxceS5Kk1OVMjiRJSkmWHEmSlJIsOZIkKSVZciRJUkqy5FTQ+PHjadu2LbVr16ZTp05Mnz79J18/bdo0OnXqRO3atdl999258847qylpcivPeX7yySfp3bs3TZo0IScnh27duvHSSy9VY9rkVt6/08XeeustMjMzOeCAA6o2YIoo73kuKiriqquuok2bNmRnZ7PHHntw//33V1Pa5FXe8/zoo4+y//77U7duXZo1a8Y555zD6tWrqyltcnrjjTc44YQTaN68OWlpaTz11FM/+55q/ywMVG7//Oc/g1q1agX33HNPMG/evODSSy8N6tWrF3zxxRelvv6zzz4L6tatG1x66aXBvHnzgnvuuSeoVatWMGnSpGpOnlzKe54vvfTS4Oabbw7+85//BPPnzw+GDx8e1KpVK/jvf/9bzcmTT3nPdbFvv/022H333YM+ffoE+++/f/WETWIVOc8nnnhi0KVLlyAvLy9YtGhR8O677wZvvfVWNaZOPuU9z9OnTw/S09OD22+/Pfjss8+C6dOnB/vuu2/Qv3//ak6eXKZMmRJcddVVwb///e8ACCZPnvyTr4/is9CSUwGHHHJIMGjQoK2O7b333sGVV15Z6uuvuOKKYO+9997q2AUXXBB07dq1yjKmgvKe59Lss88+wYgRIyo7Wsqp6LkeMGBAcPXVVwfXXXedJacMynueX3jhhSA3NzdYvXp1dcRLGeU9z3/729+C3XfffatjY8aMCVq2bFllGVNNWUpOFJ+FXq4qp02bNjFr1iz69Omz1fE+ffowY8aMUt/z9ttvb/P6Y445hpkzZ7J58+Yqy5rMKnKefywej1NYWEijRo2qImLKqOi5fuCBB1i4cCHXXXddVUdMCRU5z8888wydO3fmlltuoUWLFrRv357LLruMDRs2VEfkpFSR89y9e3eWLl3KlClTCIKAr776ikmTJnHcccdVR+QaI4rPwhq9QWdFrFq1ilgsRtOmTbc63rRpU1asWFHqe1asWFHq67ds2cKqVato1qxZleVNVhU5zz926623sm7dOk477bSqiJgyKnKuFyxYwJVXXsn06dPJzPSfkbKoyHn+7LPPePPNN6lduzaTJ09m1apVXHjhhaxZs8Z1OdtRkfPcvXt3Hn30UQYMGMDGjRvZsmULJ554InfccUd1RK4xovgsdCangtLS0rYaB0GwzbGfe31px7W18p7nYo8//jjXX389EydOZJdddqmqeCmlrOc6FotxxhlnMGLECNq3b19d8VJGef5Ox+Nx0tLSePTRRznkkEM49thjGTVqFA8++KCzOT+jPOd53rx5DB48mGuvvZZZs2bx4osvsmjRIgYNGlQdUWuU6v4s9P+ClVPjxo3JyMjY5v8RrFy5cpuGWmzXXXct9fWZmZnsvPPOVZY1mVXkPBebOHEi5557Lk888QRHH310VcZMCeU914WFhcycOZPZs2dz8cUXA+GHcRAEZGZm8vLLL3PUUUdVS/ZkUpG/082aNaNFixbk5uaWHOvQoQNBELB06VLatWtXpZmTUUXO88iRI+nRoweXX345APvttx/16tWjZ8+e3Hjjjc62V5IoPgudySmnrKwsOnXqRF5e3lbH8/Ly6N69e6nv6dat2zavf/nll+ncuTO1atWqsqzJrCLnGcIZnLPPPpvHHnvM6+llVN5znZOTwwcffMCcOXNKvgYNGsRee+3FnDlz6NKlS3VFTyoV+Tvdo0cPli1bxtq1a0uOzZ8/n/T0dFq2bFmleZNVRc7z+vXrSU/f+uMwIyMD+H6mQTsuks/CKlvSnMKKb0+87777gnnz5gVDhgwJ6tWrF3z++edBEATBlVdeGZx55pklry++bW7o0KHBvHnzgvvuu89byMugvOf5scceCzIzM4Nx48YFy5cvL/n69ttvo/oVkkZ5z/WPeXdV2ZT3PBcWFgYtW7YMTjnllGDu3LnBtGnTgnbt2gXnnXdeVL9CUijveX7ggQeCzMzMYPz48cHChQuDN998M+jcuXNwyCGHRPUrJIXCwsJg9uzZwezZswMgGDVqVDB79uySW/UT4bPQklNB48aNC9q0aRNkZWUFBx10UDBt2rSS7w0cODA4/PDDt3r966+/Hhx44IFBVlZWsNtuuwUTJkyo5sTJqTzn+fDDDw+Abb4GDhxY/cGTUHn/Tv+QJafsynueP/roo+Doo48O6tSpE7Rs2TIYNmxYsH79+mpOnXzKe57HjBkT7LPPPkGdOnWCZs2aBb/+9a+DpUuXVnPq5PLaa6/95L+5ifBZmBYEzsVJkqTU45ocSZKUkiw5kiQpJVlyJElSSrLkSJKklGTJkSRJKcmSI0mSUpIlR5IkpSRLjiRJSkmWHEmSlJIsOZKSXjweZ++99+aPf/zjVsdfeuklsrKyeOKJJyJKJilKlhxJSS89PZ3hw4czYcIEvvnmGwDef/99Tj31VG666SZOPfXUiBNKioJ7V0lKCVu2bKF9+/YMHDiQc889l65du9KvXz/GjRsXdTRJEbHkSEoZd911F3/6059o0aIFu+22G5MnTyYjIyPqWJIiYsmRlDLWrl1LkyZN2HPPPXnnnXeoV6/eVt9/5ZVX+OCDDxg6dGhECSVVJ0uOpJRx9tlnM3HiRBo2bMiiRYuoXbt21JEkRciFx5JSwjXXXMPzzz/PO++8w5YtW7jvvvu2eU3fvn356KOPIkgnKQqWHElJ79577+XWW2/l2WefZf/992fIkCHccsstbN68eavXLViwgHbt2kWUUlJ1s+RISmovvPACF110EY888ghdu3YF4JJLLqGgoICHH3645HX5+fnUr1+fzMzMqKJKqmaWHElJa9asWZx66qnccsstnHzyySXHc3JyuOSSS/jrX/9KLBYD4MMPP2TfffeNKqqkCLjwWFKNcNddd7FmzRqGDx8edRRJ1cSZHEk1wty5c+nYsWPUMSRVI2dyJElSSnImR5IkpSRLjiRJSkmWHEmSlJIsOZIkKSVZciRJUkqy5EiSpJRkyZEkSSnJkiNJklKSJUeSJKUkS44kSUpJlhxJkpSSLDmSJCkl/T9bxwpX7Qf61AAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x_range = np.linspace(0,1,100)\n", "plt.plot(xi,yi,'x');\n", "plt.plot(x_range,result.x[0] * x_range ** result.x[1] + result.x[2])\n", "plt.xlabel('$X_i$')\n", "plt.ylabel('$Y_i$');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To solve uncontrained optimisation problem like:\n", "\n", "$$\n", " \\mathop {\\min }\\limits_x f\\left( x \\right)\n", "$$\n", "\n", "with $x$ the design variable of length $n$, you can make use of `scipy.optimize.minimize(fun, x0, ...)` which outputs `scipy.optimize.OptimizeResult` with the problem's solution, some characteristics of the problem and the performance of the optimization algorithm." ] }, { "cell_type": "markdown", "metadata": {}, "source": [] } ], "metadata": { "kernelspec": { "display_name": "base", "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.9.17" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }