manual_teachbooks/book/examples/Unconstrained_optimization_stripped.ipynb
2024-02-26 22:57:13 +01:00

540 lines
107 KiB
Text

{
"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": [
"<Figure size 640x480 with 1 Axes>"
]
},
"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<n<100$. It shows a clear minimum which coincides with the minimum found by `scipy.optimize.minimize`"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAGwCAYAAABmTltaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABU20lEQVR4nO3deVxU5eIG8GcWGPaRRUAEBNwVV3AB9zSXXPJmmlouLZaZmVvXrPtrd2lftLSs1NJSuy6pdc3dNHEHRBFXBGQREZhhH5h5f38gU6gpozOcmeH5fj7z+cSZF3w4t+s8veec95UJIQSIiIiIqMbkUgcgIiIisjUsUEREREQmYoEiIiIiMhELFBEREZGJWKCIiIiITMQCRURERGQiFigiIiIiEymlDmCPDAYDMjIy4O7uDplMJnUcIiIiqgEhBAoKChAQEAC5/M5zTCxQFpCRkYGgoCCpYxAREdE9SEtLQ2Bg4B3HsEBZgLu7O4DK/wE8PDwkTkNEREQ1odVqERQUZPwcvxMWKAuoumzn4eHBAkVERGRjanL7DW8iJyIiIjIRCxQRERGRiVigiIiIiEzEAkVERERkIhYoIiIiIhOxQBERERGZiAWKiIiIyEQsUEREREQmYoEiIiIiMhELFBEREZGJWKCIiIiITMQCRURERGQiFigiIiKyKUcv56KsQi9pBhYoIiIishkZ+SUYuTQGke/uRFFZhWQ5WKCIiIjIZuxKygYANPNzh6tKKVkOFigiIiKyGbvOXAUA9G3pK2kOFigiIiKyCUVlFTh44ToA4MGWfpJmYYEiIiIim7D/fA50egOCvVzQxNdN0iwsUERERGQT/n75TiaTSZqFBYqIiIisnt4gsPvGDeRSX74DWKCIiIjIBsSl5eN6kQ7uTkp0CvWSOg4LFBEREVm/qst3vZv7wkEhfX2RPgERERHRXey8UaD6Sbx8QRUWKCIiIrJqabnFOHe1EAq5DL2bsUARERER3VXV7FOnEE+oXRwkTlOJBYqIiIis2l+X76R/+q4KCxQRERFZLW1pOQ5fygUA9GWBIiIiIrq7P85dQ4VBoHF9V4T6uEodx4gFioiIiKzWzsQbl+9aWc/sE8ACRURERFaqXG8wrj5uTfc/ASxQREREZKWOJudCW1oBb1dHdAz2lDpONSxQREREZJW2J/61ebBCLu3mwTdjgSIiIiKrI4TAjhsFqn8rf4nT3IoFioiIiKzO6Qwt0vNL4OygQPemPlLHuQULFBEREVmdqst3PZv5wMlBIXGaW7FAERERkdWpunz3oBVevgNYoIiIiMjKpOUW40ymFnIZ0LeFdWwefDMWKCIiIrIqVbNPnUK84OnqKHGa22OBIiIiIquyPTELANC/tXVevgNYoIiIiMiK5BXpcPRyHgCgv5Vt3/J3LFBERERkNXYnZUNvEGjh744gLxep4/wjFigiIiKyGsbLd1Y8+wSwQBEREZGVKC3X449zOQCs+/4ngAWKiIiIrMSB8zkoKdcjQO2E1gEeUse5I5spUAsWLECnTp3g7u4OX19fDB8+HGfPnq02RgiBN998EwEBAXB2dkbv3r1x+vTpamPKysrw4osvwsfHB66urhg2bBiuXLlSbUxeXh7GjRsHtVoNtVqNcePGIT8/39K/IhERUZ227fRfT9/JZNa1efDNbKZA7du3Dy+88AIOHTqEHTt2oKKiAv3790dRUZFxzPvvv4+PP/4YixcvxtGjR+Hv748HH3wQBQUFxjHTp0/Hxo0bsWbNGhw4cACFhYUYMmQI9Hq9cczYsWMRFxeHbdu2Ydu2bYiLi8O4ceNq9fclIiKqSyr0Buw8U7n+0wArv3wHABA2Kjs7WwAQ+/btE0IIYTAYhL+/v1i4cKFxTGlpqVCr1WLp0qVCCCHy8/OFg4ODWLNmjXFMenq6kMvlYtu2bUIIIRITEwUAcejQIeOYmJgYAUAkJSXVKJtGoxEAhEajue/fk4iIqC748/w10WjOVtHh7e2ivEIvSQZTPr9tZgbqZhqNBgDg5eUFAEhOTkZWVhb69+9vHKNSqdCrVy8cPHgQAHD8+HGUl5dXGxMQEIDw8HDjmJiYGKjVanTp0sU4pmvXrlCr1cYxNysrK4NWq632IiIiopqrunz3YEs/KBXWX0+sP+FtCCEwc+ZMdO/eHeHh4QCArKzKE+/nV/2xRz8/P+N7WVlZcHR0hKen5x3H+Preuu+Or6+vcczNFixYYLxfSq1WIygo6P5+QSIiojrEYBD4/UaBGhBu3csXVLHJAjV16lScPHkSP/300y3v3XzTmRDirjei3TzmduPv9HPmzp0LjUZjfKWlpdXk1yAiIiIAcVfycVVbBjeVEtGNfaSOUyM2V6BefPFFbN68GXv27EFgYKDxuL9/5Q1nN88SZWdnG2el/P39odPpkJeXd8cxV69eveXPvXbt2i2zW1VUKhU8PDyqvYiIiKhmqmaf+rTwhZODQuI0NWMzBUoIgalTp2LDhg3YvXs3QkNDq70fGhoKf39/7Nixw3hMp9Nh3759iI6OBgBERETAwcGh2pjMzEycOnXKOCYqKgoajQZHjhwxjjl8+DA0Go1xDBEREZmHEAK/n6osUANt4em7G5RSB6ipF154AT/++CN++eUXuLu7G2ea1Go1nJ2dIZPJMH36dMyfPx9NmzZF06ZNMX/+fLi4uGDs2LHGsU8//TRmzZoFb29veHl5Yfbs2WjTpg369esHAGjZsiUGDhyISZMm4auvvgIAPPvssxgyZAiaN28uzS9PRERkp85eLcDl68VwVMrRu3l9qePUmM0UqCVLlgAAevfuXe348uXLMXHiRADAv//9b5SUlGDKlCnIy8tDly5dsH37dri7uxvHf/LJJ1AqlRg1ahRKSkrQt29frFixAgrFX1OGq1evxrRp04xP6w0bNgyLFy+27C9IRERUB227MfvUs6kPXFU2U0sgE0IIqUPYG61WC7VaDY1Gw/uhiIiI7mDQZ/txJlOLDx5ti5GR0j7Fbsrnt83cA0VERET2JeV6Ec5kaqGQy9CvpW0sX1CFBYqIiIgkUfX0XdcwL3i6OkqcxjQsUERERCSJbTb49F0VFigiIiKqdVmaUpxIzQcA9GeBIiIiIrq7/53KBABENvKEn4eTxGlMxwJFREREte5/CZWX7wa1aSBxknvDAkVERES1KltbiqMpuQCAgeG2d/kOYIEiIiKiWvb76SwIAbQPqoeG9ZyljnNPWKCIiIioVv124/LdYBu9fAewQBEREVEtyiksw+Hk6wBs9/IdwAJFREREtWj76aswCKBtoBpBXi5Sx7lnLFBERERUa6qWLxgUbruX7wAWKCIiIqoleUU6HLxYeflukA1fvgNYoIiIiKiW7Ei8Cr1BoFUDD4T4uEod576wQBEREVGt+DWh8vLdQ21se/YJYIEiIiKiWqApLsefF3IA2O7q43/HAkVEREQWt+PMVVQYBJr7uaNxfTep49w3FigiIiKyuK0nMwAAD9nB7BPAAkVEREQWllekw4HzlZfvBrdlgSIiIiK6q+2JWagwCLTwd0cTX9u/fAewQBEREZGFbT1Z+fTd0HYBEicxHxYoIiIispjrhWXGxTNtefPgm7FAERERkcVsO50FvUEgvKHtL575dyxQREREZDFb4ysv3w1paz+X7wAWKCIiIrKQ7IJSHE62v8t3AAsUERERWci2U1kwCKBdUD0EeblIHcesWKCIiIjIIqou3w21k7Wf/o4FioiIiMwuS1OKoym5AOxn9fG/Y4EiIiIis/stIRNCABGNPBFQz1nqOGbHAkVERERmV7X3nb3dPF6FBYqIiIjM6kpeMU6k5kMms5+9727GAkVERERmteXGzeNdQr3g5+EkcRrLYIEiIiIis/olLh0AMKxdQ4mTWA4LFBEREZnN+asFSMoqgFIuw6Bwf6njWAwLFBEREZnN5vjKm8d7NasPT1dHidNYDgsUERERmYUQwlighrW3r73vbsYCRURERGZx8ooGKdeL4eQgR7+WflLHsSgWKCIiIjKLqtmnfi394KpSSpzGsligiIiI6L7pDQJbqi7ftbPvy3cACxQRERGZweHk68guKIOHkxK9mteXOo7FsUARERHRfauafRoU3gAqpULiNJbHAkVERET3RVdhwG8JWQDs/+m7KixQREREdF/2n78GTUk56rur0DXMW+o4tYIFioiIiO7LL3GVl+8Gt2kAhVwmcZrawQJFRERE96ywrALbEysv3/2rg/3ufXczFigiIiK6Z7+fykJpuQFhPq5oG6iWOk6tYYEiIiKie7YpLh0AMLxDQ8hkdePyHcACRURERPfoqrYUf17IAQAMb193Lt8BLFBERER0j7bEZ8AggIhGngj2dpE6Tq1igSIiIqJ7sjH2r8t3dQ0LFBEREZns3NUCnM7QQimXYUibBlLHqXUsUERERGSyTTdmn3o394Wnq6PEaWofCxQRERGZxGAQxsUz69LaT3/HAkVEREQmOXo5F+n5JXBXKdG3pa/UcSTBAkVEREQm2XRj9mlQG384OSgkTiMNFigiIiKqsbIKPX49WVmg6uLTd1VsqkD98ccfGDp0KAICAiCTybBp06Zq7wsh8OabbyIgIADOzs7o3bs3Tp8+XW1MWVkZXnzxRfj4+MDV1RXDhg3DlStXqo3Jy8vDuHHjoFaroVarMW7cOOTn51v4tyMiIrJ+u85kQ1taAX8PJ3QJ9ZY6jmRsqkAVFRWhXbt2WLx48W3ff//99/Hxxx9j8eLFOHr0KPz9/fHggw+ioKDAOGb69OnYuHEj1qxZgwMHDqCwsBBDhgyBXq83jhk7dizi4uKwbds2bNu2DXFxcRg3bpzFfz8iIiJrt/545aTDvzo2hEJed7ZuuYWwUQDExo0bjV8bDAbh7+8vFi5caDxWWloq1Gq1WLp0qRBCiPz8fOHg4CDWrFljHJOeni7kcrnYtm2bEEKIxMREAUAcOnTIOCYmJkYAEElJSTXKptFoBACh0Wju51ckIiKyKtnaUhE291fRaM5Wcf5qgdRxzM6Uz2+bmoG6k+TkZGRlZaF///7GYyqVCr169cLBgwcBAMePH0d5eXm1MQEBAQgPDzeOiYmJgVqtRpcuXYxjunbtCrVabRxzs7KyMmi12movIiIie7M5PgN6g0C7oHpo4usmdRxJ2U2BysrKAgD4+flVO+7n52d8LysrC46OjvD09LzjGF/fWx/J9PX1NY652YIFC4z3S6nVagQFBd3370NERGRtqi7fPdqx7t48XsVuClQVmaz69VghxC3HbnbzmNuNv9PPmTt3LjQajfGVlpZ2D8mJiIis15lMLRIztXBQyDCkbYDUcSRnNwXK398fAG6ZJcrOzjbOSvn7+0On0yEvL++OY65evXrLz7927dots1tVVCoVPDw8qr2IiIjsSdXsU98WfnVy65ab2U2BCg0Nhb+/P3bs2GE8ptPpsG/fPkRHRwMAIiIi4ODgUG1MZmYmTp06ZRwTFRUFjUaDI0eOGMccPnwYGo3GOIaIiKguqdAbjItnjogIlDiNdVBKHcAUhYWFuHDhgvHr5ORkxMXFwcvLC8HBwZg+fTrmz5+Ppk2bomnTppg/fz5cXFwwduxYAIBarcbTTz+NWbNmwdvbG15eXpg9ezbatGmDfv36AQBatmyJgQMHYtKkSfjqq68AAM8++yyGDBmC5s2b1/4vTUREJLH953OQU1gGb1dH9G5eX+o4VsGmCtSxY8fQp08f49czZ84EAEyYMAErVqzAv//9b5SUlGDKlCnIy8tDly5dsH37dri7uxu/55NPPoFSqcSoUaNQUlKCvn37YsWKFVAo/lqKfvXq1Zg2bZrxab1hw4b949pTRERE9u6/Jyov3w1rHwAHhd1cvLovMiGEkDqEvdFqtVCr1dBoNLwfioiIbJqmuByd5u+ErsKArS92R3hDtdSRLMaUz2/WSCIiIvpHWxMyoKswoLmfO1oHcFKgCgsUERER/aP/3nj6bkREw7suC1SXsEARERHRbV3ILkBsaj4UchmGd+DimX/HAkVERES39fOxytmnPs194evuJHEa68ICRURERLco1xuw/kQ6AGBUJNd+uhkLFBEREd1i79lryCksg4+bI/q0uHWP2LqOBYqIiIhuse5Y5b6uj3QM5NpPt8EzQkRERNVcKyjDnqRsAMBIbt1yWyxQREREVM2m2HRUGATaB9VDUz/3u39DHcQCRUREREZCCOPlu1GRQRKnsV4mF6gHHngA+fn5txzXarV44IEHzJGJiIiIJBKXlo/z2YVwcpBjSLsGUsexWiYXqL1790Kn091yvLS0FPv37zdLKCIiIpLGzzdWHh8U3gAeTg4Sp7FeypoOPHnypPGfExMTkZWVZfxar9dj27ZtaNiQq5QSERHZqhKdHlviMgAAI7n20x3VuEC1b98eMpkMMpnstpfqnJ2dsWjRIrOGIyIiotrzW0ImCsoqEOTljK6h3lLHsWo1LlDJyckQQiAsLAxHjhxB/fr1je85OjrC19cXCoXCIiGJiIjI8tYcTQUAPBYZBLmcGwffSY0LVKNGjQAABoPBYmGIiIhIGheyC3D0ch4UchlG8um7u6pRgdq8eTMGDRoEBwcHbN68+Y5jhw0bZpZgREREVHvWHKlcuqBPc1/4eXDj4LupUYEaPnw4srKy4Ovri+HDh//jOJlMBr1eb65sREREVAvKKvRYf6Ly6bsxnTn7VBM1KlB/v2zHS3hERET2Zfvpq8grLoe/hxN6Nat/92+gmq0D5eXlhZycHADAU089hYKCAouGIiIiotpTdfP4qMhAKLlxcI3U6CzpdDpotVoAwMqVK1FaWmrRUERERFQ7Uq4X4c8L1yGTAaM68fJdTdXoEl5UVBSGDx+OiIgICCEwbdo0ODs733bsd999Z9aAREREZDlrj1bePN6jaX0EerpInMZ21KhArVq1Cp988gkuXrwImUwGjUbDWSgiIiIbV643GLduGcPZJ5PUqED5+flh4cKFAIDQ0FD88MMP8PbmCqVERES2bE9SNq4VlMHb1RF9W/pJHcem1HghzSrJycmWyEFERES17KcjlTePPxoRCEclbx43xT2drX379mHo0KFo0qQJmjZtimHDhmH//v3mzkZEREQWkpZbjL3nrgEARncOljiN7TG5QK1atQr9+vWDi4sLpk2bhqlTp8LZ2Rl9+/bFjz/+aImMREREZGZrjqZCCKB7Ex+E+rhKHcfmyIQQwpRvaNmyJZ599lnMmDGj2vGPP/4Yy5Ytw5kzZ8wa0BZptVqo1WpoNBp4eHhIHYeIiKgaXYUB0Qt3I6ewDEse74hBbRpIHckqmPL5bfIM1KVLlzB06NBbjg8bNoz3RxEREdmA7YlZyCksg6+7Cv1a8ebxe2FygQoKCsKuXbtuOb5r1y4EBfERSCIiImu3+lDlzeOjOwXBgSuP3xOTn8KbNWsWpk2bhri4OERHR0Mmk+HAgQNYsWIFPvvsM0tkJCIiIjO5kF2ImEvXIZcBj/Hm8XtmcoF6/vnn4e/vj48++gjr1q0DUHlf1Nq1a/Hwww+bPSARERGZz4+HK2efHmjhi4b1br+rCN2dSQWqoqIC8+bNw1NPPYUDBw5YKhMRERFZQGm5Hv89Xrl1y+NdGkmcxraZdOFTqVTigw8+gF6vt1QeIiIispCtJzOhLa1AoKczejarL3Ucm2bynWP9+vXD3r17LRCFiIiILGnVoRQAwJjOwVDIZRKnsW0m3wM1aNAgzJ07F6dOnUJERARcXasvvjVs2DCzhSMiIiLzOJWuQVxaPpRyGUZF8qn5+3VPN5EDlQtn3kwmk/HyHhERkRWqmn0aEO6P+u4qidPYPpMLlMFgsEQOIiIispD8Yh02xaUDACZGh0gbxk5w9SwiIiI79/OxKygtN6BlAw9ENvKUOo5duKcCtWvXLgwZMgSNGzdGkyZNMGTIEOzcudPc2YiIiOg+6Q0CP9y4fDc+qhFkMt48bg4mF6jFixdj4MCBcHd3x0svvYRp06bBw8MDDz30EBYvXmyJjERERHSP9p3LRmpuMTyclHi4fYDUceyGyfdALViwAJ988gmmTp1qPDZt2jR069YN8+bNq3aciIiIpPV9TOXs06jIILg4mvyxT//A5BkorVaLgQMH3nK8f//+0Gq1ZglFRERE9+9yThH2nr0GmQx4oitXHjcnkwvUsGHDsHHjxluO//LLLxg6dKhZQhEREdH9q7r3qVez+gjxcb3LaDKFyXN5LVu2xLx587B3715ERUUBAA4dOoQ///wTs2bNwueff24cO23aNPMlJSIiohor1lVg3bHKfe8mRIVIG8YOyYQQwpRvCA0NrdkPlslw6dKlewpl67RaLdRqNTQaDTw8PKSOQ0REddCPh1Px6sYEBHu5YO/s3pBz65a7MuXz2+QZqOTk5HsORkRERJYnhMD3MZcBVC5dwPJkflxIk4iIyM7EXLqOpKwCODsoMDKC+95ZAgsUERGRnfnuwGUAwIiIhlC7OEgbxk6xQBEREdmRlOtF2JV0FQAwMbpm9y2T6VigiIiI7MiKg5chROXSBU183aSOY7dYoIiIiOxEQWk5fj52BQDwVHfOPlnSPRWo/fv344knnkBUVBTS09MBAD/88AMOHDhg1nBERERUcz8fu4LCsgo08XVDz6Y+UsexayYXqPXr12PAgAFwdnZGbGwsysrKAAAFBQWYP3++2QMSERHR3ekNAitvLF0wMToEMhmXLrAkkwvUu+++i6VLl2LZsmVwcPjrzv7o6GicOHHCrOGIiIioZnYnZSPlejHUzg54pGNDqePYPZML1NmzZ9GzZ89bjnt4eCA/P98cmYiIiMhEy/+sXOh6TOdguDiavE42mcjkAtWgQQNcuHDhluMHDhxAWFiYWUJZiy+//BKhoaFwcnJCREQE9u/fL3UkIiKiW5zJ1OLgxetQyGUYH9VI6jh1gskF6rnnnsNLL72Ew4cPQyaTISMjA6tXr8bs2bMxZcoUS2SUxNq1azF9+nS89tpriI2NRY8ePTBo0CCkpqZKHY2IiKiabw9Uzj4NbO2PgHrOEqepG0zeTBgAXnvtNXzyyScoLS0FAKhUKsyePRvvvPOO2QNKpUuXLujYsSOWLFliPNayZUsMHz4cCxYsuOP3WnIz4YvXCgEAjetzbQ8iIgKuakvR/b3dKNcLbJwSjQ7BnlJHslkW3UwYAObNm4fXXnsNiYmJMBgMaNWqFdzc7OcDXafT4fjx43jllVeqHe/fvz8OHjx4y/iysjLj04hA5f8AlnAqXYMnvj0MV0clNk6Jhq+Hk0X+HCIish0rD15GuV6gU4gny1MtMvkS3sqVK1FUVAQXFxdERkaic+fOdlWeACAnJwd6vR5+fn7Vjvv5+SErK+uW8QsWLIBarTa+goIss3FjA7UT6jk7ID2/BBOWH0VBablF/hwiIrINRWUVWHUoBQDwTA/7ug/Z2plcoGbPng1fX1+MHj0aW7duRUVFhSVyWYWb19AQQtx2XY25c+dCo9EYX2lpaRbJ4+2mwvdPdYGPmyPOZGoxedVx6CoMFvmziIjI+v18LA3a0gqE+riiX0u/u38DmY3JBSozMxNr166FQqHA6NGj0aBBA0yZMuW2l7ZslY+PDxQKxS2zTdnZ2bfMSgGV94B5eHhUe1lKsLcLlk/sDFdHBf68cB2zf46HwWDybWxERGTj9AaBb28sXfBU91Ao5Fw4szaZXKCUSiWGDBmC1atXIzs7G59++ilSUlLQp08fNG7c2BIZa52joyMiIiKwY8eOasd37NiB6OhoiVL9pU2gGkueiIBSLsPm+Aws+N8ZqSMREVEt+/10FtJyS+Dp4oBHOwZKHafOua+VtlxcXDBgwADk5eUhJSUFZ87Yzwf5zJkzMW7cOERGRiIqKgpff/01UlNTMXnyZKmjAQB6NquP9x9ti5nr4rFsfzL8PJx4/ZuIqI4QQuDrPy4BAMZ1bQRnR4XEieqeeypQxcXF2LhxI1avXo2dO3ciKCgIY8aMwc8//2zufJJ57LHHcP36dbz99tvIzMxEeHg4fvvtNzRqZD0LlD3SMRDZBWVY+L8kvPvrGfh6OGFYuwCpYxERkYUdT8lDXFo+HJVyjIsKkTpOnWRygRozZgy2bNkCFxcXjBw5Env37rWKy1qWMGXKFKtfHPS5nmHI0pRixcHLmLUuDj6ujohuwh24iYjsWdXs0yMdGqK+u0riNHWTyQVKJpNh7dq1GDBgAJRK7rUjNZlMhteHtMK1gjL8mpCJZ384jrXPdUXrALXU0YiIyAIuXivEjjNXAQDP9AiVOE3dZfJN5D/++CMGDx7M8mRF5HIZPhrVDl1CvVBYVoGJy48iLbdY6lhERGQBX++7BCGAfi390MTXXeo4dVaNWtDnn3+OZ599Fk5OTvj888/vOHbatGlmCUamcXJQ4OvxkXjsqxgkZRVgwvIj+O/kaHi5OkodjYiIzCRLU4oNsVcAAM/3to8n321VjfbCCw0NxbFjx+Dt7Y3Q0H+eLpTJZLh06ZJZA9oiS+6FdzdZmlKMWHIQ6fkl6BBcDz8+05VPZxAR2Yn5v53B139cQudQL6x7LkrqOHbH7HvhJScn3/afyfr4q52w8qlOGLEkBrGp+Zj64wl8NS4CSoXJV2uJiMiKaIrLsfrGti3P9+Lsk9RM/lR9++23UVx86/01JSUlePvtt80Siu5PE193fDcxEiqlHLuSsvHaxlOowUQjERFZsR8OXUaRTo8W/u7o3by+1HHqPJML1FtvvYXCwsJbjhcXF+Ott94ySyi6fxGNvLBoTAfIZcDaY2n4ZMc5qSMREdE9Ki3XY/mflwFU3vt0u31ZqXaZXKD+aUPd+Ph4eHl5mSUUmUf/1v54d3gbAMDnuy8Yd+wmIiLb8vOxNFwv0iHQ0xmD2zSQOg7BhHWgPD09IZPJIJPJ0KxZs2olSq/Xo7Cw0Gq2OaG/jO0SjKvaUny26zxe/+UUfNxUGBjuL3UsIiKqoQq9AV/dWDjz2Z5hvKfVStS4QH366acQQuCpp57CW2+9BbX6r4UaHR0dERISgqgoPhFgjab3a4rsgjL8dCQV09bEYtXTXdA5lLOFRES24NeETFzJK4G3qyNGRgRJHYduqHGBmjBhAoDKJQ2io6Ph4OBgsVBkXjKZDO883Bo5hWXYkXgVz6w8ip8nR6O5PxdgIyKyZgaDwBd7LgAAnuwWwmVprIjJ84C9evUylqeSkhJotdpqL7JOSoUci8Z0QEQjT2hLKzDhuyPIyC+ROhYREd3B9sSrOHe1EO5OSoyPDpE6Dv2NyQWquLgYU6dOha+vL9zc3ODp6VntRdbLyUGBbydEoqmvG7K0pRj/3RHkF+ukjkVERLchhMDiPecBABOjQ+DhxCs/1sTkAvXyyy9j9+7d+PLLL6FSqfDNN9/grbfeQkBAAL7//ntLZCQzqufiiJVPdUYDtRMuZBfi6ZXHUKLTSx2LiIhusvfcNZxK18LFUYEnu3HTYGtjcoHasmULvvzySzz66KNQKpXo0aMH/vOf/2D+/PlYvXq1JTKSmQXUc8bKpzrDw0mJ4yl5ePGnE6jQG6SORURENwghsGhX5ezTE10bcV9TK2RygcrNzTXuh+fh4YHc3FwAQPfu3fHHH3+YNx1ZTDM/d3w7sRNUSjl2nuFq5URE1iTm0nWcSM2Ho1KOZ3pw9skamVygwsLCcPnyZQBAq1atsG7dOgCVM1P16tUzZzaysE4hXvj8b6uVf8zVyomIrMLi3ZVP3o3pFARfdyeJ09DtmFygnnzyScTHxwMA5s6da7wXasaMGXj55ZfNHpAsa8DfVitftPsCvo+5LG0gIqI67nhKLg5evA4HhQzPctNgq1XjdaCqzJgxw/jPffr0QVJSEo4dO4bGjRujXbt2Zg1HtWNsl2BkF5Ti053n8cbm0/B2VWFwW24VQEQkharZpxEdA9GwnrPEaeifmFygbhYcHIzg4GBzZCEJvdS3Ka4VlGH14VTMWBsHT1cHRDf2kToWEVGdEp+Wjz1nr0EuAyZz9smqmVygPv/889sel8lkcHJyQpMmTdCzZ08oFFwt1ZbIZDK8/XA4rhfqsO10Fp79/jjWPNsV4Q3Vd/9mIiIyi093Vt6LOrxDQ4T4uEqchu5EJkx89Co0NBTXrl1DcXExPD09IYRAfn4+XFxc4ObmhuzsbISFhWHPnj0ICqqbe/ZotVqo1WpoNBp4eHhIHcckpeV6TPjuCA4n58LHTYUNz0cj2NtF6lhERHYvNjUP//ryIBRyGXbN7MUCJQFTPr9Nvol8/vz56NSpE86fP4/r168jNzcX586dQ5cuXfDZZ58hNTUV/v7+1e6VItvh5KDA1+Mj0cLfHTmFZRj/3WHkFJZJHYuIyO59urNy3ad/cfbJJpg8A9W4cWOsX78e7du3r3Y8NjYWI0aMwKVLl3Dw4EGMGDECmZmZ5sxqM2x5BqrKVW0pRiw5iCt5JQhv6IGfJnWFO7cRICKyiOMpeRixpHL2afesXmjkzQIlBYvOQGVmZqKiouKW4xUVFcjKygIABAQEoKCgwNQfTVbEz8MJ3z/VGd6ujjiVrsVzPxxHWQW3fCEisoSqe58e6dCQ5clGmFyg+vTpg+eeew6xsbHGY7GxsXj++efxwAMPAAASEhKMq5WT7Qqr74blT3aCq6MCBy9ex8y18dAbuFo5EZE5HU/Jxf7zOVDKZXjxgaZSx6EaMrlAffvtt/Dy8kJERARUKhVUKhUiIyPh5eWFb7/9FgDg5uaGjz76yOxhqfa1DayHr8ZFwkEhw68JmXhz82lu+UJEZEaf7Ki892lEx0A+tGNDTL4HqkpSUhLOnTsHIQRatGiB5s2bmzubzbKHe6ButvVkBl78KRZCANP7NcX0fs2kjkREZPOOXs7FyKUxUMpl2DO7N4K8WKCkZMrn9z0vpBkWFgaZTIbGjRtDqbzv9TjJyg1pG4DcIh1e/+U0Pt15Ht5uKozr2kjqWERENksIgY+3V977NDIykOXJxph8Ca+4uBhPP/00XFxc0Lp1a6SmpgIApk2bhoULF5o9IFmP8VEhmNa38vr867+cwtaTGRInIiKyXX9euI6YS9fhqJBjKu99sjkmF6i5c+ciPj4ee/fuhZPTXztE9+vXD2vXrjVrOLI+M/o1xeNdgiEEMGNtHA6cz5E6EhGRzRFC4IPfkwBU7kfKPe9sj8kFatOmTVi8eDG6d+8OmUxmPN6qVStcvHjRrOHI+lRt+TK4TQOU6wWe/eEY4tPypY5FRGRTtideRfwVDZwdFHihTxOp49A9MLlAXbt2Db6+vrccLyoqqlaoyH4p5DJ8/Fg7dGvijWKdHk+uOIqL1wqljkVEZBP0BoGPtp8FADzVPQT13VUSJ6J7YXKB6tSpE3799Vfj11WladmyZYiKijJfMrJqKqUCX42LRNtANXKLdBj3zWFk5JdIHYuIyOptjk/HuauF8HBS4tkejaWOQ/fI5MfnFixYgIEDByIxMREVFRX47LPPcPr0acTExGDfvn2WyEhWyk2lxPKJnTDyqxhculaEcd8exs+To+Hl6ih1NCIiq1SuNxjXfXquV2OoXbhFlq0yeQYqOjoaf/75J4qLi9G4cWNs374dfn5+iImJQUREhCUykhXzdlPhh6e7oIHaCRevFeHJ5UdQWHbrVj9ERASsO5aG1Nxi+Lg54sluIVLHoftwzwtp0j+zx4U07+ZCdgFGLo1BXnE5ujfxwbcTI6FSKqSORURkNUrL9ej1wR5c1ZbhzaGtMLEbtzyzNhbdTJjodpr4umPFk53h4qjAgQs5mLE2jvvmERH9zfI/L+OqtgwN6zljTJdgqePQfapxgZLL5VAoFHd8cUXyuq1dUD0sGx8JR4UcvyVk4bWNCdw3j4gIQF6RDl/uvQAAmPlgM87Q24EaN56NGzf+43sHDx7EokWL+GFJ6NbEB5+PaY8pq09gzdE0qF0cMHdQS6ljERFJavGeCygorUDLBh4Y3qGh1HHIDGpcoB5++OFbjiUlJWHu3LnYsmULHn/8cbzzzjtmDUe2aWB4Ayx8pC3+vf4kvtp3CfWcHfF8bz6qS0R1U1puMX6ISQEAvDKoBRRyrploD+7pHqiMjAxMmjQJbdu2RUVFBeLi4rBy5UoEB/OaLlUa1SkIrz1UOfP03rYk/HQkVeJERETS+Gj7Wej0BnRr4o2eTX2kjkNmYlKB0mg0mDNnDpo0aYLTp09j165d2LJlC8LDwy2Vj2zYpJ5heKFP5czTqxsTuPkwEdU5p9I12BRX+XffKwNbcscOO1LjAvX+++8jLCwMW7duxU8//YSDBw+iR48elsxGdmB2/+bVNh/eczZb6khERLXmvW2VGwYPaxeANoFqidOQOdV4HSi5XA5nZ2f069cPCsU/Pz2wYcMGs4WzVXVxHag70RsEZqyNw+b4DDg5yPH9U13QOdRL6lhERBa1//w1jPv2CBwUMuye1RtBXi5SR6K7MOXzu8Y3kY8fP55Tj3RPFHIZPhrVDoVlFdidlI2nVxzFT892RXhD/tcYEdknvUFg/m+Vs0/juoawPNkhrkRuAZyBur3Scj3Gf3cER5Jz4eXqiHXPRaGJr5vUsYiIzG7NkVS8siEBHk5K7Hu5Dzy5R6hN4ErkZJWcHBT4dkIk2jRUI7dIhye+OYy03GKpYxERmVVhWQU+3H4OADCtb1OWJzvFAkW1yt3JASuf6owmvm7I0pbi8W8O46q2VOpYRERm8+WeC8gpLEOojyvGR4VIHYcshAWKap2XqyNWPd0FQV7OSM0txhPfHEZukU7qWERE9y0ttxjfHEgGAMwd1AKOSn7M2iv+L0uS8Fc7YfXTXeHnocL57EJM+O4ICkrLpY5FRHRf3tuWBF2FAVFh3niwlZ/UcciCWKBIMsHeLlj1dBd4uToiIV2Dp1ccQ4lOL3UsIqJ7cjwlF1tPZkImA/4zhItm2jsWKJJUUz93fP9UZ7irlDhyORfPrTqOsgqWKCKyLQaDwNtbzwAARkUEoXUAl2mxdyxQJLnwhmosf7ITnB0U+OPcNbz4YyzK9QapYxER1djG2HTEp+XD1VGBWQOaSR2HagELFFmFyBAvLBsfCUelHNsTr2LWunjoDVyijIisX0FpORb8r3LRzBceaAJfdyeJE1FtYIEiq9G9qQ+WPN4RSrkMm+Mz8OqGBBhYoojIyn2+67xx2YKnu4dKHYdqCQsUWZW+Lf3w6ej2kMuAtcfS8PbWRHCxfCKyVheyC7D8z8sAgNeHtoJK+c97xZJ9sZkCNW/ePERHR8PFxQX16tW77ZjU1FQMHToUrq6u8PHxwbRp06DTVV9fKCEhAb169YKzszMaNmyIt99++5YP6H379iEiIgJOTk4ICwvD0qVLLfVr0W0MaRuA9x9tBwBYcfAy3tt2liWKiKyOEAJvbk5EhUGgX0tf9GnuK3UkqkU2U6B0Oh1GjhyJ559//rbv6/V6DB48GEVFRThw4ADWrFmD9evXY9asWcYxWq0WDz74IAICAnD06FEsWrQIH374IT7++GPjmOTkZDz00EPo0aMHYmNj8eqrr2LatGlYv369xX9H+sujEYF4Z3g4AGDpvov4bNd5iRMREVX3++ksHLiQA0elHP83pJXUcaiWKaUOUFNvvfUWAGDFihW3fX/79u1ITExEWloaAgICAAAfffQRJk6ciHnz5sHDwwOrV69GaWkpVqxYAZVKhfDwcJw7dw4ff/wxZs6cCZlMhqVLlyI4OBiffvopAKBly5Y4duwYPvzwQ4wYMeK2f3ZZWRnKysqMX2u1WvP94nXYuK6NUFaux7u/nsGnO89DpVTg+d6NpY5FRIQSnR7v3Fi24LmeYWjk7SpxIqptNjMDdTcxMTEIDw83licAGDBgAMrKynD8+HHjmF69ekGlUlUbk5GRgcuXLxvH9O/fv9rPHjBgAI4dO4by8tuvlL1gwQKo1WrjKygoyMy/Xd31TI8wvDygOYDKFX6/u7FFAhGRlJbsu4j0/BIEqJ0wpXcTqeOQBOymQGVlZcHPr/qy+Z6ennB0dERWVtY/jqn6+m5jKioqkJOTc9s/e+7cudBoNMZXWlqaWX4nqvRCnyaY9kDlX1Bvb03E6sMpEiciorosOacIS/ddBAC8NrgVnB1543hdJGmBevPNNyGTye74OnbsWI1/3u2WzRdCVDt+85iqm5NNHfN3KpUKHh4e1V5kXjMebIbneoYBAF7beArrjrKkElHtE0Lg9V9OQVdhQI+mPniojb/UkUgikt4DNXXqVIwePfqOY0JCQmr0s/z9/XH48OFqx/Ly8lBeXm6cUfL39zfONFXJzs4GgLuOUSqV8Pb2rlEWMj+ZTIZXBrVAWYUBKw5expwNJ6FUyPBIx0CpoxFRHbL1ZCb2n6+8cfydh8O5310dJmmB8vHxgY+Pj1l+VlRUFObNm4fMzEw0aNAAQOWN5SqVChEREcYxr776KnQ6HRwdHY1jAgICjEUtKioKW7Zsqfazt2/fjsjISDg4OJglK90bmUyGN4a2QoXBgFWHUjH753goFXIMaxdw928mIrpP2tJyvL01EQDwQu8mCPHhjeN1mc3cA5Wamoq4uDikpqZCr9cjLi4OcXFxKCwsBAD0798frVq1wrhx4xAbG4tdu3Zh9uzZmDRpkvGS2tixY6FSqTBx4kScOnUKGzduxPz5841P4AHA5MmTkZKSgpkzZ+LMmTP47rvv8O2332L27NmS/e70F5lMhreHhWN0pyAYBDBjbRx+S8iUOhYR1QEf/X4W1wrKEObjism9w6SOQ1ITNmLChAkCwC2vPXv2GMekpKSIwYMHC2dnZ+Hl5SWmTp0qSktLq/2ckydPih49egiVSiX8/f3Fm2++KQwGQ7Uxe/fuFR06dBCOjo4iJCRELFmyxKSsGo1GABAajeaef1+6M73eIGatixON5mwVjef+KradypQ6EhHZsfi0PBHyylbRaM5WceD8NanjkIWY8vktE4JLPJubVquFWq2GRqPhDeUWpDcIzFoXh01xGVDKZfjy8Y7o35o3dBKReekNAsO/+BMJ6RoMbx+AT0d3kDoSWYgpn982cwmP6GYKuQwfjmyHYe0CUGEQeOHHE9iZeFXqWERkZ1YevIyEdA3cnZR4bTBXHKdKLFBk05QKOT4e1Q5D2jZAuV5gyuoT2J3EEkVE5pGWW4wPfj8LAJgzsAXqu6vu8h1UV7BAkc1TKuT49LH2GNymAXR6Ayb/cAJ7krKljkVENk4IgVc3JqCkXI/OoV4Y2zlY6khkRVigyC4oFXJ8Oro9BoX7Q6c34LkfjrNEEdF92XAi3bjm08JH2kAu55pP9BcWKLIbDgo5Ph/ToVqJ2nWGl/OIyHQ5hWV459fKNZ9e6tsUYfXdJE5E1oYFiuxKVYl6qE1liZq86jhvLCcik725+TTyi8vRqoEHnu3JNZ/oVixQZHccFHJ8NroDBrepvLH8+dXHsYMliohqaGfiVWw9mQm5DHhvRFs4KPhRSbfivxVklypLVHsMNj6ddxy/n866+zcSUZ2mKSnHfzadAgBM6hGGNoFqiRORtWKBIrulVMjx2WPtMbRdAMr1Ai+sPoFfT3LbFyL6Z29vSUSWthShPq6Y3q+Z1HHIirFAkV1TKuT4ZFQ7/KtDQ1QYBF786QR+iUuXOhYRWaEdiVex/sQVyGXAhyPbwtlRIXUksmIsUGT3lAo5PhzZDiMjAo0bEK8/fkXqWERkRfKKdJi7IQFA5aW7iEZeEicia8cCRXWCQi7DeyPaYkznYBgEMPu/8Vh7NFXqWERkJV7ffBo5hWVo4uuGGQ/y0h3dHQsU1RlyuQzz/xWOCVGNIAQwZ30Cvo+5LHUsIpLYbwmZ2BKfAYVcho9GtoOTAy/d0d2xQFGdIpPJ8Oaw1nimeygA4PVfTuPrPy5KnIqIpJJTWGZ86m5yrzC0C6onbSCyGSxQVOfIZDK8NrglXnygCQBg/m9J+GzneQghJE5GRLVJCIFX1icgt0iHFv7umNa3qdSRyIawQFGdJJPJMKt/c7w8oDkA4JOd5/DetrMsUUR1yE9H0rDzzFU4KGT4aFQ7qJS8dEc1xwJFddoLfZrg9SGtAABL913EG5tPw2BgiSKyd5euFeKdrZV73b08oDlaB3DBTDINCxTVeU91D8X8f7WBTAZ8H5OC2f+NR4XeIHUsIrKQcr0B09fGoaRcj+jG3nimO/e6I9OxQBEBGNslGJ8+1h4KuQwbTqRj6o+xKKvQSx2LiCzg053ncPKKBmpnB3w0qh3kcpnUkcgGsUAR3fBw+4ZY8nhHOCrk2HY6C8+sPIZiXYXUsYjIjI4k5+LLvZVP3s7/Vxs0UDtLnIhsFQsU0d/0b+2P7yZ2grODAvvP52D8t0egKSmXOhYRmYGmpBwz1sZBCGBEx0AMbttA6khkw1igiG7SvakPVj3TGe5OShxLycPorw/hWkGZ1LGI6D4IITDnvyeRnl+CYC8XvDmsldSRyMaxQBHdRkQjL6x5tit83FQ4k6nFyKUHkZZbLHUsIrpHPxxKwbbTWXBQyLB4bAe4OzlIHYlsHAsU0T9oHaDGfydHIdDTGZevF+PRpQdx7mqB1LGIyESn0jV4d+sZAMDcQS3RNrCetIHILrBAEd1BiI8r/js5Gs383HBVW4aRS2NwIjVP6lhEVEMFpeWY+uMJ6PQG9Gvphye7hUgdiewECxTRXfirnbDuuSi0D6oHTUk5Hl92GHvOZksdi4juQgiBVzeewuXrxQhQO+HDkW0hk3HJAjIPFiiiGqjn4ojVz3RBz2b1UVKux6SVx7DhxBWpYxHRHaw5moYt8RlQyGVYNLYD6rk4Sh2J7AgLFFENuaqU+GZ8JIa3D0CFQWDmungs++OS1LGI6Dbi0/Lxxi+nAQCz+zdHRCMviRORvWGBIjKBo1KOj0e1x9PdQwEA8347g/m/neH+eURWJLdIhymrK+97erCVH57rya1ayPxYoIhMJJfL8J/BLfHKoBYAgK//uIQZ6+Kgq+D+eURS0xsEXloTi/T8EoR4u3CrFrIYFiiieyCTyTC5V2N8OLIdlHIZfonLwMTlR6At5arlRFL6dOc57D+fA2cHBZaOi4AH13siC2GBIroPj0YE4ruJneDqqMDBi9cxamkMsjSlUsciqpN2Jl7Fot0XAAALR7RBC38PiRORPWOBIrpPPZvVx9rnolDfXYWkrAI88uWfXHCTqJZdvFaIGeviAAATo0PwcPuG0gYiu8cCRWQG4Q3V2PB8NMLquyJDU4oRSw7i4IUcqWMR1Qma4nJMWnkMBaUViGzkiVcfail1JKoDWKCIzCTIywXrJ0ejU4gnCkorMP67I/jvca4VRWRJFXoDpv50ApdyihCgdsKSJyLgqORHG1ke/y0jMiNPV0f88HQXDG1XuVbU7J/j8fH2sxCCyxwQWcL835KMN40vmxCJ+u4qqSNRHcECRWRmTg4KfPZYe7zQpzEA4PPdFzBzXTzKKvQSJyOyL+uOpuG7P5MBAB+PaofWAWqJE1FdwgJFZAFyuQwvD2iB90a0gUIuw8bYdDzxzWHkFumkjkZkF45ezsVrmxIAANP7NcWgNg0kTkR1DQsUkQU91ikYK57sBHcnJY5ezsPwL/7EhWw+oUd0Py7nFOHZ74+hXC8wKNwf0x5oKnUkqoNYoIgsrEfT+tg4JRrBXi5IzS3Gv748iP3nr0kdi8gm5RbpMHH5EeQVl6NtoJorjZNkWKCIakETX3dseqGb8Qm9icuP4oeYy1LHIrIppeV6TPr+GC5fL0bDes74ZkIkXByVUseiOooFiqiWeLk6YtUzXfBIx4bQGwT+75fTeG1jAsr13EOP6G4MBoFZ6+JxPCUP7k5KrHiyE3zdnaSORXUYCxRRLVIpFfhoZDvMGdgCMhmw+nAqby4nqoH3tiXh14RMOChk+GpcBJr6uUsdieo4FiiiWiaTyfB878b4Znwk3FRKHE7OxbDFB3AmUyt1NCKr9N2BZHz1xyUAwHsj2iK6sY/EiYhYoIgk07elHzZOiUYjbxdcySvBiCUHse1UptSxiKzKxtgreHtrIgBgdv9meKRjoMSJiCqxQBFJqKmfO355oRu6NfFGsU6PyatO4IPfk6A3cOVyot1JVzH755MAgCe7heCFPk0kTkT0FxYoIonVc3HEyic7Y1KPUADAF3su4qkVR6EpLpc4GZF0jiTn4vlVJ6A3CPyrQ0P83+BWkMm4XAFZDxYoIiugVMjx2uBW+Gx0ezg5yLHv3DUM++IAkrJ4XxTVPYkZWjy98ijKKgx4oIUv3n+0Ldd6IqvDAkVkRR5u3xAbnu+GQE9npFwvxr++OIhNselSxyKqNReyCzD+u8MoKK1ApxBPfDG2IxwU/Kgi68N/K4msTKsAD2yZ2h09mvqgpFyP6Wvj8H+bTnEzYrJ7l64VYsyyw8gp1KFVAw98M6ETnB0VUsciui0WKCIr5OnqiBVPdsa0vpV7fP1wKAWjvjqE9PwSiZMRWUbK9SKMXXYY1wrK0MLfHaue6QK1s4PUsYj+EQsUkZVSyGWY+WAzLJ/YCWpnB8Sn5WPI5/ux92y21NGIzCottxhjlx1GlrYUTXzdsOqZLvBydZQ6FtEdsUARWbk+LXyx9cXuaNNQjbzickxcfhTvbUtCBbeAITuQnl+Csd9Uzq6G+bjix2e6wMdNJXUsortigSKyAUFeLvh5chTGdW0EAFiy9yJGf30IGbykRzYs5XoRRi2NQVpuCRp5u+DHSV3h68H97cg2sEAR2QgnBwXeGR6OL8Z2hLtKiWMpeRj8+X7sSeIlPbI9F7ILMHJpjHHm6adJXeGvZnki28ECRWRjBrdtgK3TuiO8oQfyisvx5IqjeGvLaT6lRzYjMUOLx746hOyCMjT3c8ea57oioJ6z1LGITGITBery5ct4+umnERoaCmdnZzRu3BhvvPEGdLrqO9inpqZi6NChcHV1hY+PD6ZNm3bLmISEBPTq1QvOzs5o2LAh3n77bQhRfduMffv2ISIiAk5OTggLC8PSpUst/jsSmaKRtyvWPx+NidEhAIDlf17G8C8O4kJ2gbTBiO4iPi0fY5YdwvUiHcIbemDNs13h686ZJ7I9SqkD1ERSUhIMBgO++uorNGnSBKdOncKkSZNQVFSEDz/8EACg1+sxePBg1K9fHwcOHMD169cxYcIECCGwaNEiAIBWq8WDDz6IPn364OjRozh37hwmTpwIV1dXzJo1CwCQnJyMhx56CJMmTcKqVavw559/YsqUKahfvz5GjBgh2TkguplKqcCbw1qjZzMfzP75JM5kajFk0QG8PqQ1xnQO4rYXZHX2n7+GyT8cR5FOj47B9bD8yc5cqoBslkzcPP1iIz744AMsWbIEly5dAgD873//w5AhQ5CWloaAgAAAwJo1azBx4kRkZ2fDw8MDS5Yswdy5c3H16lWoVJVPeSxcuBCLFi3ClStXIJPJMGfOHGzevBlnzpwx/lmTJ09GfHw8YmJibpulrKwMZWVlxq+1Wi2CgoKg0Wjg4eFhqVNAZJStLcWsn+Ox/3wOAODBVn5Y8EgbPs1EVmNTbDpm/xyPCoNAdGNvLBsfCVeVTfw3PNUhWq0WarW6Rp/fNnEJ73Y0Gg28vLyMX8fExCA8PNxYngBgwIABKCsrw/Hjx41jevXqZSxPVWMyMjJw+fJl45j+/ftX+7MGDBiAY8eOobz89pu7LliwAGq12vgKCgoy169JVCO+Hk5Y+WRnvPpQCzgoZNiReBUDP/0DOxKvSh2NCN/sv4Tpa+NQYRAY2i4Ay5/sxPJENs8mC9TFixexaNEiTJ482XgsKysLfn5+1cZ5enrC0dERWVlZ/zim6uu7jamoqEBOTs5t88ydOxcajcb4SktLu79fkOgeyOUyPNuzMTa90A3N/dyRU6jDpO+P4ZX1J1FYViF1PKqDDAaBeb8m4t1fK2f0n+oWis8eaw+VktuzkO2TtEC9+eabkMlkd3wdO3as2vdkZGRg4MCBGDlyJJ555plq793ung8hRLXjN4+puoJp6pi/U6lU8PDwqPYikkrrADV+mdoNz/YMg0wGrDmahkGf/YGYi9eljkZ1SIlOjxd/isWy/ckAgLmDWuD/hrSEXM5788g+SDqHOnXqVIwePfqOY0JCQoz/nJGRgT59+iAqKgpff/11tXH+/v44fPhwtWN5eXkoLy83zij5+/sbZ5qqZGdXrqFztzFKpRLe3t41/+WIJOTkoMCrD7XEAy18MWtdPNJySzBm2SFMiGqEfw9swcsnZFFXtaWY9P0xnLyigYNChoWPtMWIiECpYxGZlaR/i/r4+MDHx6dGY9PT09GnTx9ERERg+fLlkMurT55FRUVh3rx5yMzMRIMGDQAA27dvh0qlQkREhHHMq6++Cp1OB0dHR+OYgIAAY1GLiorCli1bqv3s7du3IzIyEg4OfFqEbEvXMG9sm94D839Lwk9HUrEyJgW7z2bj/RHtENWY/0FA5ncqXYNnVh5DlrYUni4OWPpEBLqE8d81sj828RReRkYGevXqheDgYHz//fdQKP66fu7v7w+gchmD9u3bw8/PDx988AFyc3MxceJEDB8+3LiMgUajQfPmzfHAAw/g1Vdfxfnz5zFx4kS8/vrr1ZYxCA8Px3PPPYdJkyYhJiYGkydPxk8//VTjZQxMuYufqLYcOJ+DOetPIv3G9i+PdwnGnEEt4OHE/zAg89h2KhMz1sajpFyPJr5u+HZCJBp5u0odi6jGTPn8tokCtWLFCjz55JO3fe/v8VNTUzFlyhTs3r0bzs7OGDt2LD788MNqT90lJCTghRdewJEjR+Dp6YnJkyfj9ddfr3Z/0759+zBjxgycPn0aAQEBmDNnTrUb1u+GBYqsVWFZBeb/dgY/Hk4FAPh5qPDWsHAMDPeXOBnZMr1B4OMdZ/HFnosAgB5NffDF4x1Zzsnm2F2BsjUsUGTtDl7MwWsbTyE5pwgAMKC1H94aFs69yMhkuUU6TPspFgcuVD6l/GS3ELz2UEsoFTb5kDfVcSxQEmOBIltQWq7Hot3n8dW+S6gwCLirlJjxYDOMj2rEDz+qkbi0fExZdRwZmlI4OyiwcEQbPNy+odSxiO4ZC5TEWKDIliRlafHK+gTEpeUDAFo28MC7w8MR0chT2mBktYQQ+OFQCt7degY6vQFhPq5Y8kQEmvu7Sx2N6L6wQEmMBYpsjd4gsOZoKt7fdhaaksoV90dFBmLOwBbw5nYw9De5RTr8+7/x2HmmcgmYAa398OHIdnDn/U5kB1igJMYCRbbqemEZ3tuWhHXHrgAAPJyUeKlf5WU9B17Wq/P+vJCDGWvjkF1QBkeFHHMfaoGJ0SHcuJrsBguUxFigyNYdT8nF/206jcRMLQAgrL4r/m9wK/Rp4StxMpJCWYUeH+84h6//uAQhgMb1XbFoTEe0CuDfb2RfWKAkxgJF9kBvEFh3LA0f/n4W14t0AIBezerj1Yda8l6XOiQ+LR+zf47H+exCAMCYzsF4fUgrODtyPzuyPyxQEmOBInuiLS3H4t0XsPzPZJTrBeQyYETHQMzs3wwN1M5SxyMLKS3X47Nd5/HVvoswCMDHzRHvDm/DNcPIrrFASYwFiuxRck4R3t+WhP+dqtwrUqWU48luoXi+d2OonXkDsT05npKHV9afNM46Pdw+AG8MbQ0vV0eJkxFZFguUxFigyJ6dSM3Dwt+ScORyLoDKG82f7RmGid1C4cZNim1aXpEO721LwpqjaQAAHzcV5v0rHANac9aJ6gYWKImxQJG9E0Jg15lsvLctyThL4enigOd6Ncb4qEZwcWSRsiWGG/e7LdyWhPziymUsRkYE4tWHWsKTs05Uh7BASYwFiuoKvUFg68kMfLrzvHFbGB83RzzTIwyPdwnm2kA24HhKHt79NRGxqfkAgBb+7nh3eDgiQ7ykDUYkARYoibFAUV1ToTdgU1wGPt91Hqm5xQAqL+1NjA7BxG6hvHfGCl3OKcL7vyfht4TKe9pcHRWY8WAzTIwO4VY+VGexQEmMBYrqqnK9AZti07Fk30VculY5I+XsoMDozkF4qlsogrxcJE5IuUU6LNp9HqsOpaBcLyCTVV6um9W/Ofw8uJk01W0sUBJjgaK6Tm8Q2H46C1/svYBT6ZWLccplQP9W/ni6RygiG3ly9epallukw7L9l/D9wcso0ukBVK7rNfehFmjhz7+niAAWKMmxQBFVEkLgj/M5+Gb/Jew/n2M83i5QjfFRIRjctgGcHLggoyVdLyzDsv3J+D7mMopvFKfWAR54ZVAL9GhaX+J0RNaFBUpiLFBEtzp3tQDfHUjGhth06CoMAAC1swMejQjE2C7BaFzfTeKE9uVyThGW/5mMn49fMRan8IYeeKlvM/Rr6csZQKLbYIGSGAsU0T+7XliGNUfT8OPhVKTnlxiPdwn1wqMRgRjUpgHXk7pHQggcTs7FN/uTsSvpKqr+dm/TUI3p/ZrigRYsTkR3wgIlMRYoorvTGwT+OHcNqw+nYHdSNgw3/iZydlBgYLg/RnQMRFRjbyjk/MC/G01xOTbFpeOnI6lIyiowHn+ghS+e6R6KqMbeLE5ENcACJTEWKCLTpOeXYOOJK1h/It24nhQA1HdXYWBrfzzUpgE6h3qxTP2NwVA527T2aCp+O5VlvCzq5CDHiI6BeLJbKJr48rIokSlYoCTGAkV0b4QQiE3Lx4YTV7AlPhOaknLjez5ujhjQ2h99W/oiurFPnbz5XAiBE6n5+PVkJn5LyESWttT4Xgt/d4zpHIyH2wegngvX3SK6FyxQEmOBIrp/ugoDDl7MwW8Jmfj99NVqZcrJQY7oxj7o08IXvZrWR5CXs91eoiqr0ONIci72JF3D76ezqt035q5SYmj7AIzuFIQ2DdV2ew6IagsLlMRYoIjMq1xvwMGL17H9dBb2JGUjQ1Na7f2G9ZzRNcwbUY0rXwFqJ8uUCYMeSDkIFF4F3PyARtGA3LwzYUIIXMopQszF69h79hoOXswxPkUHVK4Y3q+VH4a0DUDPZj5QKeveTByRpbBASYwFishyhBA4e7UAu85kY+/ZbMSm5qPCUP2vMV93FdoH1UP74HpoH1QP4Q3V8LjfffkSNwPb5gDajL+OeQQAA98DWg275x9botMjKUuL4yl5OHY5D8dScpFTqKs2pr67Cr2b1Ufflr7o3dy3Tl6+JKoNLFASY4Eiqj1FZRU4lpKHmIvXEXPpOk6la6A33PrXWoDaCc383dHc3x3NfN0R4uOCIC8X1HdT3X22KnEzsG48gJt/7o3vG/X9XUtUabkeabnFuHy9GJeuFSIxU4vTGVpculaIm+M6KuVoH1QPvZrVR69m9dGqgQfkvIGeyOJYoCTGAkUknWJdBU6laxGXloe4tHzEpebfcsnv75wc5AjydEGDes7wcXNEfTcVvN0c4e2qgpuTEq5KGTr/0gsOxZm4XYURkEHn4o9Dw/aiuFxAW1qOnEIdrhWU4VphGa4VlOFKbjEytaX4p79tvV0d0TZQjU6hXugc4oU2gWpemiOSAAuUxFigiKyLprgc57ILkJRVgHNZBTifXYC03BJkakpumf25WVd5ItY4vnvXP2O07j84ZGh1xzHuKiUa+bggxNsVLRt4oFWAB1o38EB99xrMghGRxZny+c3lfonI7qldHNApxAudQryqHddVGJCRX4LU3GJkaUtxvVCHnMIy5BSW4XqhDkW6CrQvLAWK7/5ndPAsg97dE24qJXzcVKjvroKPmwo+7ioEejojxNsVni4OLEpEdoIFiojqLEelHCE+rgjxcf3nQckGYOWHd/1Zcx7tBYRGmzEdEVkzudQBiIisWqPoyqftbnsHFCqPezSsHEdEdQYLFBHRncgVlUsVALi1RN34euBCs68HRUTWjQWKiOhuWg2rXKrAo0H14x4BNVrCgIjsD++BIiKqiVbDgBaDLb4SORHZBhYoIqKakiuA0B5SpyAiK8BLeEREREQmYoEiIiIiMhELFBEREZGJWKCIiIiITMQCRURERGQiFigiIiIiE7FAEREREZmIBYqIiIjIRCxQRERERCbiSuQWIIQAAGi1WomTEBERUU1VfW5XfY7fCQuUBRQUFAAAgoKCJE5CREREpiooKIBarb7jGJmoSc0ikxgMBmRkZMDd3R0ymcysP1ur1SIoKAhpaWnw8PAw68+mv/A81w6e59rB81x7eK5rh6XOsxACBQUFCAgIgFx+57ucOANlAXK5HIGBgRb9Mzw8PPh/zlrA81w7eJ5rB89z7eG5rh2WOM93m3mqwpvIiYiIiEzEAkVERERkIhYoG6NSqfDGG29ApVJJHcWu8TzXDp7n2sHzXHt4rmuHNZxn3kROREREZCLOQBERERGZiAWKiIiIyEQsUEREREQmYoEiIiIiMhELlA358ssvERoaCicnJ0RERGD//v1SR7JpCxYsQKdOneDu7g5fX18MHz4cZ8+erTZGCIE333wTAQEBcHZ2Ru/evXH69GmJEtuHBQsWQCaTYfr06cZjPM/mk56ejieeeALe3t5wcXFB+/btcfz4ceP7PNf3r6KiAv/5z38QGhoKZ2dnhIWF4e2334bBYDCO4Xk23R9//IGhQ4ciICAAMpkMmzZtqvZ+Tc5pWVkZXnzxRfj4+MDV1RXDhg3DlStXLBNYkE1Ys2aNcHBwEMuWLROJiYnipZdeEq6uriIlJUXqaDZrwIABYvny5eLUqVMiLi5ODB48WAQHB4vCwkLjmIULFwp3d3exfv16kZCQIB577DHRoEEDodVqJUxuu44cOSJCQkJE27ZtxUsvvWQ8zvNsHrm5uaJRo0Zi4sSJ4vDhwyI5OVns3LlTXLhwwTiG5/r+vfvuu8Lb21ts3bpVJCcni59//lm4ubmJTz/91DiG59l0v/32m3jttdfE+vXrBQCxcePGau/X5JxOnjxZNGzYUOzYsUOcOHFC9OnTR7Rr105UVFSYPS8LlI3o3LmzmDx5crVjLVq0EK+88opEiexPdna2ACD27dsnhBDCYDAIf39/sXDhQuOY0tJSoVarxdKlS6WKabMKCgpE06ZNxY4dO0SvXr2MBYrn2XzmzJkjunfv/o/v81ybx+DBg8VTTz1V7dgjjzwinnjiCSEEz7M53FyganJO8/PzhYODg1izZo1xTHp6upDL5WLbtm1mz8hLeDZAp9Ph+PHj6N+/f7Xj/fv3x8GDByVKZX80Gg0AwMvLCwCQnJyMrKysauddpVKhV69ePO/34IUXXsDgwYPRr1+/asd5ns1n8+bNiIyMxMiRI+Hr64sOHTpg2bJlxvd5rs2je/fu2LVrF86dOwcAiI+Px4EDB/DQQw8B4Hm2hJqc0+PHj6O8vLzamICAAISHh1vkvHMzYRuQk5MDvV4PPz+/asf9/PyQlZUlUSr7IoTAzJkz0b17d4SHhwOA8dze7rynpKTUekZbtmbNGpw4cQJHjx695T2eZ/O5dOkSlixZgpkzZ+LVV1/FkSNHMG3aNKhUKowfP57n2kzmzJkDjUaDFi1aQKFQQK/XY968eRgzZgwA/jttCTU5p1lZWXB0dISnp+ctYyzxWckCZUNkMlm1r4UQtxyjezN16lScPHkSBw4cuOU9nvf7k5aWhpdeegnbt2+Hk5PTP47jeb5/BoMBkZGRmD9/PgCgQ4cOOH36NJYsWYLx48cbx/Fc35+1a9di1apV+PHHH9G6dWvExcVh+vTpCAgIwIQJE4zjeJ7N717OqaXOOy/h2QAfHx8oFIpbGnR2dvYtbZxM9+KLL2Lz5s3Ys2cPAgMDjcf9/f0BgOf9Ph0/fhzZ2dmIiIiAUqmEUqnEvn378Pnnn0OpVBrPJc/z/WvQoAFatWpV7VjLli2RmpoKgP9Om8vLL7+MV155BaNHj0abNm0wbtw4zJgxAwsWLADA82wJNTmn/v7+0Ol0yMvL+8cx5sQCZQMcHR0RERGBHTt2VDu+Y8cOREdHS5TK9gkhMHXqVGzYsAG7d+9GaGhotfdDQ0Ph7+9f7bzrdDrs27eP590Effv2RUJCAuLi4oyvyMhIPP7444iLi0NYWBjPs5l069btlqU4zp07h0aNGgHgv9PmUlxcDLm8+senQqEwLmPA82x+NTmnERERcHBwqDYmMzMTp06dssx5N/tt6WQRVcsYfPvttyIxMVFMnz5duLq6isuXL0sdzWY9//zzQq1Wi71794rMzEzjq7i42Dhm4cKFQq1Wiw0bNoiEhAQxZswYPopsBn9/Ck8InmdzOXLkiFAqlWLevHni/PnzYvXq1cLFxUWsWrXKOIbn+v5NmDBBNGzY0LiMwYYNG4SPj4/497//bRzD82y6goICERsbK2JjYwUA8fHHH4vY2Fjjcj01OaeTJ08WgYGBYufOneLEiRPigQce4DIGJMQXX3whGjVqJBwdHUXHjh2Nj9vTvQFw29fy5cuNYwwGg3jjjTeEv7+/UKlUomfPniIhIUG60Hbi5gLF82w+W7ZsEeHh4UKlUokWLVqIr7/+utr7PNf3T6vVipdeekkEBwcLJycnERYWJl577TVRVlZmHMPzbLo9e/bc9u/kCRMmCCFqdk5LSkrE1KlThZeXl3B2dhZDhgwRqampFskrE0II889rEREREdkv3gNFREREZCIWKCIiIiITsUARERERmYgFioiIiMhELFBEREREJmKBIiIiIjIRCxQRERGRiVigiIiIiEzEAkVERERkIhYoIiIiIhOxQBERERGZiAWKiKgGLly4AJlMhl9//RV9+/aFi4sLmjdvjsOHD0sdjYgkwAJFRFQD8fHxkMlk+Oijj/Cf//wH8fHxCA4OxiuvvCJ1NCKSAAsUEVENxMfHQ61WY+3atejTpw+aNm2K4cOH49q1a1JHIyIJsEAREdVAfHw8hg4divr16xuPXbp0CU2aNJEwFRFJhQWKiKgG4uPjERUVVe1YbGws2rdvL00gIpIUCxQR0V1oNBqkpKSgQ4cO1Y7HxcWxQBHVUSxQRER3ER8fD4VCgXbt2hmPpaSkIC8vjwWKqI5igSIiuov4+Hi0aNECzs7OxmOxsbGoV68eQkJCpAtGRJKRCSGE1CGIiIiIbAlnoIiIiIhMxAJFREREZCIWKCIiIiITsUARERERmYgFioiIiMhELFBEREREJmKBIiIiIjIRCxQRERGRiVigiIiIiEzEAkVERERkIhYoIiIiIhP9P68VRxOobWnXAAAAAElFTkSuQmCC",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"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": [
"<Figure size 640x480 with 1 Axes>"
]
},
"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": [
"<Figure size 640x480 with 1 Axes>"
]
},
"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
}