{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"# Minimization problem, application to a uniformly-loaded simple beam"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"In this example, we will look at a minimization problem linked to the deflection of a uniformly-loaded simple beam. The beam is made of a cylinder of length $L$ with an inner diameter $d$ and thickness $t$.\n",
"\n",
"
\n",
"\n",
"\n",
"The deflection of this beam is expressed as:\n",
"\n",
"$$ u = \\frac{5qL^4}{284EI} $$\n",
"\n",
"where $q$ is the load applied to the beam, E its Young's modulus and I is the area moment of inertia.\n",
"\n",
"In this minimization problem, we want to minimize the deflection $u$ with the lighest beam as possible using an unconstrained and unbounded optimization.\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:49:51.485Z",
"iopub.status.busy": "2021-04-05T16:49:51.478Z",
"iopub.status.idle": "2021-04-05T16:49:51.885Z",
"shell.execute_reply": "2021-04-05T16:49:51.890Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [],
"source": [
"# We import here the python modules we need for this example\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy import optimize"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"We define three python functions to compute:\n",
"\n",
"1. the are moment of inertia $I$, as function of the diameter $d$ and thickness $t$ of the cylinder (function I(t,d))\n",
"\n",
"2. the mid-point deflection $u$ as function of the applied load $q$, the length of the beam $L$, its Young's modulus $E$ and the area moment of inertia $I$ (function mid_point(q,L,E,I)).\n",
"\n",
"3. the mass of the beam given its diameter $d$, thickness $t$ and density $\\rho$ (function mass(rho,L,d,t))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:49:54.126Z",
"iopub.status.busy": "2021-04-05T16:49:54.121Z",
"iopub.status.idle": "2021-04-05T16:49:54.135Z",
"shell.execute_reply": "2021-04-05T16:49:54.139Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [],
"source": [
"def I(t,d):\n",
" return np.pi*((d+2.0*t)**4.0-d**4.0)/64"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:49:55.098Z",
"iopub.status.busy": "2021-04-05T16:49:55.093Z",
"iopub.status.idle": "2021-04-05T16:49:55.108Z",
"shell.execute_reply": "2021-04-05T16:49:55.112Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [],
"source": [
"def mid_point(q,L,E,I):\n",
" return (5*q*L**4.0)/(384.0*E*I)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:49:56.085Z",
"iopub.status.busy": "2021-04-05T16:49:56.078Z",
"iopub.status.idle": "2021-04-05T16:49:56.095Z",
"shell.execute_reply": "2021-04-05T16:49:56.101Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [],
"source": [
"def mass(rho,L,d,t):\n",
" V_tot = L*rho*np.pi*(d/2.0+t)**2.0\n",
" V_in = L*rho*np.pi*(d/2.0)**2.0\n",
" return (V_tot-V_in)*1000.0"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"We now give values to the parameters of our problem."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:49:57.822Z",
"iopub.status.busy": "2021-04-05T16:49:57.815Z",
"iopub.status.idle": "2021-04-05T16:49:57.830Z",
"shell.execute_reply": "2021-04-05T16:49:57.835Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [],
"source": [
"q = 100.0 #N/mm\n",
"L = 1000.0 #mm\n",
"E = 70000.0 #MPa\n",
"d = 50.0 #mm\n",
"rho = 2.7e-9"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"We can now plot the deflection of the beam $u$ (red curve) and the mass of the beam (blue curve). While the mass of the beam evolves with the square of the thickness $t$, in the range of interest it results in an almost linear response. The deflection of the beam $u$ is linked to the thickness $t$ through the area moment of inertia $I$ and shows a non-linear response."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:49:59.972Z",
"iopub.status.busy": "2021-04-05T16:49:59.966Z",
"iopub.status.idle": "2021-04-05T16:50:00.144Z",
"shell.execute_reply": "2021-04-05T16:50:00.158Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": true
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8vElEQVR4nO3dd3gVZfbA8e8JoQVUBERBUFBYEBEQUHFBBBHXgmXXhmvBirqsYlkV/dnXgl2xF1Tsuq4KAiKIxI6KioAggigCgvSS0Mn5/XEmm0sSyCXJ3LnlfJ7nfe7cuXNnDkPg5H3nLaKqOOecc+kmK+oAnHPOuTB4gnPOOZeWPME555xLS57gnHPOpSVPcM4559KSJzjnnHNpKTvqACoiKytLa9asGXUYSaOgoICsLP+dJZbfky35/SjJ70lJFb0na9asUVWN/KamdIKrWbMm+fn5UYeRNHJzc+nevXvUYSQVvydb8vtRkt+Tkip6T0RkbeVFU36RZ1jnnHMuDJ7gnHPOpSVPcM4559KSJzjnnHNpyROcc865tOQJzjnnXFryBOeccy4teYJzzjmXljI3wU2cCFdcAStXRh2Jc84lhaVLoV8/yMtL6TlA/idzE9xPP8EDD8Dvv0cdiXPORW7KFDjgABg6FKZN2yHqcCpFeqTp8mjQwF4XL4Z99ok2Fueci9B//wt9+8KOO8JHH8G6dctDv6aI/AqsBjYDm1S1U2VfI3NrcIUJbtGiaONwzrmIFBTAjTfCSSfBfvvZk5vOnRMaQg9VbR9GcoNMrsHtsou9Ll4cbRzOOReBlSvhjDNgxAg491x47DGoXj3qqCpX5tbg6te3V6/BOecyzIwZcNBBMHo0PPIIPPNMpSe3bBGZGFP6lXKMAmNE5JutfF7xIMI4aUqoWhXq1vUE55zLKCNGwOmnW0L74AM49NBQLhPPM7WuqjpfRBoAY0XkR1X9uDKDyNwaHFgzpTdROucyQEEB3HYbHHccNG9uz9tCSm5xUdX5wesi4G3gwMq+RmYnuAYNvAbnnEt7q1fDySfDDTdY7e3TT2GPPaKLR0RqicgOhdvAEcDUyr5OaAlORFqKyKSYskpELhORuiIyVkRmBq87B8eLiAwWkVkiMllEOoQV2/94gnPOpblZs6xn5DvvwH33wQsvQM2aUUfFrsCnIvI98BUwUlVHV/ZFQnsGp6ozgPYAIlIFmI9VQwcC41R1kIgMDN5fAxwFtAjKQcDjwWt4vInSOZfG3nsP/v53qFIFxoyBnj2jjsio6mygXdjXSVQTZU/gZ1WdAxwPDA32DwVOCLaPB15QMwGoIyINQ42qQQObm2bTplAv45xziaQKd94JxxwDTZva87ZkSW6JlKgE1wd4NdjeVVUXBNsLsaoqwO7A3JjvzAv2hadBA/tJWLo01Ms451yi5OXZ87brroM+feCzzyzJZaLQhwmISDXgOODa4p+pqoqIbuf5+gH9ALKzs8nNzS13bLssWsS+wNejRpHfrFm5z5Ms8vLyKnQ/0pHfky35/Sgpne7J/Pk1uf76Nvz2Ww4XX/wzJ588j6++2v7zpM09UdVQC9b0OCbm/QygYbDdEJgRbD8JnFbacVsrOTk5WiHjx6uC6rhxFTtPkhg/fnzUISQdvydb8vtRUrrck1GjVOvUUa1bV3Xs2Iqdq6L3BMjXkHNLPCURTZSnUdQ8CTAc6Bts9wWGxew/K+hN2RlYqUVNmeHw6bqccylOFW6/3Z637bmnPW87/PCoo0oOoTZRBuMbegEXxuweBLwhIucBc4BTgv2jgKOBWcAa4JwwYwN8wmXnXEpbvdpWAXj7bRvf9tRTkJMTdVTJI9QEp6r5QL1i+5ZivSqLH6tA/zDjKaFuXcjK8gTnnEs5M2bAX/9qS1vefz9cdhmIRB1VcsncuSjBBofUq+dNlM65lDJ8OJx5JlSrZvNJdu8edUTJKbOn6gKfzcQ5lzIKCuCmm+D446FFC/jmG09u25LZNTjwBOecSwkrVtj6bSNH2nO3xx9Piim3kprX4Hy6Ludckps6FQ44AN5/Hx59FJ57zpNbPLwG5zU451wSe+MNW3F7hx0gNxe6dIk6otThNbgGDazuv2FD1JE459z/bNoE//oXnHoqtGtnz9s8uW0fT3CFg72XLIk2DuecCyxaBEccYcvb9O8P48dDo0ZRR5V6PMH5YG/nXBL56ivo2BE+/9yetT3yiA0HcNvPE1xhgvOOJs65iA0ZAoccYkN0P/8czj476ohSmye4wiZKr8E55yKyfj306wfnnw+HHmrP2zp0iDqq1OcJzpsonXMRmjvXam1PPw3XXmurcNerV/b3XNl8mECdOpCd7U2UzrmE+/BD6yW5fj289ZbNLekqj9fgRKyZ0mtwzrkEUYW774Zevey/n6+/9uQWBq/BgQ/2ds4lzKpVcM45VmM7+WR49lmoXTvqqNKT1+DAEpw3UTrnQjZ9Ohx0EAwbBvfeC6+/7sktTF6DA2sj+PnnqKNwzqWx//zHptzKyfElbhLFa3DgNTjnXGgKp9w65RTYbz/49ltPboniCQ6sBrd6NaxdG3Ukzrk0snAhHH540ZRbubmw++5RR5U5PMGBz2binKt0n31mg7W/+gpefNGn3IqCJzjwBOecqzSqMHiwNUPm5MCECbZQqUs8T3Dg03U55ypFXh6cfjoMGABHHQUTJ0LbtlFHlbk8wYFP1+Wcq7AZM6BzZ+v6f8cd8M47NlGSi44PEwDYbTd7XbAg2jiccynpv/+1wdvVq8OYMdCzZ9QROQi5BicidUTkTRH5UUSmi8jBIlJXRMaKyMzgdefgWBGRwSIyS0Qmi0ji5tKuVQvq1rVZT51zLk4bN9oQgJNOgtatbQiAJ7fkEXYT5UPAaFVtBbQDpgMDgXGq2gIYF7wHOApoEZR+wOMhx7alPfaA335L6CWdc6lrwYIthwB8/DE0aRJ1VC5WaAlORHYCugFDAFR1g6quAI4HhgaHDQVOCLaPB15QMwGoIyINw4qvBE9wzrk4ffwx7L+/dSLxIQDJK8waXDNgMfCciHwnIs+ISC1gV1UtfNi1ENg12N4diG0jnBfsSwxPcM65MqjaHJKHHQY77QRffulDAJJZmJ1MsoEOwCWq+qWIPERRcyQAqqoiottzUhHphzVhkp2dTW5ubqUE22TjRvZesYJPRo5kc61alXLORMvLy6u0+5Eu/J5sye9HSfHek7y8Ktx9dys++WQXunVbzNVX/8iSJZtJx9uZNj8nqhpKAXYDfo15fwgwEpgBNAz2NQRmBNtPAqfFHP+/47ZWcnJytNK89poqqE6dWnnnTLDx48dHHULS8XuyJb8fJcVzT77/XrV5c9UqVVTvu0+1oCD8uKJU0Z8TIF9Dyi3bU0JrolTVhcBcEWkZ7OoJTAOGA32DfX2BYcH2cOCsoDdlZ2ClFjVlhm+PPezVmymdczFeeMHGt+Xn21ySV1xh6yS75Bf2OLhLgJdFpBowGzgHe+73hoicB8wBTgmOHQUcDcwC1gTHJo4nOOdcjHXr4NJL4emnoUcPePVV2HXXsr/nkkeoCU5VJwGdSvmoxEiRoFrbP8x4tmm33SA72xOcc45ffrGxbd9+C9deC7feav89uMojIlWAicB8Ve0dxjX8r6xQlSq2joUnOOcy2rvvwllnFW33DuW/XgcMwMZG7xjWBXwuylg+VMC5jLVpEwwcCMcdB3vtZbU3T27hEJHGwDHAM2FexxNcLE9wzmWkwllJ7roLLrrI1nJr1izqqFJatohMjCn9in3+IHA1UBBqEGGePOXssYdNBb55szVZOufS3qRJdejTB1avth6TZ54ZdURpYZOqltb/AhHpDSxS1W9EpHuYQXgNLtYee1g7xcKFUUfinAtZQQHceSdceWU76tSxlbc9uSVEF+A4EfkVeA04TEReCuNCnuBi+VAB5zLCsmX2rO2666B790V8/TXsu2/UUWUGVb1WVRuralOgD/ChqoYy4Zk3UcaKTXAHHxxtLM65UHz1FZxyCvz+u02S3Lr1dHbYwQe4pSOvwcXyGpxzaUvVElrXrvb+s89smRuflSQ6qpob1hg48AS3pR13tCnCPcE5l1ZWrYJTT4VLLoG//MWGABxwQNRRubB5givOhwo4l1YmT4ZOneCtt2wYwLBhULdu1FG5RPBncMV5gnMuLajCc89ZM+TOO8OHH0K3blFH5RLJa3DFeYJzLuXl58PZZ8N550GXLvDdd57cMpEnuOL22MP6EOflRR2Jc64cpk+HAw+EF1+Em2+G99/3VQAylSe44gp7Us6dG20czrnt9uKL9rxtyRIYMwZuusknJcpknuCK86ECzqWctWvhggtsFYADDrAmycMPjzoqFzVPcMV5gnMupcyYAQcdBM88YzOTfPABNGoUdVQuGXgvyuIaNYKsLE9wzqWAV1+Ffv2genUYNQqOOirqiFwy8RpccdnZvvCpc0lu3Tpb1ubvf4d27WDSJE9uriRPcKXZc09bs945l3RmzoTOneHJJ+Gqq2D8eGjcOOqoXDLyBFeali2tYd85l1Refx06dLBOziNGwN13Q9WqUUflkpUnuNK0bAmLFsHy5VFH4pzDmiQvvhj69IG2ba1J8phjoo7KJTtPcKVp1cpevRbnXOR++smaJJ94Aq6+GnJzoUmTqKNyqcATXGkKE9yPP0Ybh3MZ7rXXoGNHmDcPRo60yZK9SdLFK9QEJyK/isgUEZkkIhODfXVFZKyIzAxedw72i4gMFpFZIjJZRDqEGds2NWtm/4o8wTkXibVrrfv/aadZL8nvvoOjj446KpdqElGD66Gq7VW1U/B+IDBOVVsA44L3AEcBLYLSD3g8AbGVLjsbmjf3JkrnIvDjjzZw++mnYeBA6yXpTZKuPKJoojweGBpsDwVOiNn/gpoJQB0RaRhBfKZVK6/BOZdghXNJLlgA770Hd97pTZKu/MKeyUSBMSKiwJOq+hSwq6ouCD5fCBTO8707EDvD8bxg34KYfYhIP6yGR3Z2Nrm5uaEE3qxmTZrMmsUnH3yAZqfGhC95eXmh3Y9U5fdkS8l6P9auzWLw4BaMHt2Qtm1XcP3106hRYwOJCDVZ70mU0uaeqGpoBdg9eG0AfA90A1YUO2Z58DoC6BqzfxzQaVvnz8nJ0dA8/7wqqM6YEd41Ktn48eOjDiHp+D3ZUjLejylTVPfZR1VE9YYbVDduTOz1k/GeRK2i9wTI1xBzS7wl1CZKVZ0fvC4C3gYOBP4obHoMXhcFh88HYlvaGwf7ouE9KZ0LlapNkHzAAbYE45gxcOut9gjcucoQWoITkVoiskPhNnAEMBUYDvQNDusLDAu2hwNnBb0pOwMrtagpM/FatrRX72jiXKVbvRpOP92WuOnaFb7/3pe3cZUvzN+VdgXeFpHC67yiqqNF5GvgDRE5D5gDnBIcPwo4GpgFrAHOCTG2stWpY8sAew3OuUr17bdw6qkwezbcfrv1lMzyEbkuBKElOFWdDbQrZf9SoGcp+xXoH1Y85eJzUjpXaVTh4YdtguQGDeCjj6z25lxY/PembfGhAs5VimXL4K9/hQED4IgjbC5JT24ubJ7gtqVVK1i6FJYsiToS51LWp59C+/a2IOn998Pw4VCvXtRRuUzgCW5bvKOJc+W2ebM9Y+veHapVg88/h8svB3ss71z4PMFtiw8VcK5cFiyAv/wFrr8eTj7ZOpZ06lT295yrTD7iZFv23BOqV/canHPbYfRoOOssyMuDIUPgnHO81uai4TW4balSBVq08Bqcc3HYsMHWazvqKBthM3EinHuuJzcXHa/BlaVVKxuF6pzbqtmzbWmbr76CCy+EBx6AmjWjjsplOq/BlaVVK/vXu3591JE4l5Refx32399a8t9801be9uTmkoEnuLK0a2fdwaZMiToS55JKfj6cfz706QOtW9vYthNPjDoq54p4gitLx472+s030cbhXBKZPNl6RT77LFx3HXz8MTRtGnVUzm3JE1xZmjaFnXe2fs7OZThVePRROPBAWLkSxo61sW6+KKlLRt7JpCwi0KGD1+Bcxlu6FM47D4YNg6OPhuefh112iToql85EpAHQBWgErMVWpJmoqgXxfN9rcPHo2NGewW3YEHUkzkXio4/scfSoUdZDcsQIT24uPCLSQ0TeB0YCRwENgdbA9cAUEblFRHYs6zxeg4tHhw6W3H74wbqLOZchNm2yRUhvuw2aN4cJE+yfg3MVISI1gI+B6lgeelNVb4o55GjgAlX9rZTvZgO9gV7Af7d1HU9w8YjtaOIJzmWIOXPg73+3OSTPOQcGD4bataOOyqWJ9cBhqponIlWBT0XkPVWdAKCqV23ti6q6CXgnnot4govH3nvDTjt5RxOXMd54A/r1g4ICeOUVG8TtXGUJ1v/MC95WDYoWP05Erijl6yuBb1R1UlnX8Wdw8RCxmpt3NHFprnBs26mnwj772CQ+ntxcGESkiohMAhYBY1X1y1IO6wRcBOwelAuBI4GnReTqsq5RZg0uaCvtDRzClj1ZRqrqD/H9UdJAx47wyCOwcaP3iXZp6bvvLJn99JONbbv5Zv9Rd+WWLSITY94/papPxR6gqpuB9iJSB3hbRNqo6tRi52kMdFDVPAARuQnreNIN+Aa4e5tBbOtDEbkFS265wJdYpq0B/AkYFCS/K1V18rbOkxY6dLDpuqZPh7Zto47GuUpTUAAPPQQDB0L9+vDBB3DYYVFH5VLcJlWNa4EkVV0hIuOxmlnxBNcAe15XaCOwq6quFZEy508sqwb3VbGeLbHuD8Yo7FHWRdJCbEcTT3AuTfzxB5x9ti1xc9xxtrxN/fpRR+XSnYjsAmwMkltNrEfkXaUc+jLwpYgMC94fC7wiIrWAaWVdZ5vP4FR1ZBmfL1LVids6Jm20aGFdyPw5nEsT771nv6vl5trsJO+848nNJUxDYLyITAa+xp7BjSh+kKr+G+gHrAjKRap6q6rmq+rpZV0krl6UItIJ+D9gz+A7YtfWzKnKZGVZRxPvSelS3Pr1cM011iy5337w4Yew775RR+UySfBYq8wxVyJynqoOASbG7BukqgPjuU68wwReBq4CpgBxTZGSljp2hCeftNGv2T7CwqWeX3/N4bLLrHfkJZfA3XdDjRpRR+XcVp0oIutU9WUAEXkU6wcSl3j/l16sqsPLE52IVMGy73xV7S0izYDXgHpYL5gzVXWDiFQHXgA6AkuBU1X11/JcMzQdO8LatbbCd5s2UUfjXNxUbZ22yy7ryI47wrvvQu/eUUflXJlOBIaLSAHWCWWFqp4X75fjHQd3k4g8IyKnicjfCkuc3x0ATI95fxfwgKo2B5YDhcGeBywP9j9A6Q8co9Up6BT0ZWnDNZxLTkuWwF//Cv/4B7Rrt5IpUzy5ueQmInVFpC5QEzgfuBpYDdwS7I9LvAnuHKA9lkGPDUqZ/0REpDFwDPBM8F6Aw4A3g0OGAicE28cH7wk+7xkcnzxatoQGDWzmWedSwAcfWEeS996D+++HQYMms9tuUUflXJm+wVr+vgHGA3WwXFK4Py7xNlEeoKottzNAgAexzLtD8L4eVsXcFLyfh41OJ3idCzbXmIisDI5fEntCEemH9aohOzub3NzccoRVfq332Ycdx4xhQoKvG4+8vLyE349kl6n3ZMMGYciQvXjjjSbsuWc+jz46nebN8zL2fmyL35OSor4nqtqssk5UZgGeA1rHc2zMd3oDjwXb3YERQH1gVswxTYCpwfZUoHHMZz8D9bd1jZycHE24wYNVQfWXXxJ/7TKMHz8+6hCSTibek2nTVNu3tx/Tiy9Wzc8v+iwT70dZ/J6UVNF7AuTrduSL4gXoWsbnOwJtyjpPvDW4zsAkEfkFG1UezzCBLsBxInI01utlR+AhoI6IZKvV4hoD84Pj5wcJb16wHMJOWGeT5HLoofb60Ue22rdzSULVOvlecQXUqgXDh8Oxx0YdlXPlcqKI3A2MxpolF2N5pDnQAxuydmVZJ4k3wR25vdGp6rXAtQAi0h34l6qeLiL/AU7CelL2BQpHqA8P3n8RfP5h8JtAcmnTBurWtQTXt2/U0TgHwOLFttr2u+/CEUfYatsNG0YdlXPlo6qXB51JTgROxgaGr8U6LD6pqp/Gc564EpyqzhGRnbEaVux35mxX1OYa4DURuQ34DhgS7B8CvCgis4BlQJ9ynDt8WVlwyCHe0cQljdGjbb22Zctste1LL7UfU+dSmaouA54OSrnEO5PJv4GzsedihbUqxXpElklVc7EJm1HV2cCBpRyzDsvUye/QQ2HYMJg3Dxo3jjoal6HWrbMJkh96yGYief99nybVuVjxNlGeAuytqhvCDCZlxD6HO73M6dCcq3RTpthq21On2owkd90FNWtGHZVzySXehoyp2DgEB9Cuna3w7c2ULsEKCuDBB+GAA+y523vvweDBntycK028Nbg7ge9EZCoxa/Oo6nGhRJXsqlSBrl09wbmE+v13W9pm7FjrHTlkCOyyS9RRORceETkZGK2qq0XkeqADcJuqxjXrfbwJbig2dVZmT7Yc69BDYeRIWLDAu6u50L31FlxwgU2F+sQT0K8fJNk8P86F4QZV/Y+IdAUOB+4BHgcOiufL8TZRrlHVwao6XlU/KizlDDg9FD6H++STaONwaW31auv+f+KJ0KwZfPcdXHihJzeXMTYHr8cAT6mtUVot3i/Hm+A+EZE7ReRgEelQWLY30rTSoYMtgPrhh1FH4tLUhAm2BOHzz8N118Hnn9t0qM5lkPki8iRwKjAqWHUm7kEw8TZRFi5M1zlmX9zDBNJSdjb07AmjRtkUEv4rtaskGzfC7bfDbbfZKJTcXBt66VwGOgWbaOReVV0hIg2xtUnjEu9A7x7lDC69HXusjYebPNl6VjpXQTNnwhlnwFdfwZlnwsMPW4dd5zJUQ2Ckqq4PZsRqi60bGpdtVvVE5AwR2eoxIrJ38PAvMx1zjL2++260cbiUpwpPPQXt21uSe/11eOEFT24u4/0X2CwizYGnsNm0Xon3y2XV4OphwwO+oeSEl4diS9kMLEfQ6WG33WxA0ogRcP31UUfjUtSiRdaRZMQIa/V+/nmfIMe5QIHa8ml/Ax5W1YdF5Lt4v7zNGpyqPoSNO3gV2AXoGbyfD5ypqieq6szyx54Gjj3W2pP++CPqSFwKevdd2G8/G9v24IMwZownN+dibBSR04CzsCXXAKrG++Uye6Oo6mZVHauqN6vqhap6mao+qaq/lTPg9HLssda+NGpU1JG4FJKXZ2PZjjsOGjWCb76BAQN8kmTnijkHOBi4XVV/EZFmwIvxftn/OVVUu3b2K7c/h3Nx+uILe9b2zDNwzTU2HGDffaOOyrnko6rTVPVSVX01eP+Lqt4V7/c9wVWUCPTubW1L69ZFHY1LYhs3wg032CxvmzZZ9/9Bg6B69agjcy45iUgLEXlTRKaJyOzCEu/3PcFVhmOPhfx8n5vSbdWPP8LBB9vYtjPOgO+/h27doo7KuaT3HDY11yZsJe8XgJfi/XK868HdWNp+Vb013gultR49bDr3d9+Fv/wl6mhcEikogEcfhauvhlq14M03bdot51xcaqrqOBERVZ0D3Bz06i81JxUXbw0uP6ZsBo4CmpYj2PRUsyb06gXDh9v/aM4B8+fDkUfaCts9etgabp7cnNsu64Ox2DNF5J8i8legdrxfjivBqep9MeV2oDuwV7nCTVennAJz58Jnn0UdiUsCr70GbdrYj8MTT9jCE77ohHPbbQCQA1wKdATOBPrG++V456IsLgfw0TqxTjjB2qBeesknDsxgy5bBP/5hM5F07myzkbRoEXVUzqUmVf062MzDhgxsl3ifwU3BJlcGqIIN+vbnb7Fq1bIk98YbtsSyd43LOO+/D+eeazOT3HabDQHILu+vkM5lMBEZvq3P411sO95/fr1jtjcBf6jqpji/mznOOANefhnee8+SncsI+flw1VXw+OPQurX1NeqQ2YtJOVdRBwNzsVm0vgTKtVxLvM/g5sSU+Z7ctuLww6FBA2umdBmhcND2E0/AFVfYjCSe3JyrsN2A64A2wENAL2DJ9i627ePgKlN2NvTpY7/Cr1gRdTQuRBs2wP/9nw3a3rjR1r297z6oUSPqyJxLfcEUkaNVtS+2DuksIFdE/rk95wktwYlIDRH5SkS+F5EfROSWYH8zEflSRGaJyOsiUi3YXz14Pyv4vGlYsYXqjDPsf78334w6EheSKVPgwAPhjjugb19bDrB796ijci69BDnhb9jA7v7AYODt7TlHmDW49cBhqtoOaA8cKSKdgbuAB1S1ObAcOC84/jxgebD/geC41NOpE/zpT/YszqWVzZvh7rvtr3jBAlvr9tlnYccdo47MufQiIi8AX2Cr19yiqgeo6r9Vdf72nCe0BKcmL3hbNSgKHAYUVm+GAicE28cH7wk+7yki5XqwGCkRq8Xl5sKcOVFH4yrJrFlw6KHWM7J3b5g61VYCcM6F4gygBTYO7nMRWRWU1SKyKt6TiKqWfVQ5iUgVbKHU5sCjwD3AhKCWhog0Ad5T1TYiMhU4UlXnBZ/9DBykqkuKnbMf0A8gOzu749ixY0OLv7yqL1xI59NP57c+ffjlggsSdt28vDxq1457kH9GqOg9UYV3323E44/vTZUqyoABMzn88D9IwV+9AP8ZKY3fk5Iqek969OixRlVrVWJI5aOqoRegDjAe6ArMitnfBJgabE8FGsd89jNQf1vnzcnJ0aR1wgmq9eurrl2bsEuOHz8+YddKFRW5J3Pnqh5xhCrY69y5lRdXVPxnpCS/JyVV9J4A+ZqA3FJWSUgvSlVdESS4g4E6IlI4/q4xtjo4wWsTgODznYCliYgvFP/8JyxZYgO/XUpRtZEebdrAp5/a+LbRo32lbecqi4g0EZHxwTI4P4jIgDCuE2Yvyl1EpE6wXRMbxzAdS3QnBYf1BYYF28MpmmPsJODD4DeB1HTYYdCqFTzySNSRuO2waJFNiHzmmZbgJk+Giy4iZZsknUtSm4ArVbU1Ngygv4i0ruyLhFmDawiMF5HJwNfAWFUdAVwDXCEis4B6wJDg+CFAvWD/FcDAEGMLn4jV4r7+Gr76KupoXBzeestW1h450npLfvQR7L131FE5l35UdYGqfhtsr8YqP7tX9nVCmylPVScD+5eyfzZwYCn71wEnhxVPJM46C6691mpxL7wQdTRuK5YtsyVtXn4ZOnaEoUMt0TnnwheMed4fm5KrUvlMJmHaYQc4+2ybWn7RoqijcaUYNcqaIl9/HW6+2abe8uTmXIVli8jEmNKvtINEpDbwX+AyVY27+3+8PMGFrX9/m9nk8cejjsTFWLkSzjsPjjkG6tWDL7+Em26CqlWjjsy5tLBJVTvFlKeKHyAiVbHk9rKqvhVGEJ7gwtaypY0IfughWFXpv6C4cvjgA9hvP3j+eRg4ECZO9AmSnUukYBKPIcB0Vb0/rOt4gkuEG26A5cvh0UejjiSj5eVZhbpXL8jJsdW277zTl+5zLgJdsNW5DxORSUE5urIv4ssxJkKnTnD00Tbd/CWXgM+akHAffQTnnAO//gqXXw633w41a0YdlXOZSVU/pZxrvG0Pr8Elyg03wNKl/iwuwdaty+Kyy2y2/6wsS3T33+/JzblM4AkuUTp3hiOOgHvugTVroo4mI3z2GZx/ficeesiGJH7/PRxySNRROecSxRNcIt14IyxeDE8+GXUkaW3tWrjySktmmzcLH34IDz8MtaKf+tU5l0Ce4BKpSxebwuvOO62fuqt0X3wB++9vzZAXXQTPPjuRHj2ijso5FwVPcIl29902CfMdd0QdSVpZuxauugq6drXtsWPhscegZs3NUYfmnIuIJ7hE69jRpvB68EGYPTvqaNLChAlWa7v3Xjj/fJgyBQ4/POqonHNR8wQXhTvugOxsWx7alVthra1LF+u38/779nhzxx2jjsw5lww8wUWhUSNLbm++CZ98EnU0Kal4rW3qVOuk6pxzhTzBReVf/7IVNC+/HAoKoo4mZaxda7fuz3+2WtuYMV5rc86VzhNcVHJy4K674JtvrDeEK9Pnn0P79jYhTL9+Vmvr1SvqqJxzycoTXJROOw3+8hdbM27OnKijSVpr1sAVV1gPyfXrrYfkE094rc05t22e4KIkYu1rqjZoSzXqiJLOJ59Au3bwwANw8cXeQ9I5Fz9PcFHbc08b+D16NLz0UtTRJI28PJuXuls32LwZPvzQFmPYYYeoI3POpQpPcMngH/+wXhOXXQZ//BF1NJH78ENo29YS2qWXwuTJ+Gwkzrnt5gkuGVSpAs88A/n5tqZLhvaqXLXKWmp79rRhgh9/bOvE+upCzrny8ASXLPbZx7oHvveezXKSYUaPhjZt4OmnrUPJpEnWqcQ558rLE1wy+cc/4IQTYOBAmDgx6mgSYvlyq7QedZTV1D77zPJ8Tk7UkTnnUp0nuGQiAkOGwG67QZ8+1maXxoYNg9at4cUX4brr4LvvbNk855yrDKElOBFpIiLjRWSaiPwgIgOC/XVFZKyIzAxedw72i4gMFpFZIjJZRDqEFVtSq1sXXnkFfvkFLrggLYcOLF5s+fuEE2DXXeHrr+H226F69agjc86lkzBrcJuAK1W1NdAZ6C8irYGBwDhVbQGMC94DHAW0CEo/4PEQY0tuXbva0IE33kirZXVULXfvsw+8/Tb8+9+W3PbfP+rInHPpKLQEp6oLVPXbYHs1MB3YHTgeGBocNhQ4Idg+HnhBzQSgjog0DCu+pHfVVXD66XD99daWl+Lmz4fjj7c/UvPm1hx5/fVQtWrUkTnn0pVoAprARKQp8DHQBvhNVesE+wVYrqp1RGQEMEhVPw0+Gwdco6oTi52rH1bDIzs7u+PYsWNDjz8qWevX0/6yy6j16698++ij5O+11zaPz8vLo3aS9alXhZEjG/LEE3uzaZNw7rm/cOKJ86hSJTHXT8Z7EiW/HyX5PSmpovekR48ea1S1ViWGVD6qGmoBagPfAH8L3q8o9vny4HUE0DVm/zig07bOnZOTo2lv/nzVhg1V99zTtrdh/PjxCQkpXrNmqfbooQr2OmtW4mNItnsSNb8fJfk9Kami9wTI15BzSzwl1F6UIlIV+C/wsqq+Fez+o7DpMXhdFOyfDzSJ+XrjYF9ma9QIhg+HpUttYubly6OOqEybN8P998N++9liCU8+CePGwd57Rx2Zcy6ThNmLUoAhwHRVvT/mo+FA32C7LzAsZv9ZQW/KzsBKVV0QVnwppVMneOcd+OknOPZYm14/SU2darOOXXmlzUjyww+2tI1I1JE55zJNmDW4LsCZwGEiMikoRwODgF4iMhM4PHgPMAqYDcwCngb+EWJsqadnT+uC+MUXcPLJsHFj1BFtYf16uOkm6NABZs+GV1+1imfjxlFH5pzLVNlhnVits8jWfm/vWcrxCvQPK560cOKJthBav36W5F5/PSkGj02YAOedB9OmwRln2NI29etHHZVzLtP5TCap5oIL4JFHbOjAiSfCunWRhZKXBwMGWJPk6tUwapTNSuLJzTmXDDzBpaL+/a3nxsiRNrhs7dqEh/D++zY58sMP2xSaP/xg80k651yy8ASXqvr1g2efhbFjrXflsmUJueySJXDWWXDkkVCzpq24/cgjvhCpcy75eIJLZeecY705vvwSunaleoiLpcZOs/XqqzYLyXffQZcuoV3SOecqxBNcqjv1VGsv/P13OvTvD99/X+mXmDMHjjnGptnaay/49lubR7JGjUq/lHPOVRpPcOmge3f49FM0K8uqVG+9VeZX4rF5MwweDPvua6trP/AAfP65DeB2zrlk5wkuXbRpw7ePPWY9P048EW68EQoKyn26qVNtUYMBA+CQQ+z9ZZeRsDkknXOuojzBpZEN9etDbq49m/v3v23Bte2c2mvdOnu+tv/+MGsWvPSSdf9v2jSMiJ1zLjye4NJNjRq2KvjDD8N771mm+uKLuL760UfQrp0tPvr3v8P06fbczafZcs6lIk9w6UgE/vlP+OwzyMqyNsZBg7baZLl8uY0f797dZgAbMwaGDvUB28651OYJLp0deKD15f/b3+Daa+Gww2yiyIAq/Oc/1vX/uedsjdWpU6FXrwhjds65SuIJLt3ttJPNWTlkiPXvb9sWnniC3+Yoxx0Hp5xiEyJ//TXcfTfk5EQdsHMu3YnIsyKySESmhnkdT3CZQATOPRemTmVz5y48ePGPtN57HR+OK+D++22y5P33jzpI51wGeR44MuyLeILLIN8t3YPOK0dzOQ/STT5h2qaWXL7qFrI3r486NOdcBlHVj4HQ5xf0BJcB8vPt+doBB8DcucJrr8HIOfux5986ws0328jtkSPtoZxzzqUJT3BpbtQom4nk3nutlXL6dJvdSxo1hNdeg9GjrQmzd29bDmDatKhDds6lvmwRmRhT+kURhCe4NLVwoSWyY46xjiMffwxPPQU771zswL/8BaZM4X8P49q2hQsvhPnzI4nbOZcWNqlqp5jyVBRBeIJLMwUFtuh3q1bwzjtw6602UuCQQ7bxpWrV4PLLbeqSiy+2MQPNm8PVV8PSpYkK3TnnKpUnuDTyyy+16NrVclSHDjB5MtxwA1SvHucJ6te3GVBmzICTT7Z2zaZN4brrbCE455yrBCLyKvAF0FJE5onIeWFcxxNcGlizxsZxX3BBR376yWYhGTcOWrYs5wmbNYMXXrAMecwxNgtK06bWU8WbLp1zFaSqp6lqQ1WtqqqNVXVIGNfxBJfi3n/fFhAYNAh69fqDH3+0FbcrZf7INm2sI8oPP8Dxx9tzumbN4OyzbcoT55xLYp7gUtSCBdCnDxx5pD1CGz8errlmRjjzR+6zD7z8sj2ju+gim99rv/3giCPg3Xdt4TjnnEsynuBSzObN8NhjRZ1IbrnFFvHu3j0BF2/WzFZA/e03W3Jg2jQ47jj405/gnntg0aIEBOGcc/EJLcGVNteYiNQVkbEiMjN43TnYLyIyWERmichkEekQVlypbNIk+POfoX9/G7Q9ZYqtaxp3J5LKUq+edTz55Rd44w2bzPLqq+311FPhgw8qtNiqc85VhjBrcM9Tcq6xgcA4VW0BjAveAxwFtAhKP+DxEONKOXl5cOWV0KkT/PqrLUI6diy0aBFxYFWrWm/Ljz6y2lz//hZYr17WKeX662HmzIiDdM5lqtAS3FbmGjseGBpsDwVOiNn/gpoJQB0RaRhWbKlCFd5+2x6B3X8/nH8+/Phjki5Cus8+8MAD8Pvv1jFl333hzjut+fKgg+Chh2z0uXPOJYhoiPMPikhTYISqtgner1DVOsG2AMtVtY6IjAAGqeqnwWfjgGtUdWIp5+yH1fLIzs7uOHbs2NDij9LChTUYPLg5X3xRn732yuOKK35i331XbfM7eXl51K5dO0ERlq3a4sXsOm4cDcaNY4dZs9CsLFa0bcviQw9lSbdubKhbN/QYku2eRM3vR0l+T0qq6D3p0aPHGlWtVYkhlY+qhlaApsDUmPcrin2+PHgdAXSN2T8O6FTW+XNycjTdrF+vOmiQas2aqrVqqd57r+rGjfF9d/z48aHGViHTpqnecINqq1aqoCqiesghqvfdpzprVmiXTep7EgG/HyX5PSmpovcEyNcQc0u8JdG9KP8obHoMXgu73c0HmsQc1zjYl1E++cTWZRs40KaInDbNnr1lZ0cdWSXYZx+bN2zaNBtDd+ONsHKl/QGbN7cxdwMHwqefwqZNUUfrnEsDiU5ww4G+wXZfYFjM/rOC3pSdgZWquiDBsUVm8WI45xzo1s2Wthk+3J697bFH1JGFQMSez918s41vmD3bHjA2aAD33WeTZjZoYL0xn3vOZ05xzpVbaHWDYK6x7kB9EZkH3AQMAt4I5h2bA5wSHD4KOBqYBawBzgkrrmRSUABDhsA118Dq1VaBuf56qBV9y3XiNGtmEz1ffrnV6MaMsbXpxoyxIQgArVtDz55w+OFw6KGw007RxuycSwmhJThVPW0rH/Us5VgF+ocVSzKaNMkmRZ4wwf7Pfuwx+388o+20kw07OPlk60I6ZYrNRTZuHDzzjE0EnZVl7bjdu1vp0qWUNYCccy7EBOdKt2qVPX56+GEbLz10KJx5ZhJ2+4+aiK1N17atTfK8fj188QXk5lp5+GFr0gR7fnfIIZbs/vxnG4PnN9S5jOcJLkFUrcXt8sttONiFF8Idd3jlI27VqxfV2gDWroUvv7ROKZ9+aqPfHw/mB9htN+jcGQ46iDrVqkHHjrDDDlFF7pyLiCe4BPjpJ5vk44MPbJ22d96BAw+MOqoUV7Pmlglv82brnfn551a+/BLeeYf2AP/6l03e2amTJbuOHaFdO096zqU5T3AhWrPGaml33w05OfDIIzYZf5UqUUeWhqpUsaTVrp093ARYtozJzzxD23XrYOJEm0bsxRftMxGbZaV9+6LSrp3V/rx507m04AkuJMOHw6WXwpw59oztnntg112jjirD1K3LsgMP3HKphd9/h2+/LSoTJsDrrxd9Xq+eLQW03372bK9NGxvW4D03nUs5nuAq2ezZMGAAjBhhvSJzc62XpEsSjRpZ6d27aN+KFbZ6+fffW8/NyZPh2WdtUGKh3Xe3weqtW9trq1ZWdt3Va3zOJSlPcJVk3Tqrpd1xh7WW3XOPJbqqVaOOzJWpTh0bZd+tW9G+ggJb927qVCvTp9ssLEOGbJn4dtrJmjr/9Cdb3iG21KmT6D+Jcy6GJ7hK8N57cMkl8PPPNgHHvffa0mguhWVl2XCDpk23rO0VFNjsKjNm2NIOP/5ovYg+/RReecW6yxaqVw/23ttK8+Y2qH2vvaw0auQPY50LmSe4CpgzBy67zHpFtmxpfRgOPzzqqFyosrKgSRMrxf+y1661NuqZM638/LOVwud8sYvAVq1qc7EVJtE99yx63WMPaxL16r9zFeIJrhzWrbNa2h132OOXQYNsfFu1alFH5iJVs6Z1SNl335KfbdxoTZ6zZ1v59Vcrv/xiD2z/+GPL47OyrJbXpIk1BxQm1caNrey+u/X49CTo3FZ5gttOo0dbc+SsWXDSSTZPcJMmZX/PZbiqVYuaK0uzdq0lwDlzYO7cLbcnTYJ337XfrGKJWCeX3Xe3ZLj77tCwoW03bFhUdtklTZakcG77+E99nH791Zojhw2z/gRjxkCvXlFH5dJGzZrWzt2yZemfq8KyZTBvnpX587csc+bYVGZLlpT8blaWJbnddqNttWrWC3S33WzVhl13LSq77AL163ut0KUNT3BliO0dmZUFd95pzZHVq0cdmcsoItZppV49G5C+NevXW3PnggU25m/hQtsOXrNnzbKxKwsXwoYNpZ+jbl1LdsVL/fpFr/XrWyz169vyFz5UwiUhT3DbMGKEdfWfPRtOOcWeu3lzpEtq1atbJ5WtLCb4bW4u3bt3txrhqlWWDP/4AxYtsvLHH7ZAYWGZMQM++8xqhrGdZGJVq1aUfAtL3bpFr7Fl552LSu3anhhdqDzBlWLWLGuOHDnSWnM++MCWI3MubYjYGL7CcXxlKSiA5cth6VJLdosXF20vXbpl+fFHa05dutQ612xNdraNFaxTpyjpFW4X7t9pp6LX4qV2bWtWcW4rPMHFyM+3Jsh77rFfSu+5x6bb8t6RLuNlZRXVzuJJiGC1xPx8S4zLlhWV5cuLyooVW27PmVO0vbUm1EIisOOOVgqTXuH72LLDDkWvhSXmfdaGDRar1ybTjic47Gf7rbfs2drcuXD66TZBcqNGUUfmXAoTsVpW7drla9tft84S3YoVttp74faqVfa+sMS+X7zYxh6uXAmrV9uM52XoBlabLEx+tWsXvdaubc8YC7eL76tVq+R2bPEOO5HK+AQ3bZrV0saNs7U1X37Z1s50zkWsRg3r7bnbbuU/x6ZNkJdnyW7VKiurVxe9X72a2ZMns1eDBrYvL6/o+Lw8a4ItfJ+fv+U0bfGoWtWWEqlVa8vX4tullZo1S76WVnJyPJFuRcYmuJUr4ZZbbGHo2rVtKZsLL/ThQs6lldjnfFvxW24ue8WuOLEtBQVWK8zPL0qGhYmvcF9+ftExsdux+9auteRZuG/tWivFxzrGq0qVLZNejRpbfy2+XUqplpNTvjiSTMb+d/7KK/Dgg3D++XD77db72Tnntikrq6iZMoz1rzZvtiS3Zo2VtWtLvhbfji3r1pW+b+nSou3CY9avt9fNm0uEkXPvvZX/Z4tAxia4Cy6Agw+2dS6dcy4pVKlS9PwuUTZtKkp4QQJcNXNm4q4fooxNcNnZntycc26LDjaBgnnzIgyo8vggEuecc2kpqRKciBwpIjNEZJaIDIw6Huecc6kraRKciFQBHgWOAloDp4lI62ijcs45l6qSJsEBBwKzVHW2qm4AXgOOjzgm55xzKUpUNeoYABCRk4AjVfX84P2ZwEGq+s9ix/UD+gFkZ2d3HDt2bMJjTVZ5eXnUrl076jCSit+TLfn9KMnvSUkVvSc9evRYo6pb7QoqIkcCDwFVgGdUdVC5L7YNKdeLUlWfAp4CqFWrlnaPd4BmBsgtnCne/Y/fky35/SjJ70lJYd6TmMdRvYB5wNciMlxVp1X2tZKpiXI+EDthXeNgn3POufSRsMdRyZTgvgZaiEgzEakG9AGGRxyTc8657ZctIhNjSr+Yz3YH5sa8nxfsq/wgwjhpeajqJhH5J/A+1i77rKr+EHFYzjnntt8mVe0UdRBJk+AAVHUUMCre49esWaMisjbEkFJNNrAp6iCSjN+TLfn9KMnvSUkVvSc1t/FZwh5HJU0vSldxIjIxGX5rSiZ+T7bk96MkvyclhXlPRCQb+AnoiSW2r4G/h9Fil1Q1OOecc+ktkY+jPME555xLqO19HFVeydSL0lXcU1EHkIT8nmzJ70dJfk9KSot74s/gnHPOpSWvwTnnnEtLnuBSnIg0EZHxIjJNRH4QkQFRx5QsRKSKiHwnIiOijiUZiEgdEXlTRH4UkekicnDUMUVJRC4P/s1MFZFXRaRG1DElmog8KyKLRGRqzL66IjJWRGYGrztHGWNFeIJLfZuAK1W1NdAZ6O/LDP3PAGB61EEkkYeA0araCmhHBt8bEdkduBTopKptsN58faKNKhLPA0cW2zcQGKeqLYBxwfuU5AkuxanqAlX9Nthejf2nFcq0N6lERBoDxwDPRB1LMhCRnYBuwBAAVd2gqisiDSp62UDNYFxWDvB7xPEknKp+DCwrtvt4YGiwPRQ4IZExVSZPcGlERJoC+wNfRhxKMngQuBooiDiOZNEMWAw8FzTbPiMiW13OJN2p6nzgXuA3YAGwUlXHRBtV0thVVRcE2wuBXaMMpiI8waUJEakN/Be4TFVXRR1PlESkN7BIVb+JOpYkkg10AB5X1f2BfFK46amigudKx2OJvxFQS0TOiDaq5KPWzT5lu9p7gksDIlIVS24vq+pbUceTBLoAx4nIr9hSHIeJyEvRhhS5ecA8VS2s3b+JJbxMdTjwi6ouVtWNwFvAnyOOKVn8ISINAYLXRRHHU26e4FKciAj2XGW6qt4fdTzJQFWvVdXGqtoU6zjwoapm9G/nqroQmCsiLYNdPYFKX2AyhfwGdBaRnODfUE8yuNNNMcOBvsF2X2BYhLFUiCe41NcFOBOrpUwKytFRB+WS0iXAyyIyGWgP3BFtONEJarJvAt8CU7D/C9Ni9o7tISKvAl8ALUVknoicBwwCeonITKymOyjKGCvCZzJxzjmXlrwG55xzLi15gnPOOZeWPME555xLS57gnHPOpSVPcM4559KSJzjnnHNpyROcyxgiUi9mrOBCEZkfbK8QkVIHPYvIrSJy+DbOebaIPBJe1GUTkRNE5MYQzrufiDxf2ed1LlGyow7AuURR1aXYAGdE5GYgT1XvDSapLnXNOFWt9MQRgquB4yr7pKo6RUQai8geqvpbZZ/fubB5Dc45U0VEng4WwBwjIjUBROR5ETkp2D5ARD4Xke9F5CsR2SH2BCJyjIh8ISL1g+8NDo6fXXiO4LirRORrEZksIrcE+2qJyMjg3FNF5NRg/6BgMdvJInJv8aBF5E/AelVdEhPv4yIyIbhu92BRy+mxtTERyRORe4I/7wcicqCI5AbfiU2W75KZ66S5NOAJzjnTAnhUVfcFVgAnxn4oItWA14EBqtoOm8Jobcznf8Vm5z+6MNkADYGuQG+C6Y5E5IjgWgditcmOItINW3Tyd1VtFyzAOVpE6gF/BfZV1bbAbaXE3QWbbirWzsDBwOXYvIIPAPsC+4lI++CYWtgcnfsCq4Nz9wqud2vMuSYCh2ztpjmXzDzBOWd+UdVJwfY3QNNin7cEFqjq1wCqukpVNwWfHQZcAxyjqstjvvOOqhao6jSK1tQ6IijfYYmpFZbwpmDz/90lIoeo6kpgJbAOGCIifwPWlBJ3Q2ydt1jvBsucTAH+UNUpqloA/BDz59oAjA62pwAfBbPqTyn2Z1+ELSfjXMrxBOecWR+zvZntez79M7AD8KdtnFNiXu9U1fZBaa6qQ1T1J2z5minAbSJyY5BAD8QmBe5NUUKKtRaosZXrFhSLoSDmz7VRiyai/d9xQSKM/bPXIKam6lwq8QTnXHxmAA1F5AAAEdlBRAoTwRysSfMFEdm3jPO8D5wbLFCLiOwuIg1EpBGwRlVfAu4BOgTH7KSqo7DmxnalnG860Lyif7ht+BMwNcTzOxca70XpXBxUdUPQ8ePhoAPKWuw5XOHnP4rI6cB/ROTYbZxnjIjsA3xhy5CRB5yBJal7RKQA2AhcjNUKh4lIDazmd0Upp/wYuE9EJKZGVpl6ACNDOK9zofPlcpxLcSLyEPbc7YNKPm914COga8zzRudShjdROpf67gByQjjvHsBAT24uVXkNzjnnXFryGpxzzrm05AnOOedcWvIE55xzLi15gnPOOZeWPME555xLS/8P0TIN1z7icGoAAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"t = np.linspace(0.5,10.0,100) # We span a range of thicknesses\n",
"#\n",
"fig, ax = plt.subplots(constrained_layout=True)\n",
"ax2=ax.twinx()\n",
"ax.plot(t,mid_point(q,L,E,I(t,d)),'r',label='Deflection')\n",
"ax.set_ylabel('u (mm)')\n",
"ax2.plot(t,mass(rho,L,d,t),'b',label='Mass')\n",
"ax2.set_ylabel('Mass (kg)')\n",
"ax.set_xlabel('Thickness (mm)')\n",
"ax.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"Next, we define a cost function $\\Phi$ for the minimization problem. We will start with the simplest one:\n",
"$$\\Phi = u+M$$\n",
"M being the mass of the beam, minimizing this function $\\Phi$ should give us back, _on principle_, the lowest deflection and mass.\n",
"\n",
"The minimization will be carried out using the Nelder-Mead algorithm in the first part of this notebook."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:50:02.440Z",
"iopub.status.busy": "2021-04-05T16:50:02.435Z",
"iopub.status.idle": "2021-04-05T16:50:02.449Z",
"shell.execute_reply": "2021-04-05T16:50:02.454Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [],
"source": [
"def cost_function(x,data):\n",
" [q,L,E,d,rho] = data\n",
" u = mid_point(q,L,E,I(x,d))\n",
" m = mass(rho,L,d,x)\n",
" return u+m"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T18:30:43.243Z",
"iopub.status.busy": "2021-04-05T18:30:43.235Z",
"iopub.status.idle": "2021-04-05T18:30:43.255Z",
"shell.execute_reply": "2021-04-05T18:30:43.267Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [],
"source": [
"# We define here some variables we need for the minimization problem\n",
"options = {'disp': True, 'maxiter':1000} # display the results\n",
"constraints = [] # We do not have constraints thus an empty list\n",
"bounds = [] # We do not have bounds thus an empty list\n",
"data = [q,L,E,d,rho]\n",
"arguments = (data, ) # This list of parameters is needed by our different functions"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T18:30:44.907Z",
"iopub.status.busy": "2021-04-05T18:30:44.899Z",
"iopub.status.idle": "2021-04-05T18:30:44.935Z",
"shell.execute_reply": "2021-04-05T18:30:44.947Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 18.170286\n",
" Iterations: 23\n",
" Function evaluations: 46\n",
"The optimum thickness of the beam is 18.73593750000003\n",
"The resulting deflection is 7.2465060773519765 mm\n",
"The resulting mass of the beam is 10.923779528176095 kg\n"
]
}
],
"source": [
"#We choose a thickness of 3 mm as a starting point\n",
"x0 = 3.0\n",
"res = optimize.minimize(cost_function, x0, method='Nelder-Mead',options=options,args=arguments)\n",
"x_opt = res.x[0]\n",
"print('The optimum thickness of the beam is {}'.format(x_opt))\n",
"print('The resulting deflection is {} mm'.format(mid_point(q,L,E,I(x_opt,d))))\n",
"print('The resulting mass of the beam is {} kg'.format(mass(rho,L,d,x_opt)))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"The Nelder-Mead algorithm gives us a thickness of 18.73 mm as an optimum. While this thickness definitely gives a small deflection of the beam, the resulting mass is rather large (i.e. more than 10 kg).\n",
"The reason for this result is linked to the cost function definition which simply translate the sum of the deflection and the mass. For small values of the thickness, the mass is one or two orders of magnitude smaller than the deflection and thus gives a small contribution to the cost function.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"To remedy this we can try to define a new cost function where the deflection and the mass are normalized. We choose a deflection of 700 mm and a mass of 5 kg to normalize our quantities."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:50:11.111Z",
"iopub.status.busy": "2021-04-05T16:50:11.105Z",
"iopub.status.idle": "2021-04-05T16:50:11.120Z",
"shell.execute_reply": "2021-04-05T16:50:11.125Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [],
"source": [
"def cost_function_normalized(x,data):\n",
" [q,L,E,d,rho] = data\n",
" max_u = 700.0\n",
" max_M = 5.0\n",
" u = mid_point(q,L,E,I(x,d))\n",
" m = mass(rho,L,d,x)\n",
" return u/max_u+m/max_M"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:50:13.164Z",
"iopub.status.busy": "2021-04-05T16:50:13.158Z",
"iopub.status.idle": "2021-04-05T16:50:13.178Z",
"shell.execute_reply": "2021-04-05T16:50:13.187Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.408901\n",
" Iterations: 15\n",
" Function evaluations: 30\n",
"The optimum thickness of the beam is 18.73593750000003\n",
"The resulting deflection is 7.2465060773519765 mm\n",
"The resulting mass of the beam is 10.923779528176095 kg\n"
]
}
],
"source": [
"#We choose a thickness of 3 mm as a starting point\n",
"x0 = 3.0\n",
"res = optimize.minimize(cost_function_normalized, x0, method='Nelder-Mead',options=options,args=arguments)\n",
"x_opt_normalized = res.x[0]\n",
"print('The optimum thickness of the beam is {}'.format(x_opt))\n",
"print('The resulting deflection is {} mm'.format(mid_point(q,L,E,I(x_opt,d))))\n",
"print('The resulting mass of the beam is {} kg'.format(mass(rho,L,d,x_opt)))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"We now get a thickness of around 2.4 mm with a mass of about 1.06 kg. To understand even further this result we can plot the cost function we just used as function of the thickness of the beam. The following graphs shows the cost functions and their associated optimums. When normalizing the deflection and the mass we obtain a cost function with a clear minimum compared to the original sum. It is worth mentionning that the cost function does not reach a zero value in this particular problem."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:50:16.822Z",
"iopub.status.busy": "2021-04-05T16:50:16.817Z",
"iopub.status.idle": "2021-04-05T16:50:17.186Z",
"shell.execute_reply": "2021-04-05T16:50:17.191Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": true
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABJVklEQVR4nO2dd3gU5fbHPwdICIRQVQgiTYNyRQVB0J+ocO0FUcCKCtjbvRaUoqiIqIgVxYYFVBS9gAiKV+UKEUFRUZEiIE0wSJFOaIHk/P6Y2biEJLub7O5sOZ/nmWdn3nln5ju7e+bMW4+oKoZhGIYRa1TwWoBhGIZhFIc5KMMwDCMmMQdlGIZhxCTmoAzDMIyYxByUYRiGEZOYgzIMwzBiEnNQhmEYRkxiDioERCRHRDpG+Zr1RSQnmtcsDyJSS0RURFYWSW8gIjtFZFOR9ML7C/VYwzASG3NQARCR6iLynIisBQ4FJovI5yJyfJQknAd8FqVrlYqIDBSRgQGytQRWA7VEJMMv/TEgB/ilSH7/+wv1WMMICS9eMsuKvbCZgyoVERHgE6At0AHn4XkFMAuYISLHhOEag0XkOb/tBiKyQ0R8v815wKdB5o0FWgI/AwuAowFcZ/5/wEx3nz+F91eGYw0jIDHwklmcpki87CUcsfRgi0X+DzgFuFRVF7lpW1X1IeAb4F/+mUXkExHZUsLySQnXaMX+f7TjgAWqWiAiKcCpwJRAect1l+GlFTAHmAu0cNOeBvoD/3D3AVDC/QV1rGEEQ4y8ZJaVliT5C5s5qNJpBGxQ1eLagH4GGvsnqOoFqlqzhOWCEq7REueB7OM4v+1TgV9UdXsQeWOFlvztZI4WkQuBysAE4Bj2N6ri7i/YYw0jGDx9ySyn9qR/YavktYAYZyVQR0QyVXVNkX3HAr+X5+QicghQF+cNycdxwAx33b96L1DeiOAaZXt3M81Nu9PdnuHveEWkMtAcx3DqA92Ac4BeQDOcF6KFfqf3v79QjzWMYAj0knmcf0IpL5Kl0RJ40G87XC+OLYGJQD0OfGF7i2R4YVNVW0pYAAG+Br4CmuLU+56P8wazA2hRJP9/gdwSlv8Wc/4OwEK/7UrAeuA0d3sh0DyYvG7abGA48BvOg/0NYAlwvrv/EWAqztveyUAq8D2O47kYeDnA9zEQGFjK/tbAVvd7qwko8L67rzswu0h+//sL+thI36ctibO4v38BkOlu5wAd3PXPgRHlPP8hwD4gzS/tA+BfxeT9BNjiLrvdxbf9SZG8lYG9wOE4JcBprr2ciFN62g1U8vr7jfjv57WAWF+A6sBzwFr3obnT/WO3DMO5O+KU0irhlBCeda9RC2gCLA8mr7v/IPchXQ1oA8wHquKUfp5186S7nycAT7vrNwEjgA8D/eEJ7KCuB7722+4AHOSuPwm87rev6P0FdWw07tOWxFnw+CWzFF2BbCmkl71EXayKLwCqug24E7jTHa9ztapOC9PpZ+BUBSwCVuG8JeWo6mYR6c7fvdtKzevuPxb4QFVzRaQ2zp95p7v+u1tF+LSINMBxgG+7x30DvAQcqqr7ynk/LfGrF1fVbL99rYCP/LbPL3J/wR4bC/dpxAmqqiJyPjAI5z9QFxiL47ROVtX5RfKfG+IlBKgqIpVwSmpPAgdT/iq+VsBcdTzSFrdrvE9rS5Kg/QlAXI9sxBgi8ikwXFU/DZiZwnahP1R1vIjcC/yqqpNF5AEc59YFmKyqn4nISOAdHCP6APgSyFPVpyJxLyXoDen+/I67kzi6TyO2CPdLptsT9UOc9lPfi+ONqnpYgOMGAqjqwBL2D3d267+K2fc/4CNVHV4u8fGA10U4W4pfgD5AlRDyvwEc7q6PxikpAIwD6gCX47yBjcRpj6mLU4I5BkgBfsCtGovF+4vX+7TFFlvKvlgJyjAMw4hJbByUYRiGEZOYgzIMwzBiEnNQhmEYRkySsN3MDzroIG3cuDEAO3bsID093VtBfpiewMSaptL0/PjjjxtU9eAoS4oo/vYD8fV7eIHpCUyZbMjrXhqRWlq3bq0+pk2bprGE6QlMrGkqTQ8JOGjS334C3b8XmJ7SiTU9qmWzIaviMwzDMGISc1CGYRhGTGIOyjAMw4hJEraTRHHs3buXnJwcdu/e7amOGjVqsHBh7ESOCEVPWloaDRo0ICUlJcKqDCN2CPTsiGebjhY1atRgxYoVIT0/kspB5eTkkJGRQePGjXECbXrD9u3bycjICJwxSgSrR1XZuHEjOTk5NGnSJArKDCM2CPTsiFebjibbtm0jLy8vpOdHUlXx7d69mzp16njqnOIZEaFOnTqel0ANI9rYs6P8lOX5kVQOCrA/WDmx789IVuy/X35C/Q6TzkEZRqg8/fTTrF692msZMcGqVat44YUXWLNmjddSjCTAHFSUycnJ4fLLLycrK4vDDz+cO+64g7y8vBLzb9myhZdeeqlw+88//6Rbt27RkGoA3377Lffccw8zZ870WkpMsHjxYv7973+zfPlyr6UkJTk5OXTu3Dlpnh/moKKIqtKlSxcuuOAClixZwm+//UZubi73339/iccU/YPVr1+fcePGRUOuATz++OPUrl2bTp06eS0lJkhNTQVgz549HitJPnzPj4suuihpnh/moErj3XehcWOoUMH5fPfdcp1u6tSppKWlcdVVVwFQsWJFnn32Wd58801eeuklOnfuTIcOHcjKyuLhhx8GoF+/fixbtoyWLVty77338vvvv9OiRQsARo0axUUXXcSZZ55J48aNGT58OM888wytWrXixBNPZNOmTQB06NCB2bNnA7BhwwZ8c6z5ju/cuXOpxycr8+bN4+OPP+aOO+6gSpUqXsuJCXwOqrS3doOwPzvg7+dHr169gNKfH48//jgQnedHoOPLgzmoknj3XbjxRli5ElSdzxtvLNcfbcGCBbRu3Xq/tOrVq9OwYUP27dvH999/z/jx45k7dy5jx45l9uzZDBkyhMMPP5w5c+bw5JNPHnDO+fPn8+GHH/LDDz9w//33U7VqVX7++WdOOukk3n777YCa5s+fz+jRo8t8fCIzZMgQ0tPTuf32272WEjOYgwqCCDw7ILTnx0cffRS150d5jg+EOaiSuP9+2Llz/7SdO530CHHmmWdSp04dqlSpQpcuXZgxY0bAYzp27EhGRgYHH3wwNWrUKKyKOuaYY/j9998jfnyismzZMt5//31uueUWateu7bWcmKFy5cqAOahS8eDZAfs/Pzp16pQQzw9zUCWxalVo6UHwj3/8gx9//HG/tG3btrFq1SoqVap0QBfMYLpk+h4YABUqVCjcrlChAvv27QOgUqVKFBQUABwwBiGY45ORoUOHkpKSwt133+21lJjCSlBBEIFnByTn88McVEk0bBhaehCcfvrp7Ny5k/feew+A/Px8evfuTc+ePalatSpTpkxh06ZN7Nq1i48++oiTTz6ZjIwMtm/fXuZrAjRu3Ljwjx1PDaResXr1akaNGsW1115LZmam13JiCnNQQRCBZwf8/fzwVZ2V9vyYPHlyQjw/zEGVxKOPQtWq+6dVreqklxERYcKECXz00UdkZWXRrFkz0tLSeOyxxwBo27YtXbt25dhjj6Vr1660adOGOnXqcPLJJ9OiRQvuvffeMl33nnvu4eWXX6ZVq1Zs2LChzPqThaeeeor8/Hz69OnjtZRiEZEjRWSO37JNRO4UkdoiMkVElriftdz8IiLPi8hSEZkrIseX9drmoIIgAs8O+Pv5MXbs2IDPjwsvvDAxnh/FBYlKhKW4gIW//vprwKBa+zF6tGqjRqoizufo0aEdXwLbtm07IG3kyJF62223heX8oVKcntII+XssA14FXFu/fr1WqVJFr7nmmqD14GHAQqAisBZoBAwF+rnp/YAn3PXzgP8CApwIfBfovCUFLPzrr78U0Oeffz7YrzQiRPv/Eeg/f4ANRejZURJFnx+h2nQ08Gkq7rssyYaSarLYkOne3VmMpOG5555j9+7d9OvXz2spwXI6sExVV4pIZ6CDm/4WkA30BToDb7sPglkiUlNEMlU15OkgrAQVJPbsCAvmoGKEnj170rNnT69lJDVbtmxh+PDhdOvWjebNm3stJ1guB8a463X9nM5aoK67fijwh98xOW7afg5KRG4EbgSoW7cu2dnZhftyc3PJzs4udEyLFi3ab3+08emJFjVq1Ci1LSc/P7/cbT3loWvXrnTt2rVQg9d6isOnaffu3UH/dlFxUCJyGPA2jsEoMEJVhxXJ0x3nbU+A7cAtqvqLu+93Ny0f2KeqbaKh20guhg8fzrZt27jvvvu8lhIUIpIKXAj0L7pPVVVENJTzqeoIYARAmzZttEOHDoX7srOz6dChQ2FvrgYNGuC/P9r49ESLhQsXlhq+ItbCW8SaHvhbU1paGq1atQrqmGiVoPYBvVX1JxHJAH4UkSmq+qtfnhXAaaq6WUTOxTGUdn77O6qqtfAbESE3N5dnn32W888/n5YtW3otJ1jOBX5S1XXu9jpf1Z2IZALr3fTVwGF+xzVw00KmQoUKVKpUyar4jKgQlV58qrpGVX9y17cDC3GqGPzzfKOqm93NWThGZBhR4dVXX2XTpk0MGDDAaymhcAV/V+8BTAJ6uOs9gIl+6de4vflOBLaWpf3JR2pqqjkoIypEvQ1KRBoDrYDvSsl2HU6vIx8KfOFWWbzqVkUUd+5i69B99dWB6pGjRazVD4eqJ5Q65LISzTaGvLw8HnvsMVq3bl3ivUW7zSMQIpIOnAnc5Jc8BPiPiFwHrAQuddM/xenJtxTYCfQqz7VTU1NtslgjKkTVQYlINWA8cKeqbishT0ccB9XeL7m9qq4WkUOAKSKySFWnFz22pDp0X311oHrkaBFr9cOh6gmlDrmsRLON4cUXX2TTpk2MHz++xGtGu80jEKq6A6hTJG0jTq++onkVuC1c17YSlBEtojZQV0RScJzTu6r6YQl5jgVeBzq7xgaAqq52P9cDE4C2kVdsJAN5eXk88cQTtG/fntNOO81rOXGBOSgjWkTFQYkzKdQbwEJVfaaEPA2BD4GrVfU3v/R0t2OFr1rjLGB+5FVHnuuuu26/WC1G9Hnrrbf4448/GDBggIX0DpLKlSubg/KYZHl2RKsEdTJwNfBPv+lZzhORm0XkZjfPgzhVFi+5+2e76XWBGSLyC/A9MFlVP4uS7ogwefJkmjZtyqRJkxg8eDAnnHACCxYs8FpW0rF3714ef/xx2rRpw1lnneW1nLjBSlDeMXHixKR6dkSlDUpVZ+CMbyotz/XA9cWkLweOi5C0qLNs2TLuuusuvvrqK55//nnatGlD9erV6datG/Pnz6dixYpeS0wa3nvvPVasWMGwYcOs9BQC5qC8YdmyZdxyyy1kZ2cnzbMjaWeSuPPOO5kzZ05Yz9myZUuee+65UvN88cUXnH/++TRr1qwwrUuXLjzwwAMsWbKEo446KqyajOLJz8/nscceo2XLllxwwQVey4krkt1BFffsyM/PL5eDCPbZcdFFFyXVsyNpHVQs4nS2MqLBBx98wG+//ca4ceOs9BQi1s089kjUZ0fSOqhAbyuR4owzzmDQoEEsW7asMG3ixIns3bt3vzcjI3IUFBQwePBgjj76aC6++GKv5cQdyV6CKu7ZEY2hI2eccQaPPPIIvXv3LkxL9GdH0joor8jKyuKpp57i9NNPZ9euXUycOJHatWszfvz4hKxDjkXGjx/PwoULef/996lQwUKihUrlypXZsWOH1zKSjqysLIYPH55Uzw6zTg+48MILWbFiBeeddx59+vThl19+4ZhjjvFaVlJQUFDAI488wpFHHkm3bt28lhOXJHsJyku6dOmSVM8OK0F5hIgwcuRIr2UkHRMnTmTevHmMHj06Yd86I405KG9JpmeHlaCMpEFVGTRoEFlZWVx22WVey4lbzEEZ0cJKUEbS8PHHHzNnzhxGjRpFpUr21y8r1ovPiBZJV4JK1O6Y0SJevz9V5eGHH+bwww+nu4XiLhfJWoKK1/9+LBHqd5hUDiotLY2NGzfaH62MqCobN24kLS3NaykhM3nyZH766Sfuv/9+Kz2Vk2R0UPbsKD9leX4klaU2aNCAnJwc/vrrL0917N69O6Ye8qHoSUtLo0GD+Iol6Ss9NW3alKuuusprOXFPMk4WG+jZEc82HS12795NzZo1Q3p+JJWDSklJoUmTJl7LIDs7O+LxlEIh1vSEm08//ZTZs2fzxhtvkJKS4rWcuCcZS1CBnh2xZkOxpgfKpimpqviM5ENVGThwIE2aNOHqq6/2Wk5CkJqayt69e626y4g4SVWCMpIPKz2Fn9TUVMAJ9li5cmWP1RiJjJWgjITFSk+Rwd9BGUYkMQdlJCyTJ09m9uzZ3H///QlXehKRmiIyTkQWichCETlJRGqLyBQRWeJ+1nLziog8LyJLRWSuiBxfnmubgzKiRbRCvh8mItNE5FcRWSAidxSTp0QjEpEertEtEZEe0dBsxDf+padrrrnGazmRYBjwmaoehRPQcyHQD/hSVbOAL91tgHOBLHe5EXi5PBf2VeuZgzIiTbTaoPYBvVX1JxHJAH4UkSmq+qtfHn8jaodjRO1EpDbwENAGUPfYSaq6OUrajTjkk08+4ccff0zIticRqQGcCvQEUNU8IE9EOgMd3GxvAdlAX6Az8LY6vRpmuaWvTFVdU5brWwnKiBbRCvm+Bljjrm8XkYXAoYC/gyrWiHAMboqqbgIQkSnAOcCYaGg34g9V5aGHHqJp06aJ2vbUBPgLGCkixwE/AncAdf2czlqgrrt+KPCH3/E5btp+DkpEbsQpYVG3bl2ys7ML9+Xm5hZuL126FICvv/6aFStWhO+uQsBfTyxgegJTFk1R78UnIo2BVsB3RXaVZEQlpRd37mINLNZ+LNMTmPJomjFjBj///DN9+/Zl5syZnuuJAJWA44F/qep3IjKMv6vzAFBVFZGQ+oGr6ghgBECbNm20Q4cOhfuys7PxbfsGq7Zq1Yqjjz66zDdRHvz1xAKmJzBl0RRVByUi1YDxwJ2qui3c5y/JwGLtxzI9gSmrpoKCAu68806ysrIYPHhw2KY1irHvKAfIUVXfS944HAe1zld159Y+rHf3rwYO8zu+gZtWJnxVfDZhrBFpotaLT0RScJzTu6r6YTFZSjKisBqXkdhMmDCBX375hQceeCBh59xT1bXAHyJypJt0Ok51+STA14moBzDRXZ8EXON2RDoR2FrW9iewNigjekTFgkVEgDeAhar6TAnZJgG3i8j7OJ0ktrpvgp8Dj/m6zAJnAf0jLtqIOwoKCnjooYc48sgjufLKK72WE2n+BbwrIqnAcqAXzgvnf0TkOmAlcKmb91PgPGApsNPNW2bMQRnRIlqvmCcDVwPzRGSOm3Yf0BBAVV+hBCNS1U0i8gjwg3vcIF+HCcPwZ+zYsSxYsIAxY8YkfLRcVZ2D07O1KKcXk1eB28J1betmbkSLaPXimwFIgDwlGpGqvgm8GQFpRoKQn5/PwIEDOfroo7nkkku8lpPQWAnKiBaJWUlvJB3vvfceixYtYty4cQlfevIac1BGtLCpjoy4Z+/evQwaNIjjjjuOiy++2Gs5CY/14jOihZWgjLjn7bffZunSpUycOJEKFeydK9JYCcqIFmbNRlyzZ88eBg0axAknnECnTp28lpMUmIMyooWVoIy45o033mDVqlWMGDECZzSDEWmsF58RLawEZcQtu3btYvDgwbRv356zzjrLazlJg5WgjGhhJSgjbnn55ZdZs2YNY8aMsdJTFDEHZUQLK0EZccn27dt5/PHHOeOMMzjttNO8lpNUmIMyooU5KCMuGTZsGBs2bGDw4MFeS0k6fHMcWjdzI9KYgzLijs2bN/PUU0/RqVMn2rVr57WcpENESE1NtRKUEXHMQRlxx5NPPsnWrVt55JFHvJaStJiDMqJByA5KRNJFxOaSMTxh3bp1DBs2jMsvv5zjjjvOazlhIR5tqnLlyuagjIgT0EGJSAURuVJEJovIemARsFZEfhWRJ0XkiMjLNAyHxx9/nD179vDwww97LaXMJIJNWQnKiAbBlKCmAYfjxGCqp6qHqerBQHtgFvCEiFwVQY2GAcCqVat4+eWX6dGjB82aNfNaTnmIe5syB2VEg2DGQZ2hqnv9E0SknhvVczww3o2WaxgRZdCgQQA89NBDHispN3FvU6mpqezevdtrGUaCE7AEVdSQXD4NIo9hhI3FixczatQobrnlFho2bOi1nHKRCDZVtWpVdu7c6bUMI8Ep60wSIQ3bF5E3gQuA9araopj99wLd/TQ1Bw52o+n+DmwH8oF9qlpcFFEjwXnggQdIS0vjvvvu81pKpIirqTCqVavGjh07vJZhJDhl7Wb+Woj5RwHnlLRTVZ9U1Zaq2hKnXv6rImHdO7r7zTklIT/99BNjx47l7rvv5pBDDvFaTqQIyaZE5HcRmScic0RktptWW0SmiMgS97OWmy4i8ryILBWRuSJyfHnFpqenm4MyIk4wvfgOeLNT1ZcC5SmSfzqwqbQ8flwBjAkyr5EE3HfffdSuXZvevXt7LSUshMOmXIq+uPUDvlTVLOBLdxvgXCDLXW4EXi6rdh/moIxoEEwV3zQRGQ9MVNVVvkQRScXpddQDp1fSqPKKEZGqOCWt2/2SFfhCRBR4VVVHlHL8jTgGSN26dcnOzgYgNze3cD0WMD2B8WmaM2cOn3/+OTfffDM///yz53rCRKRsqjPQwV1/C8gG+rrpb6uqArNEpKaIZKrqmrLegDkoIxoE46DOAa4FxohIE2ALUAWn9PUF8JyqhuvJ0QmYWaR6r72qrhaRQ4ApIrLILZEdgOu8RgC0adNGO3ToAEB2dja+9VjA9AQmOzub0047jf79+3PooYfyzDPPUKVKFU/1hPE7CodNFffiVtfP6awF6rrrhwJ/+B2b46aZgzJimoAOSlV3Ay8BL7ldXw8CdqnqlgjouZwi1Xuqutr9XC8iE4C2QLEOykgsJk2axKxZsxgxYoSnzinchMmmDnhxK3INdZ1X0JRUAwEHliA3bdrEtm3bPCt5x1qp3/QEpiyaQurF53Z9LfNbV2mISA3gNOAqv7R0oIKqbnfXzwIGReL6RmyRn5/PfffdR7NmzejVq5fXciJGWW2qhBe3db6qOxHJBNa72VcDh/kd3sBNK3rOYmsg4MAS5NSpUxk3bhynnXaaJ7G4Yq3Ub3oCUxZNUZksVkTGAN8CR4pIjohcJyI3i8jNftkuBr5QVf96g7rADBH5BfgemKyqn0VDs+EtU6ZM4ddff+XRRx8tDO9gOLhz92X41nFe3OYDk3Dar3A/J7rrk4Br3N58JwJby9P+BE4Vn6qya9eu8pzGMEolKpavqlcEkWcURRqFVXU5kBgzghpBs3v3bkaOHMkJJ5xA165dvZYTi9QFJrgll0rAe6r6mYj8APxHRK4DVgKXuvk/Bc4DlgI7gXIXSdPT0wHYsWMHVatWLe/pDKNYgnZQbrfXLJwBhb+5PYIMI+y89NJLrF+/PuFDuZfVpkp6cVPVjcDpxaQrcFv51O6Pv4M6+OCDw3lqwygkKAclIv8AxgEFvuNEpIuq/hoxZUZSsnXrVh599FHatGnDP//5T6/lRIx4tyl/B2UYkSLYNqiRwH2q2sKdqqgv8G7kZBnJytChQ9m0aRM33nij11IiTVzblDkoIxoELEGJyC04Yyaqi8g1vmSgjojcWnQEvGGUlT///JNnn32WK664gqysLK/lRIxEsCmfg8rNzfVYiZHIBFPF1xpIAU4okp7q7jOMsPDwww+zb98+Bg8ezKpVqwIfEL/EvU1ZCcqIBsEM1L1eRJYBw1R1KYCINAUuUNXrIi3QSA4WLVrEG2+8wW233UbTpk0T2kElgk2ZgzKiQbC9+O4GvnIHBAJcxP7z5RlGuejfvz9VqlRhwIABXkuJFnFtU+agjGgQlINS1YkiMg9nQCA4c4UtjZwsI5mYOXMmH330EY888kjSdFmOd5syB2VEg6DHQbljL16JoBYjCVFV+vbtS2ZmJnfddZfXcqJKPNuUOSgjGtgcMoanfPTRR8ycOZNXX3218KFnxD5paWmIiDkoI6JEZS4+wyiOvXv30q9fP5o3b861117rtRwjBETEQm4YEcdKUIZnvP766/z2229MnDjRJoSNQ8xBGZEmlLn4KgNdgcb+x6mqhb8wQmb79u0MHDiQU089lU6dOnktxxPi3abMQRmRJpTX1onAVuBHYE9k5BjJwtChQ1m/fj0ff/xxQk8IG4C4tilzUEakCcVBNVDVcyKmxEgaVq9ezdNPP81ll11G27ZtvZbjJXFtU+agjEgTSieJb0TkmIgpMZKGBx98kPz8fB5//HGvpXhNXNuUOSgj0oTioNoDP4rIYhGZKyLzRGRuMAeKyJsisl5E5pewv4OIbBWROe7yoN++c9xrLhWRfiHoNWKQuXPnMnLkSP71r3/RpEkTr+V4TZltKhYwB2VEmlCq+M4tx3VGAcOBt0vJ87WqXuCfICIVgReBM4Ec4AcRmRQvMXOMA7n33nupWbMm9913n9dSYoHy2JTnmIMyIk3QJShVXQnUBDq5S003LZhjpwObyqCvLbBUVZerah7wPtC5DOcxYoDPP/+cL774ggcffJDatWt7LcdzymNTsYA5KCPShNLN/A7gBuBDN2m0iIxQ1RfCpOUkEfkF+BO4R1UX4MTM+cMvTw7QrhSNNwI3AtStW5fs7GzAiVnjW48FklFPfn4+N998M/Xr1+foo48OeL1k+I6iYFMRxRyUEWlCqeK7DminqjsAROQJ4FsgHMb0E9BIVXNF5DzgIyDkiHWqOgIYAdCmTRvt0KEDANnZ2fjWY4Fk1DNixAh+//13xo0bx5lnnhkTmkIhQnrKZVNuFfhsYLWqXiAiTXBqGergdF2/WlXz3PFWb+PEmtoIXKaqv5dXvM9BqWoyDxUwIkgonSQEyPfbznfTyo2qblPVXHf9UyBFRA4CVgOH+WVt4KYZccT27dt54IEHaN++PV26dPFaTixRXpu6A1jot/0E8KyqHgFsxnGAuJ+b3fRn3XzlJj09nYKCAvbsibshXEacEEoJaiTwnRu/RnDi17wZDhEiUg9Yp6oqIm1xHOdGYAuQ5b4ZrgYuB64MxzWN6DFkyBDWr1/PJ598Ym/a+1NmmxKRBsD5wKPA3eJ8sf/kb/t4CxgIvIzTbjvQTR8HDBcRUVUtj/hq1aoBTvVnWlpaeU5lGMUSSriNZ0QkGzjZTeqhqnOCOVZExgAdgINEJAd4CCfkNar6CtANuEVE9gG7gMtd49knIrcDnwMVgTfdtikjTli5ciVPP/00V111FSecUDTCeXJTHpsCngP6ABnudh1gi6ruc7dzcNpwwa8tV1X3ichWN/8G/xOW1IYLxbfBrVmzBnA6vxx66KFEk2RooywPsaYHyqYpoIMSkRmq2l5EtgOKXxWE+xJWPdA5VPWKAPuH43RDL27fp8Cnga5hxCb9+/dHRHjssce8lhIzlNemROQCYL2q/igiHcKlq6Q2XCi+DW779u0ANGvWLOovH0nSRllmYk0PlE1TQAelqu3dz4xAeQ3Dn1mzZjFmzBgGDBjAYYcdFviAJCEMNnUycKHboSgNqA4MA2qKSCW3FOXfXutry80RkUpADZwq9HJRq1YtADZv3lzeUxlGsQTdScLtYRQwzTDAiZR71113Ua9ePfr27eu1nJikrDalqv1VtYGqNsZpl52qqt2BaTjV5QA9cCajBZjkbuPun1re9iegcCzbpk1lGeJoGIEJpRdfcX2D43okvBE53n//fWbNmsWjjz5a2JhuHEC4baovToeJpThtTG+46W8Addz0u4GwTBlmJSgj0gTTBnULcCvQtMg8YRnAzEgJM+KXnTt30rdvX1q1akXPnj29lhNzhNOmVDUbyHbXl+PMvlI0z27gkjLKLRFzUEakCaYX33vAf4HH2f/Na7uqWtneOICnnnqKP/74g9GjR1OhQiiF9KQhIWwqLS2NtLQ0c1BGxAimk8RWnKBqpfbEMwyAnJwcnnjiCbp168app57qtZyYJJFsqnbt2tYGZUSMUDpJvCUiNf22a4lIWAbqGolD//79yc/PZ+jQoV5LiXkSwaZq1aplJSgjYoRS/3Ksqm7xbajqZqBV2BUZccusWbMYPXo0d999t8V6Co64tylzUEYkCcVBVRCRWr4NEalNaFMlGQlMQUEBd9xxB5mZmfTv399rOfFC3NuUOSgjkoRiDE8Ds0TkPzgj37sBNj2AAcDo0aP5/vvvGTVqFBkZNqY7SOLepmrXrs2cOXO8lmEkKKHMxfe2iMzGmZBSgS4W2dYAZ8qbfv360bZtW66++mqv5cQNiWBTVoIyIkkoAQsrAy1xplWpBHQTEVR1UIS0GXHCY489xpo1a5gwYYJ1Kw+BRLCpWrVqkZuby969e0lJSfFajpFghFLFNxGna+yPgAWAMQBYunQpzzzzDNdccw3t2pUY7Ngonri3Kd90R5s3b+aQQw7xWI2RaITioBqo6jkRU2LEJXfffTepqakMGTLEaynxSNzblP9sEuagjHATSn3MNyJyTMSUGHHHZ599xscff8yAAQPIzMz0Wk48Evc2ZdMdGZEklBJUe6CniKzAqY4QQFX12IgoM2KavLw87rjjDrKysrjzzju9lhOvxL1NmYMyIkkoDqrMsyy7o+N9QdZaFLO/O85MzAJsB25R1V/cfb+7afnAPlVtU1YdRvgYNmwYv/32G5MnT6Zy5cpey4lX4j4agIXcMCJJKN3MV5bjOqNwIua+XcL+FcBpqrpZRM7Fierp3+LeUVU3FH+oEW3+/PNPBg0axAUXXMB5553ntZy4pZw2FRNYCcqIJKF0M3+wuPRgusSq6nQRaVzK/m/8NmfhRAM1YpS+ffuSl5fHs88+67WUuKY8NhUr+ByUlaCMSBBKFd8Ov/U0nCq7heGVA8B1OKEIfCjwhYgo8KqqjijpQBG5EbgRoG7dumRnZwOQm5tbuB4LxLOeuXPnMnr0aLp3705OTg45OTmea4oGEdITLZuKGCkpKdSqVYt169Z5LcVIQEKp4nvaf1tEngI+D6cYEemI46Da+yW3V9XVInIIMEVEFqnq9BI0jsCpHqRNmzbaoUMHALKzs/GtxwLxqmffvn3ccccdHHbYYbz66qukp6d7rilaREJPNGwqGmRmZrJmzRqvZRgJSHkmpqxKGKviRORY4HXgXFXd6EtX1dXu53oRmYATMbRYB2VElldeeYW5c+cyduzYiDqnJCasNhUtzEEZkSKUNqh5ONVtABWBg4Gw1JWLSEPgQ+BqVf3NLz0dqKCq2931s8J1TSM01q1bx4ABAzjzzDPp2rWr13ISgkjaVDTJzMzk66+/9lqGkYAEdFAi8o6qXg28AUxwk/cB61R1XzAXEZExQAfgIBHJAR4CUgBU9RXgQaAO8JKIwN/dyesCE9y0SsB7qvpZ0HdnhI2+ffuyc+dOXnjhBdzfwygj5bUpEUnDqUWojGMX41T1IRFpAryPY0s/4rzw5blz/r0NtAY2Apep6u/hup/MzEzWrl2Lqtp/wwgrwZSgWotIfaAX8BbOWCWA6u7ElgG776hqqaGtVfV64Ppi0pcDxwWh0YggM2fO5K233qJfv34ceeSRXstJBMprU3uAf6pqroikADNE5L/A3cCzqvq+iLyC0577svu5WVWPEJHLgSeAy8J1M5mZmezZs4ctW7YU9uozjHAQjIN6BfgSaIrzVub/iqRuupGg7Nu3j1tvvZXDDjuMAQMGeC0nUSiXTamqArnuZoq7KE7Yjivd9LeAgTgOqrO7DjAOGC4i4p6n3PimuVqzZo05KCOsBHRQqvo88LyIvKyqt0RBkxFDvPjii8ydO5dx48ZZx4gwEQ6bEpGKOM7tCOBFYBmwxa+KMAc41F0/FPjDvfY+EdmKUw24ocg5ix2mAaV3s1+7di0A//3vf1m/fn1ZbidkkmQYQpmJNT1QNk2hdDM355RkrFmzhgceeIBzzjmHLl26eC0n4SiPTalqPtBSRGritGMdFQY9xQ7TgNK72devX5+77rqLunXrRm1oQDIMQygPsaYHyqbJossZJdK7d2/y8vKsY0QMo6pbgGnASUBNEfG9dDYAVrvrq4HDANz9NXA6S4QF/yo+wwgn5qCMYvnf//7HmDFj6NevH0cccYTXcgw/RORgt+SEiFQBzsSZgWIa0M3N1gMnICLAJHcbd//UcLU/AWRkZJCenm4Oygg7QTsoEXkimDQj/tmzZw+33XYbhx9+OP369fNaTsJSDpvKBKaJyFzgB2CKqn6CExHgbhFZitPG9Iab/w2gjpt+NxD2H7VevXrmoIywE8pMEmfiGIA/5xaTZsQ5Q4cO5bfffuOzzz4jLS3NazmJTJlsSlXnAq2KSV+OM9NK0fTdwCVllxkYm03CiATBDNS9BbgVONx9Y/ORAXxT/FFGvLJ06VIeffRRLr30Us4++2yv5SQkiWhTmZmZ/PLLL17LMBKMYEpQ7+HMLv44+1cNbA9mkK4RP6gqt956K6mpqRZKI7IknE0deuihfPrppzabhBFWArZBqepWd1qUD4FNbpC1q4HXReSAaoZYZu3atVx44YV8+eWXXkuJST744AOmTJnCY489Rv369b2Wk7Akkk35aNKkCTt27GDDBosraoSPUHrxPeBO2toeOAOn4fWVyMiKDAUFBXz88ccsW7bMaykxx5YtW7jrrrto06YNt9xiQ96iRNzblI+mTZ3JL5YvX+6xEiORCMVB5buf5wMjVHUykBp+SZGjWrVqAGzfvt1jJbFH//79Wb9+Pa+++ioVK1b0Wk6yEPc25cMclBEJQnFQq0XkVZxJJj91Z0iOq3FUvql6cnNzA+RMLhYsWMCrr77KHXfcwfHHH++1nGQi7m3KR+PGjQFzUEZ4CcUYLsWJ9nm2O3q9NnBvJERFiooVK1K1alUrQfmRl5fH008/TYMGDXj44Ye9lpNsxL1N+ahatSqZmZnmoIywEspcfDtxGnV922uAuBv4kJGRYQ7Kj6eeeooVK1YwadIkMjIyvJaTVCSKTflo2rSpOSgjrIQU8l1EjgNOcTe/VtW4G/iQkZFhVXwuS5YsYdCgQZx66ql06tTJazlJSSLYlI+mTZvy1VdfeS3DSCBCmeroDuBd4BB3GS0i/wrh+DdFZL2IzC9hv4jI8yKyVETmisjxfvt6iMgSd+lR3PHBUq1aNStB4Yx5uvnmm0lLS+Pf//6313KSkvLaVKzRtGlT/vjjD/Ly8ryWYiQIoZSgrgPaqeoOKJwz7FvghSCPHwUMxwk9XRznAlnu0g4n0Fo7EamNEyK+DU5Qth9FZJKqbg5BeyFWgnIYNWoUU6dO5ZVXXqFOnTpey0lWymtTMUXTpk1RVVauXElWVpbXcowEIJROEsLf3WJx14MeMq6q04HSRsl3Bt5Wh1k4oQMygbNxJsPc5DqlKcA5IejeDytBOQOWe/fuzSmnnMINN9zgtZxkplw2FWv4uprbOEMjXIRSghoJfCciE9zti/h7tuRwUBj108UXEbSk9AMoKSKofyTHXbt2sX79ek+jTXod7XLgwIHk5uZy/fXXM336dM/1FEesaYqQnkjbVFTxlZoWL17MOeeU+R3SMAoJZrLYI4C6qvqMiGQD7d1d/+bvgGgxQUkRQf0jOb7zzjssWbLE02iTXka7nDhxIl999RWDBw/mmmuu8VxPScSapnDqiSebCoW6dety8MEHM39+sc3MhhEywVTxPQdsA1DVn1T1eVV9Htjs7gsXhVE/XXwRQUtKLxPJ3M18y5Yt3HrrrRx77LH06dPHaznJzHNEx6aiTosWLZg3b57XMowEIRgHVVdVD/jHuWmNw6hlEnCN25vvRGCrOy7kc+AsEaklIrWAs9y0MlGtWjVyc3MJY0DRuKFPnz6sXbuWN954g5SUFK/lJDPRsqmo06JFCxYsWEBBQYHXUowEIJg2qJql7KsS7IVEZAzQAThIRHJweualAKjqK8CnwHnAUmAn0Mvdt0lEHsGJHAowqDwhCTIyMigoKGDXrl1UrVq1rKeJO6ZOncprr73GvffeS5s2bbyWk+zULGVf0DYVixxzzDHk5uaycuVKmjRp4rUcI84JxkHNFpEbVPU1/0QRuR74MdgLqeoVAfYrcFsJ+94E3gz2WqXhmy0hNzc3aRzUjh07uP766zniiCNsOqPYICw2FYu0aNECgPnz55uDMspNMA7qTmCCiHTnb+NpgzPr8sUR0hUx/Gc0P+SQQzxWEx3uu+8+VqxYwfTp06lSJa5f0BOFOymHTYnIYTjjCevijA0coarD3DGDH+BUE/4OXKqqm8WJIDgMp4ZiJ9BTVX8K5w35OProowGYN2+ezU5ilJuADkpV1wH/JyIdgRZu8mRVnRpRZRHCvwSVDMyYMYMXXniB2267jVNOOSXwAUbECYNN7QN6q+pPIpKBM3h9CtAT+FJVh4hIP5xovX0pYRB82G7Ij+rVq9OoUSPrKGGEhVAmi50GTIuglqiQTDGhdu3axbXXXkujRo0YMmSI13KMIpTVpvwnlXUDHi7EGRvYGaedF+AtIBvHQRUOggdmiUhNEcl0zxN2jj32WH7++edInNpIMuIy9kx58JWgksFBDRgwgCVLlvD6668XOmYjsRCRxkAr4Duc3oE+p7MWpwoQQhjsHg7atWvH4sWL2by5TLORGUYhIc1mngj4HtSJXsU3c+ZMnn32WW6++WZOP/10r+UYEUBEqgHjgTtVdZvT1OSgqioiIY2lKGkmFghtJg1fO+drr71G27ZtQ5EQNEky00iZiTU9UDZNSeegkqEEtXPnTnr16kXDhg0ZOnSo13KMCCAiKTjO6V1V9cWUWuerunPnsVzvpgc12L2kmVggtJk0WrduzT333MPOnTsjNhtIIs80Eg5iTQ+UTVPSVvElcgmqf//+LFmyhJEjR1oQwgTE7ZX3BrBQVZ/x2zUJ8IWj6QFM9EsvbhB8RMjIyKBFixbMmjUrUpcwkoSkc1CJ3kkiOzub559/nttvv52OHTt6LceIDCcDVwP/FJE57nIeMAQ4U0SWAGe42+AMgl+OMwj+NeDWSAs86aST+O6772xGCaNcJF0VX2pqKqmpqQlZgtq+fTu9evXiiCOOsF57CYyqzqDksBwHNDiWNgg+Upx44omMGDGCRYsW8Y9//COalzYSiKQrQUHixoTq3bs3q1at4q233iI9Pd1rOUYS42tr+PLLL70VYsQ1SemgEnFG808//ZTXXnuN3r1783//939eyzGSnCZNmnD44YczZcoUr6UYcUxSOijfjOaJwoYNG7juuuto0aIFjzzyiNdyDAOAs846i2nTppGXl+e1FCNOSUoHlUglKFXlpptuYuPGjYwePZrKlSt7LckwAMdB5ebmWm8+o8wkpYNKpBLU22+/zYcffsjgwYM57rjjvJZjGIV07NiRihUr8vnnZQ7fZiQ5SemgatWqxcaNG72WUW5WrFjBv/71L0455RR69+7ttRzD2I8aNWrQvn17Jk6cGDizYRRDUjqoevXqsW7dOq9llIt9+/Zx9dVXIyK8/fbbVKxY0WtJhnEA3bp1Y8GCBSxcuNBrKUYcEjUHJSLniMhiEVnqhgIouv9Zv0GHv4nIFr99+X77JpVXS7169di6dSu7d+8u76k8Y8iQIcycOZPhw4fTuHFjr+UYRrF06dIFEWHcuHFeSzHikKg4KBGpCLyIE5fmH8AVIrLf6D1VvUtVW6pqS+AF4EO/3bt8+1T1wvLqqVvXmeQ5XktR3333HQMHDuSKK67gqquu8lqOYZRI/fr1Ofnkkxk7dqzXUow4JFolqLbAUlVdrqp5wPs4MWpK4gpgTKTE1KtXD4C1a9dG6hIRY/v27Vx55ZUceuihvPTSS/jPYG0Yscill17KvHnzmDt3rtdSjDgjWg4q6Hg0ItIIaAL4RxdNE5HZIjJLRC4qr5h4LkHddttt/P7777z77rvUrFnTazmGEZArr7yS1NRUXnvtNa+lGHFGLM7FdzkwTlXz/dIaqepqEWkKTBWReaq6rOiBJcWzKRqH5K+//gLgq6++onr16hG6jZIpa6yWKVOm8M4779CjRw/27dsXtngviRI7JpLEmp54ok6dOnTt2pXRo0czdOjQwnhRhhGIaDmooOLRuFxOkYktVXW1+7lcRLJxIoge4KBKimdTNA6Jb2R7jRo1PImZUpa4KEuXLuX555+nffv2vP7661SqFL6fLlFix0SSWNMTb9xwww2MGTOGsWPHcs0113gtx4gTolXF9wOQJSJNRCQVxwkd0BtPRI4CagHf+qXVEpHK7vpBOKEGfi2PmNTUVGrXrh03VXx79uzhsssuIyUlhXfffTeszskwokGHDh1o3rw5zzzzDM7k6oYRmKg4KFXdB9wOfA4sBP6jqgtEZJCI+PfKuxx4X/f/BzcHZovIL8A0YIiqlstBgdNRIl46SfTt25effvqJkSNH0rBhQ6/lGEbIiAj33nsvv/zyC1988YXXcow4IWqv4qr6KU7gNP+0B4tsDyzmuG+AY8KtJ14G63700UcMGzaMf//733TuXFrHR8OIbbp3786AAQMYMmQIZ599ttdyjDggKWeSAKcTRayXoFasWEGvXr1o3bo1Q4cO9VqOYZSL1NRU+vTpQ3Z2toXhMIIiaR1UrFfx7dmzh0svvRRV5T//+Y/NUm4kBDfffDONGjWiX79+Fg7eCEhSO6gdO3bE7Kzmd999N7Nnz2bUqFE0bdrUazmGERYqV67M4MGD+emnnxg9erTXcowYJ2kdVCwP1n3vvfd46aWXuOeee7jooou8lmPEICLypoisF5H5fmm1RWSKiCxxP2u56SIiz7vzYM4VkeO9U+4M3D3ppJO455572LRpk5dSjBgnaR1UZmYmAKtXlzQcyxvmzZvHDTfcwCmnnMLjjz/utRwjdhkFnFMkrR/wpapmAV+62+DMgZnlLjcCL0dJY7FUqFCBV155hU2bNnHPPfd4KcWIcZLWQWVlZQGwZMkSj5X8zZYtW+jSpQs1atTgP//5j413MkpEVacDRYsfnYG33PW3gIv80t9Wh1lATRHJjIrQEjj22GPp27cvI0eOZMKECV5KMWKYpH0CNmzYkMqVK7No0SKvpQBQUFDAVVddxe+//860adMKJ7Q1jBCoq6pr3PW1QF13vaS5MNf4pZU4VRhEZqqnDh06MHbsWHr27MnevXs55JBDgj421qaeMj2BKYumpHVQFStWpFmzZjHjoAYOHMjkyZN58cUXad++vddyjDhHVVVEQpqyoaSpwiByUz19/PHHnHDCCTz55JNMnz496Hn6Ym3qKdMTmLJoStoqPoAjjzySxYsXey2D8ePH88gjj9CrVy9uueUWr+UY8cs6X9Wd+7neTQ9lLsyocuSRRzJ69Ghmz57NddddZ13Pjf1Iagd11FFHsXz58sLJY73gl19+oUePHpx44om8/PLLFt/JKA+TgB7ueg9gol/6NW5vvhOBrX5VgZ5z4YUXMmTIEMaMGUPv3r1trj6jkKSt4gPn7S0/P59ly5bRvHnzqF9/3bp1dOrUiZo1a/Lhhx/aYFwjaERkDNABOEhEcoCHgCHAf0TkOmAlcKmb/VPgPGApsBPoFXXBAejTpw9r167lueeeIyUlhSeeeMJe1ozkdlBHHXUUAIsWLYq6g9q1axedO3dmw4YNzJgxo7Dbu2EEg6peUcKu04vJqxQJYRNriAjPPPMMeXl5PPnkk+zYsYPnn3+eihUrei3N8JCkdlDNmjUDiHo7VEFBAb169eL7779n/PjxHH+8p+MmDSMmEBGGDx9Oeno6Tz75JKtXr+add94hIyPDa2mGRyR1G1T16tWpX78+v/5a7ugdIfH666/zwQcfMGTIEC6++OKoXtswYhkRYejQobzwwgt88skntGvXjoULF3oty/CIpHZQAG3btuWbb76J2vVefvllxowZw0033cS9994btesaRjxx++23M2XKFP766y9at27NK6+8Yp0nkpCkd1Cnnnoqy5Yti8qURx9++CG33XYbJ510EsOHD7dGYMMohY4dOzJ37lzat2/PLbfcwplnnhlTM78Ykccc1KmnAjB9+vSIXmfq1KlcccUVtGvXjgcffNCmMTKMIMjMzOTzzz/nlVde4YcffqBFixb069cvZqMQGOElag5KRM4RkcXujMr9itnfU0T+EpE57nK9374e7gzNS0SkR9Fjy8Nxxx1HRkZGRB3U999/T+fOncnKymLy5MmkpaVF7FqGkWiICDfddBOLFy/m8ssv54knnqB79+488cQT5qgSnKg4KBGpCLyIM6vyP4ArROQfxWT9QFVbusvr7rG1ccZ4tAPaAg/5wgiEg0qVKtG+fXu++uqrcJ1yP+bMmcPZZ5/NIYccwhdffEHt2rUjch3DSHTq1avHW2+9xezZsznqqKPo168fDRs25IEHHuDPP//0Wp4RAaJVgmoLLFXV5aqaB7yPM8NyMJwNTFHVTaq6GZjCgWEGysVpp53GwoULycnJCedpmTNnDqeffjoZGRlMnTqV+vXrh/X8hpGMtG7dmieeeIJZs2Zx2mmn8eijj9KoUSMuueQSvvjiC/Lz872WaISJaDWEFDebcrti8nUVkVOB34C7VPWPEo49tLiLlDQbc6BZdBs0aADAoEGDuPLKK4O5n4AsWrSIvn37kpaWxpAhQ1ixYgUrVqwISk+0iTU9EHuaYk2PAe3atWPChAksW7aMl156iVGjRjFu3Djq16/PZZddxiWXXEK7du2oUCHpm9rjF1WN+AJ0A173274aGF4kTx2gsrt+EzDVXb8HGOCX7wHgnkDXbN26tfqYNm2aBqJ9+/Z65JFHakFBQcC8gcjOztaMjAxt0qSJLlu27ID9weiJJrGmRzX2NJWmB5itUbCjaC7+9hPo/r2gOD27d+/WsWPH6oUXXqgpKSkKaGZmpt5www06YcIE3bp1a1T1eEms6VEtmw1F69Ui4GzKqrpRVfe4m68DrYM9Nhz07NmTxYsX8/3335frPGPHjuWss86iQYMGfP311zRt2jRMCg3DKI3KlSvTrVs3Jk6cyPr16xk9ejQnn3wy77//PhdffDG1a9fmpJNOon///vz3v/9l8+bNXks2AhAtB/UDkCUiTUQkFbgcZ4blQopE+LwQ8A0f/xw4S0RquZ0jznLTwsoll1xCeno6Q4cOLdPxqspjjz3GpZdeSps2bZgxYwaHHlpsTaRhGBGmZs2adO/enbFjx7JhwwamTZtGv379EBGeeuopzjvvPGrXrs1RRx3FNddcw7Bhw5g+fTpbt271WrrhR1TaoFR1n4jcjuNYKgJvquoCERmEU7SbBPxbRC4E9uGEsu7pHrtJRB7BcXIAg1S1aKjrclO9enX69evHAw88EHJgra1bt3Lttdfy4Ycf0r17d15//XXrSm4YMUJqaiodOnQotOkdO3bw/fff88033/Ddd98xZcoU3nnnncL8DRs25Oijj6Z58+YcddRRNGvWjKysLDIzM21wfZSJ2mhRVf0UZ9p//7QH/db7A/1LOPZN4M2ICgR69+7N66+/zu233863334b1CSVX331FT179uSPP/7gqaee4u6777Y/sWHEMOnp6XTs2JGOHTsCTu3HmjVr+Pnnn5k7dy7z589nwYIFTJs2jd27dxceV6VKFZo0aUKTJk1o1KgRDRs25LDDDqNBgwasXr2anTt3UrVqVa9uKyGx6Qz8qFKlCq+++irnn38+Xbt25eOPP6byuHFw//2wahU0bAiPPgrdu/Pnn39y//33M2rUKA4//HC+/vprTjrpJK9vwTCMEBER6tevT/369Tn//PML0/Pz8/njjz9YvHgxS5cuZfny5SxfvpwVK1Ywc+ZMtmzZcsC5qlevTr169ahbty6HHHIIBx98MAcffDB16tQpXGrVqlW41KhRI37jwL37brHPxnBiDqoIZ599Nq+99hrXXnstxx9+OE9u2MAZe/aQCuxeuZJvr7uOd197jXe+/RZVpU+fPjz44IOkp6d7Ld0wjDBSsWJFGjduTOPGjTn77LMP2L99+3ZycnLIycnhyy+/pEaNGqxdu5a1a9eybt06fv31V/766y82btxY6kS3lStXpnr16oVLRkYG1apVo1q1aqSnp5Oenk7VqlVJT0+nSpUq+y1paWlUrly58LNy5cqkpqayatUqli1bRmpqKikpKaSkpFCpUqXCz0qVKpWv+/2778KNN8LOnc72ypXONoTVSZmDKoZevXpRr149burUifPz86kMVAc2AgV79lB1+nR63Xgjffr0sV56hpGkZGRk0Lx5c5o3b05KSkqJ7db5+fls3ryZTZs2FS6bN29m69atbNmyha1bt7J161a2b9/Otm3b2L59O+vWrWPZsmXs2LGjcMnLywv7PVSqVImKFStSsWJFKlSoUPjpv4hI4Sc4JU5Zswby8/FvzNCdO6FHD+jfH1WloKDASXfX77nnnpDa9sEcVImce+65LM7P50tgOrAdOBhoBZylSvorr3iqzzCM+KBixYocdNBBHHTQQeU6T35+Prt27dpv2bNnT+Gnb8nLy+Pnn38mKyuLvLw89u7dW7js27eP/Px89u7dS35+fuF2fn4+BQUFFBQUFK6rKvn5+Qc4G1WFkSMprkwo+flwxhkA+zm1ChUqUKtW6DPUmYMqhSqNGnHBypVcUHRHo0ZeyDEMI4mpWLFiYdVfIGrXrh1yaSUkpk51qvWK0qgRvFl8f7ayzMRic4CUxqOPQtFeOVWrOumGEWcEiihgGEETpWejOajS6N4dRoxw3gpEnM8RI8LeU8UwIk0IEQUMIzBRejZaFV8gunc3h2QkAoURBQBExBdR4FdPVRnxSxSejeagDCM5CBhRoKRoABB7s7mbntKJNT1QNk3moAzDAEBVRwAjANq0aaP+jeyhTv8VaUxP6cSaHiibJmuDMozkICpRAQwjnJiDMozkIGBEAcOINayKzzCSgJIiCngsyzBKRUqbIyqeEZG/AN9IsoOADR7KKYrpCUysaSpNTyNVPTiaYiJNEfuB+Po9vMD0BCZkG0pYB+WPiMxW1TZe6/BhegITa5piTU+0ibX7Nz2lE2t6oGyarA3KMAzDiEnMQRmGYRgxSbI4qBFeCyiC6QlMrGmKNT3RJtbu3/SUTqzpgTJoSoo2KMMwDCP+SJYSlGEYhhFnmIMyDMMwYpKEd1CxFgNHRH4XkXkiMkdEZntw/TdFZL2IzPdLqy0iU0RkifsZeujL8OoZKCKr3e9ojoicF0U9h4nINBH5VUQWiMgdbrpn35GXxJr9gNlQkHoSwoYS2kHFcAycjqra0qNxCqOAc4qk9QO+VNUs4Et320s9AM+631FLVf00inr2Ab1V9R/AicBt7n/Gy+/IE2LYfsBsKJAeSAAbSmgHhV8MHFXNA3wxcJIWVZ0ObCqS3Bl4y11/C7jIYz2eoaprVPUnd307sBAnVIVn35GHmP0Ug9lQ6YTThhLdQRUXA+dQj7T4UOALEfnRjb8TC9RV1TXu+lqgrpdiXG4Xkblu9YUn1Wki0hhoBXxHbH5HkSYW7QfMhoIl7m0o0R1ULNJeVY/HqTa5TURO9VqQP+qMO/B67MHLwOFAS2AN8HS0BYhINWA8cKeqbvPfFyPfUTJjNhSYhLChRHdQMRcDR1VXu5/rgQk41Shes05EMgHcz/VeilHVdaqar6oFwGtE+TsSkRQcw3pXVT90k2PqO4oSMWc/YDYUDIliQ4nuoGIqBo6IpItIhm8dOAuYX/pRUWES0MNd7wFM9FCL78/r42Ki+B2JiABvAAtV9Rm/XTH1HUWJmLIfMBsKloSxIVVN6AU4D/gNWAbc77GWpsAv7rLACz3AGJwi/16cNoXrgDo4vWqWAP8Danus5x1gHjDX/VNnRlFPe5yqh7nAHHc5z8vvyMslluzH1WM2FJyehLAhm+rIMAzDiEkSvYrPMAzDiFPMQRmGYRgxiTkowzAMIyYxB2UYhmHEJOagDMMwjJjEHFQUEZE6frMLr/WbbXiLiPxawjGDROSMUs7ZU0SGR051YETkIhF5MALnPUZERoX7vEb8YjYU8nnj2oYqeS0gmVDVjThTjyAiA4FcVX3Kna/qkxKOCfufNgL0AS4M90lVdZ6INBCRhqq6KtznN+IPs6HQiHcbshJU7FBRRF5z46d8ISJVAERklIh0c9dPEJFvROQXEfneN6Leh4icLyLfishB7nHPu/mX+87h5rtXRH5wJ5J82E1LF5HJ7rnni8hlbvoQceK6zBWRp4qKFpFmwB5V3eCn92URmeVet4M7WeVC/zc5EckVkSfd+/2fiLQVkWz3GH9D/RhnBgPDCITZUILZkDmo2CELeFFVjwa2AF39d4oz1cwHwB2qehxwBrDLb//FOPFVzvP90YFMnFHdFwBD3Hxnuddqi/Mm2lqcyTbPAf5U1eNUtQXwmYjUwZkm5WhVPRYYXIzuk4GfiqTVAk4C7sIZxf4scDRwjIi0dPOkA1Pd+93unvtM93qD/M41GzilpC/NMPwwG0owGzIHFTusUNU57vqPQOMi+48E1qjqDwCquk1V97n7/gn0Bc5X1c1+x3ykqgWq+it/T21/lrv8jGMUR+EY2zzgTBF5QkROUdWtwFZgN/CGiHQBdhajOxP4q0jax+pMUTIPWKeq89SZtHKB333lAZ+56/OAr1R1r7vuf+/rgfrFXNcwimI2lGA2ZA4qdtjjt55PaO2Dy4AMoFkp5xS/z8f170ibR6jqG6r6G3A8zp97sIg86BpvW2AczhvkZxzILiCthOsWFNFQ4Hdfe/XvebYK87lG6H/vafi95RpGKZgNkVg2ZA4qflgMZIrICQAikiEivj/hSpzqjLdF5OgA5/kcuFacWC2IyKEicoiI1Ad2qupo4EngeDdPDXXCRd8FHFfM+RYCR5T35kqhGbExW7UR/5gNxRnWiy9OUNU8t9H1BbfxdxdOHbpv/yIR6Q6MFZFOpZznCxFpDnwrIgC5wFU4BvKkiBTgzIp8C84b5UQRScN5a7y7mFNOB54WEfF7mwsnHYHJETivkWSYDcUfNpu5UW5EZBhOnfn/wnzeysBXOBFU9wXKbxjxitlQ8VgVnxEOHgOqRuC8DYF+8WhYhhEiZkPFYCUowzAMIyaxEpRhGIYRk5iDMgzDMGISc1CGYRhGTGIOyjAMw4hJzEEZhmEYMcn/A01V5KcTabBPAAAAAElFTkSuQmCC",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"t = np.linspace(0.5,20.0,100)\n",
"fig,axs = plt.subplots(1,2)\n",
"axs[0].scatter(x_opt_normalized,cost_function_normalized(x_opt_normalized,data),c='r',label='Optimum')\n",
"axs[0].plot(t,cost_function_normalized(t,data),'k',label='$\\Phi$')\n",
"axs[0].set_xlabel('Thickness (mm)')\n",
"axs[0].set_ylabel('Cost function $\\Phi$ (-)')\n",
"axs[0].legend()\n",
"axs[0].grid()\n",
"axs[0].set_title('$\\Phi = u/u_{max} +M/M_{max} $')\n",
"\n",
"axs[1].scatter(x_opt,cost_function(x_opt,data),c='r',label='Optimum')\n",
"axs[1].plot(t,cost_function(t,data),'k',label='$\\Phi$')\n",
"axs[1].set_xlabel('Thickness (mm)')\n",
"axs[1].set_ylabel('Cost function $\\Phi$ (-)')\n",
"axs[1].legend()\n",
"axs[1].grid()\n",
"axs[1].set_title('$\\Phi = u +M $')\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"In the next figure we plot the optimum obtained with the normalized cost function on the evolution of the deflection $u$ and mass of the beam $M$ as a function of the thickness $t$ of the beam. The _optimum_ value we obtain is a point where the deflection is not necessarily minimum but beyond that point the increase in mass becomes more critical than the gain in rigidity."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:50:20.684Z",
"iopub.status.busy": "2021-04-05T16:50:20.678Z",
"iopub.status.idle": "2021-04-05T16:50:20.799Z",
"shell.execute_reply": "2021-04-05T16:50:20.814Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": true
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA9A0lEQVR4nO3dd3gVZfbA8e8JoQVUBERAEFBYEFEQUHFFBRFXxbZrw0XFioVVLKuiP9vuWrAr9oKKLrZ1VRAQQSR2VFQEBFFEEBDpLaGT8/vjTDaXJJBLkrlzy/k8z/vcuXPnzpwM5eR95y2iqjjnnHPpJivqAJxzzrkweIJzzjmXljzBOeecS0ue4JxzzqUlT3DOOefSkic455xzaSk76gAqIisrS2vWrBl1GEmjoKCArCz/nSWW35Ot+f0oye9JSRW9J2vXrlVVjfympnSCq1mzJvn5+VGHkTRyc3Pp1q1b1GEkFb8nW/P7UZLfk5Iqek9EZF3lRVN+kWdY55xzLgye4JxzzqUlT3DOOefSkic455xzackTnHPOubTkCc4551xa8gTnnHMuLXmCc845l5YyN8FNmgRXXw2rVkUdiXPOJYVly6BfP8jLS+k5QP4ncxPcjz/Cgw/Cb79FHYlzzkVu6lQ48EAYOhSmT98p6nAqRXqk6fJo0MBelyyBffaJNhbnnIvQf/8LffvCzjvDhx/C+vUrQr+miMwB1gBbgM2q2rmyr5G5NbjCBLd4cbRxOOdcRAoK4JZb4NRTYb/97MlNly4JDaG7qnYII7lBJtfgdtvNXpcsiTYO55yLwKpVcNZZMHIknH8+PP44VK8edVSVK3NrcPXr26vX4JxzGWbmTDj4YBgzBh59FJ59ttKTW7aITIop/Uo5RoGxIvL1Nj6veBBhnDQlVK0Kdet6gnPOZZSRI6FPH0to778PRxwRymXieabWVVUXiEgDYJyI/KCqH1VmEJlbgwNrpvQmSudcBigogNtvhxNPhJYt7XlbSMktLqq6IHhdDLwFHFTZ18jsBNeggdfgnHNpb80aOO00uPlmq7198gnsuWd08YhILRHZqXAbOBqYVtnXCS3BiUhrEZkcU1aLyJUiUldExonIT8HrrsHxIiKDRWSWiEwRkY5hxfY/nuCcc2lu1izrGfn223D//fDii1CzZtRRsTvwiYh8B3wJjFLVMZV9kdCewanqTKADgIhUARZg1dCBwHhVHSQiA4P31wPHAq2CcjDwRPAaHm+idM6lsXffhb/+FapUgbFjoUePqCMyqjobaB/2dRLVRNkD+FlV5wInAUOD/UOBk4Ptk4AX1UwE6ohIo1CjatDA5qbZvDnUyzjnXCKpwl13Qa9e0Ly5PW9LluSWSIlKcL2BV4Lt3VV1YbD9O1ZVBdgDmBfznfnBvvA0aGB/E5YtC/UyzjmXKHl59rztxhuhd2/49FNLcpko9GECIlINOBG4ofhnqqoiojt4vn5AP4Ds7Gxyc3PLHdtuixezL/DV6NHkt2hR7vMki7y8vArdj3Tk92Rrfj9KSqd7smBBTW66qR2//prDpZf+zGmnzefLL3f8PGlzT1Q11II1PY6NeT8TaBRsNwJmBttPAWeWdty2Sk5OjlbIhAmqoDp+fMXOkyQmTJgQdQhJx+/J1vx+lJQu92T0aNU6dVTr1lUdN65i56roPQHyNeTcEk9JRBPlmRQ1TwKMAPoG232B4TH7zwl6U3YBVmlRU2Y4fLou51yKU4U77rDnbc2a2fO2o46KOqrkEGoTZTC+oSdwcczuQcDrInIBMBc4Pdg/GjgOmAWsBc4LMzbAJ1x2zqW0NWtsFYC33rLxbU8/DTk5UUeVPEJNcKqaD9Qrtm8Z1quy+LEK9A8znhLq1oWsLE9wzrmUM3Mm/PnPtrTlAw/AlVeCSNRRJZfMnYsSbHBIvXreROmcSykjRsDZZ0O1ajafZLduUUeUnDJ7qi7w2UyccymjoABuvRVOOglatYKvv/bktj2ZXYMDT3DOuZSwcqWt3zZqlD13e+KJpJhyK6l5Dc6n63LOJblp0+DAA+G99+Cxx+D55z25xcNrcF6Dc84lsddftxW3d9oJcnPh0EOjjih1eA2uQQOr+2/cGHUkzjn3P5s3w9//DmecAe3b2/M2T247xhNc4WDvpUujjcM55wKLF8PRR9vyNv37w4QJ0Lhx1FGlHk9wPtjbOZdEvvwSOnWCzz6zZ22PPmrDAdyO8wRXmOC8o4lzLmJDhsBhh9kQ3c8+g3PPjTqi1OYJrrCJ0mtwzrmIbNgA/frBhRfCEUfY87aOHaOOKvV5gvMmSudchObNs1rbM8/ADTfYKtz16pX9PVc2HyZQpw5kZ3sTpXMu4T74wHpJbtgAb75pc0u6yuM1OBFrpvQanHMuQVThnnugZ0/77+errzy5hcFrcOCDvZ1zCbN6NZx3ntXYTjsNnnsOateOOqr05DU4sATnTZTOuZDNmAEHHwzDh8N998Frr3lyC5PX4MDaCH7+OeoonHNp7D//sSm3cnJ8iZtE8RoceA3OOReawim3Tj8d9tsPvvnGk1uieIIDq8GtWQPr1kUdiXMujfz+Oxx1VNGUW7m5sMceUUeVOTzBgc9m4pyrdJ9+aoO1v/wSXnrJp9yKgic48ATnnKs0qjB4sDVD5uTAxIm2UKlLPE9w4NN1OecqRV4e9OkDAwbAscfCpEmw//5RR5W5PMGBT9flnKuwmTOhSxfr+n/nnfD22zZRkouODxMAaNjQXhcujDYO51xK+u9/bfB29eowdiz06BF1RA5CrsGJSB0ReUNEfhCRGSJyiIjUFZFxIvJT8LprcKyIyGARmSUiU0QkcXNp16oFdevarKfOORenTZtsCMCpp0LbtjYEwJNb8gi7ifJhYIyqtgHaAzOAgcB4VW0FjA/eAxwLtApKP+CJkGPb2p57wq+/JvSSzrnUtXDh1kMAPvoImjaNOioXK7QEJyK7AIcDQwBUdaOqrgROAoYGhw0FTg62TwJeVDMRqCMijcKKrwRPcM65OH30ERxwgHUi8SEAySvMGlwLYAnwvIh8KyLPikgtYHdVLXzY9Tuwe7C9BxDbRjg/2JcYnuCcc2VQtTkkjzwSdtkFvvjChwAkszA7mWQDHYHLVfULEXmYouZIAFRVRUR35KQi0g9rwiQ7O5vc3NxKCbbppk3svXIlH48axZZatSrlnImWl5dXafcjXfg92Zrfj5LivSd5eVW45542fPzxbhx++BKuu+4Hli7dQjrezrT5e6KqoRSgITAn5v1hwChgJtAo2NcImBlsPwWcGXP8/47bVsnJydFK8+qrqqA6bVrlnTPBJkyYEHUIScfvydb8fpQUzz357jvVli1Vq1RRvf9+1YKC8OOKUkX/ngD5GlJu2ZESWhOlqv4OzBOR1sGuHsB0YATQN9jXFxgebI8Azgl6U3YBVmlRU2b49tzTXr2Z0jkX48UXbXxbfr7NJXn11bZOskt+YY+DuxwYJiLVgNnAedhzv9dF5AJgLnB6cOxo4DhgFrA2ODZxPME552KsXw9XXAHPPAPdu8Mrr8Duu5f9PZc8Qk1wqjoZ6FzKRyVGigTV2v5hxrNdDRtCdrYnOOccv/xiY9u++QZuuAH++U/778FVHhGpAkwCFqjq8WFcw//IClWpYutYeIJzLqO98w6cc07R9vGh/NfrgAHY2Oidw7qAz0UZy4cKOJexNm+GgQPhxBNhr72s9ubJLRwi0gToBTwb5nU8wcXyBOdcRiqcleTuu+GSS2wttxYtoo4qpWWLyKSY0q/Y5w8B1wEFoQYR5slTzp572lTgW7ZYk6VzLu1NnlyH3r1hzRrrMXn22VFHlBY2q2pp/S8QkeOBxar6tYh0CzMIr8HF2nNPa6f4/feoI3HOhaygAO66C665pj116tjK257cEuJQ4EQRmQO8ChwpIv8O40Ke4GL5UAHnMsLy5fas7cYboVu3xXz1Fey7b9RRZQZVvUFVm6hqc6A38IGqhjLhmTdRxopNcIccEm0szrlQfPklnH46/PabTZLctu0MdtrJB7ilI6/BxfIanHNpS9USWteu9v7TT22ZG5+VJDqqmhvWGDjwBLe1nXe2KcI9wTmXVlavhjPOgMsvhz/9yYYAHHhg1FG5sHmCK86HCjiXVqZMgc6d4c03bRjA8OFQt27UUblE8GdwxXmCcy4tqMLzz1sz5K67wgcfwOGHRx2VSySvwRXnCc65lJefD+eeCxdcAIceCt9+68ktE3mCK27PPa0PcV5e1JE458phxgw46CB46SW47TZ47z1fBSBTeYIrrrAn5bx50cbhnNthL71kz9uWLoWxY+HWW31SokzmCa44HyrgXMpZtw4uushWATjwQGuSPOqoqKNyUfMEV5wnOOdSysyZcPDB8OyzNjPJ++9D48ZRR+WSgfeiLK5xY8jK8gTnXAp45RXo1w+qV4fRo+HYY6OOyCUTr8EVl53tC586l+TWr7dlbf76V2jfHiZP9uTmSvIEV5pmzWzNeudc0vnpJ+jSBZ56Cq69FiZMgCZNoo7KJSNPcKVp3doa9p1zSeW116BjR+vkPHIk3HMPVK0adVQuWXmCK03r1rB4MaxYEXUkzjmsSfLSS6F3b9h/f2uS7NUr6qhcsvMEV5o2bezVa3HORe7HH61J8skn4brrIDcXmjaNOiqXCjzBlaYwwf3wQ7RxOJfhXn0VOnWC+fNh1CibLNmbJF28Qk1wIjJHRKaKyGQRmRTsqysi40Tkp+B112C/iMhgEZklIlNEpGOYsW1Xixb2r8gTnHORWLfOuv+feab1kvz2WzjuuKijcqkmETW47qraQVU7B+8HAuNVtRUwPngPcCzQKij9gCcSEFvpsrOhZUtvonQuAj/8YAO3n3kGBg60XpLeJOnKI4omypOAocH2UODkmP0vqpkI1BGRRhHEZ9q08RqccwlWOJfkwoXw7rtw113eJOnKL+yZTBQYKyIKPKWqTwO7q+rC4PPfgcJ5vvcAYmc4nh/sWxizDxHph9XwyM7OJjc3N5TAW9SsSdNZs/j4/ffR7NSY8CUvLy+0+5Gq/J5sLVnvx7p1WQwe3IoxYxqx//4ruemm6dSosZFEhJqs9yRKaXNPVDW0AuwRvDYAvgMOB1YWO2ZF8DoS6BqzfzzQeXvnz8nJ0dC88IIqqM6cGd41KtmECROiDiHp+D3ZWjLej6lTVffZR1VE9eabVTdtSuz1k/GeRK2i9wTI1xBzS7wl1CZKVV0QvC4G3gIOAhYVNj0Gr4uDwxcAsS3tTYJ90fCelM6FStUmSD7wQFuCcexY+Oc/7RG4c5UhtAQnIrVEZKfCbeBoYBowAugbHNYXGB5sjwDOCXpTdgFWaVFTZuK1bm2v3tHEuUq3Zg306WNL3HTtCt9958vbuMoX5u9KuwNviUjhdV5W1TEi8hXwuohcAMwFTg+OHw0cB8wC1gLnhRhb2erUsWWAvQbnXKX65hs44wyYPRvuuMN6Smb5iFwXgtASnKrOBtqXsn8Z0KOU/Qr0DyuecvE5KZ2rNKrwyCM2QXKDBvDhh1Z7cy4s/nvT9vhQAecqxfLl8Oc/w4ABcPTRNpekJzcXNk9w29OmDSxbBkuXRh2Jcynrk0+gQwdbkPSBB2DECKhXL+qoXCbwBLc93tHEuXLbssWesXXrBtWqwWefwVVXgT2Wdy58nuC2x4cKOFcuCxfCn/4EN90Ep51mHUs6dy77e85VJh9xsj3NmkH16l6Dc24HjBkD55wDeXkwZAicd57X2lw0vAa3PVWqQKtWXoNzLg4bN9p6bcceayNsJk2C88/35Oai4zW4srRpY6NQnXPbNHu2LW3z5Zdw8cXw4INQs2bUUblM5zW4srRpY/96N2yIOhLnktJrr8EBB1hL/htv2MrbntxcMvAEV5b27a072NSpUUfiXFLJz4cLL4TevaFtWxvbdsopUUflXBFPcGXp1Mlev/462jicSyJTplivyOeegxtvhI8+gubNo47Kua15gitL8+aw667Wz9m5DKcKjz0GBx0Eq1bBuHE21s0XJXXJyDuZlEUEOnb0GpzLeMuWwQUXwPDhcNxx8MILsNtuUUfl0pmINAAOBRoD67AVaSapakE83/caXDw6dbJncBs3Rh2Jc5H48EN7HD16tPWQHDnSk5sLj4h0F5H3gFHAsUAjoC1wEzBVRP4hIjuXdR6vwcWjY0dLbt9/b93FnMsQmzfbIqS33w4tW8LEifbPwbmKEJEawEdAdSwPvaGqt8Ycchxwkar+Wsp3s4HjgZ7Af7d3HU9w8YjtaOIJzmWIuXPhr3+1OSTPOw8GD4bataOOyqWJDcCRqponIlWBT0TkXVWdCKCq127ri6q6GXg7not4govH3nvDLrt4RxOXMV5/Hfr1g4ICePllG8TtXGUJ1v/MC95WDYoWP05Eri7l66uAr1V1clnX8Wdw8RCxmpt3NHFprnBs2xlnwD772CQ+ntxcGESkiohMBhYD41T1i1IO6wxcAuwRlIuBY4BnROS6sq5RZg0uaCs9HjiMrXuyjFLV7+P7UdJAp07w6KOwaZP3iXZp6dtvLZn9+KONbbvtNv+r7sotW0Qmxbx/WlWfjj1AVbcAHUSkDvCWiLRT1WnFztME6KiqeQAicivW8eRw4Gvgnu0Gsb0PReQfWHLLBb7AMm0N4A/AoCD5XaOqU7Z3nrTQsaNN1zVjBuy/f9TROFdpCgrg4Ydh4ECoXx/efx+OPDLqqFyK26yqcS2QpKorRWQCVjMrnuAaYM/rCm0CdlfVdSJS5vyJZdXgvizWsyXWA8EYhT3LukhaiO1o4gnOpYlFi+Dcc22JmxNPtOVt6tePOiqX7kRkN2BTkNxqYj0i7y7l0GHAFyIyPHh/AvCyiNQCppd1ne0+g1PVUWV8vlhVJ23vmLTRqpV1IfPncC5NvPuu/a6Wm2uzk7z9tic3lzCNgAkiMgX4CnsGN7L4Qar6L6AfsDIol6jqP1U1X1X7lHWRuHpRikhn4P+AZsF3xK6tmVOVycqyjibek9KluA0b4PrrrVlyv/3ggw9g332jjsplkuCxVpljrkTkAlUdAkyK2TdIVQfGc514hwkMA64FpgJxTZGSljp1gqeestGv2T7CwqWeOXNyuPJK6x15+eVwzz1Qo0bUUTm3TaeIyHpVHQYgIo9h/UDiEu//0ktUdUR5ohORKlj2XaCqx4tIC+BVoB7WC+ZsVd0oItWBF4FOwDLgDFWdU55rhqZTJ1i3zlb4btcu6mici5uqrdN25ZWd2HlneOcdOP74qKNyrkynACNEpADrhLJSVS+I98vxjoO7VUSeFZEzReQvhSXO7w4AZsS8vxt4UFVbAiuAwmAvAFYE+x+k9AeO0eocdAr6orThGs4lp6VL4c9/hssug/btVzF1qic3l9xEpK6I1AVqAhcC1wFrgH8E++MSb4I7D+iAZdATglLmPxERaQL0Ap4N3gtwJPBGcMhQ4ORg+6TgPcHnPYLjk0fr1tCggc0861wKeP9960jy7rvwwAMwaNAUGjaMOirnyvQ11vL3NTABqIPlksL9cYm3ifJAVW29gwECPIRl3p2C9/WwKubm4P18bHQ6wes8sLnGRGRVcPzS2BOKSD+sVw3Z2dnk5uaWI6zya7vPPuw8diwTE3zdeOTl5SX8fiS7TL0nGzcKQ4bsxeuvN6VZs3wee2wGLVvmZez92B6/JyVFfU9UtUVlnajMAjwPtI3n2JjvHA88Hmx3A0YC9YFZMcc0BaYF29OAJjGf/QzU3941cnJyNOEGD1YF1V9+Sfy1yzBhwoSoQ0g6mXhPpk9X7dDB/ppeeqlqfn7RZ5l4P8ri96Skit4TIF93IF8UL0DXMj7fGWhX1nnircF1ASaLyC/YqPJ4hgkcCpwoIsdhvV52Bh4G6ohItlotrgmwIDh+QZDw5gfLIeyCdTZJLkccYa8ffmirfTuXJFStk+/VV0OtWjBiBJxwQtRROVcup4jIPcAYrFlyCZZHWgLdsSFr15R1kngT3DE7Gp2q3gDcACAi3YC/q2ofEfkPcCrWk7IvUDhCfUTw/vPg8w+C3wSSS7t2ULeuJbi+faOOxjkAliyx1bbfeQeOPtpW227UKOqonCsfVb0q6ExyCnAaNjB8HdZh8SlV/SSe88SV4FR1rojsitWwYr8zd4eiNtcDr4rI7cC3wJBg/xDgJRGZBSwHepfj3OHLyoLDDvOOJi5pjBlj67UtX26rbV9xhf01dS6Vqepy4JmglEu8M5n8CzgXey5WWKtSrEdkmVQ1F5uwGVWdDRxUyjHrsUyd/I44AoYPh/nzoUmTqKNxGWr9epsg+eGHbSaS997zaVKdixVvE+XpwN6qujHMYFJG7HO4PmVOh+ZcpZs61VbbnjbNZiS5+26oWTPqqJxLLvE2ZEzDxiE4gPbtbYVvb6Z0CVZQAA89BAceaM/d3n0XBg/25OZcaeKtwd0FfCsi04hZm0dVTwwlqmRXpQp07eoJziXUb7/Z0jbjxlnvyCFDYLfdoo7KufCIyGnAGFVdIyI3AR2B21U1rlnv401wQ7GpszJ7suVYRxwBo0bBwoXeXc2F7s034aKLbCrUJ5+Efv0gyeb5cS4MN6vqf0SkK3AUcC/wBHBwPF+Ot4lyraoOVtUJqvphYSlnwOmh8Dncxx9HG4dLa2vWWPf/U06BFi3g22/h4os9ubmMsSV47QU8rbZGabV4vxxvgvtYRO4SkUNEpGNh2dFI00rHjrYA6gcfRB2JS1MTJ9oShC+8ADfeCJ99ZtOhOpdBFojIU8AZwOhg1Zm4B8HE20RZuDBdl5h9cQ8TSEvZ2dCjB4webVNI+K/UrpJs2gR33AG3326jUHJzbeilcxnodGyikftUdaWINMLWJo1LvAO9u5czuPR2wgk2Hm7KFOtZ6VwF/fQTnHUWfPklnH02PPKIddh1LkM1Akap6oZgRqz9sXVD47Ldqp6InCUi2zxGRPYOHv5lpl697PWdd6KNw6U8VXj6aejQwZLca6/Biy96cnMZ77/AFhFpCTyNzab1crxfLqsGVw8bHvA1JSe8PAJbymZgOYJODw0b2oCkkSPhppuijsalqMWLrSPJyJHW6v3CCz5BjnOBArXl0/4CPKKqj4jIt/F+ebs1OFV9GBt38AqwG9AjeL8AOFtVT1HVn8ofexo44QRrT1q0KOpIXAp65x3Ybz8b2/bQQzB2rCc352JsEpEzgXOwJdcAqsb75TJ7o6jqFlUdp6q3qerFqnqlqj6lqr+WM+D0csIJ1r40enTUkbgUkpdnY9lOPBEaN4avv4YBA3ySZOeKOQ84BLhDVX8RkRbAS/F+2f85VVT79vYrtz+Hc3H6/HN71vbss3D99TYcYN99o47KueSjqtNV9QpVfSV4/4uq3h3v9z3BVZQIHH+8tS2tXx91NC6JbdoEN99ss7xt3mzd/wcNgurVo47MueQkIq1E5A0RmS4iswtLvN/3BFcZTjgB8vN9bkq3TT/8AIccYmPbzjoLvvsODj886qicS3rPY1NzbcZW8n4R+He8X453PbhbStuvqv+M90JprXt3m879nXfgT3+KOhqXRAoK4LHH4LrroFYteOMNm3bLOReXmqo6XkREVecCtwW9+kvNScXFW4PLjylbgGOB5uUINj3VrAk9e8KIEfY/mnPAggVwzDG2wnb37raGmyc353bIhmAs9k8i8jcR+TNQO94vx5XgVPX+mHIH0A3Yq1zhpqvTT4d58+DTT6OOxCWBV1+Fdu3sr8OTT9rCE77ohHM7bACQA1wBdALOBvrG++V456IsLgfw0TqxTj7Z2qD+/W+fODCDLV8Ol11mM5F06WKzkbRqFXVUzqUmVf0q2MzDhgzskHifwU3FJlcGqIIN+vbnb7Fq1bIk9/rrtsSyd43LOO+9B+efbzOT3H67DQHILu+vkM5lMBEZsb3P411sO95/fsfHbG8GFqnq5ji/mznOOguGDYN337Vk5zJCfj5cey088QS0bWt9jTpm9mJSzlXUIcA8bBatL4ByLdcS7zO4uTFlgSe3bTjqKGjQwJopXUYoHLT95JNw9dU2I4knN+cqrCFwI9AOeBjoCSzd0cW2fRxcZcrOht697Vf4lSujjsaFaONG+L//s0HbmzbZurf33w81akQdmXOpL5gicoyq9sXWIZ0F5IrI33bkPKElOBGpISJfish3IvK9iPwj2N9CRL4QkVki8pqIVAv2Vw/ezwo+bx5WbKE66yz73++NN6KOxIVk6lQ46CC4807o29eWA+zWLeqonEsvQU74Czawuz8wGHhrR84RZg1uA3CkqrYHOgDHiEgX4G7gQVVtCawALgiOvwBYEex/MDgu9XTuDH/4gz2Lc2llyxa45x77I1640Na6fe452HnnqCNzLr2IyIvA59jqNf9Q1QNV9V+qumBHzhNaglOTF7ytGhQFjgQKqzdDgZOD7ZOC9wSf9xCRcj1YjJSI1eJyc2Hu3KijcZVk1iw44gjrGXn88TBtmq0E4JwLxVlAK2wc3Gcisjooa0RkdbwnEVUt+6hyEpEq2EKpLYHHgHuBiUEtDRFpCryrqu1EZBpwjKrODz77GThYVZcWO2c/oB9AdnZ2p3HjxoUWf3lV//13uvTpw6+9e/PLRRcl7Lp5eXnUrh33IP+MUNF7ogrvvNOYJ57YmypVlAEDfuKooxaRgr96Af53pDR+T0qq6D3p3r37WlWtVYkhlY+qhl6AOsAEoCswK2Z/U2BasD0NaBLz2c9A/e2dNycnR5PWySer1q+vum5dwi45YcKEhF0rVVTknsybp3r00apgr/PmVV5cUfG/IyX5PSmpovcEyNcE5JaySkJ6UarqyiDBHQLUEZHC8XdNsNXBCV6bAgSf7wIsS0R8ofjb32DpUhv47VKKqo30aNcOPvnExreNGeMrbTtXWUSkqYhMCJbB+V5EBoRxnTB7Ue4mInWC7ZrYOIYZWKI7NTisLzA82B5B0RxjpwIfBL8JpKYjj4Q2beDRR6OOxO2AxYttQuSzz7YEN2UKXHIJKdsk6VyS2gxco6ptsWEA/UWkbWVfJMwaXCNggohMAb4CxqnqSOB64GoRmQXUA4YExw8B6gX7rwYGhhhb+ESsFvfVV/Dll1FH4+Lw5pu2svaoUdZb8sMPYe+9o47KufSjqgtV9Ztgew1W+dmjsq8T2kx5qjoFOKCU/bOBg0rZvx44Lax4InHOOXDDDVaLe/HFqKNx27B8uS1pM2wYdOoEQ4daonPOhS8Y83wANiVXpfKZTMK0005w7rk2tfzixVFH40oxerQ1Rb72Gtx2m0295cnNuQrLFpFJMaVfaQeJSG3gv8CVqhp39/94eYILW//+NrPJE09EHYmLsWoVXHAB9OoF9erBF1/ArbdC1apRR+ZcWtisqp1jytPFDxCRqlhyG6aqb4YRhCe4sLVubSOCH34YVlf6LyiuHN5/H/bbD154AQYOhEmTfIJk5xIpmMRjCDBDVR8I6zqe4BLh5pthxQp47LGoI8loeXlWoe7ZE3JybLXtu+7ypfuci8Ch2OrcR4rI5KAcV9kX8eUYE6FzZzjuOJtu/vLLwWdNSLgPP4TzzoM5c+Cqq+COO6Bmzaijci4zqeonlHONtx3hNbhEuflmWLbMn8Ul2Pr1WVx5pc32n5Vlie6BBzy5OZcJPMElSpcucPTRcO+9sHZt1NFkhE8/hQsv7MzDD9uQxO++g8MOizoq51yieIJLpFtugSVL4Kmnoo4kra1bB9dcY8lsyxbhgw/gkUegVvRTvzrnEsgTXCIdeqhN4XXXXdZP3VW6zz+HAw6wZshLLoHnnptE9+5RR+Wci4InuES75x6bhPnOO6OOJK2sWwfXXgtdu9r2uHHw+ONQs+aWqENzzkXEE1yidepkU3g99BDMnh11NGlh4kSrtd13H1x4IUydCkcdFXVUzrmoeYKLwp13Qna2LQ/tyq2w1nboodZv57337PHmzjtHHZlzLhl4gotC48aW3N54Az7+OOpoUlLxWtu0adZJ1TnnCnmCi8rf/24raF51FRQURB1Nyli3zm7dH/9otbaxY73W5pwrnSe4qOTkwN13w9dfW28IV6bPPoMOHWxCmH79rNbWs2fUUTnnkpUnuCideSb86U+2ZtzcuVFHk7TWroWrr7Yekhs2WA/JJ5/0Wptzbvs8wUVJxNrXVG3QlmrUESWdjz+G9u3hwQfh0ku9h6RzLn6e4KLWrJkN/B4zBv7976ijSRp5eTYv9eGHw5Yt8MEHthjDTjtFHZlzLlV4gksGl11mvSauvBIWLYo6msh98AHsv78ltCuugClT8NlInHM7zBNcMqhSBZ59FvLzbU2XDO1VuXq1tdT26GHDBD/6yNaJ9dWFnHPl4QkuWeyzj3UPfPddm+Ukw4wZA+3awTPPWIeSyZOtU4lzzpWXJ7hkctllcPLJMHAgTJoUdTQJsWKFVVqPPdZqap9+ank+JyfqyJxzqc4TXDIRgSFDoGFD6N3b2uzS2PDh0LYtvPQS3HgjfPutLZvnnHOVIbQEJyJNRWSCiEwXke9FZECwv66IjBORn4LXXYP9IiKDRWSWiEwRkY5hxZbU6taFl1+GX36Biy5Ky6EDS5ZY/j75ZNh9d/jqK7jjDqhePerInHPpJMwa3GbgGlVtC3QB+otIW2AgMF5VWwHjg/cAxwKtgtIPeCLE2JJb1642dOD119NqWR1Vy9377ANvvQX/+pcltwMOiDoy51w6Ci3BqepCVf0m2F4DzAD2AE4ChgaHDQVODrZPAl5UMxGoIyKNwoov6V17LfTpAzfdZG15KW7BAjjpJPuRWra05sibboKqVaOOzDmXrkQT0AQmIs2Bj4B2wK+qWifYL8AKVa0jIiOBQar6SfDZeOB6VZ1U7Fz9sBoe2dnZncaNGxd6/FHJ2rCBDldeSa05c/jmscfI32uv7R6fl5dH7STrU68Ko0Y14skn92bzZuH883/hlFPmU6VKYq6fjPckSn4/SvJ7UlJF70n37t3XqmqtSgypfFQ11ALUBr4G/hK8X1ns8xXB60iga8z+8UDn7Z07JydH096CBaqNGqk2a2bb2zFhwoSEhBSvWbNUu3dXBXudNSvxMSTbPYma34+S/J6UVNF7AuRryLklnhJqL0oRqQr8Fximqm8GuxcVNj0Gr4uD/QuApjFfbxLsy2yNG8OIEbBsmU3MvGJF1BGVacsWeOAB2G8/Wyzhqadg/HjYe++oI3POZZIwe1EKMASYoaoPxHw0AugbbPcFhsfsPyfoTdkFWKWqC8OKL6V07gxvvw0//ggnnGDT6yepadNs1rFrrrEZSb7/3pa2EYk6MudcpgmzBncocDZwpIhMDspxwCCgp4j8BBwVvAcYDcwGZgHPAJeFGFvq6dHDuiB+/jmcdhps2hR1RFvZsAFuvRU6doTZs+GVV6zi2aRJ1JE55zJVdlgnVusssq3f23uUcrwC/cOKJy2ccoothNavnyW5115LisFjEyfCBRfA9Olw1lm2tE39+lFH5ZzLdD6TSaq56CJ49FEbOnDKKbB+fWSh5OXBgAHWJLlmDYwebbOSeHJzziUDT3CpqH9/67kxapQNLlu3LuEhvPeeTY78yCM2heb339t8ks45lyw8waWqfv3guedg3DjrXbl8eUIuu3QpnHMOHHMM1KxpK24/+qgvROqcSz6e4FLZeedZb44vvoCuXake4mKpsdNsvfKKzULy7bdw6KGhXdI55yrEE1yqO+MMay/87Tc69u8P331X6ZeYOxd69bJptvbaC775xuaRrFGj0i/lnHOVxhNcOujWDT75BM3KsirVm2+W+ZV4bNkCgwfDvvva6toPPgiffWYDuJ1zLtl5gksX7drxzeOPW8+PU06BW26BgoJyn27aNFvUYMAAOOwwe3/llSRsDknnnKsoT3BpZGP9+pCba8/m/vUvW3BtB6f2Wr/enq8dcADMmgX//rd1/2/ePIyInXMuPJ7g0k2NGrYq+COPwLvvWqb6/PNSDx02bBjNmzcnKyuL5s2bc9NN42jf3hYf/etfYcYMe+7m02w551KRJ7h0JAJ/+xt8+ilkZVkb46BBWzVZDhs2jH79+jF37lxUd2Hu3P/jjjt6smLFGsaOhaFDfcC2cy61eYJLZwcdZH35//IXuOEGOPJImygS+L//+z/Wrl0LnIqtRXsecA81ahxIz54Rxuycc5XEE1y622UXm7NyyBDr37///vDkk8ydq9gCDv8B5gMHAtczf/6PkYbrnEt/IvKciCwWkWlhXscTXCYQgfPPh2nT2NLlUB669AeE74EjgauALsBkAOrWrRtdnM65TPECcEzYFwltNQGXfL5dtif9Vo1hEkI7RjOLS1nPr1GH5ZzLMKr6kYg0D/s6XoPLAPn5cO21cOCBMG+e8OqrsJReXMCvNMPWNCqclGR5gua0dM65sHmCS3OjR9tMJPfdZ62UM2bY7F7VmzXjMaA11sVkHfAucGTDhpHG65xLC9kiMimm9IsiCE9waer33y2R9eoFOTk21dbTT8Ouu9rnd9xxBzk5OYwF9qPoSdzYRYvg4othwYLognfOpbrNqto5pjwdRRCe4NJMQYEt+t2mDbz9NvzznzZS4LDDtj6uT58+PP300zRr1ozNIrzVrBnvP/EEWZddBs8/Dy1bwnXXwbJlkfwczjlXUZ7g0sgvv9Sia1e49FLo2BGmTIGbb4bq1Us/vk+fPsyZM4eCggLmzJnDqZdcYjOgzJwJp51m7ZrNm8ONN9pCcM45VwlE5BXgc6C1iMwXkQvCuI4nuDSwdq2N477ook78+KPNQjJ+PLRuXc4TtmgBL75oGbJXL5sFpXlz66niTZfOuQpS1TNVtZGqVlXVJqo6JIzreIJLce+9ZwsIDBoEPXsu4ocfbMXtSpk/sl07ePVV+P57OOkkeOABS37nnmvLCzjnXBLzBJeiFi6E3r3hmGOgWjWYMAGuv35mOPNH7rMPDBtmywtccgn85z+2KNzRR8M779jCcc45l2Q8waWYLVvg8ceLOpH84x+2iHe3bgm4eIsWtgLqr7/akgPTp8OJJ8If/gD33guLFycgCOeci09oCa60ucZEpK6IjBORn4LXXYP9IiKDRWSWiEwRkY5hxZXKJk+GP/4R+ve3QdtTp9q6ptvqRBKaevWs48kvv8Drr0OTJtbjskkTG5vw/vsVWmzVOecqQ5g1uBcoOdfYQGC8qrYCxgfvAY4FWgWlH/BEiHGlnLw8uOYa6NwZ5syxRUjHjYNWrSIOrGpV62354YdWm+vf3wLr2dM6pdx0E/z0U8RBOucyVWgJTlU/AorP+3QSMDTYHgqcHLP/RTUTgToi0iis2FKFKrz1lj0Ce+ABuPBC+OGHJF2EdJ994MEH4bffrGPKvvvCXXdZ8+XBB8PDD9voc+ecSxBR1fBObpNpjlTVdsH7lapaJ9gWYIWq1hGRkcAgVf0k+Gw8cL2qTirlnP2wWh7Z2dmdxo0bF1r8Ufr99xoMHtySzz+vz1575XH11T+y776rt/udvLw8ateunaAIy1ZtyRJ2Hz+eBuPHs9OsWWhWFiv3358lRxzB0sMPZ2MCVi5ItnsSNb8fJfk9Kami96R79+5rVbVWJYZUPqoaWgGaA9Ni3q8s9vmK4HUk0DVm/3igc1nnz8nJ0XSzYYPqoEGqNWuq1qqlet99qps2xffdCRMmhBpbhUyfrnrzzapt2qiCqojqYYep3n+/6qxZoV02qe9JBPx+lOT3pKSK3hMgX0PMLfGWRPeiXFTY9Bi8Fna7WwA0jTmuSbAvo3z8MRxwAAwcCH/6kz3WuuYayE6HRY322cfmDZs+3cbQ3XILrFplP2DLljbmbuBA+OQT2Lw56midc2kg0QluBNA32O4LDI/Zf07Qm7ILsEpVFyY4tsgsWQLnnQeHH25L24wYYc/e9twz6shCIGLP5267zcY3zJ5tDxgbNID777dJMxs0sN6Yzz/vM6c458ottLpBMNdYN6C+iMwHbgUGAa8H847NBU4PDh8NHAfMAtYC54UVVzIpKIAhQ+D662HNGqvA3HQT1Iq+5TpxWrSAq66ysmoVjB0Lo0bZ6+uv2zFt20KPHnDUUXDEEbDLLtHG7JxLCaElOFU9cxsf9SjlWAX6hxVLMpo82SZFnjjR/s9+/HH7fzyj7bKLDTs47TTrQjp1qs1FNn48PPusTQSdlWXtuN26WTn00KI1gJxzLkY6PN1JKatX2+OnRx6x8dJDh8LZZydht/+oicD++1u59lrYsAE+/xxyc6088og1aYI9vzvsMEt2f/yjjcHzG+pcxvMElyCq1uJ21VU2HOzii+HOO73yEbfq1YtqbQDr1sEXX1inlE8+sdHvTwTzAzRsCF26wMEHU6daNejUCXbaKarInXMR8QSXAD/+aJN8vP++rdP29ttw0EFRR5XiatbcOuFt2WK9Mz/7zMoXX8Dbb9MB4O9/t8k7O3e2ZNepE7Rv70nPuTTnCS5Ea9daLe2eeyAnBx591Cbjr1Il6sjSUJUqlrTat7eHmwDLlzPl2WfZf/16mDTJphF76SX7TMRmWenQoai0b2+1P2/edC4teIILyYgRcMUVMHeuPWO7917Yffeoo8owdeuy/KCDtl5q4bff4JtvisrEifDaa0Wf16tnSwHtt58922vXzoY1eM9N51KOJ7hKNns2DBgAI0dar8jcXOsl6ZJE48ZWjj++aN/KlbZ6+XffWc/NKVPguedsUGKhPfawwept29prmzZWdt/da3zOJSlPcJVk/Xqrpd15p7WW3XuvJbqqVaOOzJWpTh0bZX/44UX7Cgps3btp06zMmGGzsAwZsnXi22UXa+r8wx9seYfYUqdOon8S51wMT3CV4N134fLL4eefbQKO++6zpdFcCsvKsuEGzZtvXdsrKLDZVWbOtKUdfvjBehF98gm8/LJ1ly1Urx7svbeVli1tUPtee1lp3NgfxjoXMk9wFTB3Llx5pfWKbN3a+jAcdVTUUblQZWVB06ZWiv9hr1tnbdQ//WTl55+tFD7ni10EtmpVm4utMIk2a1b0uuee1iTq1X/nKsQTXDmsX2+1tDvvtMcvgwbZ+LZq1aKOzEWqZk3rkLLvviU/27TJmjxnz7YyZ46VX36xB7aLFm19fFaW1fKaNrXmgMKk2qSJlT32sB6fngSd2yZPcDtozBhrjpw1C0491eYJbtq07O+5DFe1alFzZWnWrbMEOHcuzJu39fbkyfDOO/abVSwR6+Syxx6WDPfYAxo1su1GjYrKbrulyZIUzu0Y/1sfpzlzrDly+HDrTzB2LPTsGXVULm3UrGnt3K1bl/65KixfDvPnW1mwYOsyd65NZbZ0acnvZmVZkmvYkP2rVbNeoA0b2qoNu+9eVHbbDerX91qhSxue4MoQ2zsyKwvuusuaI6tXjzoyl1FErNNKvXo2IH1bNmyw5s6FC23M3++/23bwmj1rlo1d+f132Lix9HPUrWvJrnipX7/otX59i6V+fVv+wodKuCTkCW47Ro60rv6zZ8Ppp9tzN2+OdEmtenXrpLKNxQS/yc2lW7duViNcvdqS4aJFsHixlUWLbIHCwjJzJnz6qdUMYzvJxKpWrSj5Fpa6dYteY8uuuxaV2rU9MbpQeYIrxaxZ1hw5apS15rz/vi1H5lzaELExfIXj+MpSUAArVsCyZZbsliwp2l62bOvyww/WnLpsmXWu2ZbsbBsrWKdOUdIr3C7cv8suRa/FS+3a1qzi3DZ4gouRn29NkPfea7+U3nuvTbflvSNdxsvKKqqdxZMQwWqJ+fmWGJcvLyorVhSVlSu33p47t2h7W02ohURg552tFCa9wvexZaedil4LS8z7rI0bLVavTaYdT3DY3+0337Rna/PmQZ8+NkFy48ZRR+ZcChOxWlbt2uVr21+/3hLdypW22nvh9urV9r6wxL5fssTGHq5aBWvW2IznZTgcrDZZmPxq1y56rV3bnjEWbhffV6tWye3Y4h12IpXxCW76dKuljR9va2sOG2ZrZzrnIlajhvX2bNiw/OfYvBny8izZrV5tZc2aovdr1jB7yhT2atDA9uXlFR2fl2dNsIXv8/O3nqYtHlWr2lIitWpt/Vp8u7RSs2bJ19JKTo4n0m3I2AS3ahX84x+2MHTt2raUzcUX+3Ah59JK7HO+bfg1N5e9Ylec2J6CAqsV5ucXJcPCxFe4Lz+/6JjY7dh969ZZ8izct26dleJjHeNVpcrWSa9GjW2/Ft8upVTLySlfHEkmY/87f/lleOghuPBCuOMO6/3snHPblZVV1EwZxvpXW7ZYklu71sq6dSVfi2/HlvXrS9+3bFnRduExGzbY65YtJcLIue++yv/ZIpCxCe6ii+CQQ2ydS+ecSwpVqhQ9v0uUzZuLEl6QAFf/9FPirh+ijE1w2dme3JxzbqsONoGC+fMjDKjy+CAS55xzaSmpEpyIHCMiM0VklogMjDoe55xzqStpEpyIVAEeA44F2gJnikjbaKNyzjmXqpImwQEHAbNUdbaqbgReBU6KOCbnnHMpSlQ16hgAEJFTgWNU9cLg/dnAwar6t2LH9QP6AWRnZ3caN25cwmNNVnl5edSuXTvqMJKK35Ot+f0oye9JSRW9J927d1+rqtvsCioixwAPA1WAZ1V1ULkvth0p14tSVZ8GngaoVauWdot3gGYGyC2cKd79j9+Trfn9KMnvSUlh3pOYx1E9gfnAVyIyQlWnV/a1kqmJcgEQO2Fdk2Cfc8659JGwx1HJlOC+AlqJSAsRqQb0BkZEHJNzzrkdly0ik2JKv5jP9gDmxbyfH+yr/CDCOGl5qOpmEfkb8B7WLvucqn4fcVjOOed23GZV7Rx1EEmT4ABUdTQwOt7j165dqyKyLsSQUk02sDnqIJKM35Ot+f0oye9JSRW9JzW381nCHkclTS9KV3EiMikZfmtKJn5Ptub3oyS/JyWFeU9EJBv4EeiBJbavgL+G0WKXVDU455xz6S2Rj6M8wTnnnEuoHX0cVV7J1IvSVdzTUQeQhPyebM3vR0l+T0pKi3viz+Ccc86lJa/BOeecS0ue4FKciDQVkQkiMl1EvheRAVHHlCxEpIqIfCsiI6OOJRmISB0ReUNEfhCRGSJySNQxRUlErgr+zUwTkVdEpEbUMSWaiDwnIotFZFrMvroiMk5Efgped40yxorwBJf6NgPXqGpboAvQ35cZ+p8BwIyog0giDwNjVLUN0J4MvjcisgdwBdBZVdthvfl6RxtVJF4Ajim2byAwXlVbAeOD9ynJE1yKU9WFqvpNsL0G+08rlGlvUomINAF6Ac9GHUsyEJFdgMOBIQCqulFVV0YaVPSygZrBuKwc4LeI40k4Vf0IWF5s90nA0GB7KHByImOqTJ7g0oiINAcOAL6IOJRk8BBwHVAQcRzJogWwBHg+aLZ9VkS2uZxJulPVBcB9wK/AQmCVqo6NNqqksbuqLgy2fwd2jzKYivAElyZEpDbwX+BKVV0ddTxREpHjgcWq+nXUsSSRbKAj8ISqHgDkk8JNTxUVPFc6CUv8jYFaInJWtFElH7Vu9inb1d4TXBoQkapYchumqm9GHU8SOBQ4UUTmYEtxHCki/442pMjNB+aramHt/g0s4WWqo4BfVHWJqm4C3gT+GHFMyWKRiDQCCF4XRxxPuXmCS3EiIthzlRmq+kDU8SQDVb1BVZuoanOs48AHqprRv52r6u/APBFpHezqAVT6ApMp5Fegi4jkBP+GepDBnW6KGQH0Dbb7AsMjjKVCPMGlvkOBs7FayuSgHBd1UC4pXQ4ME5EpQAfgzmjDiU5Qk30D+AaYiv1fmBazd+wIEXkF+BxoLSLzReQCYBDQU0R+wmq6g6KMsSJ8JhPnnHNpyWtwzjnn0pInOOecc2nJE5xzzrm05AnOOedcWvIE55xzLi15gnPOOZeWPMG5jCEi9WLGCv4uIguC7ZUiUuqgZxH5p4gctZ1znisij4YXddlE5GQRuSWE8+4nIi9U9nmdS5TsqANwLlFUdRk2wBkRuQ3IU9X7gkmqS10zTlUrPXGE4DrgxMo+qapOFZEmIrKnqv5a2ed3Lmxeg3POVBGRZ4IFMMeKSE0AEXlBRE4Ntg8Ukc9E5DsR+VJEdoo9gYj0EpHPRaR+8L3BwfGzC88RHHetiHwlIlNE5B/BvloiMio49zQROSPYPyhYzHaKiNxXPGgR+QOwQVWXxsT7hIhMDK7bLVjUckZsbUxE8kTk3uDnfV9EDhKR3OA7scnyHTJznTSXBjzBOWdaAY+p6r7ASuCU2A9FpBrwGjBAVdtjUxiti/n8z9js/McVJhugEdAVOJ5guiMROTq41kFYbbKTiByOLTr5m6q2DxbgHCMi9YA/A/uq6v7A7aXEfSg23VSsXYFDgKuweQUfBPYF9hORDsExtbA5OvcF1gTn7hlc758x55oEHLatm+ZcMvME55z5RVUnB9tfA82Lfd4aWKiqXwGo6mpV3Rx8diRwPdBLVVfEfOdtVS1Q1ekUral1dFC+xRJTGyzhTcXm/7tbRA5T1VXAKmA9MERE/gKsLSXuRtg6b7HeCZY5mQosUtWpqloAfB/zc20ExgTbU4EPg1n1pxb72Rdjy8k4l3I8wTlnNsRsb2HHnk//DOwE/GE755SY17tUtUNQWqrqEFX9EVu+Zipwu4jcEiTQg7BJgY+nKCHFWgfU2MZ1C4rFUBDzc23Soolo/3dckAhjf/YaxNRUnUslnuCci89MoJGIHAggIjuJSGEimIs1ab4oIvuWcZ73gPODBWoRkT1EpIGINAbWquq/gXuBjsExu6jqaKy5sX0p55sBtKzoD7cdfwCmhXh+50LjvSidi4Oqbgw6fjwSdEBZhz2HK/z8BxHpA/xHRE7YznnGisg+wOe2DBl5wFlYkrpXRAqATcClWK1wuIjUwGp+V5dyyo+A+0VEYmpklak7MCqE8zoXOl8ux7kUJyIPY8/d3q/k81YHPgS6xjxvdC5leBOlc6nvTiAnhPPuCQz05OZSldfgnHPOpSWvwTnnnEtLnuCcc86lJU9wzjnn0pInOOecc2nJE5xzzrm09P9G1B5dwol4VQAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"t = np.linspace(0.5,10.0,100)\n",
"#\n",
"fig, ax = plt.subplots(constrained_layout=True)\n",
"ax2=ax.twinx()\n",
"ax.plot(t,mid_point(q,L,E,I(t,d)),'r')\n",
"ax.scatter(x_opt_normalized,mid_point(q,L,E,I(x_opt_normalized,d)),c='k',label='Optimum')\n",
"ax.set_ylabel('u (mm)')\n",
"ax.set_xlabel('Thickness (mm)')\n",
"ax2.plot(t,mass(rho,L,d,t),'b',label='Mass')\n",
"ax2.scatter(x_opt_normalized,mass(rho,L,d,x_opt_normalized),c='k')\n",
"ax2.set_ylabel('Mass (kg)')\n",
"ax.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"In the next cell we can evaluate the effect of changing the starting point of the optimization problem. Considering the minimization problem at hand we do not have problems of local minima but changing the starting point can affect the convergence of the algorithm. For various starting points we extract the optimum value, the number of iterations _niter_ as well as the number of function evaluations _nfev_."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:50:22.654Z",
"iopub.status.busy": "2021-04-05T16:50:22.647Z",
"iopub.status.idle": "2021-04-05T16:50:22.669Z",
"shell.execute_reply": "2021-04-05T16:50:22.674Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"optimum: 2.400488 nit: 18 nfev: 36\n",
"optimum: 2.400488 nit: 14 nfev: 28\n",
"optimum: 2.400439 nit: 15 nfev: 30\n",
"optimum: 2.400452 nit: 18 nfev: 36\n",
"optimum: 2.400439 nit: 18 nfev: 36\n"
]
}
],
"source": [
"x0s = [1.0,2.0,3.0,5.0,6.0]\n",
"niter,nfev = [], []\n",
"for x0 in x0s:\n",
" res = optimize.minimize(cost_function_normalized, x0, method='Nelder-Mead', options={'disp': False},args=arguments)\n",
" x_opt_tmp = res.x[0]\n",
" print('optimum: {:0.6f} nit: {:2d} nfev: {:2d}'.format(x_opt_tmp,res.nit,res.nfev))\n",
" niter.append(res.nit)\n",
" nfev.append(res.nfev)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"A first observation we can make is that the optimum value does not change significantly in the range of starting points we are looking at. The variations are of the order of 1$e^{-5}$ which is quite negligible in this case. On the other hand, the number of iterations and the number of function evaluations is changing and fewer iterations are needed when the starting point is close to the optimum. This is a usual observation in many local optimization algorithms, the closer we are from the solution, the faster we will reach it."
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"To conclude this notebook, we will evaluate the effect of the algorithm used to minimize the normalized cost function. To this end, we make a loop over three different optimization algorithms, the Nelder-Mead, the conjugate gradient method (CG) and the Broyden, Fletcher, Goldfarb, and Shanno (BFGS)."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": true,
"execution": {
"iopub.execute_input": "2021-04-05T16:50:28.205Z",
"iopub.status.busy": "2021-04-05T16:50:28.196Z",
"iopub.status.idle": "2021-04-05T16:50:28.222Z",
"shell.execute_reply": "2021-04-05T16:50:28.230Z"
},
"jupyter": {
"outputs_hidden": false,
"source_hidden": false
},
"nteract": {
"transient": {
"deleting": false
}
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"algorithm: Nelder-Mead optimum: 2.400439 niter: 15 nfev: 30\n",
"algorithm: CG optimum: 2.400416 niter: 3 nfev: 38\n",
"algorithm: BFGS optimum: 2.400446 niter: 4 nfev: 12\n"
]
}
],
"source": [
"optimizer = ['Nelder-Mead','CG','BFGS']\n",
"for algorithm in optimizer:\n",
" res = optimize.minimize(cost_function_normalized, 3.0, method=algorithm, options={'disp': False},args=arguments)\n",
" x_opt_tmp = res.x[0]\n",
" print('algorithm: {:11s} optimum: {:0.6f} niter: {:2d} nfev: {:2d}'.format(algorithm,x_opt_tmp,res.nit,res.nfev))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nteract": {
"transient": {
"deleting": false
}
}
},
"source": [
"Once again the optimum value shows only small variations but we see large differences in the number of iterations _niter_ and the number of function evaluations _nfev_. The Nelder-Mead algorithm has the highest number of iterations _niter_ compared to the two other algorithms employed here. This high number is a direct result of its formulation. Both the CG and BFGS algorithms employ the first derivative of the cost function, its gradient, which accelerates significantly the convergence speed.\n",
"\n",
"Since we have not supplied a closed form expression for the gradient of the cost function, a numerical differentiation is used to its values. This numerical approach is responsible for the higher number of function evaluations required by the CG and BFGS algorithms.\n",
"\n",
"This highlights an important aspect of numerical optimization, changing the optimization algorithm (as long as we use a local optimizer) might not have a large impact on the optimum value. The impact of the algorithm will be on the convergence speed, _how fast will we reach an optimum?_ and the numerical cost _how many function evaluation do we need?_"
]
}
],
"metadata": {
"kernel_info": {
"name": "python3"
},
"kernelspec": {
"argv": [
"/Library/Frameworks/Python.framework/Versions/3.9/bin/python3.9",
"-m",
"ipykernel_launcher",
"-f",
"{connection_file}"
],
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.10"
},
"nteract": {
"version": "0.28.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}