{ "cells": [ { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "uWElxhs6QiCg" }, "source": [ "# Electricity Market Optimization\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Home Activity

\n", " You are expected to read this entire notebook before class.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Learning Objectives\n", "\n", "After studying this notebook, completing the activities, and asking questions in class, you should be able to:\n", "* Solve more complex optimization problems using pyomo.\n", "* Create mathematical models on paper.\n", "* Answer questions using/interpreting the pyomo model results." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''\n", "This cell installs Pyomo and Ipopt on Google Colab. To run this notebook\n", "locally (e.g., with anaconda), you first need to install Pyomo and Ipopt.\n", "'''\n", "\n", "import sys\n", "if \"google.colab\" in sys.modules:\n", " !wget \"https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py\"\n", " import helper\n", " helper.install_idaes()\n", " helper.install_ipopt()" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [] }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "from pyomo.environ import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Optimal Operation of Battery Energy Storage" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Grid-scale battery energy storage systems (BESS) are expected to play a critical role in supporting wide-spread adoption of renewable electricity. In an effort to *future proof* their electric grid, California has mandated over **1 GW** of BESS power capacity be brought online by 2020. For context, the [peak electricity demand in California](https://www.caiso.com/documents/californiaisopeakloadhistory.pdf) for 2021 was 44 GW.\n", "\n", "For policy-makers, technology developers, and investors, there is a critical need to understand the true value of energy storage systems. In California and many regions throughout the United States, electricity is purchased and sold in a wholesale market with time-varying prices (units of \\$ / MWh). In principle, a smart battery operator wants to **buy low** (charge) and **sell high** (discharge). This is known as energy arbitrage. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualize Price Data\n", "\n", "The text file `Prices_DAM_ALTA2G_7_B1.csv` contains an entire year of wholesale energy prices for Chino, CA. Let's import and inspect the data using Pandas. Our text file contains only one column and no header. We will manually specify \"Price\" as the name for the single column." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Class Activity

\n", " Run the code below.\n", "
" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Price
036.757
134.924
233.389
332.035
433.694
\n", "
" ], "text/plain": [ " Price\n", "0 36.757\n", "1 34.924\n", "2 33.389\n", "3 32.035\n", "4 33.694" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = pd.read_csv('https://raw.githubusercontent.com/ndcbe/data-and-computing/main/notebooks/data/Prices_DAM_ALTA2G_7_B1.csv',names=['Price'])\n", "data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's compute some summary statistics." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Price
count8760.000000
mean32.516994
std9.723477
min-2.128700
25%26.510000
50%30.797500
75%37.544750
max116.340000
\n", "
" ], "text/plain": [ " Price\n", "count 8760.000000\n", "mean 32.516994\n", "std 9.723477\n", "min -2.128700\n", "25% 26.510000\n", "50% 30.797500\n", "75% 37.544750\n", "max 116.340000" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's make a histogram." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEGCAYAAACdJRn3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXE0lEQVR4nO3dfbRldX3f8fdHUAGNDshI7QzNoE6k+ACSEQkmjRHl0Qht1JJFypSyMm1KU2yTKsQuMRrWwpoGpSsSiRAHqyI+MlWjjohJqvIwI8ijlFFQZgTn6vAgakDw2z/27+phnDv3bLhnzj1z36+1zrp7//Zv7/37sS/3M/vh/HaqCkmS+njcuBsgSZo8hockqTfDQ5LUm+EhSerN8JAk9bbruBswCnvvvXctW7Zs3M2QpImyfv3671XV4mHq7pThsWzZMtatWzfuZkjSREnyrWHretlKktSb4SFJ6m2k4ZHk9iTXJ7k2ybpWtleStUlubT/3bOVJcm6SDUmuS3LwwHZWtvq3Jlk5yjZLkma3I848fquqDqqqFW3+dOCyqloOXNbmAY4GlrfPKuA86MIGOBN4MXAIcOZ04EiSxmMcl62OA1a36dXA8QPlF1XnCmBRkmcARwJrq2pLVd0NrAWO2sFtliQNGHV4FPC5JOuTrGpl+1TVnW36LmCfNr0EuGNg3Y2tbKbyR0iyKsm6JOumpqbmsg+SpK2M+lHdX6+qTUmeDqxN8vXBhVVVSeZkWN+qOh84H2DFihUOFSxJIzTSM4+q2tR+bgY+TnfP4rvtchTt5+ZWfROw78DqS1vZTOWSpDEZWXgkeVKSX5qeBo4AbgDWANNPTK0ELm3Ta4CT2lNXhwL3tstbnwWOSLJnu1F+RCuTJI3JKC9b7QN8PMn0fj5QVZ9JcjVwSZJTgG8Br231Pw0cA2wAfgScDFBVW5K8Fbi61XtLVW0ZYbsXnGWnf2ps+7797GPHtm9Jj97IwqOqvgkcuI3y7wOHb6O8gFNn2NaFwIVz3UZJ0qPjN8wlSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb2NPDyS7JLkmiSfbPP7JbkyyYYkH0ryhFb+xDa/oS1fNrCNM1r5LUmOHHWbJUnbtyPOPE4Dbh6YfxtwTlU9G7gbOKWVnwLc3crPafVIcgBwAvBc4CjgXUl22QHtliTNYKThkWQpcCzwnjYf4GXAR1qV1cDxbfq4Nk9bfnirfxxwcVU9UFW3ARuAQ0bZbknS9o36zOMdwOuBn7b5pwH3VNVDbX4jsKRNLwHuAGjL7231f1a+jXV+JsmqJOuSrJuamprjbkiSBo0sPJK8EthcVetHtY9BVXV+Va2oqhWLFy/eEbuUpAVr1xFu+yXAq5IcA+wGPAV4J7Aoya7t7GIpsKnV3wTsC2xMsivwVOD7A+XTBteRJI3ByM48quqMqlpaVcvobnh/oapOBC4HXt2qrQQubdNr2jxt+Reqqlr5Ce1prP2A5cBVo2q3JGl2ozzzmMkbgIuT/BlwDXBBK78AeF+SDcAWusChqm5McglwE/AQcGpVPbzjmy1JmrZDwqOqvgh8sU1/k208LVVV/wi8Zob1zwLOGl0LJUl9+A1zSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb0NFR5Jnj/qhkiSJsewZx7vSnJVkv+Y5KkjbZEkad4bKjyq6jeAE4F9gfVJPpDkFdtbJ8luLXC+luTGJH/ayvdLcmWSDUk+lOQJrfyJbX5DW75sYFtntPJbkhz5aDsrSZobQ9/zqKpbgf8OvAH4TeDcJF9P8q9mWOUB4GVVdSBwEHBUkkOBtwHnVNWzgbuBU1r9U4C7W/k5rR5JDgBOAJ4LHEV3FrRLr15KkubUsPc8XpDkHOBm4GXAb1fVP2/T52xrnerc32Yf3z7V1vlIK18NHN+mj2vztOWHJ0krv7iqHqiq24ANwCFD91CSNOeGPfP4X8BXgQOr6tSq+ipAVX2H7mxkm5LskuRaYDOwFvgGcE9VPdSqbASWtOklwB1tuw8B9wJPGyzfxjqD+1qVZF2SdVNTU0N2S5L0aAwbHscCH6iqHwMkeVySPQCq6n0zrVRVD1fVQcBSurOF/R9bc2dWVedX1YqqWrF48eJR7UaSxPDh8Xlg94H5PVrZUKrqHuBy4NeARUl2bYuWApva9Ca6G/K05U8Fvj9Yvo11JEljMGx47DZw/4I2vcf2VkiyOMmiNr078Aq6eyaXA69u1VYCl7bpNW2etvwLVVWt/IT2NNZ+wHLgqiHbLUkagV1nrwLAD5McPH2vI8mvAj+eZZ1nAKvbk1GPAy6pqk8muQm4OMmfAdcAF7T6FwDvS7IB2EL3hBVVdWOSS4CbgIeAU6vq4eG7KEmaa8OGx+uADyf5DhDgnwD/ensrVNV1wAu3Uf5NtvG0VFX9I/CaGbZ1FnDWkG2VJI3YUOFRVVcn2R94Tiu6pap+MrpmSZLms2HPPABeBCxr6xychKq6aCStkiTNa0OFR5L3Ac8CrgWm7zcUYHhI0gI07JnHCuCA9vSTJGmBG/ZR3RvobpJLkjT0mcfewE1JrqIb8BCAqnrVSFolSZrXhg2PN4+yEZKkyTLso7p/l+SXgeVV9fk2rpXDokvSAjXskOy/TzdM+rtb0RLgEyNqkyRpnhv2hvmpwEuA++BnL4Z6+qgaJUma34YNjweq6sHpmTbqrY/tStICNWx4/F2SPwF2b+8u/zDwf0bXLEnSfDZseJwOTAHXA/8e+DTbeYOgJGnnNuzTVj8F/rp9JEkL3LBjW93GNu5xVNUz57xFkqR5r8/YVtN2o3vvxl5z3xxJ0iQY9rLV97cqekeS9cCb5r5JWkiWnf6psez39rOPHct+pZ3FsJetDh6YfRzdmUifd4FIknYiwwbA/xyYfgi4HXjtnLdGkjQRhr1s9VujbogkaXIMe9nqv25veVX9xdw0R5I0Cfo8bfUiYE2b/23gKuDWUTRKkjS/DRseS4GDq+oHAEneDHyqqn5vVA2TJM1fww5Psg/w4MD8g61MkrQADXvmcRFwVZKPt/njgdUjaZEkad4b9mmrs5L8LfAbrejkqrpmdM2SJM1nw162AtgDuK+q3glsTLLfiNokSZrnhn0N7ZnAG4AzWtHjgf89qkZJkua3Yc88/iXwKuCHAFX1HeCXRtUoSdL8Nmx4PFhVRRuWPcmTRtckSdJ8N2x4XJLk3cCiJL8PfB5fDCVJC9asT1slCfAhYH/gPuA5wJuqau2I2yZJmqdmDY+qqiSfrqrnAwaGJGnoy1ZfTfKiPhtOsm+Sy5PclOTGJKe18r2SrE1ya/u5ZytPknOTbEhy3eA7RJKsbPVvTbKyTzskSXNv2PB4MXBFkm+0P+zXJ7lulnUeAv6oqg4ADgVOTXIAcDpwWVUtBy5r8wBHA8vbZxVwHnRhA5zZ2nAIcOZ04EiSxmO7l62S/LOq+jZwZN8NV9WdwJ1t+gdJbgaWAMcBL23VVgNfpPsOyXHARe2priuSLEryjFZ3bVVtaW1aCxwFfLBvmyRJc2O2ex6foBtN91tJPlpVv/NodpJkGfBC4EpgnxYsAHfx8wEWlwB3DKy2sZXNVC5JGpPZLltlYPqZj2YHSZ4MfBR4XVXdN7hs8Lsjj1WSVUnWJVk3NTU1F5uUJM1gtvCoGaaHkuTxdMHx/qr6WCv+brscRfu5uZVvAvYdWH1pK5up/JENrTq/qlZU1YrFixf3baokqYfZwuPAJPcl+QHwgjZ9X5IfJLlveyu274dcANy81Wtq1wDTT0ytBC4dKD+pPXV1KHBvu7z1WeCIJHu2G+VHtDJJ0phs955HVe3yGLb9EuDfANcnubaV/QlwNt031k8BvgW8ti37NHAMsAH4EXBya8OWJG8Frm713jJ981ySNB7Dvgyqt6r6vzzynsmgw7dRv4BTZ9jWhcCFc9c6SdJj0ed9HpIkAYaHJOlRMDwkSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6MzwkSb0ZHpKk3gwPSVJvhockqTfDQ5LUm+EhSerN8JAk9WZ4SJJ6G1l4JLkwyeYkNwyU7ZVkbZJb2889W3mSnJtkQ5Lrkhw8sM7KVv/WJCtH1V5J0vBGeebxXuCorcpOBy6rquXAZW0e4GhgefusAs6DLmyAM4EXA4cAZ04HjiRpfEYWHlX198CWrYqPA1a36dXA8QPlF1XnCmBRkmcARwJrq2pLVd0NrOUXA0mStIPt6Hse+1TVnW36LmCfNr0EuGOg3sZWNlP5L0iyKsm6JOumpqbmttWSpEcY2w3zqiqg5nB751fViqpasXjx4rnarCRpG3Z0eHy3XY6i/dzcyjcB+w7UW9rKZiqXJI3Rjg6PNcD0E1MrgUsHyk9qT10dCtzbLm99FjgiyZ7tRvkRrUySNEa7jmrDST4IvBTYO8lGuqemzgYuSXIK8C3gta36p4FjgA3Aj4CTAapqS5K3Ale3em+pqq1vwkuSdrCRhUdV/e4Miw7fRt0CTp1hOxcCF85h0yRJj5HfMJck9WZ4SJJ6MzwkSb0ZHpKk3gwPSVJvhockqbeRPaorzWfLTv/U2PZ9+9nHjm3f0lzxzEOS1JvhIUnqzfCQJPVmeEiSejM8JEm9GR6SpN4MD0lSb4aHJKk3w0OS1JvhIUnqzfCQJPVmeEiSenNgxHlknIP1SVIfnnlIknozPCRJvRkekqTeDA9JUm+GhySpN8NDktSb4SFJ6s3wkCT1ZnhIknozPCRJvRkekqTeHNtqGxxjSpK2z/CQdrBx/ePk9rOPHct+tXOamMtWSY5KckuSDUlOH3d7JGkhm4jwSLIL8JfA0cABwO8mOWC8rZKkhWtSLlsdAmyoqm8CJLkYOA64aaytkibIQryX56W60ZmU8FgC3DEwvxF48WCFJKuAVW32/iS3zHEb9ga+N8fbHJedqS+wc/VnZ+oLjLk/educbm5nOjYz9eWXh93ApITHrKrqfOD8UW0/ybqqWjGq7e9IO1NfYOfqz87UF9i5+mNfHmki7nkAm4B9B+aXtjJJ0hhMSnhcDSxPsl+SJwAnAGvG3CZJWrAm4rJVVT2U5D8BnwV2AS6sqht3cDNGdklsDHamvsDO1Z+dqS+wc/XHvgxIVc1FQyRJC8ikXLaSJM0jhockqTfDYwiTPDRKkn2TXJ7kpiQ3Jjmtle+VZG2SW9vPPcfd1mEl2SXJNUk+2eb3S3JlOz4fag9VTIQki5J8JMnXk9yc5Ncm9dgk+S/td+yGJB9MstskHZskFybZnOSGgbJtHot0zm39ui7JweNr+S+aoS9vb79n1yX5eJJFA8vOaH25JcmRw+zD8JjFTjA0ykPAH1XVAcChwKmt/acDl1XVcuCyNj8pTgNuHph/G3BOVT0buBs4ZSytenTeCXymqvYHDqTr18QdmyRLgP8MrKiq59E92HICk3Vs3gsctVXZTMfiaGB5+6wCzttBbRzWe/nFvqwFnldVLwD+H3AGQPt7cALw3LbOu9rfve0yPGb3s6FRqupBYHpolIlQVXdW1Vfb9A/o/jgtoevD6lZtNXD8WBrYU5KlwLHAe9p8gJcBH2lVJqkvTwX+BXABQFU9WFX3MKHHhu7pzd2T7ArsAdzJBB2bqvp7YMtWxTMdi+OAi6pzBbAoyTN2SEOHsK2+VNXnquqhNnsF3ffloOvLxVX1QFXdBmyg+7u3XYbH7LY1NMqSMbXlMUmyDHghcCWwT1Xd2RbdBewzrnb19A7g9cBP2/zTgHsG/qeYpOOzHzAF/E27DPeeJE9iAo9NVW0C/hz4Nl1o3AusZ3KPzbSZjsWk/134d8DftulH1RfDY4FI8mTgo8Drquq+wWXVPa8975/ZTvJKYHNVrR93W+bIrsDBwHlV9ULgh2x1iWqCjs2edP+C3Q/4p8CT+MXLJhNtUo7FbJK8ke5y9vsfy3YMj9lN/NAoSR5PFxzvr6qPteLvTp9mt5+bx9W+Hl4CvCrJ7XSXD19Gd89gUbtUApN1fDYCG6vqyjb/EbowmcRj83LgtqqaqqqfAB+jO16TemymzXQsJvLvQpJ/C7wSOLF+/iW/R9UXw2N2Ez00SrsncAFwc1X9xcCiNcDKNr0SuHRHt62vqjqjqpZW1TK64/CFqjoRuBx4das2EX0BqKq7gDuSPKcVHU73moGJOzZ0l6sOTbJH+52b7stEHpsBMx2LNcBJ7amrQ4F7By5vzUtJjqK75PuqqvrRwKI1wAlJnphkP7qHAK6adYNV5WeWD3AM3dMJ3wDeOO729Gz7r9Odal8HXNs+x9DdK7gMuBX4PLDXuNvas18vBT7Zpp/Zftk3AB8Gnjju9vXox0HAunZ8PgHsOanHBvhT4OvADcD7gCdO0rEBPkh3v+YndGeFp8x0LIDQPYX5DeB6uqfMxt6HWfqyge7exvTfgb8aqP/G1pdbgKOH2YfDk0iSevOylSSpN8NDktSb4SFJ6s3wkCT1ZnhIknozPCRJvRkeWjCSPJzk2jZk+IeT7DFDvS/Pwb7enGRTkrdsXT5D/dOTnNjWqyTPHlj2ula2IslpSd4xsOzdST4/MP+HbajwZYPDcW+1r7cnuSvJHz/WfmrhMjy0kPy4qg6qbsjwB4H/MLhwehiNqjpsjvZ3TlW9qW37yUkuAf6gvU/hf2xV90jgc236erpv0E97DXBjm/4SMNi+A4GnDgyhfRiw3fCrqv8G/FXfzkiDDA8tVP8APDvJS5P8Q5I1dMNpkOT+6UpJ3pDk+iRfS3J2K3tWks8kWd/W3X+I/Z0E3E/33oeDgIsG9vEU4AlVNdWKPkEb9j/Js+hGqP1eW3Yt8CtJdm9Duv+4lT2/LT+MLmAAdkny1+le0PS5JLsP+x9Hmo3hoQWnnWEcTfcvfOgGIzytqn5lq3pH0/0Rf3FVHQhMny2cD/xhVf0q8MfAu4bY7YPAU4Ddq+qnVTV4SenldENgTLuPbsyr59GdgXxoekF1w5tfA7yI7uVeV9K9m+Gw9kKmVNX08NrLgb+squcC9wC/M0Q7paEYHlpIdk9yLd1YUt+mvYQJuKq6l+Bs7eXA31QbRK6qtrSh7Q8DPty29W5gmJcAXQR8E1iZ5MtJXj2w7Ch+/m6FaRfTBcfxwMe3Wvbl1obDgK+0z/T84CWr26rq2ja9Hlg2RDuloew6exVpp/HjqjposKAbAJYf9tjG4+hecHTQbBUHVfcWytcn+RHdmcRnk6yrqtvp3tr2B1ut8kng7cC6qrqvtXPal+ju1+xGNzjfFN0rkqd4ZHg8MDD9MOBlK80Zzzykma0FTp5+KivJXtW9SOu2JK9pZUly4GwbSrK8DekP3Qit9wJ7JHku8PWqeniwfjvbeQNw1jY29xW6S1aLq2pzdaObTtFdYvvSNupLc87wkGZQVZ+he9fBunaJavrR1hOBU5J8je4pqGHeab8/8EXgZLp7Fp+qqpvo7r18Zob9X1zt/fNbld9NFxY3DhR/BXg68LUh2iI9Zg7JLo1A+z7H/VX151uXV9WbB+bXAifVDn6R0Eztk4blmYc0GvcDq7b+kiDd2cfPVNUrxhAcbwd+j373eqRH8MxDktSbZx6SpN4MD0lSb4aHJKk3w0OS1Nv/B9HqlN6joWdHAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.hist(data['Price'])\n", "plt.xlabel(\"Price [$/MWh]\")\n", "plt.ylabel(\"Frequency\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, let's plot the prices for the first three days in the data set." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAifUlEQVR4nO3dfZgcZZnv8e8vGUiyASYgIScIIRERl1WS4BzcCAsBDi4qC+yuq4iroK7x/V0R8BxFXUQ8Irir4kbQDUcQFUU8viBckRzfIpuEGQISRQ1GgRCCksaJJjDMff6ommSY9PRUd7q6qrt/n+uaa7qru6vvnp6++6mnnud+FBGYmVn3mFR0AGZm1lpO/GZmXcaJ38ysyzjxm5l1GSd+M7Mu01N0AFnsv//+MXfu3KLDMDNrK2vWrHk4ImaO3d4WiX/u3LmsXr266DDMzNqKpA3Vtrurx8ysy+Sa+CXNkHS9pJ9LWidpkaT9JN0i6Zfp733zjMHMzJ4s7xb/J4GbIuKZwHxgHXAesDwiDgOWp9fNzKxFckv8knqB44CrACLisYjYApwOLEvvtgw4I68YzMxsV3m2+OcBm4EvSOqXdKWk6cCsiNiY3udBYFaOMZiZ2Rh5Jv4e4CjgiohYCGxlTLdOJBXiqlaJk7RE0mpJqzdv3pxjmGZm3SXPxH8fcF9E3JZev57ki2CTpNkA6e+Hqj04IpZGRF9E9M2cucswVMtBZWWFDRdvoLKyUnQoZpaj3MbxR8SDkn4n6fCI+AVwEnB3+nM28NH09415xWDZVVZWGDh+gHg8YBJMP3I6Pb09zDprFgcuObDo8MysifKewPUW4BpJewLrgVeRHGV8RdJrgA3AS3KOwTLYsmJLkvQBhmGoMsS2e7cBOPGbdZhcE39EDAB9VW46Kc/ntfrNWDwj+UoehknTJnHENUew/vz1RYdlZjloi5INlr/eRb1MP3I6Q5UhjrjmCHoX9RYdkpnlxInfdujp7aGnt8dJ36zDuVaPmVmXceI3M+syTvxmZl3Gid/MrMs48ZuZdRknfjOzLuPEb2bWZZz4zcy6jBO/mVmXceI3M+syTvxmZl3Gid/MrMs48ZuZdRknfjOzLuPEb2bWZZz4zcy6jBO/mVmXyTXxS/qNpDslDUhanW67UNL96bYBSS/MMwYzM3uyViy9eEJEPDxm22UR8fEWPLeZmY3hrh4zszpVVlbYcPEGKisrRYfSkLxb/AHcLCmA/4iIpen2N0t6JbAaeFdEPDL2gZKWAEsA5syZk3OYZmbZVFZWGDh+gHg8YBJMP3I6Pb09zDprFgcuObDo8DLJu8V/bEQcBbwAeJOk44ArgEOBBcBG4NJqD4yIpRHRFxF9M2fOzDlMM7NstqzYkiR9gGEYqgwxODDIpms3FRtYHXJN/BFxf/r7IeAG4OiI2BQRT0TEMPA54Og8YzAza6YZi2fsyJyTpk3iiGuOYK8FexUaU71yS/ySpkvae+Qy8HzgLkmzR93t74G78orBzKzZehf1Mv3I6UyZN4X5y+fTu6i36JDqlmcf/yzgBkkjz3NtRNwk6f9IWkDS//8b4HU5xmBm1nQ9vT309Pa0ZdKHHBN/RKwH5lfZ/oq8ntPMzCbm4ZxmZl3Gid/MbALtPm5/rFbM3DUza1vVxu1vu3db243kGc0tfjOzGqqN299rwV7MOmtWsYHtBrf4zcxq2DFuf3jnuP12Hc0zwonf6lJZWWHLii3MWDyj7f/5zbIYGbc/VBnqiKQPTvxWw1BliKHKEJWVFXoX9XZEjRKzRrT7uP2xnPitqsrKClvXboVh6D+2f0eLZ2xf57Z7twE48Zu1EZ/ctaq2rNiSzK2GHUm+p7cn+Y+Z1L41SszMLX4bx4zFM5g0dRLDjw0zac+dJ7Tcx2/W/pz4rareRb3MXz5/lyTfu6jXCT9H/mK1VnDit3E5ybfWeCfPhypDTJs7jYPPPdjvhzWFE79ZSVSbKASwdWArWwe28vA3H/ZIKmsKn9w1K4lqC3wc8JIDQOkd2nS1JysfJ36zkqi2wMfISXYmeySVNY+7erqYTySWz9iJQuOdZDfbHU78XaoTKw52Kp9kt2abMPFLOgA4BjgQ+DPJGrmr08XSrU11YsVBM8tm3MQv6QTgPGA/oB94CJgKnAEcKul64NKIeLQFcVqTdWLFQTPLplaL/4XAayPit2NvkNQDnAqcDHxtvB1I+g3wR+AJYCgi+iTtB3wZmEuy2PpLIuKRBuO3BnVixUEzy2bcUT0R8Z5qST+9bSgivhER4yb9UU6IiAUR0ZdePw9YHhGHAcvT61aAnt4eps6Z6qRv1mWy9PFPAf6RpIW+4/4R8aEGn/N0YHF6eRmwAnhvg/syM7M6ZRnHfyNJsh4Cto76ySKAmyWtkbQk3TYrIjamlx8Eqp5NlLRE0mpJqzdv3pzx6czMbCJZhnMeFBGnNLj/YyPi/nRk0C2Sfj76xogISVHtgRGxFFgK0NfXV/U+ZmZWvywt/p9IenYjO4+I+9PfDwE3AEcDmyTNBkh/P9TIvs3MrDHjJn5Jd0paCxwL3C7pF5LWjtpek6TpkvYeuQw8n2QOwDeBs9O7nU3SlWRm1lEqKytsuHgDlZWVokPZRa2unlN3c9+zgBskjTzPtRFxk6RVwFckvQbYALxkN5/HzKxUyr4+da3E/w7gJ8CPR7ps6hER64H5Vbb/Hjip3v1ZOY1dkN3Mqs+ML9P61LX6+H9FMkv3x5J+I+laSW+WtFCSq3rajgXZt9+7nf5j+1m1cBX9i/t5YOkDRYdm1lJDlSG2/Xbbjm6daiW2y1QHa9wWf0R8CvgUgKQDgeelP28HDgD2aUF8VmLVFmQvU6vGrBVGGkAMQ/+x/Tu6dSZNm8QeB+xRypnxNVvuShwJnEYylv94kiOBS1sQm5Wca8WbVW8AAezdtzeHnHdI6ZI+1C7SdgtJq34A+CnwkYhY16K4rA24VrzZzgbQ8GPDTNqzPQoe1jq5ux44EjgM+D3wsKTNEfFwSyKztuBa8dbt2rEBVKuP/3UAkvYB/pqkf/9NkmYCd0XE2eM91sysm7RbAyhLyYbtwJ9IFmHZDhwE7JlnUGZmlp9aM3cvk3QbsBH4ILA38Fng8IhoqISDmZkVr1aL/17gi8BARDzRonjMzCxntRL/j9Lf89OyC08SEbfnEpGZmeWqVuJfTVJUbWQUz+jsH8CJeQVlZs1VWVlpq1Enlq9aif+dwItJTupeB9wQEYMticralmv35K/ev3HZC4ZZ69Uaznk5cLmkpwFnAsslbSCZyDXQmvCsnYw3dd0Jpnka+RuXvWCYtd6ExdbSKps3AjeTLKTyjLyDsvZUber64MAgm67dVGRYHaWRv3HZC4ZZ69Uazvk0SRekQzo/CNwB/GVEfKVl0Vlbce2e/DXyN+5d1Mv0I6czZd4U5i+f7y44q9nH/ytgLUlr/1FgDvCGkRE+EfGJ3KOzttKOU9fbTaN/457eHnp6e/yeGFA78X9w1GU32yyTdpu63o78N7bdVSvx3wPcnK6YZWZmHaLWyd05wFcl/VDShZKeq2ozucysUGNXfzKbSK3hnJcAl0jaG/gfwKuBz0paB9wEfC8iJhyuIWkyyWSw+yPiVEn/SbKgy8h/6TkeHmrWGA+htUZMWJ0zIv4I3JD+IOkI4AXA1cDfZniOtwHrePJSje+JiOvrjtbMnsTLX1oj6lo0XdIcYDgiLo2ICZO+pIOAFwFXNhifmdXgIbTWiInW3P1o2sJH0j8CPwS+LOmijPu/HDgXGB6z/SJJa9PSz1PGee4lklZLWr158+aMT2e1VFZW2HDxBvcFd5CR4Z3zPjzPY/Qts4m6ek6JiPPSy+8Ank8yvv924H21HijpVOChiFgjafGom84HHiRZzGUp8F7gQ2MfHxFL09vp6+uLsbdbfarVa9l27za3DjuAh3davWottv4BYJak9wPTgEOBl5JU6exNt6+IiB+Ms4tjgNMkvRCYCuwj6YsR8c/p7dslfQF4d5Nei9VQrV7LXgv2YtZZs4oNzMxartaong+m3TyHkJyYvToiPiRpT+D5EbFLK33M488nad2TtvjfHRH/LGl2RGxMh4aeQVL62XK2o17L8M6+4E5qJbZr2eE84naFVJvIRF09rwZeCTxGMooHkvH9F+/Gc16TLtguYAB4/W7syzIaqdcyVBnqyKTfjmWH8+h+G2945+DAoLv1ClamL+SaiT8itgJXjNn2K5J+/swiYgWwIr3sBVwK0qn1Wtq17HAe3W/Vhnf29Pa4W69gZZtvUauP/3PAv0XEnVVum07S3789Iq7JMT6zCVXrxlp//vqiw5pQHt1vI8M7hx8bZtKendel167KNt+iVov/08D/kvRskn74zSQnaQ8j6fP/POCkb4Vr126sPOJ2hdRyqvaFXGTjpNbJ3QHgJZL2AvqA2STLMK6LiF+0JjyzbNq1GyuPuD28s3zK9oWcpWTDIGn/vJlZN8hjtFWZvpAnTPxmZt2kGyY71lWrx8ys03XDZMfMLX5JfxERf8ozGDOzonX6ZEfI0OKX9DxJdwM/T6/Pl/SZ3CMzMytANyxOn6Wr5zKSuvu/B4iIO4Dj8gzKzKxIPb09TJ0ztSOTPmTs44+I343Z9EQOsZiZWQtk6eP/naTnASFpD3auqGVmZm0oS4v/9cCbgKcC9wML0utmZtagocoQ2367rZCFkbJM4HoYeHkLYjEz6wpFF23LMqpnmaQZo67vK+nzuUZl1mRedtLKpFrRtsGBQTZdu6klz5+lj//IiNgyciUiHpG0ML+QzJqrXev1W+cqumhblsQ/SdK+EfEIgKT9Mj7OrBTatV6/da6ii7ZlSeCXAislfZVk1awXAxflGpVZE7VrvX7rbEUWbctycvdqSauBkZWz/iEi7s43LLPmadd6/WZ5qbUC1z4R8WjatfMgcO2o2/aLiD+0IkCzZmjXev1meajV4r8WOBVYw87zz5B09wTwtCxPIGkysBq4PyJOlTQPuA54SrrvV0TEYw3EbmZmDRh3OGeapAUcHxFPG/UzLyIyJf3U2Jm+lwCXRcTTgUeA1zQUuZmZNaTmOP6ICODbje5c0kHAi4Ar0+siOVdwfXqXZcAZje7fzMzql6Vkw+2S/nuD+78cOBcYTq8/BdgSEUPp9ftISkHsQtISSaslrd68eXODT29mZmNlSfzPBX4q6deS1kq6U9LaiR4k6VTgoYhY00hgEbE0Ivoiom/mzJmN7MLMzKrIMo7/bxvc9zHAaZJeCEwF9gE+CcyQ1JO2+g8iKfxm1lRDlSGGKkNUVlY8ksdsjHFb/JIOkHQ58GmSCp2PRMSGkZ+JdhwR50fEQRExFzgT+H5EvBy4lWQSGMDZwI27+RpKx3VhijVSAGv7vdvpP7afVQtXMTgwWHRYZqVRq8V/Nclwy38nGdb5b8A5TXjO9wLXSfpXoB+4qgn7LA3XhSletQJYnbZYttnuqJX4Z0fE+9LL35N0e6NPEhErgBXp5fXA0Y3uq4wqKys7am64LkzxqhXAcneP2U41+/gl7UsyYQtg8ujrnrmbGNvCn3ro1B23uS5MMX3tRRfAMiu7Wom/l6SrR6O2jbT6M8/c7XRjW/gxFExfMJ1pc6dx8LkHd3XSGW+xCSD3rq8iC2CZld24iT89KWsTqFb50QknUa2vvae3Z8eJVnd9mRUjyzh+q2Gk8uOUeVOYv3y+k/4oI33tTN75pbhwxUL2WrBX0aGZdTUvqNIErvxYnfvazcrJid9y5b52s/LJstj6pZL+qhXBmJlZ/rL08a8Dlkq6TdLrJbn5ZrtlqDLEtt9u88zmAvk96G5Zll68ErhS0uHAq4C1kn4MfC4ibs07QOss4w3x9Mzm1vF7YJlG9aSraD0z/XkYuAN4p6TrcoytlFyHZ/dUG+I5ODDIpms3FRlWV/F7YBO2+CVdBvwdsBz4SET8V3rTJZJ+kWdwu2t0KYVmnGAcrw7P4MCghyhmVK2cQjfPbC6C3wPLMqpnLfA/I2JrldtKW3OnsrLCwHEDxFDziqVVq8PT09vjAmB18BDP4vk9sCyJ/w7g8GTVxB0qwIaIKG1/x5YVW4gnmlsszbN0m6MMQzy7vV5/Gd4DK06WxP8Z4CiSlr+AZwE/A3olvSEibs4xvoblcTg7Mkt3qDLkpN/GiqwhZFYGWU7uPgAsTJdBfA6wEFgPnAx8LM/gdsfI4ey8D89raimFnt4eps6Z6qTfxqqd3AR8gtO6RpYW/zMi4mcjVyLibknPjIj1Y7p/SseHs1bNePX6+xf3Fx2adbFWdj9mSfx3S7oCGBm6+dJ02xTg8dwiM8uJT262RrNH1XWyVnc/Zkn8ZwNvBN6eXv8x8G6SpH9CU6Mxy2h3k4qPBvPlJUjr0+oS5hOtwDUZ+E5EnABcWuUu465gLWkq8ANgSvo810fEByT9J3A8ycgggHMiYqD+0K1bVUsq2+7d5rkUJeIlSOvT6u7Hmok/Ip6QNCypt4Ghm9uBEyNiUNIewI8kfTe97T0RcX0jAbejbh862GzVkornUpRLtaHPniQ2vlZ3P2bp6hkE7pR0C7BjEldEvLXWgyIi2HlEsEf6E+M/In9FJGDXRWk+z6dovmZ/Njz0uX6t7H7Mkvi/nv7ULe0qWgM8Hfh0RNwm6Q3ARZLeT1IG4ryI2N7I/utR1Njtan13PuR9snqTjpNKc+XVOPECReWVpTrnMknTgDkRUVdtnoh4AlggaQZwg6RnAecDDwJ7AkuB9wIfGvtYSUuAJQBz5syp52mrKmr9V9dFqa1a0gEmXKzeSaV53DjpPlkWYvk7YAC4Kb2+QNI363mSiNgC3AqcEhEbI7Ed+ALj1PuJiKXppLG+mTNn1vN0VRW1/mteE8k6xdik89imx9g6sJWHv/Ew/cf2s2rhKvoX9/PA0geKDLOjVfts+ER5OeS1bkKWrp4LSZLzCoCIGJD0tIkeJGkm8HhEbEmPGE4mqeg5OyI2Kpn9dQZwV4Ox16XIsdseOji+sUdE+5++Pxv/Y2PyZeDWZ0t4XkM5jT4avuOkO5pbgSDDfR6PiMqYWbrDGR43G1iW9vNPAr4SEd+S9P30S0EkRxKvrzPmhjkBl8/YpAOwadkmd421mD8b5TP6aHj4sWG2rNjS0sT/M0lnAZMlHQa8FfjJRA+KiLUkdX3Gbj+x7iito41NOm59mu16NDzSMGqGLIn/LcD7SMblfwn4HvDhpkVQcp523npufZrl2wWXZVTPn0gS//ua9qxtwjNEzaxIeTWCsiy9+AyS2jxzR9+/U7psao0h9wxRM+tEWbp6vgp8FrgSeCLfcFprookrniFaDi550bnclVqMLIl/KCKuyD2SAkw0ccUzRIs33pezF7jPVx5ftmP36QqexcmS+P+vpDcCN5Cc4AUgIv6QW1QtkmVWrWeIFmu8GdfucstPHuVNqu1zqDK0S1fqn3/5Z7b9dhvTnz3dn7kcZa3HD/CeUdsCmHASV9l54kr5jVeu1vLTaHmTWt024+1z+6SkLTlpyiQOOe8Q7nnDPWy/d7sLGuYsy6ieea0IpCgeOlhu/nJuvfG+bFctXLWjfMDY92GiEXDj7XP0l4VrBrXOuIlf0rkR8bH08j9FxFdH3faRiLigFQGa+cu5tap92U40EGKiEXDjfYGPfW9d0LA1arX4zwQ+ll4+n2R0z4hTgI5M/B5BYrZrQp6oNZ5lBNxEX+A+umudWolf41yudr0jVGvVeMKW2cQDIZo1As5Hd61RK/HHOJerXe8I1Vo1Hj1iVr01Pvbo2CPg2ketxD9f0qMkrftp6WXS61Nzj6wAHkFiNr7RrXEfHbe3cRN/RExuZSBl4D5Gs2x8dNzesozj7yruY+xePrGfnY+O25sTvxn5LTjeqXx03N6c+M3wguON8NFx+5pwsXWzbuAFx62buMVvhrsurLvklvglTQV+AExJn+f6iPiApHnAdcBTgDXAKyLisbziMMvKXRfWLfLs6tkOnBgR84EFwCmS/hq4BLgsIp4OPAK8JscYSmmoMrSj2FWeKisrbLh4Q+7PY2btJbcWf0QEMJhe3SP9CeBE4Kx0+zLgQqAjF3qpplWjR7xecDl5xSkrg1z7+CVNJunOeTrwaeDXwJaIGErvch/w1HEeuwRYAjBnzpw8w2ypVo0e8XrBu6/Z4/r9ZVw/z63IR66JPyKeABZImkGygtcz63jsUmApQF9fX8fUBsqy6leznsfrBTcujyMzfxnXx3Mr8tOSUT0RsUXSrcAiYIaknrTVfxBwfytiKItWjR7xesG7J48jM38Z18dzK/KT28ldSTPTlj6SpgEnA+uAW4EXp3c7G7gxrxjKqndRL4ecf0juH/qe3h6mzpnq5NKAPMb1j3wZT5k3hfnL5/t9mUCr5lZ04yCIPFv8s4FlaT//JOArEfEtSXcD10n6V6AfuCrHGMwakteRmUsXZ9eKo+Nq511G1hfu5HMveY7qWQssrLJ9PXB0Xs9r1iwe11+8vN+Dauddenp7Ov7ci2fumlnX6tbzLk78Zta1unUQhIu0mVnbyGPWezcOgnCL38zawnjj+gGP7a+TW/xm1haqjesHGBwYZNO1mwqLqx058ZtZW6g2rn/hioUdPewyL+7qKQHXIzGbmNdMaB4n/oK539IsO8+taA539RTM/ZZmu6dV61t0Erf4C1atWmfvol76F/cXHZqN4S658nEFz8Y48RfM/ZbtwQmmnFzBszFO/CXgfsvyc4Ipp1atb9FpnPjNMnCCKScfMTfGid8sAyeY8hp7xOxzMRNz4jfLqJEuOS+u3lo+F5ONE39JudVSfhO9R15cvfV8LiYbJ/4Scqul/LJMvPPi6q2X5VyMj8Kc+EvJrZbyq/YejSzZB8n71K2LfBRponMxPgpLOPGXkEeQlF+WiXfdushH0Wqdi/FRWCK3xC/pYOBqYBZJ22hpRHxS0oXAa4HN6V0viIjv5BVHO2pkBIkPX1sr63vkxdXLxUdhiTxb/EPAuyLidkl7A2sk3ZLedllEfDzH52579YwgqXb4OtLt0G2HsK3kiXftx0dhidyKtEXExoi4Pb38R2Ad8NS8nq+bVTt8BbryELZoLhhWft241OJYLenjlzQXWAjcBhwDvFnSK4HVJEcFj1R5zBJgCcCcOXNaEWap1Ro66MPXcqg20qcbTxxa+eVellnSXsDXgLdHxKPAFcChwAJgI3BptcdFxNKI6IuIvpkzZ+YdZqmNJJTt926n/9h+Vi1cRf/ifh5Y+gCw8/B1yrwpzF8+30m/INVG+vioy8oo1xa/pD1Ikv41EfF1gIjYNOr2zwHfyjOGTpBleKdPIhZvvJE+ZmWTW4tfkoCrgHUR8YlR22ePutvfA3flFUOnqLbWqLsPymdkpM+8D8/zkVeJ+LzLrvJs8R8DvAK4U9JAuu0C4GWSFpC0YX8DvC7HGDqCC4S1D4/0KRefd6kut8QfET8CVOUmj9lvgCsQmtXP512q88zdNuRWjFk2Pu9SnRN/G3Irxiwbd5NW58TfhtyKMcvO51125cTfhtyKMbPd4cTfptyKMbNG5T5z18zMysWJ38ysyzjxm5l1GSd+M7Mu48RvZtZlnPjNzLqMImLiexVM0mZgQ4MP3x94uInhFMmvpXw65XWAX0tZ7c5rOSQidlnQpC0S/+6QtDoi+oqOoxn8WsqnU14H+LWUVR6vxV09ZmZdxonfzKzLdEPiX1p0AE3k11I+nfI6wK+lrJr+Wjq+j9/MzJ6sG1r8ZmY2ihO/mVmX6ejEL+kUSb+Q9CtJ5xUdTz0kfV7SQ5LuGrXtf0v6uaS1km6QNKPAEDORdLCkWyXdLelnkt425vZ3SQpJ+xcVY1aSpkr6L0l3pK/lg+l2SbpI0j2S1kl6a9GxZiFpsqR+Sd9Kr58k6XZJA5J+JOnpRceYhaQZkq5PPxvrJC2StJ+kWyT9Mv29b9FxTkTS4enffuTnUUlvz+Nz37F9/JImA/cAJwP3AauAl0XE3YUGlpGk44BB4OqIeFa67fnA9yNiSNIlABHx3gLDnJCk2cDsiLhd0t7AGuCMiLhb0sHAlcAzgedERKkn3EgSMD0iBiXtAfwIeBvwl8AJwDkRMSzpgIh4qMhYs5D0TqAP2CciTpV0D3B6RKyT9Ebg6Ig4p9AgM5C0DPhhRFwpaU/gL4ALgD9ExEfTRt++Zf+sjJbmr/uB5wKH0+TPfSe3+I8GfhUR6yPiMeA64PSCY8osIn4A/GHMtpsjYii9+lPgoJYHVqeI2BgRt6eX/wisA56a3nwZcC47VxAutUgMplf3SH8CeAPwoYgYTu/XDkn/IOBFJF+8IwLYJ73cCzzQ6rjqJakXOA64CiAiHouILSSf9WXp3ZYBZxQR3244Cfh1RGzI43PfyYn/qcDvRl2/j50JpxO8Gvhu0UHUQ9JcYCFwm6TTgfsj4o5io6pP2j0yADwE3BIRtwGHAi+VtFrSdyUdVmiQ2VxO8qU7PGrbvwDfkXQf8ArgowXEVa95wGbgC2m31ZWSpgOzImJjep8HgVmFRdiYM4EvVdnelM99Jyf+jiXpfcAQcE3RsWQlaS/ga8DbSWK/AHh/kTE1IiKeiIgFJK2uoyU9C5gCbEun1X8O+HyBIU5I0qnAQxGxZsxN7wBeGBEHAV8APtHy4OrXAxwFXBERC4GtwJPO50XSn90WR5UAaXfVacBXx2xv2ue+kxP//cDBo64flG5ra5LOAU4FXh5tcoIm7Q//GnBNRHydpIU8D7hD0m9I3pvbJf234qKsT9qdcCtwCsnR5NfTm24AjiworKyOAU5L//bXASdK+jYwPz2CAfgy8LyC4qvHfcB9o+K+nuSLYFN6fmnkPFPpu99GeQFwe0RsGtnQ7M99Jyf+VcBhkual36BnAt8sOKbdIukUksPz0yLiT0XHk0V6QvQqYF1EfAIgIu6MiAMiYm5EzCX58B4VEQ8WGOqEJM0cGVEhaRrJwIGfA98gObkLcDzJoILSiojzI+Kg9G9/JvB9kj7xXknPSO92Msn5mFJL/2d+J+nwdNNJwN0kn/Wz021nAzcWEF6jXsaobp48Pvc9zdhJGaVnwN8MfA+YDHw+In5WcFiZSfoSsBjYP+1z/QBwPkm3wi1JPuWnEfH6woLM5hiS/uI7075xgAsi4jvFhdSw2cCydMTFJOArEfEtST8CrpH0DpKRWP9SZJCNSD8vrwW+JmkYeISkP7kdvIXk778nsB54Fen7I+k1JCXdX1JgfJml5ydOBl43avOnaPLnvmOHc5qZWXWd3NVjZmZVOPGbmXUZJ34zsy7jxG9m1mWc+M3MuowTv3U8SU8ZVfHwQUn3p5cHJX0mp+d8u6RXppdXSJpwsex0nsBNecRjNlrHjuM3GxERvwcWAEi6EBiMiI/n9XySekjGwB9Vz2MiYrOkjZKOiYgf5xWfmVv81rUkLR5Vi/5CScsk/VDSBkn/IOljku6UdFNadgJJz5H0/yStkfS9kbIAY5xIMuV+aNS2f1JSy/8eSX+T7uscSd+U9H1geXq/bwAvz+1Fm+HEbzbaoSRJ+zTgi8CtEfFs4M/Ai9Lk/+/AiyPiOSTF2C6qsp9jSNYdGK0nIo4mKVL3gVHbj0r3d3x6fTXwN815OWbVuavHbKfvRsTjku4kKfMx0t9+JzCXZEGMZ7Fz6vxkYGOV/cxm1zo3I0Xc1qT7GnFLRIxed+Eh4MDGX4LZxJz4zXbaDpCuovX4qCqIwySfFQE/i4hFE+znz8DUavsGnuDJn7utY+43NX28WW7c1WOW3S+AmZIWQVJuWtJfVbnfOqDR9WqfAdw14b3MdoMTv1lG6RKeLwYukXQHMED1mvXfJVkOsBEnAN9u8LFmmbg6p1kOJN0AnBsRv6zzcT8gWfD8kXwiM3PiN8tFujDIrIj4QR2PmQkcExHfyC0wM5z4zcy6jvv4zcy6jBO/mVmXceI3M+syTvxmZl3Gid/MrMv8f/oTPCwXV1KZAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# determine number of hours to plot\n", "nT = 3*24\n", "\n", "# determine hour we should start counting\n", "# 0 means start counting at the first hour, i.e., midnight on January 1, 2015\n", "t = 0 + np.arange(nT+1)\n", "price_data = data[\"Price\"][t]\n", "\n", "# Make plot.\n", "plt.figure()\n", "plt.step(t,price_data,'m.-',where='post')\n", "plt.xlabel('Time (hr)')\n", "plt.ylabel('Energy Price ($/MWh)')\n", "plt.xticks(range(0,nT+1,12))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create Mathematical Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For simplicity, assume a battery energy storage system is *price-taker*, i.e., they are subject to the market price but their actions do not influence the market price. During each hour $t$, the battery operator decided to either charge (buy energy) or discharge (sell energy) at rates $c_t$ and $d_t$ (units: MW), respectively, subject to the market price $p_t$ (units: \\$/MWh)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Assume the battery has a round trip efficiency of $\\eta = 88\\%$. Let $E_t$ represent the state-of-charge at time $t$ (units: MW). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](https://ndcbe.github.io/data-and-computing/_images/battery2.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Class Activity

\n", " Using the picture above, sketch the battery system. Label $d_t$, $c_t$, $E_t$, and $p_t$ on your sketch. Verify the units are consistent.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now write a constrained optimization problem to compute the optimal market participation strategy (when to buy and sell)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{align*}\n", " \\max_{E,d,c} \\quad & \\psi := \\sum_{t = 1}^{N} p_{t} (d_{t} - c_{t}) \\\\\n", "\\mathrm{s.t.} \\quad & E_{t} = E_{t-1} + c_{t} \\sqrt{\\eta} - \\frac{d_{t}}{\\sqrt{\\eta}}, ~~ \\forall ~ t \\in \\{1,...,N\\} \\\\\n", " & 0 \\leq c_{t} \\leq c_{max}, ~~ 0 \\leq d_{t} \\leq d_{max}, \\nonumber \\\\\n", " & 0 \\leq E_{t} \\leq E_{max}, ~~\\forall ~ t \\in \\{1,...,N\\}\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Class Activity

\n", " Discuss the optimization problem above with a partner. Then, write a few sentences to explain each equation.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define Pyomo Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now write the optimization problem in Pyomo. We will create a set `TIME` to write the model compactly. For a review on creating concrete models, refer to the previous notebook (01-Pyomo-Basics)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Class Activity

\n", " Finish the Pyomo model below.\n", "
" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "tags": [] }, "outputs": [], "source": [ "# define a function to build model\n", "def build_model(price,e0 = 0):\n", " '''\n", " Create optimization model for battery operation\n", "\n", " Inputs:\n", " price: Pandas DataFrame with energy price timeseries\n", " e0: initial value for energy storage level\n", " \n", " Output:\n", " model: Pyomo optimization model\n", " '''\n", " \n", " # Create a concrete Pyomo model.\n", " model = ConcreteModel()\n", "\n", " ## Define Sets\n", "\n", " # Number of timesteps in planning horizon\n", " model.TIME = Set(initialize = price.index)\n", "\n", " ## Define Parameters\n", "\n", " # Square root of round trip efficiency\n", " model.sqrteta = Param(initialize = sqrt(0.88))\n", "\n", " # Energy in battery at t=0\n", " model.E0 = Param(initialize = e0, mutable=True)\n", "\n", " # Charging rate [MW]\n", " model.c = Var(model.TIME, initialize = 0.0, bounds=(0, 1))\n", "\n", " # Add your solution here\n", " \n", " ## Define constraints\n", " \n", " # Define Energy Balance constraints. [MWh] = [MW]*[1 hr]\n", " # Note: this model assumes 1-hour timestep in price data and control actions.\n", " def EnergyBalance(model,t):\n", " # First timestep\n", " if t == 0 :\n", " return model.E[t] == model.E0 + model.c[t]*model.sqrteta-model.d[t]/model.sqrteta \n", " \n", " # Subsequent timesteps\n", " else :\n", " # Add your solution here\n", " \n", " model.EnergyBalance_Con = Constraint(model.TIME, rule = EnergyBalance)\n", " \n", " ## Define the objective function (profit)\n", " # Receding horizon\n", " def objfun(model):\n", " return sum((-model.c[t] + model.d[t]) * price[t] for t in model.TIME)\n", " model.OBJ = Objective(rule = objfun, sense = maximize)\n", " \n", " return model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solve Optimization Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now create an instance of our Pyomo model. Notice the function `build_model` requires we pass in a Pandas DataFrame with the price data. Let's try the first day only." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ipopt 3.11.1: \n", "\n", "******************************************************************************\n", "This program contains Ipopt, a library for large-scale nonlinear optimization.\n", " Ipopt is released as open source code under the Eclipse Public License (EPL).\n", " For more information visit http://projects.coin-or.org/Ipopt\n", "******************************************************************************\n", "\n", "NOTE: You are using Ipopt by default with the MUMPS linear solver.\n", " Other linear solvers might be more efficient (see Ipopt documentation).\n", "\n", "\n", "This is Ipopt version 3.11.1, running with linear solver mumps.\n", "\n", "Number of nonzeros in equality constraint Jacobian...: 95\n", "Number of nonzeros in inequality constraint Jacobian.: 0\n", "Number of nonzeros in Lagrangian Hessian.............: 0\n", "\n", "Total number of variables............................: 72\n", " variables with only lower bounds: 0\n", " variables with lower and upper bounds: 72\n", " variables with only upper bounds: 0\n", "Total number of equality constraints.................: 24\n", "Total number of inequality constraints...............: 0\n", " inequality constraints with only lower bounds: 0\n", " inequality constraints with lower and upper bounds: 0\n", " inequality constraints with only upper bounds: 0\n", "\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 0 2.1094237e-015 1.13e-002 2.99e+001 -1.0 0.00e+000 - 0.00e+000 0.00e+000 0\n", " 1 -1.4124684e-001 1.07e-002 2.72e+001 -1.0 2.03e-001 - 8.78e-002 4.88e-002f 1\n", " 2 -2.6373063e+000 8.29e-003 2.74e+001 -1.0 8.21e-001 - 5.78e-002 2.27e-001f 1\n", " 3 -2.2360880e+001 6.20e-003 2.78e+001 -1.0 3.20e+000 - 3.56e-002 2.52e-001f 1\n", " 4 -3.8049600e+001 5.73e-003 2.51e+001 -1.0 8.23e+000 - 9.39e-002 7.63e-002f 1\n", " 5 -5.9077018e+001 4.87e-003 2.34e+001 -1.0 6.97e+000 - 8.18e-002 1.50e-001f 1\n", " 6 -7.9367583e+001 3.64e-003 2.18e+001 -1.0 4.22e+000 - 9.20e-002 2.52e-001f 1\n", " 7 -8.5414120e+001 3.02e-003 1.71e+001 -1.0 2.06e+000 - 2.10e-001 1.72e-001f 1\n", " 8 -8.9903789e+001 1.64e-003 1.17e+001 -1.0 1.23e+000 - 3.29e-001 4.57e-001f 1\n", " 9 -9.2816651e+001 9.05e-004 3.65e+000 -1.0 1.25e+000 - 6.67e-001 4.48e-001f 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 10 -9.3848256e+001 4.09e-004 1.17e+000 -1.0 9.08e-001 - 6.61e-001 5.48e-001f 1\n", " 11 -9.3998203e+001 3.54e-005 6.77e-002 -1.0 7.33e-001 - 1.00e+000 9.13e-001f 1\n", " 12 -9.8037848e+001 4.44e-016 2.55e-003 -1.7 1.37e-001 - 9.88e-001 1.00e+000f 1\n", " 13 -9.9097566e+001 4.44e-016 9.77e-015 -2.5 5.13e-002 - 1.00e+000 1.00e+000f 1\n", " 14 -9.9251854e+001 4.44e-016 2.64e-003 -3.8 2.84e-002 - 9.49e-001 1.00e+000f 1\n", " 15 -9.9252231e+001 4.44e-016 7.61e-015 -3.8 6.38e-004 - 1.00e+000 1.00e+000f 1\n", " 16 -9.9260699e+001 4.44e-016 1.02e-014 -5.7 6.12e-004 - 1.00e+000 1.00e+000f 1\n", " 17 -9.9260804e+001 4.44e-016 1.55e-005 -8.6 8.42e-006 - 1.00e+000 9.97e-001f 1\n", " 18 -9.9260804e+001 8.88e-016 2.84e-001 -8.6 2.56e-008 - 8.71e-001 1.00e+000f 1\n", " 19 -9.9260804e+001 8.88e-016 7.58e-015 -8.6 2.64e-009 - 1.00e+000 1.00e+000h 1\n", "\n", "Number of Iterations....: 19\n", "\n", " (scaled) (unscaled)\n", "Objective...............: -9.9260804110396435e+001 -9.9260804110396435e+001\n", "Dual infeasibility......: 7.5766950348697468e-015 7.5766950348697468e-015\n", "Constraint violation....: 8.8817841970012523e-016 8.8817841970012523e-016\n", "Complementarity.........: 3.0241425605247924e-009 3.0241425605247924e-009\n", "Overall NLP error.......: 3.0241425605247924e-009 3.0241425605247924e-009\n", "\n", "\n", "Number of objective function evaluations = 20\n", "Number of objective gradient evaluations = 20\n", "Number of equality constraint evaluations = 20\n", "Number of inequality constraint evaluations = 0\n", "Number of equality constraint Jacobian evaluations = 20\n", "Number of inequality constraint Jacobian evaluations = 0\n", "Number of Lagrangian Hessian evaluations = 19\n", "Total CPU secs in IPOPT (w/o function evaluations) = 0.056\n", "Total CPU secs in NLP function evaluations = 0.000\n", "\n", "EXIT: Optimal Solution Found.\n" ] } ], "source": [ "# Build the model\n", "instance = build_model(data[\"Price\"][0:24],0.0)\n", "\n", "# Specify the solver\n", "solver = SolverFactory('ipopt')\n", "\n", "# Solve the model\n", "results = solver.solve(instance, tee = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Extract Solution from Pyomo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Excellent. Ipopt terminated with the message `Optimal Solution Found`. Let's inspect the answer." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Class Activity

\n", " Run the code below.\n", "
" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Declare empty lists\n", "c_control = []\n", "d_control = []\n", "E_control = []\n", "time = []\n", "\n", "# Loop over elements of TIME set.\n", "for i in instance.TIME: \n", " # Record the time\n", " time.append(value(i))\n", " \n", " # Use value( ) function to extract the solution for each variable and append to the results lists\n", " c_control.append(value(instance.c[i]))\n", " \n", " # Adding negative sign to discharge for plotting\n", " d_control.append(-value(instance.d[i]))\n", " E_control.append(value(instance.E[i]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot Results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's plot the optimal charge and discharge profile." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAq+ElEQVR4nO3debRcVZn38e8vE4PMJEDIyCQQAgRyRWIEqwDbGG1wAFtbVGgxjaKCaLOM3a/SuFw0bzugDciLgiIimsYJ9JaASS6DhEiCYchACAFMAOESICEEQobn/WOfMkVR91bde+vUOafO81mrVk2nqp5b91Q9dfZ+9t4yM5xzzuXXoKQDcM45lyxPBM45l3OeCJxzLuc8ETjnXM55InDOuZwbknQAfTV8+HAbP3580mE451ymLFy48DkzG1HrvswlgvHjx7NgwYKkw3DOuUyR9ERP93nTkHPO5ZwnAuecyzlPBM45l3OeCJxzLuc8ETjnXM7FnggkDZb0F0m/q3HfdpJ+IWmFpPmSxscdj3POuddrxRHBucDSHu77JPCCmR0IfAe4pAXxONe25s2Diy8O52l7jEuvWMcRSBoNvAf4BnB+jU1OAS6MLt8IXCZJ5nNjO9dnd98NhQJs3gyDB8PHPgajR/f+mNWr4brrYMuWvj9m61bYbjuYPRumTGnan+ESEPeAskuBC4Cde7h/FLAKwMw2S1oL7Ak8V7mRpBnADICxY8fGFatzmTZrFmzaFC5v3gw/+hFIvT+m8idXfx7z2mvQ1eWJIOtiaxqS9F7gWTNbONDnMrOrzKzDzDpGjKg5Qtq53CsngcGDYYcdwhHC1q29n+6+O2zb18dst9221yoUEvuTXZPEeUQwFThZ0nRge2AXST81s9MrtnkSGAOsljQE2BVYE2NMzrWtJUtg//3hrLPCl3Mjv9KnTAlNO11dfX/M9OlwxBF+NNAO1IrmeEkF4Etm9t6q288BDjezsyV9GPiAmX2ot+fq6Ogwn2vIuddbtw723BPOPx8uaVHJxfnnw+WXQ3c37LJLa17T9Z+khWbWUeu+lo8jkHSRpJOjq1cDe0paQehM/nKr43GuHcyeHdr43/3u1r3maaeFPoKbbmrda7p4tGT2UTPrArqiy1+tuP1V4LRWxOBcOyuVYOedYerU1r3mW98KY8bA//4vnH56/e1devnIYucyziwkgne+E4YObd3rDhoEp54Kf/gDrF3butd1zeeJwLmMW7w41Pa3slmo7EMfCs1DN9/c+td2zeOJwLmMK5XC+bRprX/tcvPQrFmtf23XPJ4InMu4UgkOP7z+iOA4SKHT+JZbvHkoyzwROJdhL70Ed92VTLNQmVcPZZ8nAucybPbsMKI4yUTw1rfC2LGheshlkycC5zIsibLRalKoHvLmoezyROBcRplBZyecdFJry0Zr8eahbPNE4FxGJVk2Ws2bh7LNE4FzGVUuG01DIvDmoWzzROBcRpVKMHFiMmWjtZQHl3nzUPZ4InAug9JQNlrtmGNC85APLsseTwTOZVAaykarlQeX3XorvPhi0tG4vvBE4FwGlUqw007Jlo3W4tVD2eSJwLmMqZxtdNiwpKN5vXLzkFcPZYsnAucyZskSWLUqXc1CZZVzD3nzUHbEuXj99pL+LOl+SYsl/WeNbc6Q1C1pUXQ6K654nGsXaSobreVDHwr9F948lB1xHhFsBE4wsyOBScA0ScfW2O4XZjYpOv0wxnicawtpKxut9pa3wLhxXj2UJbElAgvWR1eHRieL6/Wcy4OXXoI770zv0QB49VAWxdpHIGmwpEXAs8BtZja/xmYflPSApBsljenheWZIWiBpQXd3d5whO5dqc+akr2y0ltNOC3H+9rdJR+IaEWsiMLMtZjYJGA0cI2li1SY3A+PN7AjgNuDaHp7nKjPrMLOOESNGxBmyc6mW1rLRauXmIa8eyoaWVA2Z2YvAXGBa1e1rzGxjdPWHwORWxONcFpXLRk86KX1lo9W8eShb4qwaGiFpt+jyDsA7gWVV24ysuHoysDSueJzLuqVL4a9/TX+zUFm5esibh9IvziOCkcBcSQ8A9xL6CH4n6SJJJ0fbfD4qLb0f+DxwRozxOJdpnZ3hPCuJoKMDxo/36qEsGBLXE5vZA8BRNW7/asXlmcDMuGJwrp2USnDYYTCmZklF+pSbhy69FF54AXbfPemIXE98ZLFzGZCFstFavHooGzwROJcBWSkbrVZuHvLqoXTzROBcBpTLRt/+9qQj6Zty89Btt4XmIZdOngicS7kslY3W4tVD6eeJwLmUy1rZaLXJk715KO08ETiXcmmfbbQeKRwVePNQenkicC7lslY2WotXD6WbJwLnUmz9+myWjVabPBlGjoSLL4Z585KOxlXzROBcis2ZE9YAznoiuOce6O6G5cvhxBM9GaSNJwLnUiyrZaPVurpg69ZweePGcN2lhycC51KqXDZ64onZLButVCjAdtuFy4MGhesuPTwROJdSy5bBE09kv1kIYMoUmD0bDjgglJJOmZJ0RK6SJwLnUirrZaPVpkyBM8+ERx+F559POhpXyROBcynV2QkTJsDYsUlH0jyFQmjyuuOOpCNxlTwROJdCs2eHDtUjj0w6kuZ6y1tgxx1h7tykI3GVPBE4lzLz5sH06bBlC/zqV+1VajlsWFhv2auG0iXOpSq3l/RnSfdHq5D9Z41ttpP0C0krJM2XND6ueJzLiq6uMAoXYPPm9vvSLBTggQdgzZqkI3FlcR4RbAROMLMjgUnANEnHVm3zSeAFMzsQ+A5wSYzxOJcJ73jHtsvDhrVfqWX577n99kTDcBViSwQWrI+uDo1OVrXZKcC10eUbgRMlKa6YnMuC3XcPHarve1/oK2i3UstyP0G7HelkWax9BJIGS1oEPEtYvH5+1SajgFUAZrYZWAvsWeN5ZkhaIGlBd3d3nCE7l7hy2eh3v9t+SQBg6NAwUtoTQXrEmgjMbIuZTQJGA8dImtjP57nKzDrMrGPEiBFNjdG5tCmV2q9stFqhAA8+GOYfcslrSdWQmb0IzAWmVd31JDAGQNIQYFfAu5Bcbq1fH2rs22UQWU/K/QQ+niAd4qwaGiFpt+jyDsA7gWVVm90EfCK6fCowx8yq+xGcy425c9tjttF6OjrgTW/y5qG0GBLjc48ErpU0mJBwZpnZ7yRdBCwws5uAq4HrJK0Angc+HGM8zqVeqRS+ILM+22g95X4CH1iWDrElAjN7ADiqxu1frbj8KnBaXDE4lyWVs42WZ+psZ4UCzJwJzz4Le+2VdDT55iOLnUuJhx+Gxx9v/2ahsmIxnHs/QfI8ETiXEp2d4TwvieDoo8OiO948lDxPBM6lRKkEhx4K48YlHUlr+HiC9PBE4FwKlMtGp09POpLWKhRgyZLQT+CS44nAuRTIS9lotXI/gR8VJMsTgXMpkJey0WrlfgJPBMnyROBcwvJWNlppyBA47jhPBEnzROBcwvJWNlqtWISlS+GZZ5KOJL88ETiXsHZbpL6vyvMO+VFBcjwROJewvJWNVjvqKNh5Z08ESfJE4FyCXn45rNSV16MBCP0Exx/viSBJngicS1Bey0arFQqwbBk8/XTSkeSTJwLnElQuGz3uuKQjSZavY5wsTwTOJcQszC90wgn5KxutNmkS7LKLNw8lxROBcwnJe9lopXI/gU9Al4w4VygbI2mupCWSFks6t8Y2BUlrJS2KTl+t9VzOtaO8l41WKxRg+XJ46qmkI8mfOI8INgNfNLMJwLHAOZIm1NjuTjObFJ0uijEe51KlXDY6fnzSkaSD9xMkp24ikDRa0pck/VbSvZLukHSFpPdI6vHxZva0md0XXX4JWAqMal7ozmWXl42+0aRJsOuu3jyUhF4TgaQfAdcArwGXAB8BPgP8EZgG3CXp+HovImk8YdnK+TXuniLpfkklSYf18PgZkhZIWtDd3V3v5ZxLPS8bfaPBg308QVLqrVn8LTN7qMbtDwG/kjQMGNvbE0jaCfglcJ6Zrau6+z5gnJmtlzQd+A1wUPVzmNlVwFUAHR0dVidm51LPy0ZrKxTg5pvhySdhlLcftEyvRwQ9JIHK+18zsxU93S9pKCEJXG9mv6rx+HVmtj663AkMlTS8ocidy6jybKNeNvpGPu9QMhrqLJY0VdJtkpZLWinpMUkr6zxGwNXAUjP7dg/b7BNth6RjonjW9O1PcC5bli+Hxx7zZqFajjwSdtvNE0Gr1WsaKrsa+AKwENjS4GOmAh8DHpS0KLrtK0RNSWZ2JXAq8GlJm4FXgA+bmTf9uLbmZaM9836CZDSaCNaaWakvT2xmdwGqs81lwGV9eV7nsq5UgkMO8bLRnhQKcNNNsHo1jB6ddDT5UK9q6GhJRwNzJf23pCnl26LbnXN94GWj9fk6xq1Xt2qo6npHxWUDTmhuOM61t7lzYeNGTwS9OeII2H33kAhOPz3paPKh10RgZsVWBeJcHpRKsOOOoR3c1TZokM871Gr1mobWSOqU9O+SipJ2bFVgzrUbLxttXKEAK1fCX/+adCT5UK98dD/gUmAoMBNYFY3w/a6kD8UdnHPtpFw2On160pGkX7mfwOcdao16A8rWmdmtZnahmf0DofTzx8B7gBtaEJ9zbcPLRht3+OHb+glc/HrtI5C0L/C26PSW6OaFwH8A8+INzbn24mWjjRs0CN7xDu8naJV6VUOrCfMBfQf4spm9Fn9IzrWfDRtCM8dnPpN0JNlRLMJvfgNPPAHjxiUdTXur10cwFfgZ8H5gnqRfRlNST5Xk3V3ONcjLRvvO1ydonXp9BPPM7NtmdqqZTQa+CGwErgXWtiJA59qBl4323cSJsMce3jzUCnWnmJB0CNv6CaYCuwH3AFfGGplzbcLLRvun3E/gHcbxq9dZ/BzwFKFj+A7gv3qbdtq11t13w6xZYWWnI45o/HFLl4Z212IRpkyJLTwXmTUr1MS///1JR5I9xSL8+tdwwQXh/fP9NR71jggOMDNvAkqhO+8MH5Itjc4FW2XQoPDrdPZs/3DFad48+NjHwuXLL4cPftDf777Yffdw/s1vwmWX+f4al3qJ4OvRcgE1mdnnmxuOa8Srr8KMGduSwKBB8NGPwqmn1n/sjTfCddfB1q1hqcSuLv9gxamrCzZtCpc3bfL3u6/KI4vNfH+NU71EcDZhWcpZhCaiXqeVdvF7+WV43/tg2TIYOjR8oQ8bBp/+dGMfkBEj4Prrtz2uXJnh4nHsseFc8ve7P4rFsEbBli3+/sWpXvnoSMJawe8iLDIzFPitmV1rZtfGHZx7vbVr4V3vgjlz4NprQ1nd17/et8PlKVO21bL//Of+6ypuGzaE8zPO8GaN/pgyBc49N1y+7jp//+JSr3x0jZldGc1CeiahYmiJpI/Ve2JJYyTNlbRE0mJJ59bYRpK+J2mFpAd8jYOerVkDJ50E8+fDL34BH/94+FDMnNn3D8cZZ4Tzl19uepiuSrls9Ior/Eusv848M5yvW5dsHO2s0TWLjwbOBU4HSoRpJurZDHzRzCYAxwLnSJpQtc27gYOi0wzg+w3GnSt/+1s4JH7wwTDSspG+gN5MmgS77OJleXGrLBvdfvuko8muww6D4cN9f41TvWmoL5K0EDgfuB3oMLNPmtmSek9sZk+b2X3R5ZeApcCoqs1OAX5iwT3AbpJG9ucPaVerVoVa6sceg85OeM97Bv6cvi5sazzySCgb9dHEAyOFH0Jz54bk6pqv3hHBfxCag44ELgbui5pwHpT0QKMvImk8cBQwv+quUcCqiuureWOyQNKMaPrrBd3d3Y2+bOY9+igcd1w4Irj11vDLslmKxTAt8lNPNe853ev5bKPNUyiEH0WPPZZ0JO2pXtXQfgN9AUk7Ab8EzjOzfrXymdlVhE5rOjo6cvGbYOnS0Cfw6quhc3jy5OY+f7n6oqsL/vmfm/vcLiiV4OCDYb8Bf4pc5TrG+++faChtqd4RwV/N7ImeThA6fHt6sKShhCRwvZn9qsYmTwJjKq6Pjm7LtUWLQnPQli2hMqjZSQDgyCNht928eSguGzaE99aPBprj0END6bPPOxSPeolgrqTPSRpbeaOkYZJOkHQt8IlaD4wSxNXAUjP7dg/PfxPw8ah66FhgrZk93ce/oa3Mnx9+/Wy/fRg9PHFiPK9T7ifwD1Y8urp8ttFmKvcTdHV5P0Ec6iWCacAW4AZJT0WloCuBR4CPAJea2Y97eOxUwtiDEyQtik7TJZ0t6exom05gJbAC+AGQ29na580Lo4WLRdhzz5AEDjoo3tcsFGDFCli9Ot7XySOfbbT5CoWwr65cmXQk7afXPgIzexW4ArgiauYZDrxiZi/We2Izu4s6I5HNzIBzGo62Tc2bFxLAxo3hl8+3v92ahTgq53v/6Efjf708KZW2Hdm55qjsJzjggERDaTsNjSMAMLNNUUnoizHGk0vlZgQI8wYtXtya1z3yyDCplzcPNdcjj4SKL28Waq5DDoG99vL9NQ4NJwIXn7dEq0G3ej6aQYN8PEEcOjvDuSeC5vJ+gvh4IkiB8tD5T32q9fPRFArh1+uqVXU3dQ0ql416mWPzFYvw5JNhn3XNU29k8YGSpta4faokb6VrklIJdt01zFff6vlofF3Y5vKy0XiV91dvHmquekcElwK1BoGti+5zA1Sej+ad74QhdRcObb4jjvB+gmbystF4HXww7LOPN2c2W71EsLeZPVh9Y3Tb+FgiypmHHgqHukl9cfi6sM3lZaPx8n6CeNRLBLv1ct8OTYwjt8odi9OmJRdDoRBqs8urQbn+87LR+BUKYY6sRx5JOpL2US8RLJD0qeobJZ1FY1NRuzpKpVDGue++ycVQWZ/t+s/LRlujcp4s1xz1EsF5wJmSuiR9KzrdDnySsD6BG4B16+BPf0r+i2PiRNhjD/9gDZTPNtoab34zjBzp+2sz1RtZ/AzwNklFoDzrze/NbE7skeXAH/8Imzcn/8Xh/QTNUSqFLykvG41X9foEPU976RrVl3EEVnFyTVAqhZXC0rCEYbEY5np/4omkI8mmV17xstFWKhTCOh3LlycdSXuoN45glKT5wIXA/tHpQkl/lvSGBWRc4yrLRocOTToab3cdqK6usHaEJ4LW8H6t5qp3RHAZ8H0ze4eZnR+d3hHdfkX84bWvpMtGqx12WJj11D9Y/VMqwQ47hCY2F78DDwwFFj7+pTnqJYIJtaaZNrOfAIfEElFOlDsWkywbreT9BAPT2ello63k4wmaq14iqHm/pEHA4OaHkx/lstFRKWpgKxbh8cfDyTWuXDY6fXrSkeRLoQDPPAMPP5x0JNlXLxH8TtIPJL2pfEN0+UrCojKuH9atg7vuSk+zUJn3E/SPl40mw/sJmqdeIrgAWAs8IWmhpPuAxwlzDX2ptwdKukbSs5Ie6uH+gqS1FauXfbUf8WfS7NnpKButNmECDB/uH6y+8rLRZBxwQDii9n6Cgas3jmAT8CVJ/wc4MLr5UTPb0MBz/5jQqfyTXra508ze20ig7SRNZaOVBg3y+uy+KpeN/uu/Jh1J/pT7CW67zffXgapXPvoWSfuY2SvRRHNHEdYv/p6kPXp7rJndATzfxFjbglnoWExL2Wi1QiHMOeT9BI3xstFkFYvw7LOwbFnSkWRbvaah/we8BiDpeOC/CL/w1wJXNeH1p0i6X1JJ0mE9bSRphqQFkhZ0d3c34WWTk7ay0WreT9A3XjaaLF+foDnqJYLBZlb+Vf9PwFVm9kszq2wq6q/7gHFmdiTwP8BvetrQzK4ysw4z6xgxYsQAXzZZaSsbrTZhAowY4YmgUT7baLL23x9Gj/b9daDqJgJJ5X6EE4HKOYYGtIyKma0zs/XR5U5gqKThA3nOLCiVwmIwaSobrVQ9j4vr2YoV4ZTWo7s8kEIi9vEEA1MvEdwA3C7pt8ArwJ0QlrAkNA/1m6R9pNC9I+mYKJY1A3nOtEtr2Wi1QiGsYfzYY0lHkm5eNpoOhQJ0d8OSJUlHkl31qoa+IWk2MBK41ezvOXcQ8LneHivpBqAADJe0GvgaMDR63iuBU4FPS9pMSDIfrnj+tpTWstFqlf0EXhLZs1IJDjoolDG65FTur4f12NPoelO3ecfM7qlxW905/8zsI3Xuv4xQXpob5bLRt70t6Uh6d+ihsNdeoXnoX/4l6WjS6ZVXwvszY0bSkbj99oOxY0MiOOecpKPJpr5MQ+0GoDzb6EknpbNstJLP41Lf7beHslGfViJ5lfvr1q1JR5NNnghaZPFiWL06/c1CZYVCiHflyqQjSafOTi8bTZNCAZ57zvsJ+ssTQYukvWy0ms/j0jsvG00X318HxhNBi5RKcPjhoeY5Cw4+GPbe2wfq1OJlo+kzfjyMG+f7a395ImiBctloltqTvZ+gZ142mk6FQui78X6CvvNE0AKzZ8OmTdn74igWw3QYjz6adCTp4mWj6VQowJo1oT/O9Y0nghbIStloNZ/H5Y3KZaNZS+p54PNk9Z8ngphlqWy02pvfDPvs4x+sSldeGcpGx49POhJXbfz4cPIfLn3niSBmWSsbreTzuLzevHlwwQXh8r//e7ju0sX7CfrHE0HMslY2Wq1QgKeeCuvy5l1XV5giBOC11/xIKY2KRXj++TDdu2ucJ4KYZa1stJq3u25zYDTx+qBBMGzYtvfGpUd5gJ83D/WNJ4IYvfRSNmYb7c1BB8HIkZ4IAJ55Jpx/4QuhEixtS426MJZgv/18f+2rAa0p4HqX1bLRSuV+Al/HOBzdHXggfPObSUfielMswq9/HfoJBvlP3Yb42xSjUgl23hmmTk06koEpFODpp+GLX8xvB+krr8CcOdkaFJhXhQK88AKce25+99e+8kQQkyyXjVbbdddwfumlcOKJ+fxwlWcbzfLRXV7svHM4v/zy/O6vfRVbIpB0jaRnJdXsv1fwPUkrJD0g6ei4YknCkiVhla92+OJYsSKcm+W3WqZUChPM+Wyj6bd0aTjP8/7aV3EeEfwY6K1o8t3AQdFpBvD9GGNpuXaaj6ZYhMGDw+W8VsuUZxvdYYekI3H1FAq+v/ZVbInAzO4Anu9lk1OAn1hwD7CbpJFxxdNqpRJMnJjdstFKU6bAv/1buHz11fmrlnn00TCOoh2Seh5MmRIG/AFccUX+9tf+SLKPYBSwquL66ui2N5A0Q9ICSQu6u7tbEtxAvPQS3Hlne31xnHVWOF+7Ntk4ktBOR3d58alPhfMXXkg2jqzIRGexmV1lZh1m1jFixIikw6mrHcpGq+2/fzi6yeNAnXLZaHlAmUu/0aPD/8v7BxqTZCJ4EhhTcX10dFvm/fjHoW2y3E7ZDvK6PsGrr/pso1lVnndoy5akI0m/JBPBTcDHo+qhY4G1ZvZ0gvE0xd13w003hWqFadPaq3StUIBnn4Vly5KOpHVuvz2MIfBEkD2FQmjKvP/+pCNJvzjLR28A5gEHS1ot6ZOSzpZ0drRJJ7ASWAH8APhMXLG00qxZ234xt1vpWnld2Dw1D5XLRr3yJHt8nqzGxTbFhJl9pM79BpwT1+snpTw75eDB7Ve6tt9+MGZM+GB9pi3Sdn2lUvgfetlo9owaFebKmjsXzj8/6WjSLROdxVmydGnoWP3619tvYrK89ROsXAnLl3uzUJYVCnDHHd5PUI8ngiYql41+8IMwc2Z7JYGyYhG6u8PI6XbnZaPZVyzCunWwaFHSkaSbJ4ImmjOn/cpGq+Wp3bWzMyxQf9BBSUfi+svXJ2iMJ4ImKpVgp52yP9tob8aPh7Fj2z8ReNloe9h337D2drvvrwPliaBJKmcbHTYs6Wjik5d1jL1stH0Ui6HJtlzI4d7IE0GTLF0Kf/1rPr44CgV47jlYvDjpSOLjZaPto1AI/QR/+UvSkaSXJ4ImyVPHYh76CcplozvumHQkbqDK/QTtvL8OlCeCJunsDLONjhlTf9usGz8+nNr1g+Vlo+1l5Eg45JD23V+bwRNBE7TjbKP1lMcTbN2adCTNl6eju7woFLyfoDeeCJogD2Wj1QoFWLOmPfsJSiUvG203hUL4wXbffUlHkk6eCJogD2Wj1dq1n+DVV0Niz1NSz4N23V+bxRPBAOWlbLTauHFh7qF2G6hzxx1eNtqO9t4bDj3UE0FPPBEMUJ7KRquV53tvp36CUgm2287LRttRuZ9g06akI0kfTwQDlOeOxWIRnn8eHnoo6Uiax8tG21exCOvXez9BLZ4IBqhUgsMOy0fZaLV2m8dl5Up4+OF8JvU88PEEPfNEMADr1+evbLTS2LFhyu12+WDl+eguD/baCyZMaJ8fLs0UayKQNE3Sw5JWSPpyjfvPkNQtaVF0OivOeJptzpywClmevzjaqZ+gVAqJzctG21ehAHfd5f0E1eJcqnIwcDnwbmAC8BFJE2ps+gszmxSdfhhXPHEol42+/e1JR5KcYhFeeAEeeCDpSAamsmxUSjoaF5diEV5+GRYuTDqSdInziOAYYIWZrTSz14CfA6fE+HotVS4bPfHEfJWNVmuXdlcvG82H448P51nfX5stzkQwClhVcX11dFu1D0p6QNKNkmp2uUqaIWmBpAXd3d1xxNpny5bBE0/A9OlJR5KsMWPCKNysf7DKZaPFYtKRuDjttVco7vB+gtdLurP4ZmC8mR0B3AZcW2sjM7vKzDrMrGPEiBEtDbAnnZ3h3H9Bhi/P22/P9rqwpVI4uvGy0fZXLHo/QbU4E8GTQOUv/NHRbX9nZmvMbGN09YfA5Bjjaao8l41WKxTgxRez20/w2GNeNponhQJs2AALFiQdSXrEmQjuBQ6StJ+kYcCHgZsqN5A0suLqycDSGONpmryXjVbL+jwu5bLRvDfz5UW7jX9phtgSgZltBj4L3EL4gp9lZoslXSTp5Gizz0taLOl+4PPAGXHF00xeNvp6o0aFksusfrC8bDRfhg+Hww/P7g+XOAyJ88nNrBPorLrtqxWXZwIz44whDl42+kaFAsyaFfoJBg9OOprGlctGzzzTy0bzpFCAq68OP+jyXPVXlnRnceZ42WhthQKsXQv33590JH1z552hvdiP7vKl3E9w771JR5IOngj6qFw26l8cr1fuJ8ha85CXjeZTu4x/aRZPBH3k89HUtu++8OY3Z++D5WWj+bTnnnDEEdnbX+PiiaCPSqUwcdXYsUlHkj6FQhihm5XxBI89Fo7wPKnnU6EAf/pT6CfIO08EfbB+ffii8y+O2opFWLcO/vKXpCNpjB/d5VuxGKYV+fOfk44keZ4I+mDu3PDrwevNa8tau+vPfga77w5r1iQdiUvC8ceHSrGs7K9x8kTQB1422ruRI+Hgg7PxwbryytAs8OKLYb3pefOSjsi12h57eD9BmSeCBnnZaGPK68Ju3px0JD3r7ITPfS5cNgtHef5lkE/lfoKNG+tu2tY8ETTo4Yfh8ce9PbmetPcT/OxncMopYSTx9tuHwW/Dhvli9XlVLIZBhXnvJ/BE0CCfbbQxae4nuPxyOP300LR3771hRPHXvw6zZ8OUKUlH55Jw3HHeTwCeCBrmZaON2WcfOOSQdH2wzOCii+Czn4V//Mfwv9xll/DlP3OmJ4E822MPOPLI7A2EbDZPBA3wstG+KRbT00+wdSucdx587WvwiU/AL38ZmoScKysWQ7HAq68mHUlyPBE0oFw26omgMYUCvPQS3HdfsnFs2hS+/L/3PfjCF+Caa2BIrNMsuiwqFLyfwBNBA0oleNObvGy0UWnoJ3jlFfjAB+CnP4VvfAO+9S0Y5Hu7q8H7CTwR1FVZNrrddklHkw177x36U5Jqd33xRXjXu+D3v4fvfx++8hWfYtr1bPfd4aij8t1P4ImgDi8b7Z9CIZl1YZ95JrT53nMP3HADnH12a1/fZVOhkO9+glgTgaRpkh6WtELSl2vcv52kX0T3z5c0Ps54+sPno+mfQiF0sreyn+Dxx0Pz3fLlcPPN8E//1LrXdtlWKIRBZfPnJx1JMmJLBJIGA5cD7wYmAB+RNKFqs08CL5jZgcB3gEviimfePLj44r5PJXDDDbDXXvDUU/HE1a7K/QRf+1rf3vP+/p+uvz5MF/C3v8Ftt4WmIecaVe4nuOii1uyv/Xlcf1+rEXHWUBwDrDCzlQCSfg6cAiyp2OYU4MLo8o3AZZJkZtbMQObNC80FGzeGf/a4cY3NP79hQ/iVKYU+Ah941LhHHw3v2y23wK23Nvaeb9gQFv0x69//CUJpqPcHuL5aujScz5kDU6fGv7/29XHlx0DYx5v9XRRnIhgFrKq4vhp4a0/bmNlmSWuBPYHnKjeSNAOYATC2HyO6urq2tVWbhYnjDjmk/uOWLdv2mPJ8NJ4IGtPVFXZys8bf82XLwrbQv/8ThP+z/59cX7V6f+3r4yofE8t3kZnFcgJOBX5Ycf1jwGVV2zwEjK64/igwvLfnnTx5svXV3Xeb7bCD2eDB4fzuu+N9nOvfe+f/J5eUtO+vzdjHgQXWw/eqrLmtMH8naQpwoZm9K7o+M0o8F1dsc0u0zTxJQ4C/ASOsl6A6OjpswYIFfY5n3ryQRQuFvmXS/j7O9e+98/+TS0ra99eB7uOSFppZR837YkwEQ4DlwInAk8C9wD+b2eKKbc4BDjezsyV9GPiAmX2ot+ftbyJwzrk86y0RxNZHYKHN/7PALcBg4BozWyzpIsIhyk3A1cB1klYAzwMfjise55xztcU684qZdQKdVbd9teLyq8BpccbgnHOudz6y2Dnncs4TgXPO5ZwnAuecyzlPBM45l3OxlY/GRVI38EQ/Hz6cqlHLKePxDUya40tzbODxDUSaY6s0zsxG1Lojc4lgICQt6KmONg08voFJc3xpjg08voFIc2yN8qYh55zLOU8EzjmXc3lLBFclHUAdHt/ApDm+NMcGHt9ApDm2huSqj8A559wb5e2IwDnnXBVPBM45l3O5SQSSpkl6WNIKSV9OOp5KkraX9GdJ90taLOk/k46pkqTdJN0oaZmkpdFaE6kh6VxJD0Xv3XkpiOcaSc9Keqjitv+O3r8HJP1a0m4pi+9CSU9KWhSdpqcotkmS7oniWiDpmCRii2IZI2mupCXR/nZudPtp0fWtkrJXStrTijXtdCJMg/0osD8wDLgfmJB0XBXxCdgpujwUmA8cm3RcFfFdC5wVXR4G7JZ0TBWxTSSsdLcjYTbdPwIHJhzT8cDRwEMVt/0DMCS6fAlwScriuxD4Ugr+n7ViuxV4d3R5OtCVYHwjgaOjyzsT1lyZABwKHAx0AR1Jv499PeXliOAYYIWZrTSz14CfA6ckHNPfWbA+ujo0OqWiF1/SroQP59UAZvaamb2YaFCvdygw38w2mNlm4HbgA0kGZGZ3ENbXqLzt1ig+gHuA0S0PbFssb4gvLXqIzYBdosu7Ak+1NKjKQMyeNrP7ossvAUuBUWa21MweTiqugcpLIhgFrKq4vjq6LTUkDZa0CHgWuM3M5iccUtl+QDfwI0l/kfRDSW9KOqgKDwHHSdpT0o6EX4xjEo6pnn8BSkkHUcNno6arayTtnnQwFc4D/lvSKuCbwMxkwwkkjQeOIhzBZ1peEkHqmdkWM5tE+KV4jKSJCYdUNoRwqP59MzsKeBlITR+LmS0lNLXcCvwBWARsSTKm3kj6d2AzcH3SsVT5PnAAMAl4GvhWotG83qeBL5jZGOALREenSZK0E/BL4DwzW5d0PAOVl0TwJK//lTg6ui11omaXucC0hEMpWw2srjhCuZGQGFLDzK42s8lmdjzwAqHdNnUknQG8F/ioRY3MaWFmz0Q/RrYCPyA0p6bFJ4BfRZf/l4RjkzSUkASuN7Nf1ds+C/KSCO4FDpK0n6RhhLWRb0o4pr+TNKJcRSJpB+CdwLJEg4qY2d+AVZIOjm46EViSYEhvIGmv6HwsoX/gZ8lG9EaSpgEXACeb2Yak46kmaWTF1fcTmtzS4ingHdHlE4BHkgpEkghHJEvN7NtJxdFsuRlZHJXDXUqoILrGzL6RbETbSDqCUJkzmJCcZ5nZRclGtY2kScAPCRVDK4EzzeyFRIOqIOlOYE9gE3C+mc1OOJ4bgAJheuJngK8R2rW3A9ZEm91jZmenKL4CoVnIgMeBfzWzp1MS28PAdwnNlK8CnzGzha2OLYrv7cCdwIPA1ujmrxD+t/8DjABeBBaZ2buSiLE/cpMInHPO1ZaXpiHnnHM98ETgnHM554nAOedyzhOBc87lnCcC55zLOU8ELnei6SjKs2z+rWLWzfWSrojpNc+T9PHoclcjM1RG40v+EEc8zlUaknQAzrWama0h1Mwj6UJgvZl9M67XkzSEML9QwyOyJQ0xs25JT0uaamZ/iis+5/yIwLmIpIKk30WXL5R0raQ7JT0h6QOS/q+kByX9IZpmAEmTJd0uaaGkW6pG6JadANxXMfsowGnRGhTLJR0XPdcZkm6SNAcoD4r7DfDR2P5o5/BE4FxvDiB8iZ8M/BSYa2aHA68A74mSwf8Ap5rZZOAaoNaI9alA9UjYIWZ2DGFmza9V3H509HzlKRUWAMc1589xrjZvGnKuZyUz2yTpQcL0H+X2+geB8YSFSCYCt4UpaBhMmLmz2kjCvPWVypOVLYyeq+w2M6ucj/9ZYN/+/wnO1eeJwLmebQQws62SNlXMGLqV8NkRsNjM6i3d+Qqwfa3nJkyZXfk5fLlqu+2jxzsXG28acq7/HgZGlNdwljRU0mE1tlsKHNjP13gz6ZoJ1LUhTwTO9VO07OmpwCWS7icsivO2GpuWCMt99kcR+H0/H+tcQ3z2UedaQNKvgQvMrE9z6Uu6AzglTdN+u/bjicC5FogW9tk7Wpy90ceMAKaa2W9iC8w5PBE451zueR+Bc87lnCcC55zLOU8EzjmXc54InHMu5zwROOdczv1/XwXCE4MSPucAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhcElEQVR4nO3de3xdZZ3v8c83KaUVpqVABwqUBoVREToFMh0jaCMXBc8MMIge7HAEkVP1yDiKzgCDXESxZRjBGUWgYhGlXOQmnbHIpRJFG5AUYtqC0FqRthYItwACveV3/lhr4yZNsldWsi9pvu/Xa7+y1rOeZ63fzl47v6zb8ygiMDMzG6i6agdgZmbDkxOImZnl4gRiZma5OIGYmVkuTiBmZpbLqGoHUEk777xzNDQ0VDsMM7NhZcmSJc9GxMSe5SMqgTQ0NNDW1lbtMMzMhhVJf+it3KewzMwsFycQMzPLxQnEzMxycQIxM7NcnEDMzCyXqiYQSfMkPSNpWR/LJem/JK2U1CHpwKJlJ0lakb5OqlzUZmYG1T8C+T5wZD/LjwL2SV+zgMsBJO0InAf8LTAdOE/ShLJGWimtrTB7dvKznG2s8ir12Xofsgqp6nMgEfELSQ39VDkG+EEkfc7fL2kHSZOAZuDuiHgeQNLdJIno+jKHXF6trTBjBmzcCHV1MHUqjB/ff5uuLujogO7u7G0KZs6EWbMGH7eVVqnPtlJtink/GrGqfQRSyu7A6qL5NWlZX+VbkDRLUpukts7OzrIFOiRaWpI/MJB8mbu6Srfp6krqDqQNQHs7XHddnigtj0p9tpVqU+D9aETb6p9Ej4i5wFyAxsbG2h49q7k5+Q+wuxvGjoX586Gpqf82ra1w2GGwYQOMHp2tTWFbVjmV+mwr1ab4fdmIVesJZC0wuWh+j7RsLclprOLylopFVS5NTcnpg66u7F/ipiZYtCj5D7e5OfsX3yqrUp9tpdqYUfsJZAFwmqQbSC6Yd0XEOkl3Al8vunD+AeCsagU5pMaPT14D+RI3NflLPxxU6rOtVBsb8aqaQCRdT3IksbOkNSR3Vm0DEBFXAAuBDwErgVeBT6TLnpf0VeDBdFUXFC6om5lZZVT7LqyPlVgewGf7WDYPmFeOuMzMrLRavwvLzMxqlBOImZnl4gRiZma5OIGYmVkuTiBmZpaLE4iZmeXiBGJmZrk4gZiZWS5OIGZmlosTiJmZ5eIEYmZmuTiBmJlZLk4gZmaWixOImZnl4gRiZma5OIGYmVkuVU0gko6U9JiklZLO7GX5pZLa09fjkl4sWra5aNmCigZuZmbVG5FQUj1wGXAEsAZ4UNKCiHikUCcivlBU/5+AA4pW8VpETKtQuGZm1kM1j0CmAysjYlVEbABuAI7pp/7HgOsrEpmZmZVUzQSyO7C6aH5NWrYFSVOAvYCfFRWPkdQm6X5Jx/a1EUmz0nptnZ2dQxC2mZnB8LmIfgJwc0RsLiqbEhGNwEzgm5Le1lvDiJgbEY0R0Thx4sRKxGpmNiJUM4GsBSYXze+RlvXmBHqcvoqItenPVUALb74+YmZmZVbNBPIgsI+kvSSNJkkSW9xNJekdwASgtahsgqRt0+mdgYOBR3q2NTOz8qnaXVgRsUnSacCdQD0wLyKWS7oAaIuIQjI5AbghIqKo+TuBKyV1kyTBOcV3b5mZWflVLYEARMRCYGGPsnN7zJ/fS7vFwP5lDc7MzPo1XC6im5lZjXECMTOzXJxAzMwsFycQMzPLxQnEzMxycQIxM7NcnEDMzCwXJxAzM8vFCcTMzHJxAjEzs1ycQMzMLBcnEDMzy8UJxMzMcnECMTOzXJxAzMwslwElEEnbSaovVzBmZjZ89JtAJNVJminpJ5KeAX4LrJP0iKSLJe09mI1LOlLSY5JWSjqzl+UnS+qU1J6+Ti1adpKkFenrpMHEYWZmA1dqRMJ7gXuAs4BlEdENIGlH4P3ARZJui4hrB7rh9EjmMuAIYA3woKQFvQxNe2NEnNaj7Y7AeUAjEMCStO0LA43DzMzyKZVADo+IjT0LI+J54BbgFknb5Nz2dGBlRKwCkHQDcAyQZWzzDwJ3p3Eg6W7gSOD6nLGYmdkAlboGcrGkj0rava8KvSWYjHYHVhfNr0nLevqwpA5JN0uaPMC2SJolqU1SW2dnZ85Qzcysp1IJZCVwLPArSU9Iuk7SaZIOkFSJO7j+G2iIiKnA3cA1A11BRMyNiMaIaJw4ceKQB2hmNlL1mwQi4tsRMTMiGoD3ALcCbwVuAl4c5LbXApOL5vdIy4q3/1xErE9nrwIOytrWzMzKq+RRhBJTgaNJrlHMIDky+cYgt/0gsI+kvSSNBk4AFvTY9qSi2aOBR9PpO4EPSJogaQLwgbTMzMwqpN+L6OnF6XFAO3A/8PWIeLS/NllFxCZJp5H84a8H5kXEckkXAG0RsQD4nKSjgU3A88DJadvnJX2VJAkBXFC4oG5mZpVR6i6sVcBUYB/gOeBZSZ0R8exQbDwiFgILe5SdWzR9FsktxL21nQfMG4o4zMxs4PpNIBHxKQBJ44B3k1wH+aykiSTPhfgBPjOzEarUEUjBeuBV4LV0eg9gdLmCMjOz2leqK5NLJT0ArAO+AvwFcAXw9ojYvwLxmZlZjSp1BPJ74FqgPSI2VyAeMzMbJkolkF+mP/9a0hYLI+KhIY/IzMyGhVIJpA1YBhTuuirOIgEcWo6gzMys9pVKIKcDx5NcPL8BuC0iXil7VGZmVvNKdWXyzYg4BPgnkq5DFkn6kaRplQjOzMxqV6YOEdMu128H7iLphv2vyhmUmZnVvlJdmbyVpI+qY0i6T7+BpDuT1yoQm5mZ1bBS10BWAh0kRx8vAXsCnynckRURl5Q1OjMzq1mlEsgFJHdbAWxf5ljMzGwYKdUX1vkVisPMzIaZUl2ZfDkdb6Ov5YdK+ruhD8vMzGpdqVNYS4H/kfQ68BDQCYwh6d59GnAP8PVyBmhl0tWVvFpboamp2tGY2TBU6hTW7cDtkvYBDgYmkVxMvxaYNdi7sSQdCfwnyYBSV0XEnB7LTwdOJRlQqhM4JSL+kC7bTJLgAJ6MiKMHE8uI0toKHR3Q3Q2HHAJTp8L48dnazpwJs2aVN75a19oKLS3Q3OzkayNapu7cI2IFsGIoNyypHrgMOAJYAzwoaUFEPFJU7WGgMSJelfQZ4N+B/50uey0ipg1lTCNGSwtEem9Ed3dyJJIlgbS3Jz9HcgJpbYUZM2DjRqirG1jybW+HadPKGZ1ZRWUdD6QcpgMr04cUkXQDyfMmbySQiLi3qP79wIkVjXBr1dwMY8bAhg0wejTMn5/tP+nm5nJHVvtaWpLkAQNLvpAkj5kzyxWZWcVVM4HsTvJwYsEa4G/7qf9J4I6i+TGS2khOb82JiB/31kjSLGAWwJ577jmYeLceTU2waJFPw+TR3JwceXR3w9ix2ZOv2VaomgkkM0knAo3AjKLiKRGxNn1a/meSlkbE73q2jYi5wFyAxsbG6Ll8xGpq8h++PJqaktNWXV1OHjbilerK5Fv8+UHCLUTE5wax7bUkHTQW7JGW9YzhcOBsYEZErC/a9tr05ypJLcABwBYJxGzIjR+fvJw8bIQr1ZliG7CE5NbdA0kupK8guYV3sGOiPwjsI2kvSaNJ+txaUFxB0gHAlcDREfFMUfkESdum0zuT3CFWfPHdzMzKrNRtvNcApHdAHRIRm9L5K4D7BrPhiNgk6TTgTpLbeOdFxHJJFwBtEbEAuJikC5Wb0v63CrfrvhO4UlI3SRKc0+PuLTMzK7Os10AmAOOA59P57dOyQYmIhcDCHmXnFk0f3ke7xcD+g92+mZnllzWBzAEelnQvybC27wPOL1dQZmZW+7I+SHi1pDv48222Z0TEU+ULy8zMal2pzhTfkf48ENiN5LmN1cBuaZmZmY1QpY5ATid5CO8bvSwL4NAhj8jMzIaFUndhzZJUB3w5In5VoZjMzGwYKPUcCBHRDXy7ArGYmdkwUjKBpBZJ+rAKg6GbmdmIlzWBfAq4CVgv6SVJL0t6qYxxmZlZjct6G+9flDsQMzMbXkp1plgPjI2IV9L5d/PnPrAejoiXyxyfmZnVqFJHIBcBz5CMBAhwPbCMpHPFh4AzyheamZnVslIJ5DDgb4rmX4yIv08vpg+qM0UzMxveSl1Eryv0wJs6AyAigqRDRTMzG6FKJZDRkt64gB4RdwFIGk9yGsvMzEaoUgnku8CNkt4YTFzSFJJrIVeVMzAzM6ttpboyuUTSq8AvJW2XFr9CMoDT5WWPzszMalaWrkyuiIg9gQagISKmDFXykHSkpMckrZR0Zi/Lt5V0Y7r8AUkNRcvOSssfk/TBoYjHzMyyy/okOhHx8lA+95E+Y3IZcBSwL/AxSfv2qPZJ4IWI2Bu4lOS2YtJ6JwDvAo4EvpOuryxa75jL7K99kNY75pa1DUDruC5m7/kkratbs7dZ3crs+2aXv03dH5k9rmNg76m1FWbPTn4ORI52lfqc8nxGUMHPqUJtIN8+keu7kWc/qlSbnO0qtr/m/FuUhZIbqipPUhNwfkR8MJ0/CyAiZhfVuTOt0yppFPAUMBE4s7hucb3+ttnY2BhtbW0DirP1jrm8r/VTbKqDuoCpL41hfGzbb5suradj3Ot0K3sbgK76TXRs/6ekXV0dU3eZyvhtx/ffZn0XHU930B3d1KmMbZ5fR0fX4wN7T5s3wSt/YuZSmPVwHUydCuP7306ysS7o6IDubqjL1q617o/MOGQFG8v8OeX5jKCCn1OF2kC+fSLXdyPPfpRjH8rVJme7iu2vaZsQjNkEiw6+kqajZpV+Tz1IWhIRjT3LMx+BlMHuJINTFaxJy3qtk95O3AXslLEtAJJmSWqT1NbZ2TngIFuW3MLmOkDQreQPSCld9Zvo1sDaAHRtU9Quuul6vat0m9e76I5uoMxtXu4c+HvatJn2XeG6/Um+XF2lt5NsrCupT/Z2Lds/y8YKfE55PiOo4OdUoTaQb5/I9d3Isx/l2IdytcnZrmL7a9omBBvqkr9nQyoiSr6AeuBo4HMkg0ydDpyepW0/6zweuKpo/v8A3+5RZxmwR9H874CdSbqXP7Go/HvA8aW2edBBB8VALV54ZYw9m6g/hxh7NrF44ZVlaRMRsfjJxTH2a2Oj/iv1MfZrY2Pxk4trp02e97R4ccw4pS5mnEzE2LERi0tvp9Auxo6NqK/P3G7xwiuj7lyC88r7OeX53eVtV8ttIir43cizH+XYh3K1ydmuYvtrzr9FPQFt0cvf1EynsCQtBF4HlgLdRcnnK3kT13A5hQXJaayWJbfQfNCHMx/+5WkDybnolidaaG5opmlyU221yfGemv/rQHjxRVqOmA9N2baTbKwVWlqguTlzuwMueTtdL3cyf/qcsn5OeX53edvVchuo3Hcj136UYx/K1SZnu4rtrzn/FhXr6xRW1gTSERFTc22573WOAh4n6S5lLfAgMDMilhfV+Sywf0R8WtIJwHER8VFJ7wKuA6aTjNW+CNgnIjb3t828CcTya/5+MwAtJ7dsVduyytoaP9vh9J76SiCZunMH7pD0gUifRB8KEbFJ0mnAnSSnyOZFxHJJF5AcLi0gOTX1Q0krgedJ7rwirfcj4BFgE/DZUsnDzMyGVtYEcj9wWzo++kZAJF1ijRvMxiNiIbCwR9m5RdOvAx/po+2FwIWD2b6ZmeWXNYFcAjQBSyPLOS8zM9vqZb2NdzWwzMnDzMwKsh6BrAJaJN0BrC8URsQlZYnKzMxqXtYE8vv0NZo/D2lrZmYjWKYEUnjeQ9L26fwr5QzKzMxqX6ZrIJL2k/QwsBxYLmlJ+iyGmZmNUFkvos8l6bpkSkRMAb5IMtiUmZmNUFkTyHYRcW9hJiJagO36rm5mZlu7zHdhSToH+GE6fyLJnVlmZjZCZT0COYWkE8NbgVtIesQ9pVxBmZlZ7St5BJKO9HdrRLy/AvGYmdkwkWVM9M1At6QMQ3OZmdlIkfUayCvAUkl3A38qFEbE58oSlZmZ1bysCeTW9GVmZgaUSCCSFkXEYcC+EXFGhWIyM7NhoNQRyCRJ7wGOlnQDyTggb4iIh8oWmZmZ1bRSCeRc4BxgD5IxQYoFcGiejUraEbgRaACeAD4aES/0qDMNuBwYB2wGLoyIG9Nl3wdmAF1p9ZMjoj1PLGZmlk+/CSQibgZulnRORHx1CLd7JrAoIuZIOjOd73mK7FXg4xGxQtJuwBJJd0bEi+nyf0njMzOzKsj0IOEQJw+AY4Br0ulrgGN72ebjEbEinf4j8AzJw4xmZlYDsj6JPtR2iYh16fRTwC79VZY0nWQckt8VFV8oqUPSpZK27aftLEltkto6OzsHHbiZmSXKlkAk3SNpWS+vY4rrpcPk9jlUrqRJJH1wfSIiutPis4B3AH8D7MiWp7+K1z83IhojonHiRB/AmJkNlaxdmSyPiHcMZMURcXg/63xa0qSIWJcmiGf6qDcO+AlwdkTcX7TuwtHLeklXA18aSGxmZjZ4WbsyeUzSnkO43QXASen0ScDtPStIGg3cBvyg58XyNOkgSSTXT5YNYWxmZpZB1ifRJ5CMRPhr3tyVydE5tzsH+JGkTwJ/AD4KIKkR+HREnJqWvQ/YSdLJabvC7brzJU0keS6lHfh0zjjMzCynrAnknKHcaEQ8BxzWS3kbcGo6fS1wbR/tcz1/YmZmQydTAomIn0uaAuwTEfdIegtQX97QzMyslmW6C0vS/wVuBq5Mi3YHflymmMzMbBjIehvvZ4GDgZcA0gf8/rJcQZmZWe3LmkDWR8SGwoykUfTz7IaZmW39siaQn0v6N2CspCOAm4D/Ll9YZmZW67ImkDOBTmAp8ClgIfDlcgVlZma1L+ttvO8Hro2I75YzGDMzGz6yHoF8HPiNpPslXSzp7yVNKGdgZmZW27I+B3ISQDoux/HAZcBuWdubmdnWJ1MCkHQi8F5gf+BZ4NvAfWWMy8zMalzWI4hvkozFcQVwb0Q8Ua6AzMxseMg6IuHOwCnAGJKBnH4t6YdljczMzGpa1q5MxgF7AlOABmA80N1fGzMz27plPYX1y6LXtyNiTflCMjOz4SDrXVhTASRtX95wzMxsuMh6Cms/SQ8Dy4FHJC2RtF/ejUraUdLdklakP3t9pkTSZknt6WtBUflekh6QtFLSjenohWZmVkFZHyScC5weEVMiYk/gi2lZXmcCiyJiH2BROt+b1yJiWvoqHv3wIuDSiNgbeAH45CBiMTOzHLImkO0i4t7CTES0ANsNYrvHANek09eQjGueSToO+qEk45MMuL2ZmQ2NrAlklaRzJDWkry8Dqwax3V0iYl06/RSwSx/1xkhqS7tQOTYt2wl4MSI2pfNrSAa46pWkWek62jo7OwcRspmZFct6F9YpwFeAW0nGAbkvLeuTpHuAXXtZdHbxTESEpL7GFpkSEWslvRX4maSlQFfGmAvrn0t6uq2xsdFjmJiZDZF+E4ikMcCngb1JunL/YkRszLLiiDi8n/U+LWlSRKyTNAl4po91rE1/rpLUAhwA3ALsIGlUehSyB7A2S0xmZjZ0Sp3CugZoJEkeRwEXD9F2FwAnpdMnAbf3rCBpgqRt0+mdSYbUfSQiAriXpFPHPtubmVl5lUog+0bEiRFxJckf7PcN0XbnAEdIWgEcns4jqVHSVWmddwJtkn5DkjDmRMQj6bIzgNMlrSS5JvK9IYrLzMwyKnUN5I3TVRGxKbkBavAi4jngsF7K24BT0+nFJL3/9tZ+FTB9SIIxM7NcSiWQv5b0UjotkjHRX0qnIyLGlTU6MzOrWf0mkIior1QgZmY2vGR9DsTMzOxNnEDMzCwXJxAzM8vFCcTMzHJxAjEzs1ycQMzMLBcnEDMzy8UJxMzMcnECMTOzXJxAzMwsFycQMzPLxQnEzMxycQIxM7NcnEDMzCyXqiQQSTtKulvSivTnhF7qvF9Se9HrdUnHpsu+L+n3RcumVfo9mJmNdNU6AjkTWBQR+wCL0vk3iYh7I2JaREwDDgVeBe4qqvIvheUR0V6BmM2sSrrWd/Fk15O0rm6tdihWpNSIhOVyDNCcTl8DtJCMc96X44E7IuLV8oZlZrWmdXUrHU930B3dHHL1IUzdZSrjtx2fqe3M/Wcy66BZZY5w5KrWEcguEbEunX4K2KVE/ROA63uUXSipQ9Klkrbtq6GkWZLaJLV1dnYOImQzq4aWJ1qICAC6o5uu17sytWt/qp3rll5XztBGvLIdgUi6B9i1l0VnF89EREiKftYzCdgfuLOo+CySxDMamEty9HJBb+0jYm5ah8bGxj63Y2a1qbmhmTGjxrBh8wZG149m/nHzaZrcVLrd95vLH9wIV7YEEhGH97VM0tOSJkXEujRBPNPPqj4K3BYRG4vWXTh6WS/pauBLQxK0mdWcpslNLPr4IlqeaKG5oTlT8rDKqNY1kAXAScCc9Oft/dT9GMkRxxuKko+AY4FlZYrTzGpA0+QmJ44aVK1rIHOAIyStAA5P55HUKOmqQiVJDcBk4Oc92s+XtBRYCuwMfK0SQZuZ2Z9V5QgkIp4DDuulvA04tWj+CWD3XuodWs74zMysND+JbmZmuTiBmJlZLk4gZmaWixOImZnl4gRiZma5OIGYmVkuTiBmZpaLE4iZmeXiBGJmZrk4gZiZWS5OIGZmlosTiJmZ5eIEYmZmuTiBmJlZLk4gZmaWixOImZnlUpUEIukjkpZL6pbU2E+9IyU9JmmlpDOLyveS9EBafqOk0ZWJ3MzMCqp1BLIMOA74RV8VJNUDlwFHAfsCH5O0b7r4IuDSiNgbeAH4ZHnDtby61nfxZNeTtK5uHVC71tWtzL5v9oDbmRXk2ffy7ncjdX+t1pC2jwJI6q/adGBlRKxK694AHCPpUeBQYGZa7xrgfODycsVr+bSubqXj6Q66o5tDrj6EqbtMZfy240u261rf9Ua7OtVlbtf+VDvTdp02BJHbcJdn38u7343k/bWWr4HsDqwuml+Tlu0EvBgRm3qU90rSLEltkto6OzvLFqxtqeWJFiICgO7opuv1rkztul7voju6B9xu2q7TmLn/zNIVbauXZ9/Lu9+N5P21bEcgku4Bdu1l0dkRcXu5tttTRMwF5gI0NjZGpbZr0NzQzJhRY9iweQOj60cz/7j5NE1uKtmudXUrh/3gsAG3MyvIs+/l3e9G8v5atgQSEYcPchVrgclF83ukZc8BO0galR6FFMqtxjRNbmLRxxfR8kQLzQ3Nmb9UeduZFeTZh7y/DpwKh3lV2bjUAnwpItp6WTYKeBw4jCRBPAjMjIjlkm4CbomIGyRdAXRExHdKba+xsTHa2rbYlJmZ9UPSkojY4o7Zat3G+w+S1gBNwE8k3ZmW7yZpIUB6dHEacCfwKPCjiFieruIM4HRJK0muiXyv0u/BzGykq+oRSKX5CMTMbOBq6gjEzMyGPycQMzPLxQnEzMxycQIxM7NcRtRFdEmdwB9yNt8ZeHYIwxlqtRxfLccGjm+wajm+Wo4Naj++gikRMbFn4YhKIIMhqa23uxBqRS3HV8uxgeMbrFqOr5Zjg9qPrxSfwjIzs1ycQMzMLBcnkOzmVjuAEmo5vlqODRzfYNVyfLUcG9R+fP3yNRAzM8vFRyBmZpaLE4iZmeXiBJKBpCMlPSZppaQzqx1PgaQxkn4t6TeSlkv6SrVj6knSDpJulvRbSY9KqqnBEiT9s6Rl6e/v8zUQzzxJz0haVlR2cfr765B0m6Qdaii28yWtldSevj5Ujdj6iW+apPvT2NokTa9SbJMl3SvpkXRf++e0/CPpfLek4Xc7b0T41c8LqAd+B7wVGA38Bti32nGlsQnYPp3eBngAeHe14+oR4zXAqen0aGCHasdUFNt+wDLgLSSDq90D7F3lmN4HHAgsKyr7ADAqnb4IuKiGYjufZEyfWvg8e4vvLuCodPpDQEuVYpsEHJhO/wXJWEf7Au8E3g60AI3V/h0O9OUjkNKmAysjYlVEbABuAI6pckwAROKVdHab9FUzd0VIGk/ypf4eQERsiIgXqxrUm70TeCAiXo1k/JmfA8dVM6CI+AXwfI+yu9L4AO4nGYWz4nqLrZb0EV8A49Lp8cAfKxpUIYiIdRHxUDr9MskYR7tHxKMR8Vg1YhoKTiCl7Q6sLppfk5bVBEn1ktqBZ4C7I+KBKodUbC+gE7ha0sOSrpK0XbWDKrIMeK+knSS9heQ/1Mkl2lTbKcAd1Q6ih9PS02vzJE2odjA9fB64WNJq4D+As6obDkhqAA4gOWMwrDmBDHMRsTkippH8Vzpd0n5VDqnYKJJTCpdHxAHAn4CauYYUEY+SnBK6C/gp0A5srmZM/ZF0NrAJmF/tWIpcDrwNmAasA75R1Wi29BngCxExGfgCVR69VNL2wC3A5yPipWrGMhScQEpby5v/K90jLasp6amhe4EjqxxKsTXAmqKjoptJEkrNiIjvRcRBEfE+4AWSc9M1R9LJwN8B/xjpifRaEBFPp//EdAPfJTnlW0tOAm5Np2+iivFJ2oYkecyPiFtL1R8OnEBKexDYR9JekkYDJwALqhwTAJImFu7IkTQWOAL4bVWDKhIRTwGrJb09LToMeKSKIW1B0l+mP/ckuf5xXXUj2pKkI4F/BY6OiFerHU8xSZOKZv+B5LRgLfkjMCOdPhRYUY0gJInk6OfRiLikGjGUg59EzyC9NfGbJHdkzYuIC6sbUULSVJK7nOpJ/hn4UURcUN2o3kzSNOAqkjuwVgGfiIgXqhpUEUn3ATsBG4HTI2JRleO5Hmgm6eb7aeA8kvP22wLPpdXuj4hP10hszSSnrwJ4AvhURKyrdGz9xPcY8J8kp1NfB/5fRCypQmyHAPcBS4HutPjfSD7XbwETgReB9oj4YKXjy8sJxMzMcvEpLDMzy8UJxMzMcnECMTOzXJxAzMwsFycQMzPLxQnELKO0y5NCr7NPFfVC+4qk75Rpm5+X9PF0uiVLj63p80E/LUc8ZsVGVTsAs+EiIp4jeeYBSecDr0TEf5Rre5JGkfR9lfnpfUmjIqJT0jpJB0fEr8oVn5mPQMwGSVKzpP9Jp8+XdI2k+yT9QdJxkv5d0lJJP027s0DSQZJ+LmmJpDt7PNFdcCjwUFFPvAAfSceAeVzSe9N1nSxpgaSfAYUHIX8M/GPZ3rQZTiBm5fA2kj/+RwPXAvdGxP7Aa8D/SpPIt4DjI+IgYB7QW+8GBwM9n5oeFRHTSXqZPa+o/MB0fYVuO9qA9w7N2zHrnU9hmQ29OyJio6SlJN3MFK5HLAUaSAYQ2g+4O+kiiXqSnmx7mkQybkSxQid8S9J1FdwdEcVjYTwD7Jb/LZiV5gRiNvTWA0REt6SNRb3ndpN85wQsj4hSw/u+Bozpbd0k3c4Xf3//1KPemLS9Wdn4FJZZ5T0GTCyMDy9pG0nv6qXeo8DeObfxV9Rez7i2lXECMauwdGjk44GLJP2GZCCr9/RS9Q6SIYHzeD/wk5xtzTJxb7xmNUzSbcC/RsSAxrGQ9AvgmFrqOt+2Pk4gZjUsHYxrl4j4xQDaTAQOjogfly0wM5xAzMwsJ18DMTOzXJxAzMwsFycQMzPLxQnEzMxycQIxM7Nc/j8zEYVP0ginFgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the state of charge (E)\n", "plt.figure()\n", "plt.plot(time,E_control,'b.-')\n", "plt.xlabel('Time (hr)')\n", "plt.ylabel('SOC (MWh)')\n", "plt.xticks(range(0,24,3))\n", "plt.show()\n", "\n", "# Plot the charging and discharging rates \n", "# we can do a step/stair plot as the control values are constant across each hr\n", "plt.figure()\n", "plt.step(time,c_control,'r.-',where=\"post\")\n", "plt.step(time,d_control,'g.-',where=\"post\")\n", "plt.xlabel('Time (hr)')\n", "plt.ylabel('Power from Grid (MW)')\n", "plt.xticks(range(0,24,3))\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How much money can a 4 MWh battery make in a year?" ] }, { "cell_type": "markdown", "metadata": { "nbgrader": { "grade": false, "locked": false, "solution": false } }, "source": [ "
\n", "

Class Activity

\n", " Copy the code from the cell below the heading Solve Optimization Model to the cell below. Then modify to calculate the revenue for an entire year. Hint: Do NOT modify the function `create_model`.\n", "
" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "tags": [ "remove-output" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ipopt 3.11.1: \n", "\n", "******************************************************************************\n", "This program contains Ipopt, a library for large-scale nonlinear optimization.\n", " Ipopt is released as open source code under the Eclipse Public License (EPL).\n", " For more information visit http://projects.coin-or.org/Ipopt\n", "******************************************************************************\n", "\n", "NOTE: You are using Ipopt by default with the MUMPS linear solver.\n", " Other linear solvers might be more efficient (see Ipopt documentation).\n", "\n", "\n", "This is Ipopt version 3.11.1, running with linear solver mumps.\n", "\n", "Number of nonzeros in equality constraint Jacobian...: 35039\n", "Number of nonzeros in inequality constraint Jacobian.: 0\n", "Number of nonzeros in Lagrangian Hessian.............: 0\n", "\n", "Total number of variables............................: 26280\n", " variables with only lower bounds: 0\n", " variables with lower and upper bounds: 26280\n", " variables with only upper bounds: 0\n", "Total number of equality constraints.................: 8760\n", "Total number of inequality constraints...............: 0\n", " inequality constraints with only lower bounds: 0\n", " inequality constraints with lower and upper bounds: 0\n", " inequality constraints with only upper bounds: 0\n", "\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 0 -1.4265256e-012 1.13e-002 2.01e+001 -1.0 0.00e+000 - 0.00e+000 0.00e+000 0\n", " 1 2.0748446e+002 1.02e-002 1.87e+001 -1.0 2.29e-001 - 7.14e-002 9.38e-002f 1\n", " 2 1.2878907e+002 8.90e-003 1.86e+001 -1.0 1.32e+000 - 3.65e-002 1.29e-001f 1\n", " 3 -3.0565744e+002 8.10e-003 1.85e+001 -1.0 5.90e+000 - 2.60e-002 8.93e-002f 1\n", " 4 -8.2367026e+002 7.70e-003 1.79e+001 -1.0 1.18e+001 - 3.67e-002 4.97e-002f 1\n", " 5 -1.6963163e+003 7.28e-003 1.73e+001 -1.0 1.63e+001 - 3.90e-002 5.48e-002f 1\n", " 6 -3.3783979e+003 6.73e-003 1.70e+001 -1.0 1.28e+001 - 2.51e-002 7.59e-002f 1\n", " 7 -4.8399880e+003 6.32e-003 1.61e+001 -1.0 9.92e+000 - 5.62e-002 6.11e-002f 1\n", " 8 -6.5735466e+003 5.89e-003 1.53e+001 -1.0 8.87e+000 - 5.36e-002 6.67e-002f 1\n", " 9 -9.4760273e+003 5.24e-003 1.50e+001 -1.0 7.89e+000 - 3.45e-002 1.11e-001f 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 10 -1.1338827e+004 4.80e-003 1.40e+001 -1.0 5.74e+000 - 7.04e-002 8.36e-002f 1\n", " 11 -1.3594338e+004 4.29e-003 1.30e+001 -1.0 5.39e+000 - 7.77e-002 1.08e-001f 1\n", " 12 -1.6138331e+004 3.72e-003 1.26e+001 -1.0 5.01e+000 - 7.33e-002 1.32e-001f 1\n", " 13 -1.8398194e+004 3.20e-003 1.19e+001 -1.0 4.17e+000 - 9.64e-002 1.40e-001f 1\n", " 14 -2.0238358e+004 2.75e-003 1.11e+001 -1.0 3.76e+000 - 1.01e-001 1.40e-001f 1\n", " 15 -2.2207864e+004 2.24e-003 9.97e+000 -1.0 3.49e+000 - 1.36e-001 1.87e-001f 1\n", " 16 -2.3865802e+004 1.75e-003 8.54e+000 -1.0 2.83e+000 - 1.72e-001 2.17e-001f 1\n", " 17 -2.4915328e+004 1.40e-003 6.62e+000 -1.0 2.54e+000 - 2.17e-001 2.02e-001f 1\n", " 18 -2.6025064e+004 9.75e-004 5.25e+000 -1.0 2.15e+000 - 2.41e-001 3.03e-001f 1\n", " 19 -2.6714253e+004 6.54e-004 3.17e+000 -1.0 1.64e+000 - 3.74e-001 3.29e-001f 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 20 -2.7173480e+004 3.93e-004 1.58e+000 -1.0 1.36e+000 - 4.66e-001 3.99e-001f 1\n", " 21 -2.7569525e+004 1.13e-004 3.67e-001 -1.0 9.92e-001 - 8.72e-001 7.13e-001f 1\n", " 22 -2.8821288e+004 3.41e-005 9.93e-002 -1.7 6.97e-001 - 5.82e-001 6.97e-001f 1\n", " 23 -2.9338373e+004 1.65e-005 1.23e-001 -2.5 6.46e-001 - 4.15e-001 5.17e-001f 1\n", " 24 -2.9576291e+004 8.74e-006 1.66e-001 -2.5 6.75e-001 - 3.67e-001 4.70e-001f 1\n", " 25 -2.9703152e+004 4.71e-006 3.83e-002 -2.5 6.19e-001 - 5.46e-001 4.62e-001f 1\n", " 26 -2.9791043e+004 1.97e-006 1.59e-002 -2.5 4.48e-001 - 5.84e-001 5.81e-001f 1\n", " 27 -2.9862087e+004 8.85e-007 3.44e-002 -3.8 3.90e-001 - 4.58e-001 5.51e-001f 1\n", " 28 -2.9888279e+004 4.89e-007 2.88e-002 -3.8 5.11e-001 - 3.99e-001 4.47e-001f 1\n", " 29 -2.9907078e+004 2.06e-007 2.71e-002 -3.8 6.00e-001 - 4.54e-001 5.78e-001f 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 30 -2.9915085e+004 8.63e-008 1.18e-002 -3.8 3.18e-001 - 5.75e-001 5.81e-001f 1\n", " 31 -2.9919911e+004 1.43e-008 2.67e-002 -3.8 2.16e-001 - 5.13e-001 8.34e-001f 1\n", " 32 -2.9920812e+004 9.66e-010 6.88e-003 -3.8 1.40e-001 - 1.00e+000 9.33e-001f 1\n", " 33 -2.9920878e+004 8.88e-016 1.07e-014 -3.8 7.17e-002 - 1.00e+000 1.00e+000f 1\n", " 34 -2.9923433e+004 8.88e-016 7.13e-004 -5.7 1.06e-001 - 8.16e-001 7.00e-001f 1\n", " 35 -2.9924388e+004 8.88e-016 9.82e-005 -5.7 5.79e-002 - 8.81e-001 8.73e-001f 1\n", " 36 -2.9924522e+004 8.88e-016 5.71e-005 -5.7 2.01e-001 - 1.00e+000 9.63e-001f 1\n", " 37 -2.9924527e+004 8.88e-016 9.84e-015 -5.7 3.46e-002 - 1.00e+000 1.00e+000f 1\n", " 38 -2.9924572e+004 1.33e-015 9.06e-004 -8.6 3.73e-002 - 8.56e-001 9.93e-001f 1\n", " 39 -2.9924572e+004 8.88e-016 1.72e-001 -8.6 6.70e-003 - 1.00e+000 7.81e-001f 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 40 -2.9924572e+004 1.33e-015 1.22e-014 -8.6 1.40e-003 - 1.00e+000 1.00e+000f 1\n", "\n", "Number of Iterations....: 40\n", "\n", " (scaled) (unscaled)\n", "Objective...............: -2.5721653825666763e+004 -2.9924572060780712e+004\n", "Dual infeasibility......: 1.2178002166165205e-014 1.4167887720116600e-014\n", "Constraint violation....: 1.3322676295501878e-015 1.3322676295501878e-015\n", "Complementarity.........: 4.1788718664542149e-009 4.8616995294328333e-009\n", "Overall NLP error.......: 4.1788718664542149e-009 4.8616995294328333e-009\n", "\n", "\n", "Number of objective function evaluations = 41\n", "Number of objective gradient evaluations = 41\n", "Number of equality constraint evaluations = 41\n", "Number of inequality constraint evaluations = 0\n", "Number of equality constraint Jacobian evaluations = 41\n", "Number of inequality constraint Jacobian evaluations = 0\n", "Number of Lagrangian Hessian evaluations = 40\n", "Total CPU secs in IPOPT (w/o function evaluations) = 3.910\n", "Total CPU secs in NLP function evaluations = 0.058\n", "\n", "EXIT: Optimal Solution Found.\n", "Revenue in 2015: $ 29924.57150833467\n" ] } ], "source": [ "# Add your solution here\n" ] } ], "metadata": { "colab": { "name": "L12-Optimization.ipynb", "provenance": [], "version": "0.3.2" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.13" } }, "nbformat": 4, "nbformat_minor": 4 }