pv_size_optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PV size optimization\n",
    "This notebook optimizes the size of the pv modules for each house, such that it has an optimal size of the pv module. Optimal is in this case that the cost over a specific time period minimizes. Therefore we consider a finite timespan. We simulate a single year and then multiply the import costs from the grid and rewards, that we get from the grid company, by the number of years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from src.grid import Grid\n",
    "import random\n",
    "from ipywidgets import IntProgress\n",
    "from IPython.display import display\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import datasets\n",
    "We only choose a single house here for the consumption. For the production by the pv panels we choose one from the data set. We then scale its size and assume that the production increases linear with the size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Consumption data\n",
    "df_cons = pd.read_csv('data/consumption_ausgrid.csv')\n",
    "df_cons['Unnamed: 0'] = pd.to_datetime(df_cons['Unnamed: 0'])\n",
    "df_cons = df_cons.set_index('Unnamed: 0')\n",
    "df_cons_red = df_cons.iloc[:,0]\n",
    "# Production data\n",
    "df_prod = pd.read_csv('data/normalized_generation_ausgrid_025steps.csv')\n",
    "df_prod['Unnamed: 0'] = pd.to_datetime(df_prod['Unnamed: 0'])\n",
    "df_prod = df_prod.set_index('Unnamed: 0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve data from consumption\n",
    "data_consumption = df_cons_red.values\n",
    "\n",
    "# Retrieve data from production and pv_sizes\n",
    "data_production = df_prod.values\n",
    "pv_peakpowers = [float(key) for key in df_prod.keys()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Single Household\n",
    "### Create a grid model \n",
    "This creates a standard grid model, which does not include any storage device (i.e. one storage with zero capacity). Therefore we get the pure consumption data and the price for this consumption. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'import_grid_kwh': 7679.200999999989, 'export_grid_kwh': 0.0, 'pv_production_kwh': 0.0, 'setup_cost_storage': 0, 'setup_cost_pv': 0.0, 'cost_import_grid': 1919.8002499999973, 'reward_export_grid': 0.0}\n"
     ]
    }
   ],
   "source": [
    "g = Grid(num_houses=1,num_storages=1,max_capacity=0,num_pvtypes=len(pv_peakpowers),pv_peakpower=pv_peakpowers)\n",
    "g.set_costs() # Set standard costs\n",
    "res_dict = g.simulate(data_cons=data_consumption,data_prod=data_production)\n",
    "print(res_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Objective function \n",
    "We now need to define an objective function, that we want to optimize. As outlined before we want to optimize here the cost within a certain time frame (e.g. 7 years). Setting up storages and PV systems are one time expenses, so we just leave them as they result from the model simulation. However the costs and reward from grid interactions are only valid for a single year. Hence we multiply them by the number of years, which we look at, to get an approximation for the total costs.\n",
    "We can also think of completly different objective functions, such as minimizing the interaction with the grid. Or combinations of costs and grid dependency factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective_function(res_dict,years=7):\n",
    "    total_cost = res_dict[\"setup_cost_storage\"] + res_dict[\"setup_cost_pv\"]\n",
    "    total_cost += res_dict[\"cost_import_grid\"]*years\n",
    "    total_cost -= res_dict[\"reward_export_grid\"]*years\n",
    "    return total_cost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying the objective function to the model without a PV installation, gives us simply the normal cost for energy in the house over the defined timespan."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13438.601749999982\n"
     ]
    }
   ],
   "source": [
    "total_cost_noPV = objective_function(res_dict)\n",
    "print(total_cost_noPV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying the objective function now to the grid model for all possible PV types gives us a \"cost-landscape\".\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEKCAYAAADaa8itAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4VGX2wPHvSQiEGqQLAQMYmpRAAtKCgIqgrrCCAiKgoojYCy6WtWxxLfzUteGiooBIESkWFHVXBRGQEBJClSIlFAkBQgkh7fz+mBuMEMIAU5KZ83meeTLz3jtzz6TMyXvf955XVBVjjDHGE0L8HYAxxpjAYUnFGGOMx1hSMcYY4zGWVIwxxniMJRVjjDEeY0nFGGOMx1hSMcYY4zGWVIwxxniMJRVjjDEeU8bfAfhajRo1NCoqyt9hGGNMqbJixYp9qlrzTPsFXVKJiooiISHB32EYY0ypIiLb3NnPTn8ZY4zxGEsqxhhjPMaSijHGGI8JujGVouTk5JCamkpWVpa/QzEeFh4eTmRkJGFhYf4OxZigYEkFSE1NpXLlykRFRSEi/g7HeIiqkp6eTmpqKg0bNvR3OMYEBTv9BWRlZVG9enVLKAFGRKhevbr1QI3xIUsqDksogcl+rsb4liUVY4wJcJvTjvDiV+vxxfLxllRKiD179jBo0CAaN25MixYtuPrqq/nll1/O+nWee+45j8b19ttvM3nyZAA++OADdu3adWJbVFQU+/btO+9jHDx4kLfeeuu020NDQ4mJiaFly5bccMMNZGZm0r17dxYsWPCH/V599VVGjx593vEYE0hWbDvAgPE/MWP5DnZneP9UsCWVEkBV+fOf/0z37t3ZvHkza9eu5bnnnuO3334769fydFIZNWoUw4YNA05NKp5ypqRSvnx5kpKSWL16NWXLluXtt99m8ODBTJ8+/Q/7TZ8+ncGDB3s8PmNKq6/X7OGmd5YSUT6M2aM7U7dqea8f05JKCfDdd98RFhbGqFGjTrTFxMQQHx+PqjJmzBhatmxJq1atmDFjBgC7d++mW7duJ/6DX7RoEWPHjuXYsWPExMQwZMgQjh49yjXXXEObNm1o2bLliecW2Lt3L7GxsQAkJycjImzfvh2Axo0bk5mZyTPPPMO4ceOYNWsWCQkJDBkyhJiYGI4dOwbA66+/Trt27WjVqhXr168HYP/+/fTr14/WrVvTsWNHVq1aBXDitQq0bNmSrVu3MnbsWDZv3kxMTAxjxowp9nsVHx/Ppk2bGDBgAJ9//jnHjx8HYOvWrezatYuuXbue88/BmEDy4dJtjPpwBc0urMInd3XmouoVfXJcm1J8sgcegKQkz75mTAy8+uppN69evfrEh/vJZs+eTVJSEsnJyezbt4/27dvTrVs3PvroI6666iqeeOIJ8vLyyMzMJD4+njfeeIMkJ/5PPvmEunXr8sUXXwCQkZHxh9euVasWWVlZHDp0iEWLFhEXF8eiRYvo2rUrtWrVokKFCif2HTBgAG+88Qbjxo0jLi7uRHuNGjVITEzkrbfeYty4cbz77rs8/fTTtG3blrlz5/K///2PYcOGnYipKM8//zyrV68udh+A3NxcvvzyS3r37k316tXp0KEDX331FX379mX69OkMHDjQBuZN0FNV/u/rX3jju030bFaLN25qS4Wyvvuot55KCffjjz8yePBgQkNDqV27NpdddhnLly+nffv2vP/++zzzzDOkpKRQuXLlU57bqlUrvv32W/7yl7+waNEiIiIiTtmnc+fOLF68mIULF/L444+zcOFCFi1aRHx8vFvxXX/99QDExsaydevWEzEPHToUgJ49e5Kenn5KQjsbBb2vuLg4GjRowIgRIwD+cArMTn0ZAzl5+Tzy8Sre+G4TgzvUZ8LQWJ8mFLCeyqmK6VF4yyWXXMKsWbOK3Ha62RrdunVj4cKFfPHFFwwdOpQxY8acGPso0KRJE1asWMH8+fN57LHH6NWrF0899dQf9omPj2fRokVs27aNvn378sILLyAiXHvttW7FXq5cOcA1mJ6bm3vamEWEMmXKkJ+ff6LN3etHCsZUTtavXz8eeughEhMTOXbsGO3atXPr9YwJREeO5zJ6aiILf0njwSuacN/lF/ul5249lRKgZ8+eHD9+nHfeeedE2/Lly/nhhx/o1q0bM2bMIC8vj7S0NBYuXEiHDh3Ytm0btWrV4o477mDEiBEkJiYCEBYWRk5ODgC7du2iQoUK3HzzzTzyyCMn9imsW7dufPjhh0RHRxMSEkK1atWYP38+Xbp0OWXfypUrc/jw4TO+n27dujF16lQAvv/+e2rUqEGVKlWIioo6EUNiYiK//vrrWb3uySpVqkT37t257bbbrJdigtrew1kMmrCExZv28UL/Vtx/RbTfTgVbT6UEEBHmzJnDAw88wPPPP094eDhRUVG8+uqrdOvWjSVLltCmTRtEhBdffJE6deowadIkXnrpJcLCwqhUqdKJab8jR46kdevWtGvXjmHDhjFmzBhCQkIICwtj/Pjxpxy7YMGybt26AdC1a1dSU1O54IILTtn3lltuYdSoUZQvX54lS5ac9v0888wz3HrrrbRu3ZoKFSowadIkAPr378/kyZOJiYmhffv2NGnSBIDq1avTpUsXWrZsSZ8+fXjppZfc/t4NHjyY66+//pSZYMYEi81pRxg+8WfSj2Tz7rA4ejSr5dd4xBcXw5QkcXFxevIiXevWraN58+Z+ish4m/18TaBasW0/IyYlECrCxFva06Z+Va8dS0RWqGrcmfaznooxxpRCX63ew/3TV1K3ank+uLW9z6YMn4klFWOMKWUmL9nK05+uoU1kVd4bHkf1SuX8HdIJllSMMaaUyM9XXlywgbd/2MwVzWvz+uC2lC8b6u+w/sCSijHGlALZufk8OiuZuUm7uOnSBvztuksoE1ryJvBaUjHGmBLuUFYOo6as4KfN6Yy5qimjuzcusdUjLKkYY0wJtjvjGLe+v5xNe48w7oY2DIiN9HdIxSp5facgFail77///nsiIiJo27YtzZs359lnn+Xo0aNUr179lNIt/fr1Y+bMmR6N35jSbMOew1z/1k+kHjjG+7e2L/EJBSyplAiBXvo+Pj6elStXkpCQwIcffsj69evp1asXc+fOPbFPRkYGP/74o9vlYYwJdEs2pzPg7Z/Iy1dm3NmR+Oia/g7JLZZUSoBAK31/OhUrViQ2NpbNmzefsh7KnDlz6N279x8qIxsTrD5N3sXwiT9Tu0o4c+7uwiV1Ty0GW1LZmMpJnv1sDWt3HfLoa7aoW4Wn/3TJabcHWun700lPT2fp0qX89a9/JTo6mttvv5309HSqV6/O9OnTuffee8/8zTQmgKkq7yzawnPz19OhYTXeGRpHRIUwf4d1VqynUsKVxtL3J1u0aBFt27alV69ejB07lksuuYSyZcty3XXXMWvWLPbt20dSUhK9evVy/xtjTIDJy1ee/Wwtz81fzzWtL2TybR1KXUIB66mcorgehbcEWun7k8XHx/P555+f0j548GD+8Y9/oKr07duXsLDS9wdkjCdk5eTxwPQkvlqzhzviG/JYn+aEhJTMKcNnYj2VEiDQSt+7q0ePHmzcuJE333zTSteboLX/aDY3vbOUBWv38NS1LXjimhalNqGA9VRKhEArfe+ukJAQ+vfvz8cff3zi+MYEk23pR7nl/eXsOniMt25qR59WF/o7pPNmpe+x0uiBzn6+piRK2nGQER8sJ0+V94bHEXtRNX+HVCwrfW+MMSXUN2t/495pidSsXI4Pbu1A45qV/B2Sx1hSMcYYH5qyxFW2vlW9CN4d3p6alUtO2XpPsKTiUNUSW6DNnLtgO71rSq4/lq2vxWuD21KhbOB9BAfeOzoH4eHhJy7Cs8QSOFSV9PR0wsPD/R2KCXLHc/MY8/EqPk3exc0dG/DMn0pm2XpPsKQCREZGkpqaSlpamr9DMR4WHh5OZGTJL8JnAldGZg4jpySw7Nf9PNq7KXddVnLL1nuCJRVc13Y0bNjQ32EYYwJM6oFMbnl/OdvTM/n3oBj6xtTzd0heZ0nFGGO8ICU1g9smLed4Th6TR3SgY6Pq/g7JJyypGGOMh323fi93f5TIBRXK8tHtlxJd+9TafIHKayNFIjJRRPaKyOoitj0iIioiNZzH3UUkQ0SSnNtThfbtLSIbRGSTiIwt1N5QRJaJyEYRmSEiZb31Xowxxl0fLdvO7ZMTaFijInNGdw6qhALerf31AdD75EYRqQ9cCWw/adMiVY1xbn9z9g0F3gT6AC2AwSLSwtn/BeAVVY0GDgAjvPIujDHGDarKSwvW8/icFOKjazDzzk7UqhJ8Mw+9llRUdSGwv4hNrwCPAu5cQNAB2KSqW1Q1G5gO9BXX1ImeQEFp30lAv/OP2hhjzl52bj4PzUzmze82M6h9fd4dFkfFcsE5uuDTidIich2wU1WTi9jcSUSSReRLESmoP18P2FFon1SnrTpwUFVzT2o/3XFHikiCiCTYtGFjjCdlHMth+MSfmbNyJ4/0asK/rm8VsNeguMNnqVREKgBPAEWtxJQIXKSqR0TkamAuEA0UNZlbi2kvkqpOACaAq6DkWYZujDFFSj2Qya3vL2dr+lFeHRhDv7aBP2X4THyZThsDDYFkEdkKRAKJIlJHVQ+p6hEAVZ0PhDmD+KlA/UKvEQnsAvYBVUWkzEntxhjjE6t3ZvDnt35iz6EsJt3WwRKKw2dJRVVTVLWWqkapahSuhNFOVfeISB1nnAQR6eDElQ4sB6KdmV5lgUHAp+oq6PQdMMB5+eHAPF+9F2NMcPtuw15u/M8SyoaG8MldnencuIa/QyoxvDmleBqwBGgqIqkiUtzsrAHAahFJBl4DBqlLLnAPsABYB8xU1TXOc/4CPCQim3CNsbznrfdijDEFPlq2ndsn/T5luEmQTRk+E1ukyxhj3KCqjPt6A29+t5nuTWvy5k3tgmqGly3SZYwxHnI8N49HZ61iXtIuBneoz9/7tgzqGV7FsaRijDHFKFxleMxVTRndPbCrDJ8vSyrGGHMaO/Zncsv7P7Nj/7GgqTJ8viypGGNMEZJ3HGTEpOXk5ClTRnTg0iCpMny+LKkYY8xJvln7G/dNW0n1SmWZPrIDF9eq5O+QSg1LKsYYU8ikn7by7GdraFUvgneHt6dm5XL+DqlUsaRijDFAfr7y3Px1vPvjr1zRvDavDY6hQln7iDxb9h0zxgS9Y9l5PDgjia/W7GF4p4t46k+XEBpiM7zOxVklFREJASqp6iEvxWOMMT6178hxbp+UQHLqQf56bQtu6xJlU4bPwxmv3hGRj0SkiohUBNYCG0RkjPdDM8YY79q09wh/fmsx6/ccYvyQWEZ0bWgJ5Ty5c0loC6dn0g+YDzQAhno1KmOM8bKlW9LpP/4njmXnMX1kJ3q3rOPvkAKCO0klTETCcCWVeaqag3urNhpjTIk0Z2UqQ99bRo1KZZkzugsx9av6O6SA4c6YytvAViAZWCgiFwE2pmKMKXVUldf/t4mXv/mFjo2q8Z+b44ioEObvsAJKsUnFGZj/TVXrFWrbDvTwdmDGGONJ2bn5PD4nhVkrUrm+bT2e79+asmWsKKSnFZtUVDVfRO4BZhZqUyD39M8yxpiSJSMzh1EfrmDJlnTuuzyaB6+ItgF5L3Hn9Nc3IvIIMAM4WtCoqvu9FpUxxnhIQVHI7fsz+b8b2tA/NtLfIQU0d5LKbc7Xuwu1KdDI8+EYY4znJG4/wB2TEsjNV6aMuJSOVhTS686YVFS1oS8CMcYYT5qfspsHZyRRu0o479/ansY1rSikL7h1Rb2ItARaAOEFbao62VtBGWPMuVJVJizcwr++XE/sRRcwYWgs1StZUUhfOWNSEZGnge64ksp8oA/wI2BJxRhTouTk5fPUvDVM+3k717a+kHE3tCE8LNTfYQUVd3oqA4A2wEpVvVVEagPvejcsY4w5O4eycrh7aiKLNu7j7h6NefjKpoRYUUifcyepHHOmFueKSBVgLzZIb4wpQXbsz2TEpOVsSTvKi/1bc2P7+v4OKWi5k1QSRKQq8A6wAjgC/OzVqIwxxk0rtx/gjskJZOfmM3lEBzo3ruHvkIKaO7O/Rjt33xaRr4AqqrrKu2EZY8yZFZ7hNX1ke1v2twRwp/S9iMjNIvKUqm4FDopIB++HZowxRVNV3vp+E6OnJtKyXgRzRne2hFJCuHP66y0gH+gJ/A04DHwCtPdiXMYYU6Ts3HyenJvCzIRUrmtTlxcHtLYZXiWIO0nlUlVtJyIrAVT1gIiU9XJcxhhzCqvhVfK5k1RyRCQUZw0VEamJq+dijDE+s3XfUW6btJzU/cd4ZWAb/tzWaniVRO4kldeAOUAtEfknrutWnvRqVMYYU8jPv+7nzikJAHx4+6V0aFjNzxGZ03Fn9tdUEVkBXA4I0E9V13k9MmOMAWYnpvKXT1ZR/4IKTLylPVE1Kvo7JFMMd8q0jAPeV9U3fRCPMcYAkJ+vvPzNL7zx3SY6NarO2zfH2iqNpYA7p7/WAxNEpAzwPjBNVTO8G5YxJphl5eTx8MfJfLFqNwPj6vP3fi1tlcZS4ow/JVV9V1W7AMOAKGCViHwkIraksDHG49IOH2fQhKXMT9nNY32a8Xz/VpZQShG3flLO7K9mzm0fkAw8JCLTvRibMSbIbNhzmH5vLmb9nkOMHxLLnZc1tinDpYw7YyovA38C/gc8p6oFdb9eEJEN3gzOGBM8vlu/l3unraRC2VBm3tmJ1pFV/R2SOQfujKmsBp5U1cwitlm5FmPMeVFVPvhpK3//fC3N6lThvVviuDCivL/DMufInSnFE4vZZgP2xphzlpuXz7OfrWXK0m1c0bw2/x4UQ8Vybi1Ia0oo++kZY/wi41gO93zkWlTrzm6NeLR3M0JtUa1Sz5KKMcbntqdnctuk5Wzdd5QX+rdiYPsG/g7JeMhpZ385C3OdMxGZKCJ7RWR1EdseEREVkRrOYxGR10Rkk4isEpF2hfYdLiIbndvwQu2xIpLiPOc1sSkixpQKy7fup99bi0k7fJwpIy61hBJgiptSvE9EvhWREeeYYD4Aep/cKCL1gSuB7YWa+wDRzm0kMN7ZtxrwNHAprkkBT4vIBc5zxjv7FjzvlGMZY0qW2YmpDHlnGRHlw5h7dxc6Na7u75CMhxWXVNYBr+JaR2WziMwTkUEi4ta0DFVdCOwvYtMrwKM4VY8dfYHJ6rIUqCoiFwJXAd+o6n5VPQB8A/R2tlVR1SWqqsBkoJ87cRljfC8/X3nxq/U8NDOZ2IsuYM7ozjS0Gl4BqbikkqOqn6vqECASmArcCKSKyEfncjARuQ7YqarJJ22qB+wo9DjVaSuuPbWIdmNMCZOZncvoqYm89f1mBneoz+QRHahawZZkClTFDdSfGKNQ1WPATGCmiERwDr0CEakAPAH0Ku5Yheg5tJ/u2CNxnSqjQQM7f2uMr+zJyOL2yctZs+sQT17TnBFdG9oV8gGuuJ7K1KIaVTVDVSedw7EaAw2BZBHZiqv3kygidXD1NOoX2jcS2HWG9sgi2oukqhNUNU5V42rWrHkOoRtjzlZKagZ93/yRX9OO8t7wOG6Pb2QJJQicNqmo6jhPHkhVU1S1lqpGqWoUrsTQTlX3AJ8Cw5xZYB2BDFXdDSwAeonIBc4AfS9ggbPtsIh0dGZ9DQPmeTJeY8y5m5+ymxv+8xNlQkL4ZHRnejar7e+QjI947ToVEZkGdAdqiEgq8LSqvnea3ecDVwObgEzgVgBV3S8ifweWO/v9TVULBv/vwjXDrDzwpXMzxviRqvLmd5sY9/UvtGtQlQnD4qhRqZy/wzI+JK7JU8EjLi5OExIS/B2GMQEnKyePsZ+sYm7SLvrG1OWF/q0JDwv1d1jGQ0RkharGnWm/Ykvfi0ioiDzoubCMMYEo7fBxBr+zlLlJu3ikVxNeHRhjCSVIFZtUVDUP1zUkxhhTpHW7D9H3jR9Zt/sQ44e0456e0TYgH8TcGVNZLCJvADOAowWNqprotaiMMaXCN2t/4/7pK6kSHsasUZ1pWS/C3yEZP3MnqXR2vv6tUJviutLeGBOEVJUJC7fw/FfraVUvgneGxVG7Sri/wzIlgDvrqdha9MaYE47n5vHEnNXMWpHKNa0vZNyANpQva+MnxsWd5YQjcBV17OY0/YBraq8t0GVMkEk/cpxRH65g+dYD3H95NA9cYeMn5o/cOf01EdeSwjc6j4cC7wPXeysoY0zJs37PIUZ8kMC+I8d5bXBbrmtT198hmRLInaTSWFX7F3r8rIgkeSsgY0zJ883a33hg+koqlivDzDs70ab+eS23ZAJYsVOKHcdEpGvBAxHpAhzzXkjGmJJCVRn//WZGTkmgca1KfHpPV0soplju9FRGAZOdsRWAA8DwYvY3xgSArJw8Hp+dwuyVO7m29YW8ZAPyxg2nTSoicr+q/huopKptRKQKgKoe8ll0xhi/2Hs4izunrGDl9oM8dGUT7u15sQ3IG7cUd/rrVufr6+BKJpZQjAl8a3Zl0O+NxazffZjxQ9px3+U2w8u4r7jTX+ucdU9qisiqQu0CqKq29mpkxhif+zJlNw/NTKZqhTA+HtXJrpA3Z+20SUVVBzsLaC0ArvNdSMYYX8vPV/793438+78badegKm8PjaVWZbtC3py9YgfqnQW02vgoFmOMH2Rm5/LwzGS+XL2H/u0iee76lpQrYwPy5tx4bZEuY0zJt/PgMe6YlMD6PbaGvPEMSyrGBKmErfu5c8oKsnPzee+W9vRoWsvfIZkAYEnFmCA0c/kOnpibQr2q5Xl3eHsurlXJ3yGZAFHcdSqrTretkDRVvdyD8RhjvCg3L59/zl/H+4u30vXiGrx5UzsiKoT5OywTQIrrqYQCVxezXYBPPRuOMcZbDmZmc89HK/lx0z5u6RzFk9c0p0yoO5WajHFfcUnlTlXdVtyTRWS0h+MxxnjBL78d5o7JCew+mMWLA1pzY1x9f4dkAlRxSSVBROqq6q7CjSJSF0hX1eOq+qN3wzPGnK+CCsPly5Zh2siOxF50gb9DMgGsuL7vO0A7ABG5sFB7O+BdbwZljDl/qsob/9vIyCkJNKpZic/u7WIJxXhdcUmlmap+7tz/QkT6Azhtzb0emTHmnGVm53LPtJWM+/oXrmtTl49HdeLCiPL+DssEgeJOf+WLSGNV3YxrUP4J4BMRiQbyfRKdMeaspR7IZOTkFazbc4jH+jRjZLdGdkGj8ZniksrDwLci8hVQHVARGQ/0Bob5IjhjzNlZuiWd0VMTycnLZ+Lw9vRoZhc0Gt867ekvZxA+FlgGhDu3ZUCcqi7yTXjGGHeoKpOXbOXmd5dxQYUw5t3dxRKK8YszFZTcD3wgIlnO4+k+icoY47bjuXk8NXcNMxJ2cHmzWrwyKIYq4XZBo/GP4q6ov7ZgoP50yaTwPsYY39t7KItRH64gcftB7ulxMQ9d2YSQEBs/Mf5TXE/lJRHZiWuQ/nSeAyypGOMHSTsOcueUBA4dy+WtIe24utWFZ36SMV5WXFL5DXj5DM/f6MFYjDFumrUilcfnpFCrcjlmj+5M8wur+DskY4DiV37s7sM4jDFuyMnL559frOODn7bSuXF13ripHdUqlvV3WMacYKXvjSkl9h05zt1TE1n2635GdG3IY32aWUFIU+JYUjGmFEhJzeDOKQmkH83mlYFt+HPbSH+HZEyRLKkYU8J9siKVx+akULNSOT65qzMt60X4OyRjTuuMfWcRuUFEKjv3nxSR2SLSzvuhGRPccvLyefazNTz8cTLtGlTl03u6WEIxJZ47J2T/qqqHRaQrcBUwCRjv3bCMCW7pR44z9L1lvL94K7d2iWLKiEupXqmcv8My5ozcOf2V53y9BhivqvNE5BnvhWRMcFuVepC7Pkwk7chxXr6xDde3s/ETU3q4k1R2ish/gCuAF0SkHO71cIwxZ2lmwg6enLvaNX4yqjOtIu10lyld3EkONwILgN6qehCoBow505NEZKKI7BWR1YXa/i4iq0QkSUS+dlaRRES6i0iG054kIk8Vek5vEdkgIptEZGyh9oYiskxENorIDBGxyfqm1MrOzeeJOSk8OmsV7aMu4LN7u1pCMaWSO0nlP6o6W1U3AqjqbmCoG8/7AFeZ/MJeUtXWqhqDq7zLU4W2LVLVGOf2NwARCQXeBPoALYDBItLC2f8F4BVVjQYOACPciMmYEmdPRhYDJyxh6rLtjLqsMZNu7WAXNJpSy52kcknhB84HfeyZnqSqC4H9J7UdKvSwIqBneJkOwCZV3aKq2cB0oK+4VhzqCcxy9psE9DtTTMaUND//up9rX/+RDXsO89aQdoy1CxpNKXfa314ReUxEDgOtReSQczsM7AXmnesBReSfIrIDGMIfeyqdRCRZRL4UkYJEVg/YUWifVKetOnBQVXNPajemVFBV3l/8Kze9s5Qq4WWYd3cXKwhpAkJxi3T9S1Ur4zplVcW5VVbV6qr62LkeUFWfUNX6wFTgHqc5EbhIVdsArwNznfaiKiRrMe1FEpGRIpIgIglpaWnnGroxHpGZncuDM5J49rO1dG9ai7n3dCG6dmV/h2WMR7jTz/5cRCoCiMjNIvKyiFzkgWN/BPQH12kxVT3i3J8PhIlIDVw9kPqFnhMJ7AL2AVVFpMxJ7UVS1QmqGqeqcTVr1vRA6Macmy1pR/jzmz8xL3kXD1/ZhAlDY21BLRNQ3Ekq44FMEWkDPApsAyafy8FEJLrQw+uA9U57HWecBBHp4MSVDiwHop2ZXmWBQcCnqqrAd8AA57WGcx6n5Izxha9W7+G6Nxaz93AWk2/rwL2XR9uCWibguHOdSq6qqoj0Bf6tqu+JyPAzPUlEpgHdgRoikgo8DVwtIk2BfFzJaZSz+wDgLhHJBY4Bg5zEkSsi9+Ca0hwKTFTVNc5z/gJMF5F/ACuB99x7y8b4Vm5ePi99vYH//LCFNpERvHVzLPWqlvd3WMZ4hbg+u4vZQeQH4CvgNiAeSAOSVLWV98PzvLi4OE1ISPB3GCZIpB0+zr3TElm6ZT83d2zAX69tQbkyof4Oy5izJiIrVDXuTPu501MZCNwE3Kaqe0SkAfDS+QboNxs2QPfPOxhuAAATxUlEQVTu/o7CBIEVleoyusl1ZISG8/KvX3P90rXwqr+jMsa7zjimoqp7cM3UihCRa4EsVT2nMRVjgoEC79WJZWCLQYTn5zJnzVSu37fW32EZ4xNn7KmIyI24eibf45rK+7qIjFHVWcU+saRq2hS+/97fUZgAdSgrh7/MWsWXq/dwZYvajLuhDxHlb/Z3WMacP3FvUok7p7+eANqr6l7X60pN4Ft+v5rdGAOs2ZXB3VMT2XHgGE9c3Zzb4xsibv4hGhMo3EkqIQUJxZGOVSk25gRVZcbyHTz16RqqVSjLjJEdiYuq5u+wjPELd5LKVyKyAJjmPB4IfOm9kIwpPTKzc3ly7mpmJ+4kProGrw6MscW0TFA7Y1JR1TEicj3QFdeYygRVneP1yIwp4TbtPcLoqSvYuPcID1wRzb09owm1ixlNkDttUhGRi4HaqrpYVWcDs532biLSWFU3+ypIY0qaeUk7eXx2CuXCQpl8Wwfio638jzFQ/NjIq8DhItozsdn2Jkgdy85j7CeruH96Ei3qVuGL+7paQjGmkOJOf0Wp6qqTG1U1QUSivBaRMSXUpr1HuHtqIht+O8zo7o156MomtvaJMScpLqmEF7PNCheZoPLJilSenLua8mVDmXRbBy5rYr0TY4pSXFJZLiJ3qOo7hRtFZASwwrthGVMyZGbn8vS8NXy8IpUODavx2qC21Iko7v8tY4JbcUnlAWCOiAzh9yQSB5QF/uztwIzxt42/HWb01EQ2pR3h3p4Xc//l0Xa6y5gzOG1SUdXfgM4i0gNo6TR/oar/80lkxvhJwcWMz3y2hkrlyjDltkvpGl3D32EZUyq4c53Kd7gWxDIm4B3KyuHx2Sl8vmo3XS6uzis3xlCrip3uMsZd7lxRb0xQSNpxkHunJbLrYBZjrmrKqMsa28WMxpwlSyom6OXnK+8s2sJLCzZQu0o4M+/sSOxFVrvLmHNhScUEtbTDx3n442QW/pJGn5Z1eP761kRUCPN3WMaUWpZUTNBatDGNB2ckczgrh3/0a8mQSxtYqXpjzpMlFRN0snPzGff1BiYs3MLFtSrx4e0daFanir/DMiYgWFIxQWVz2hHun76S1TsPMeTSBjx5TQvKlw31d1jGBAxLKiYoFFx78uxnaykXFsJ/hsZy1SV1/B2WMQHHkooJeAczs3lsdgpfrt5D58bVefnGGCu1YoyXWFIxAW3plnQenJFE2uHjjO3TjJHxjQixa0+M8RpLKiYgZefm8+q3vzD+h81cVK0Cs0d3pnVkVX+HZUzAs6RiAs6mvYd5YEYSq3ce4sa4SJ7+0yVULGe/6sb4gv2lmYCRn69MXrKVf325norlyvD2zbH0bmmD8cb4kiUVExB+O5TFIx8ns2jjPro3rcmLA1pTq7INxhvja5ZU3LQ9PZNjOXk0rVPZ36GYk8xP2c3jc1LIysmzK+ON8TNLKm7Iy1cGTVhCw5oVmXp7R3+HYxyHsnJ4Zt4aZq/cSZvICF4ZGEOjmpX8HZYxQc2WsXNDaIhwS5coFm9KJ3nHQX+HY4AffknjqlcWMjdpJ/ddHs2suzpbQjGmBLCk4qbBHRpQObwMb/+w2d+hBLXDWTmM/WQVwyf+TMVyZZg9ugsPXdmEMFvm15gSwU5/ualyeBjDOl3EW99vZkvaEfuv2A9+3LiPR2cls+dQFqMua8wDV0QTHmZ1u4wpSezfu7NwS+eGhIWGMGHhFn+HElSOHM/l8Tkp3PzeMsLLhjLrrs6M7dPMEooxJZAllbNQs3I5boyLZHbiTn47lOXvcILCT5v2cdUrC5n283ZGdmvE/PviadfgAn+HZYw5DUsqZ2lkfGNy8/OZ+OOv/g4loGVkusZObnp3GWXLhDBrVCcev7q59U6MKeFsTOUsNahegWta12Xqsu2M7nExEeVt6VlPUlXmp+zh6U/XcCAzmzu7NeKBK5rYmifGlBLWUzkHoy5rxJHjuXy4dJu/Qwkouw4e4/ZJCdz9USJ1Isox7+4uPHZ1c0soxpQi1lM5B5fUjaBbk5q8v3grI7o2tFMy5ykvX5myZCsvLdhAvsKT1zTnls5RlLFpwsaUOl77qxWRiSKyV0RWF2r7u4isEpEkEflaROo67SIir4nIJmd7u0LPGS4iG53b8ELtsSKS4jznNfFxXY67LmvMviPHmbUi1ZeHDTjr9xyi//ifeOaztcRGVePrB7txe3wjSyjGlFLe/Mv9AOh9UttLqtpaVWOAz4GnnPY+QLRzGwmMBxCRasDTwKVAB+BpESmY+jPe2bfgeScfy6s6NqpGm/pVmbBwC7l5+b48dEA4nJXDPz5fy7Wv/cj2/Zm8OjCGSbe2p361Cv4OzRhzHryWVFR1IbD/pLZDhR5WBNS53xeYrC5LgaoiciFwFfCNqu5X1QPAN0BvZ1sVVV2iqgpMBvp5670URUS467LGbN+fyZer9/jy0KWaqjI7MZUe437gvcW/ckNcfb596DL6ta1nRSCNCQA+H1MRkX8Cw4AMoIfTXA/YUWi3VKetuPbUItp9qleL2jSqWZHx32/m2tYX2ofiGazZlcHT89aQsO0AMfWrMvGWOFuN0ZgA4/MT16r6hKrWB6YC9zjNRX0a6zm0F0lERopIgogkpKWlnW3IpxUSIozq1pi1uw+xaOM+j71uoDmYmc1T81bzp9d/ZMu+o7zYvzWz77LlfY0JRP4cDf0I6O/cTwXqF9oWCew6Q3tkEe1FUtUJqhqnqnE1a9b0QOi/69u2LrWrlGP891Zo8mS5efl8tGw7Pf/vBz5cuo1hnaL47uHu3Ni+PiEh1qszJhD5NKmISHShh9cB6537nwLDnFlgHYEMVd0NLAB6icgFzgB9L2CBs+2wiHR0Zn0NA+b57p38rlyZUG7v2oglW9JJ3H7AHyGUOKrK12v2cNWrC3l8TgqNa1bk83vjeea6S4ioYBeLGhPIvDamIiLTgO5ADRFJxTWL62oRaQrkA9uAUc7u84GrgU1AJnArgKruF5G/A8ud/f6mqgWD/3fhmmFWHvjSufnF4Esb8J+FW7j3o5XMHt2Z2lWCdxnbFdsO8K/560jYdoBGNSvyn6Gx9GpR28abjAkS4po8FTzi4uI0ISHB46+bkprBwAlLaFCtAjNHdaJKeHD9R7457QgvfbWBr9bsoWblcjxwRTQD4+rb9SbGBAgRWaGqcWfaz66o95BWkRG8fXMst32wnDsmJTDptg5BcaX93kNZ/Pu/G5m+fAfhZUJ46Mom3B7fkApl7VfLmGBkf/ke1K1JTcbd0IYHZiTx0MwkXh/cjtAAHZDenXGMt7/fzLTlO8jPV26+tAH3Xh5NjUrl/B2aMcaPLKl4WL+29Ug7fJx/zl9HzUpreOa6SwJqPGHH/kzG/7CZWQmp5KvSv10ko3s05qLqFf0dmjGmBLCk4gV3dGvE3sNZvLPoV2pVCefuHhf7O6TztnXfUd78bhNzVu4kRIQb20cy6rLGRF5gZVWMMb+zpOIlj/VpTtrh47y0YIOzYmT9Mz+pBFqzK4N3F/3KvKSdhIWGcHPHixh1WWPqRATvDDdjzOlZUvGSkBDhxQFtSD+azWOzU6hRqSw9m9X2d1huOZ6bx5cpe5iydBsrth2gfFgot8c34vb4htSqbMnEGHN6NqXYy44cz2XwhKVs3HuYx69uzuAODQgrodNsdx48xkfLtjFj+Q72HckmqnoFbu54ETfE1reLFo0Jcu5OKbak4gNph49z37SVLNmSTuOaFXnimub0aFqrRAzg5+UrizftY8rSbfx33W8A9GxWm2GdLqLrxTWsnIoxBrCkclr+SCrgKl3y7bq9PDd/Hb/uO0rXi2vwxDXNaX5hFZ/HkpWTx+JN+1iwZg/frtvL/qPZVKtYloHt6zPk0gY2+G6MOYUlldPwV1IpkJ2bz9Rl2/j3fzeScSyHG2Pr83CvJtTycmmXw1k5fLchjQVr9vD9+r0czc6jcrky9GhWi6suqcPlzWsFxcWaxphzY0nlNPydVApkZObw+v82MmnJVsJCQ7itS0Pio2vQsl4EFcud//yJ3w5lsSo1g5TUg6zccZBlW/aTnZdPjUplubJFba66pA6dGlenXBlLJMaYM7OkcholJakU2JZ+lOe/XH9i9cgQgehalWlTP4I29avSJrIqTetUPmVwPy9fOZaTR2Z2Lsey8/h131FSUjNITs0gZedBfjt0/A+v1zW6Br1b1qFdgwsC9ip/Y4z3WFI5jZKWVAqkHznOqtQMklMPkrzjIMmpGew/mg1AuTIhXBgRTlZOPsdy8jiWnUd2Xn6Rr9OoZkXaRFalVb0IWkdG0KJuFavDZYw5b1ZQspSpXqkcPZrVokezWoBrYD/1wDGSUw+StP0gew8fp3xYKOXLOrcw1y3cuV+3ajit6kVQOciqIxtjShZLKiWUiFC/WgXqV6vAta3r+jscY4xxS8m8Cs8YY0ypZEnFGGOMx1hSMcYY4zGWVIwxxniMJRVjjDEeY0nFGGOMx1hSMcYY4zGWVIwxxnhM0JVpEZE0YNtpNtcA9vkwnJIkWN+7ve/gE6zv/Xzf90WqWvNMOwVdUimOiCS4U9smEAXre7f3HXyC9b376n3b6S9jjDEeY0nFGGOMx1hS+aMJ/g7Aj4L1vdv7Dj7B+t598r5tTMUYY4zHWE/FGGOMx1hScYhIbxHZICKbRGSsv+PxBRGpLyLficg6EVkjIvf7OyZfEpFQEVkpIp/7OxZfEpGqIjJLRNY7P/tO/o7JF0TkQef3fLWITBORcH/H5C0iMlFE9orI6kJt1UTkGxHZ6Hy9wBvHtqSC68MFeBPoA7QABotIC/9G5RO5wMOq2hzoCNwdJO+7wP3AOn8H4Qf/Br5S1WZAG4LgeyAi9YD7gDhVbQmEAoP8G5VXfQD0PqltLPBfVY0G/us89jhLKi4dgE2qukVVs4HpQF8/x+R1qrpbVROd+4dxfbjU829UviEikcA1wLv+jsWXRKQK0A14D0BVs1X1oH+j8pkyQHkRKQNUAHb5OR6vUdWFwP6TmvsCk5z7k4B+3ji2JRWXesCOQo9TCZIP1wIiEgW0BZb5NxKfeRV4FMj3dyA+1ghIA953Tv29KyIV/R2Ut6nqTmAcsB3YDWSo6tf+jcrnaqvqbnD9QwnU8sZBLKm4SBFtQTMtTkQqAZ8AD6jqIX/H420ici2wV1VX+DsWPygDtAPGq2pb4CheOg1SkjjjB32BhkBdoKKI3OzfqAKTJRWXVKB+oceRBHDXuDARCcOVUKaq6mx/x+MjXYDrRGQrrlOdPUXkQ/+G5DOpQKqqFvRIZ+FKMoHuCuBXVU1T1RxgNtDZzzH52m8iciGA83WvNw5iScVlORAtIg1FpCyuAbxP/RyT14mI4Dq3vk5VX/Z3PL6iqo+paqSqRuH6Wf9PVYPiv1ZV3QPsEJGmTtPlwFo/huQr24GOIlLB+b2/nCCYoHCST4Hhzv3hwDxvHKSMN160tFHVXBG5B1iAa1bIRFVd4+ewfKELMBRIEZEkp+1xVZ3vx5iM990LTHX+gdoC3OrneLxOVZeJyCwgEdesx5UE8JX1IjIN6A7UEJFU4GngeWCmiIzAlWRv8Mqx7Yp6Y4wxnmKnv4wxxniMJRVjjDEeY0nFGGOMx1hSMcYY4zGWVIwxxniMJRUT9EQkT0SSnOq1HzvXMnwvIledtN8DIvKWF+P4QEQGeOv1jfEFSyrGwDFVjXGq12YDo4BpnFrFdpDTHvCcoovGnDVLKsb80SLgYlzlS64VkXJwouBmXeDHwjuLSJSzLskkEVnlrFNSwdkWKyI/iMgKEVlQqETGHSKyXESSReSTgv1Pet2/Oz2XEBHZKiIviMjPzu1iZ5+LROS/znH/KyINnDVitohLVRHJF5Fuzv6LRORiEanorLex3Ckq2dfZfovTU/sMCLZii8ZDLKkY43D+O+8DpKhqOvAzv69JMQiYoUVfLdwUmKCqrYFDwGinptrrwABVjQUmAv909p+tqu1VtWAtkxEnxfEirgqyt6pqQRXlQ6raAXgDV4VlnPuTneNOBV5T1TzgF1zrAnUFVgDxTnKMVNVNwBO4StO0B3oALxWqVNwJGK6qPc/uu2eMiyUVY1xrbCQBCbjKV7zntBc+BVbcqa8dqrrYuf8hrg/zpkBL4BvntZ/EVagUoKXTa0gBhgCXFHqtvwJVVfXOkxLYtEJfC1Zq7AR85Nyf4hwXXL2tbs7tX057e1w17gB6AWOduL4HwoEGzrZvVPXkdTiMcZudNzXGGVMpon0u8LKItAPKFyxoVoSTey+KazmFNapa1FK9HwD9VDVZRG7BVaOpwHIgVkSqnfThrqe5X1Qci3CNC9UFngLGOMdY6GwXoL+qbij8ZBG5FFcpfGPOmfVUjDkNVT2C6z/5iRQ/QN+g0Drvg3GNu2wAaha0i0iYiBT0SCoDu51TZENOeq2vcBX++0JEKhdqH1jo6xLn/k/83pMawu/jPctwlXXPV9UsIAm4E1eyAVfh1Hudar2ISNti3psxZ8WSijHFm4ZrHffpxeyzDhguIquAargWwMoGBgAviEgyrg/2gvU7/orrg/8bYP3JL6aqHwPvAJ+KSHmnuZyILAPuBx502u4DbnWOO9TZhqoex7WS6VJnv0W4ElmK8/jvQBiwSkRWO4+N8QirUmzMeXBmhX3uTEf21jG2AnGqus9bxzDGU6ynYowxxmOsp2KMMcZjrKdijDHGYyypGGOM8RhLKsYYYzzGkooxxhiPsaRijDHGYyypGGOM8Zj/B3KjvecA+T/kAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "total_costs = []\n",
    "# For all possible PV types\n",
    "for i in range(data_production.shape[1]):\n",
    "    # Set pv type of house to this specific one\n",
    "    g.change_pvtype(num_house=0,pv_type=i)\n",
    "    # Simulate the grid for one year\n",
    "    res_dict = g.simulate(data_cons=data_consumption,data_prod=data_production)\n",
    "    # Objective function\n",
    "    total_costs.append(objective_function(res_dict))\n",
    "    \n",
    "# Plot the total costs agains the peakpower of the used module\n",
    "plt.axhline(total_cost_noPV,c='r',label='Costs without PV')\n",
    "plt.xlabel('PV peakpower')\n",
    "plt.ylabel('Costs [€] for 7 years')\n",
    "plt.plot(pv_peakpowers,total_costs,label='Costs with PV')\n",
    "plt.legend()\n",
    "plt.savefig('plots/single_household_7year.png',dpi=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The downside of this specific objective function is, that it depends on the timespan. If the timespan is too large (e.g. 15 years), the total cost only decreases with the peakpower of the installed pv module. Hence the largest module is the best. \n",
    "Otherwise, if the timespan is too small, the costs only increase with peakpower and no pv module is the cheapest choice.\n",
    "\n",
    "# Multiple households\n",
    "Now we are going to do the same as above just for multiple households that all have different consumption schedules. In order to cut down the dimensions we need to use later on we only take the first 20 PV modules as available. These correspond to sizes from $0.25$ kWp to $5.00$ kWp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose number of houses houses\n",
    "num_houses = 5\n",
    "data_consumption = df_cons.iloc[:,:num_houses].values\n",
    "\n",
    "# Retrieve data from production and pv_sizes\n",
    "df_prod = df_prod.iloc[:,0:20]\n",
    "data_production = df_prod.values\n",
    "pv_peakpowers = [float(key) for key in df_prod.keys()]\n",
    "\n",
    "# Create the grid model and simulate one year\n",
    "g = Grid(num_houses=num_houses,num_storages=1,max_capacity=0,num_pvtypes=len(pv_peakpowers),pv_peakpower=pv_peakpowers)\n",
    "# Set standard costs\n",
    "g.set_costs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we randomize the attached pv module to each house and simulate for one year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'import_grid_kwh': 22405.888250000233, 'export_grid_kwh': 11918.947249999948, 'pv_production_kwh': 15951.254999999815, 'setup_cost_storage': 0, 'setup_cost_pv': 14350.0, 'cost_import_grid': 5601.472062500058, 'reward_export_grid': 1191.8947249999949}\n"
     ]
    }
   ],
   "source": [
    "g.randomize()\n",
    "res_dict = g.simulate(data_cons=data_consumption,data_prod=data_production)\n",
    "print(res_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now this it the point, where we can start implementing a **genetic algorithm** to optimize the size of the pv module attached to each house. Since all houses do not interact (YET!) we could also simply optimize each house on its own and would result in the optimal solution. However the houses will interact through the storage systems later and we can not optimize each house on its own anymore. Thus we already build the genetic algorithm here. \n",
    "\n",
    "First we initialize a starting population of grids and randomize their initial states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_population = 100\n",
    "population = []\n",
    "total_costs = np.array([])\n",
    "for i in range(num_population):\n",
    "    # Create Grid, set costs, randomize and simulate\n",
    "    g = Grid(num_houses=num_houses,num_storages=1,max_capacity=0,num_pvtypes=len(pv_peakpowers),pv_peakpower=pv_peakpowers)\n",
    "    g.set_costs()\n",
    "    g.randomize()\n",
    "    res_dict = g.simulate(data_consumption,data_production)\n",
    "    # Append grid to population and total costs to array\n",
    "    population.append(g)\n",
    "    total_costs = np.append(total_costs,objective_function(res_dict)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to create crossovers and mutations. We cut down to the original population size and start again. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0e828e103ce84ebea380a7d0be5982c5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "IntProgress(value=0, max=20)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Rates of crossover and mutation\n",
    "mutation_rate = 0.5\n",
    "crossover_rate = 0.5\n",
    "\n",
    "# Storage arrays\n",
    "num_generations = 20\n",
    "mean_total_costs = np.zeros(num_generations)\n",
    "std_total_costs = np.zeros(num_generations)\n",
    "min_total_costs = np.zeros(num_generations)\n",
    "\n",
    "# Progress bar\n",
    "f = IntProgress(min=0, max=num_generations-1) # instantiate the bar\n",
    "display(f) # display the bar\n",
    "\n",
    "for i in range(num_generations): \n",
    "    # Mutations\n",
    "    for k in range(int(mutation_rate*num_population)): \n",
    "        idx = random.randint(0,len(population)-1)\n",
    "        g = population[idx].get_copy()\n",
    "        g.mutate(storage_connection=False,pv_type=True)\n",
    "        res_dict = g.simulate(data_consumption,data_production)\n",
    "        population.append(g)\n",
    "        total_costs = np.append(total_costs,objective_function(res_dict))\n",
    "    \n",
    "    # Crossover\n",
    "    for k in range(int(crossover_rate*num_population)):\n",
    "        idx1 = random.randint(0,len(population)-1)\n",
    "        idx2 = random.randint(0,len(population)-1)\n",
    "        g1 = population[idx1].get_copy()\n",
    "        g2 = population[idx2].get_copy()\n",
    "        g2 = g1.crossover(g2,storage_connection=False,pv_type=True)\n",
    "        population.append(g1)\n",
    "        population.append(g2)\n",
    "        res_dict = g1.simulate(data_consumption,data_production)\n",
    "        total_costs = np.append(total_costs,objective_function(res_dict))\n",
    "        res_dict = g2.simulate(data_consumption,data_production)\n",
    "        total_costs = np.append(total_costs,objective_function(res_dict))\n",
    "        \n",
    "    # Size reduction\n",
    "    del_idx = np.argsort(total_costs)[num_population:]\n",
    "    for idx in sorted(del_idx,reverse=True):\n",
    "        population.pop(idx)\n",
    "    total_costs = np.delete(total_costs,del_idx)\n",
    "    \n",
    "    # Calculate mean and min total costs\n",
    "    mean_total_costs[i] = np.mean(total_costs)\n",
    "    std_total_costs[i] = np.std(total_costs)\n",
    "    min_total_costs[i] = np.min(total_costs)\n",
    "    \n",
    "    # Update value of the proress bar\n",
    "    f.value = i\n",
    "\n",
    "# Store the best grid individual\n",
    "best_grid = population[np.argmin(total_costs)]\n",
    "min_cost_genetic = np.min(min_total_costs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XecFdX5x/HPs5221KXsrrogCIIg4AoqKtgVkWKJGGNssdckxhJrjCYaf0mMRJMYNRq7UaPE8kKiErvSUUABAaX3Xrc8vz9mFq/L7rKw997Z8n2/Xve1c8+cmfvMbHn2zDlzxtwdERGReEiJOgAREak/lFRERCRulFRERCRulFRERCRulFRERCRulFRERCRulFSkzjGzv5rZrdWs+7iZ3VXFejezznGKqyDcX1o89reLzxpnZj+Jw36SFrM0DEoqknRmNt/MlplZk5iyn5jZuOps7+6XuvuvExZgHWdm55lZiZltjHkNijouaRiUVCQqacA1UQcRL7XwP/2P3b1pzGtc1AFJw6CkIlG5D7jOzFpUtNLMupnZWDNbbWZfmdkPYtZ975KWmV1vZkvMbHHY4il/Saulmb1uZhvM7FMz27fcxw02s7lmttLM7jOzlHC/KWZ2i5l9Y2bLzeyfZtY8XFd22ehCM/sWeCdmf2eb2bfh/m6OiTPTzO4P41wcLmfGrL/IzOaExzzazHJj1h1nZl+a2Toz+zNgu3Oyq2G3Yw5bRB/E7iT23JvZYDObEZ73RWZ2XUy9IWY2xczWmtlHZtYrZt0NYf0N4ff+mDgfqySQkopEZQIwDriu/IrwsthY4BmgLXAW8JCZ9aig7onAz4Bjgc7AwAo+6yzgV0BLYA5wd7n1I4BCoC8wDLggLD8vfB0FdAKaAn8ut+1AYH/ghJiyw4GuwDHAbWa2f1h+M3AI0Bs4EOgH3BIex9HAb4EfAB2Ab4DnwnVtgJfCum2Ar4EBFRxnrD5hgphlZrdWoyW12zFXw6PAJe7eDDiAMPGaWV/gMeASoDXwN2B0mMC6AlcCB4fbnQDMr+bnSW3g7nrpldQXwR+JYwn+0KwDcoCfAOPC9WcC75fb5m/A7eHy48Bd4fJjwG9j6nUGHOgcU/eRmPWDgS9j3jtwYsz7y4G3w+W3gctj1nUFiggu3RWE23aKWV9Wlh9T9hkwMlz+Ghgcs+4EYH64/Cjwu5h1TcPPKgB+DHwSs86AhcBPKjm/nYCOBP809gRmADdVUrcmMZ8HfFBuf7Hn/luCxJFdrs5fgF+XK/uKIEF3BpaHPx/pUf+s6rX7L7VUJDLu/gXwGnBjuVX7AP3DSyNrzWwtcDbQvoLd5AILYt4vqKDO0pjlzQR/sGPFbvNNuM+yfX9Tbl0a0G4PP6+i/VX4We6+EVgF5FHuGD34K1zR55atn+vu89y91N0/B+4ETq+sfg1i3pXTCJL4N2b2PzM7NCzfB/h5ue/vXkCuu88BrgXuAJab2XOxlwGl9lNSkajdDlxE8MezzALgf+7eIubV1N0vq2D7JUB+zPu99iCG2G32BhaHy4sJ/gDGrisGlsWU7c403xXtr8LPCi8BtgYWERzjXjHrjN07TmfP+2CqinkT0Dgmru8lfXcf7+7DCC5hvgK8EK5aANxd7vvb2N2fDbd7xt0PDz/XgXv3MHaJgJKKRCr8z/R54OqY4teA/czsHDNLD18Hx1znj/UCcL6Z7W9mjYHb9iCMX5hZSzPbi2BE2vNh+bPAT82so5k1BX4DPO/uxXvwGWX7u8XMcsJ+ktuAp8J1z4TH0TvsCP8N8Km7zwdeB3qY2alh38jVVNxqA8DMTjKzduFyN+BW4NUExDw1jKu3mWURtC7KYsgws7PNrLm7FwHrgZJw9d+BS82svwWamNnJZtbMzLqa2dHhOdgKbInZTuoAJRWpDe4Edtyz4u4bgOOBkQT/FS8l+G81s/yG7v4m8ADwLkEn/Mfhqm278fmvAhOBKQR/wB8Nyx8DngTeA+YR/JG7ajf2W95dBAMUpgGfA5PCMtz9bYI//i8RtEz2JTh+3H0lcAZwD8ElsS7Ah1V8zjHANDPbBLwBvEyQpOId8yyC791/gdnAB+W2PQeYb2brgUuBH4XbTSBonf4ZWEPwfTsv3CYzPM6VBN/3tsAv9zB2iYAFl2dF6oewNfMFkFmDFoWI7CG1VKTOM7MR4eWWlgQtmv8ooYhEQ0lF6oNLgBUEw19LgIo69EUkCXT5S0RE4kYtFRERiZvaNglewrVp08YLCgqiDkNEpM6YOHHiSnfPqU7dBpdUCgoKmDBhQtRhiIjUGWb2za5rBXT5S0RE4kZJRURE4kZJRURE4qbB9amIiOxKUVERCxcuZOvWrVGHklRZWVnk5+eTnp6+x/tQUhERKWfhwoU0a9aMgoICgkmh6z93Z9WqVSxcuJCOHTvu8X50+UtEpJytW7fSunXrBpNQAMyM1q1b17h1pqQiIlKBhpRQysTjmJVUqmFrUQl/f28uH81ZGXUoIiK1mpJKNaSlGA+/P5fHP5ofdSgi0kCYGeecc86O98XFxeTk5DBkyJAIo9o1JZVqSEtNYXjvXN79ajlrNm2POhwRaQCaNGnCF198wZYtWwAYO3YseXl5u9gqekoq1TSiTz5FJc5r0xbvurKISBycdNJJvP766wA8++yznHXWWTvWbdq0iQsuuICDDz6YPn368OqrwROj58+fzxFHHEHfvn3p27cvH330EQDjxo1j0KBBnH766XTr1o2zzz6bRMxSryHF1dQ9N5tu7Zvx0qRFnHNoQdThiEiS/Oo/05mxeH1c99k9N5vbT+mxy3ojR47kzjvvZMiQIUybNo0LLriA999/H4C7776bo48+mscee4y1a9fSr18/jj32WNq2bcvYsWPJyspi9uzZnHXWWTvmO5w8eTLTp08nNzeXAQMG8OGHH3L44YfH9djUUtkNp/XNZ8qCtcxdsTHqUESkAejVqxfz58/n2WefZfDgwd9b99Zbb3HPPffQu3dvBg0axNatW/n2228pKirioosuomfPnpxxxhnMmDFjxzb9+vUjPz+flJQUevfuzfz58+Mes1oqu2FY71x+++ZM/j15ET8/vmvU4YhIElSnRZFIQ4cO5brrrmPcuHGsWrVqR7m789JLL9G16/f/Ft1xxx20a9eOqVOnUlpaSlZW1o51mZmZO5ZTU1MpLo7/U7fVUtkNbbOzOLxLDi9PWkRpqZ6YKSKJd8EFF3DbbbfRs2fP75WfcMIJjBo1ake/yOTJkwFYt24dHTp0ICUlhSeffJKSkpKkxquksptO65vHorVbGD9/ddShiEgDkJ+fzzXXXLNT+a233kpRURG9evXigAMO4NZbbwXg8ssv54knnuCQQw5h1qxZNGnSJKnxNrhn1BcWFnpNHtK1ZXsJhXeNZUivXO49vVccIxOR2mLmzJnsv//+UYcRiYqO3cwmunthdbZXS2U3NcpI5aSeHXjj8yVsLUpus1JEpLZTUtkDp/bJY8O2YsbOWBZ1KCIitUrCk4qZpZrZZDN7rVz5KDPbaWyumZ1uZm5mhTFlN5nZHDP7ysxOiCk/MSybY2Y3JvZIvnNIp9bkNs/i5UkLk/WRIiJ1QjJaKtcAM2MLwoTRonxFM2sGXA18GlPWHRgJ9ABOBB4KE1Uq8CBwEtAdOCusm3ApKcbwPnm8N3slKzZsS8ZHiojUCQlNKmaWD5wMPBJTlgrcB1xfwSa/Bn4HxE7oPwx4zt23ufs8YA7QL3zNcfe57r4deC6smxSn9s2jpNQZPVXTtoiIlEl0S+V+guRRGlN2JTDa3ZfEVjSzPsBe7v69y2RAHrAg5v3CsKyy8qTo3LYZvfKb6xKYiEiMhCUVMxsCLHf3iTFlucAZwKhydVOAPwI/r2hXFZR5FeUVxXKxmU0wswkrVqyo5hHs2ql98pi+eD1fLd0Qt32KiMCup74fPXo099xzT1ThVSqRLZUBwFAzm09waepoYDrQGZgTljc2szlAM+AAYFxYfggwOux7WQjsFbPffGBxFeU7cfeH3b3Q3QtzcnLidoCnHJhLWorx8mS1VkQkvnY19f3QoUO58cakjU+qtoQlFXe/yd3z3b2AoKP9HXdv6e7t3b0gLN/s7p3dfZ27t4kp/wQY6u4TgNHASDPLNLOOQBfgM2A80MXMOppZRvgZoxN1PBVp3TSTQV1zeGXyIko0bYuIxFlVU98//vjjXHnllQCcd955XH311Rx22GF06tSJF198MZJ4oQ5MKOnu083sBWAGUAxc4e4lAGZ2JTAGSAUec/fpyY7v1L75/Hfmcj7+ehWHd2mT7I8XkUR780ZY+nl899m+J5y060tXVU19X96SJUv44IMP+PLLLxk6dCinn356fGOupqQkFXcfB4yroLxpJfUHlXt/N3B3BfXeAN6IR4x76uhubWmWlcbLkxYqqYhIXFU19X15w4cPJyUlhe7du7NsWXQ3Ztf6lkptl5WeypBeubwyeRG/Hl5Mk0ydUpF6pRotikSqbOr78mKntY9yTkdN0xIHp/bNY0tRCWOmL406FBGpZyqb+r62UlKJg8J9WrJXq0a8PGlR1KGISD1T2dT3tZWmvo+TP46dxQPvzObjG4+hffOsXW8gIrWWpr7X1PeRG9EnD3d4ZYpaKyLScCmpxElBmyYctE9LXp60MNJOMhGRKCmpxNGpffOYtWwj0xevjzoUEamhhvjPYTyOWUkljob0zCUjNUUd9iJ1XFZWFqtWrWpQicXdWbVqFVlZNesT1k0VcdS8cTrH7N+W0VMX8cvB3UhLVc4WqYvy8/NZuHAh8ZyAti7IysoiPz+/RvtQUomzEX3yePOLpbw/eyVHdWsbdTgisgfS09Pp2LFj1GHUSfpXOs4GdW1Ly8bpvKTnrIhIA6SkEmcZaSkMPTCXt2YsY/3WoqjDERFJKiWVBBjRN5/txaW8+fmSXVcWEalHlFQS4MD85nTKacJLGgUmIg2MkkoCmBmn9c3ns3mrWbB6c9ThiIgkjZJKggzrnQvAK5PVWhGRhkNJJUHyWzbmkE6teHnyogZ1A5WINGxKKgl0at985q3cxJQFa6MORUQkKZRUEuikA9qTmaZpW0Sk4VBSSaBmWemc0KM9/5m2mO3FpVGHIyKScEoqCTaibx5rNxfx7lfLow5FRCThlFQS7IjObWjTNJOXNW2LiDQASioJlpaawvDeubzz5XLWbNoedTgiIgmlpJIEI/rmUVTiPPLBXA0vFpF6TUklCbp3yOaEHu148N2vOf/x8SxbvzXqkEREEkJJJQnMjL+cfRC/GtqDT+au4vg/vsfoqYujDktEJO6UVJIkJcU497AC3rj6CDq2acLVz07mimcmqZ9FROoVJZUk65TTlBcvPZRfnNCVt6Yv5fj73+PdLzXcWETqByWVCKSlpnDFUZ155YoBtGqcwfmPj+fGl6axcVtx1KGJiNSIkkqEeuQ2Z/RVA7h04L68MGEBJ97/Hp/MXRV1WCIie0xJJWKZaanceFI3XrjkUFJTjLP+/gl3vTaDrUUlUYcmIrLblFRqicKCVrxx9RGc3X9vHvlgHkNGfcC0hZrdWETqFiWVWqRJZhp3De/JPy/ox8atxYx46CP+OHYWRSWajFJE6gYllVroyP1yGHPtkZzSqwN/ens2pz70EbOXbYg6LBGRXVJSqaWaN07n/pF9eOjsvixcs5mTR33AH8bOYt3moqhDExGplJJKLTe4ZwfG/PRIjtu/HQ+8PZsB977DfWO+ZLVumhSRWijhScXMUs1sspm9Vq58lJltjHn/MzObYWbTzOxtM9snZt25ZjY7fJ0bU36QmX1uZnPM7AEzs0QfTxTaNsviwbP78uY1RzBwvxweGvc1h9/7Dr99cyYrN26LOjwRkR2S0VK5BpgZW2BmhUCLcvUmA4Xu3gt4EfhdWLcVcDvQH+gH3G5mLcNt/gJcDHQJXycm6Bhqhf07ZPPg2X1569ojOa57O/7+3lwOv/cdfv3aDJZrkkoRqQUSmlTMLB84GXgkpiwVuA+4Prauu7/r7pvDt58A+eHyCcBYd1/t7muAscCJZtYByHb3jz2YT/6fwPBEHk9t0aVdM/40sg9jfzaQwT078PhH8zn8d+9y+6tfsGTdlqjDE5EGLNEtlfsJkkfsmNgrgdHuvqSK7S4E3gyX84AFMesWhmV54XL58p2Y2cVmNsHMJqxYsWL3jqAW2zenKX/4QW/e+flARvTO4+lPv2Xg78Zx878/Z+GazbvegYhInCUsqZjZEGC5u0+MKcsFzgBGVbHdj4BCgtYMQEX9JF5F+c6F7g+7e6G7F+bk5FTzCOqOfVo34d7Te/HudYM4ozCfFyYsYNB947jhxWl8s2pT1OGJSAOSlsB9DwCGmtlgIAvIBqYD24A5YZ96YzOb4+6dAczsWOBmYKC7l/VALwQGxew3HxgXlueXK2/QDynZq1Vj7h7RkyuO6szf/vc1z45fwIuTFjK8dx5XHLUvnXKaRh2iiNRzlozH25rZIOA6dx9SrnyjuzcNl/sQdNCf6O6zY+q0AiYCfcOiScBB7r7azMYDVwGfAm8Ao9z9japiKSws9AkTJsTnwGq55eu38rf35vL0p9+wvbiUA/KaU9C6CQVtmtCxTWMKWjehY5smtGicEXWoIlKLmdlEdy+sTt1EtlR2131AU+BfYSvmW3cfGiaPXwPjw3p3uvvqcPky4HGgEUEfzJvIDm2zs7h1SHcuHbgv//x4PlMWrGXygjW8Nm0xpTH/S7RonM4+rZvQsXXjMOE02ZF8mjdKjyx+Eal7ktJSqU0aUkulMtuKS1iwegvzV25i/qpNzAu/zl+5mcXrthD7I9GqSQYFYbLJb9mY9tlZtG+eSbvsLNpnZ9GqSQb19PYgEQnV1ZaKJElmWiqd2zalc9ud+1i2FpXw7erNQaKJSTofzVnFsg2LKP8/SEZqCm2zM2mfnUW75kGiaZ+dtaOsffMs2mVnkZWemqSjE5EoKanI92Slp7Jfu2bs167ZTuuKSkpZsWEby9ZvZdn6rSxdt5Wl67ftWJ65eD3vfrmczdt3fhZMi8bp9NmrBcN653Fc93Y0ydSPnkh9pN9sqbb01BRyWzQit0WjSuu4Oxu2FbNs3VaWhslm2fqtLFq7lfdmreDa56eQlZ7Csfu3Y+iBuQzsmkNmmloxIvWFkorElZmRnZVOdlY6Xcq1dkpLnUnfruHVKYt5/fMlvDZtCdlZaQzu2YGhB+bSv1NrUlPUPyNSl6mjXiJRVFLKh3NWMnrqYsZ8sZRN20to2yyTIb1yGdY7l175zTUAQKSW2J2OeiUVidzWohLenrmc0VMX8e6XK9heUkpB68YMPTCXob1z6dx25/4dEUkeJZUqKKnUbuu2FDFm+lJGT1nMR1+vpNShe4dshvXO5bju7ejYpolaMCJJpqRSBSWVumP5hq28Pm0Jr05ZzJQFawHIaZZJ/46t6N+pNYd0bEXntk2VZEQSTEmlCkoqddOC1Zv5YM5KPp27ik/nrWbJuuD5Ma2bZNCvY6sdiaZru2akqLNfJK7imlTM7FSC55lsMLMbCebg+o27T6l5qMmnpFL3uTsLVm/hk3mr+HTuaj6dt4qFa4LnyDRvlM7BBa04pFMr+ndsTffcbI0oE6mheN9Rf4e7v2xmhwGnAH8A/gocUoMYRfaYmbF368bs3boxPyjcC4BFa7cErZgwyfx35jIAmmWmUVjQkv6dWnPYvq05ILe5WjIiCVSdpFJ2e/QQ4CF3f8nMbklgTCK7La9FI07tm8+pfYOnISxdt5VP5wWXyj6du4p3vwoeztaqSQaHd27DkfvlcGSXNrTNzooybJF6pzpJZYmZPUjw/PdCM8sgOc+2F9lj7ZtnMax3HsN6Bw8DXbFhGx/OWcl7s1bw3uzg/hiAbu2bhQkmh8KClpqjTKSGqtOn0hQYDExz9y/Dpzce6O51cpp59alIaakzc+l63p8dJJkJ89ewvaSUrPQUDunUmiO65DBwvzbsm6ORZSIQ/476x939vF2V1RVKKlLe5u3FfDJ3Fe/NCpLM3JXBI5hzm2cFrZj9chiwbxuaN9azZaRhindHfa9yO08BDt6TwERqo8YZaRzdrR1Hd2sHBMOXy1oxr3++hOfGLyDFoH12FlkZqWSlpdIoI5Ws9BSy0lJjylJi1qWSmZYSLIdlGakppKUa6akppKUYaakppKcaaSnh17A8vaxeSvC1bFkDDKQuqDSpmNkNwI1AMzMre9KiAQ48moTYRCKxV6vG/LD/3vyw/94Ul5QyZcFa3pu9ksVrt7ClqIRtRSVsKSpha1EpazcXhWWlYVmwLhG3fzXJSKVddhY5zTJpm51Fu2aZtM3O/K6sWRbtsjNpmpmmy3YSmUovf1nwU5kK/JYguQDg7js/LKMO0eUvSTR3Z3tJKVuLStkak2i2F5dSVOIUl5RSXOoUlZRSXOIUl4blZV+/t/xd3fVbilm+YSvL129j+YatLFu/jS1FO/86NkpPpV12kGTaxnxtlpVGqhkpKUaqGakphhmkhu+rLE8xUszUsmqg4nL5y4NsU2xmLwGZ7r7ZzM4ysz7AKHdfEKd4ReoVMyMzLZXMtFSaN0pcP4y7s3FbMcvDB6eVPUAtSDrB8ozF63l3/XI2VfDgtERLMUgrSzQpRt99WvK703vRtpmGcddn1emonwYcCPQEngYeB4a6+8CER5cAaqlIQ7RxWzGbtxVT4k5JqVNaynfLMV8rLC91Sh2KS3evZVVc4hSF22wpKuHlSQvJzkrnobP7UljQKupTIrsh3h31xe7uZjYM+JO7P2JmZ9csRBFJpqaZaTSN+BHO5xyyD5c+NZGRD3/CLwfvz/kDCtT3Uw9V5ybGTWb2C+Ac4PVw9JfGVorIbtm/QzajrzycQV3bcudrM7jq2cls2lYcdVgSZ9VJKmcSjPq6xN2XAPkE83+JiOyW5o3Seficg/jFCV154/MlDH/wQ+Ys3xh1WBJHu0wq7r4YeAzINLMTgc3u/o+ERyYi9VJKinHFUZ158sL+rN60nWF//oA3Pl8SdVgSJ7tMKmZ2GjCJ4PLXj4EJZjYi0YGJSP02oHMbXrv6cLq0a8blT0/i7tdnUFRSGnVYUkPV6bm7DTjY3ZcBmFk74C3g34kMTETqvw7NG/HCJYdy1+sz+Pv785i6YB1/PruPhh3XYdXpU0kpSyihFdXcTkRklzLSUrhz2AHcf2Zvpi1ay8kPfMD4+at3vaHUStVJDm+Z2Rtm9iMz+xEwmqClIiISN8P75PHKFQNokpHKyIc/4ZH359LQHndeH1QnqVxHcMNjP6A/8ERYJiISV93aZzP6qsM5pltb7np9Jlc+O5mNGnZcp1SaVMysk5kd6oEX3P1qd78KWAIUJC1CEWlQsrPS+ds5B3HjSd14c8ew4w1RhyXVVFVL5U/AlgrKtwP3JyYcEZFg/rRLB+7LUxf2Z82m7Qz784e8MH4B24s1Oqy2qyqpdHT3KeUL3f0zoGPiQhIRCRzWuQ2vX30E3Tpkc/1L0zj83ncY9fZsVm3cFnVoUomqkkpmFesaxzsQEZGKtG+exb8uOZR/nHcwXds34/djZ3HoPe/wi39NZcbi9VGHJ+VUdZ/KJDM7v/zd82Z2HjA5oVGJiMRISTGO6taWo7q1Zc7yDfzjw/m8PGkR/5q4kP4dW3H+gAKO696eVD3DJXJVPaSrA/AKsAGYGBYXAs2AYeE8YHWOpr4XqR/WbS7i+Qnf8sRH37Bo7RbyWjTi3MP24czCvWneWHPextPuTH1fneepHAccEL6d7u51+h4VJRWR+qW4pJT/zlzGYx/O57N5q2mUnsppB+Vx3mEd6dy2adTh1QtxTSpxCCYVmAAscvchMeWjgPPdvWn4PhP4J3AQsAo4093nh+tuAi4ESoCr3X1MWH4iwSi1VOARd79nV/EoqYjUX9MXr+PxD+fz6tTFbC8u5cj9cjh/QAEDu+To8cY1EO+HdNXUNcBMILuswMwKgRbl6l0IrHH3zmY2ErgXONPMugMjgR5ALvBfM9sv3OZB4DhgITDezEa7+4yEHo2I1Fo9cptz3xkHcsNJ3Xj202958pNvOP8f4+nUpgnHdm9HWoqRYoZZMGw5xSAl/GpheUpMuZlhBI9GzkpPpVWTDFo3zaBVk0xaN82gWWaaHjRWTkKTipnlAycDdwM/C8tSgfuAHwKxsx0PA+4Il18E/mzBd2sY8Jy7bwPmmdkcgrv7Aea4+9xwv8+FdZVURBq4Nk0zueqYLlwycF/e/GIJ//hwPo9/OB8neDRyqTvxuEiTnmq0ahImmR0JJ4PWYVmrJhm0CcvyWjYiMy215h9ayyW6pXI/cD1B536ZK4HR7r6kXIbPAxYAuHuxma0DWofln8TUWxiWUVY/prx/RUGY2cXAxQB77733nh6LiNQxGWkpDOudx7DeeRWud/8uyZQlmtivpf5dnS1FJazeuJ1Vm7axauN2Vm/azqpN21m9aRurN21n5cbtfLt6M6s3ba9wapnObZsy+soBNM6I9rHOiVbp0ZnZGqCiXG6Au3urqnZsZkOA5e4+0cwGhWW5wBnAoEr2W55XUV7RPTYV/u/h7g8DD0PQp1JV3CLScJgZqQapFf6Z2Vlei0bVqre1qIQ1m7ezamOQeL5evpE7X5vBH8fO4uaTu9ck5FqvqpTZpob7HgAMNbPBQBZBn8p0YBswJ2ylNDazOe7emaClsRew0MzSgObA6pjyMvnA4nC5snIRkchkpafSoXkjOjQPktDA/XKYs2Ijj34wj1MOzKVXfvku5fqj0jvq3b0k9kXwR75dzKtK7n6Tu+e7ewFBR/s77t7S3du7e0FYvjlMKBBMqX9uuHx6WN/D8pFmlmlmHYEuwGfAeKCLmXU0s4zwM0bv9hkQEUmCG0/qRpummdzw0uf1+gmXu7y4Z2YnA38kaAmsIujPmAV0i3MsjwJPhh3xqwmSBO4+3cxeIOiALwauCJMcZnYlMIZgSPFj7j49zjF958MHoLQoYbtPPIMeI6CVpm0TiUJ2Vjq/Hn4Alzw5kb+/P5fLB3Xe9UZ1UHVufpxCMGz3LXfvE94MeZq7X5qMAONtj+9TubsDFG2Of0DJ1GpfuOxDSK/edWEuyWVoAAAWA0lEQVQRib/LnprI218uZ8y1R9KxTZOow6mWeN+nUuzuK8wsxczM3cea2d01jLHuuWF+1BHUzPwP4KlT4X/3wrF3RB2NSIP1q6E9+GDOSm58aRrPXnRIvbspszpJZZ2ZNQE+AP5pZsuB+ntBsDJpVU3aXAd0Pgb6/Ci4jNdjBHQ4MOqIRBqkttlZ3Dx4f258+XNemLCAkf3q120O1Xmc8HBgK3AtMA5YBAypagOppY6/Cxq3htFXQYke0SoSlTMP3otDOrXi7jdmsnz91qjDiavqJJWbwhFgRe7+qLv/gfDueKljGrWEwffBkqnwyYNRRyPSYJkZvz21F9uKS7l9dOLGF0WhOknlxArKTo53IJIk3YdB15Ph3d/Aqq+jjkakwerYpgnXHtuFN79YypjpS6MOJ24qTSpmdomZTQa6mtmkmNdsNL9W3WUGJ/8fpGbAa9cSlwmQRGSPXHREJ/bvkM1tr37B+q11+ZaF71TVUnmBYEqVN8KvZa8B7n5WEmKTRMnOhePuhHnvweQno45GpMFKT03h3tN6smLDNu5988uow4mLqu6oX+Puc9z9DKARwb0qxwE5yQpOEqjvubDPAHjrFthQf5reInVNr/wWXHh4R57+9Fs+m7c66nBqbJd9KmZ2BUGrZe/w9YKZXZ7owCTBUlLglAegaCu88YuooxFp0H563H7s1aoRN748ja1FJVGHUyPV6ai/BOjn7r90918STC9fJ++ml3LadIZBN8LM0TDzP1FHI9JgNc5I4zcjejJ3xSYefHdO1OHUSHWSigGxPUhFVDwdvdRFh10F7XvC69fBlrVRRyPSYB3RJYfT+ubzl3Ff8+XS9VGHs8eqGv1Vdrf9k8AnZnaLmd0CfAQ8kYzgJAlS02HoKNi0HMbeFnU0Ig3aLSfvT/NG6dzw0ueUlNbNkZlVtVQ+A3D33xE8NXEzsAW41N3/LwmxSbLk9oFDr4BJT8C896OORqTBatkkg9uH9mDqgrU88dH8qMPZI1UllR2XuNx9vLv/wd1/7+7jkxCXJNugX0LLjvCfq6FoS9TRiDRYp/TqwFFdc/i/t75iweq6NzN6VUklx8x+VtkraRFKcmQ0hlP+BKvnwrh7oo5GpMEyM+4a0ROAm1/5gl09nqS2qSqppAJNgWaVvKS+6TQwmMn4o1HB/GAiEom8Fo24/oSuvDdrBa9OqVtPSa/0IV1mNsnd+yY5noTb44d0NRRb1sCf+0F2B/jJO5BanacjiEi8lZQ6p//1I+av3MR/fzaQ1k2je/zG7jykq1p9KtKAaCZjkVohNcW497RebNxWzF2vz4w6nGqr6t/QY5IWhdQu3YdBtyHBTMbdhkDrfaOOSKRB2q9dMy4b1JkH3p7NCT3ac3BByz3el5nRqklGHKOr5HPqWidQTenyVzWtXwIP9gueEHnuf4LZjUUk6bYVlzD4T+/z9YpNNdpPm6aZTLjl2D3aNt7PqJeGKLtDMJPxa9cGMxn3/XHUEYk0SJlpqTz9k0MYO2MpNWkCZKWnxi2mqiipSOX6ngufvwhjboEux0Oz9lFHJNIgtW+exTmHFkQdRrVUZ+4vaahSUmDoA1CsmYxFpHrUUpGqtd43mMn47V/Bvy+DrOZ7vq99j4L9TohfbCJS6yipyK4ddhV8+wl8+fqe76O0CD79C/T8AZx0LzRuFb/4RKTWUFKRXUtNh7NfqNk+irfD+7+H9/8P5v0PhtwP3QbHJz4RqTXUpyLJkZYBR90EF70DTXLgubPg5Ythc91/fKqIfEdJRZKrw4Fw0btw5PXByLKHDoWv3ow6KhGJEyUVSb60DDj65qDV0rg1PDsS/n1pMO+YiNRpSioSndzecPE4OPIXMO2FoNUya0zUUYlIDSipSLTSMuDoW+Cit4PJLJ/5QTB0ecvaqCMTkT2gpCK1Q26foNVyxHUw7Xl46BCY9VbUUYnIblJSkdojLROOuRV+8l/IagHPnAGvXK5Wi0gdoqQitU9eX7jkf3DEz2Hqc0Ffy+yxUUclItWgmx+ldkrLhGNuC57n8spl8PTp0LIArAb/B3UbAsf+KpjTTEQSIuFJxcxSgQnAIncfYmaPAoUET5acBZzn7hvNbG/gCaAFkArc6O5vhPu4CbgQKAGudvcxYfmJwJ/C+o+4+z2JPh5Jsry+cMl78NEoWPHlnu9nyxr46AHYshpOeQBSkjMNuEhDk4yWyjXATCA7fP9Td18PYGZ/AK4E7gFuAV5w97+YWXfgDaAgXB4J9ABygf+a2X7hvh4EjgMWAuPNbLS7z0jCMUkypWXCkdfVbB/uMO638L97oaQYhj+kxCKSAAlNKmaWD5wM3A38DCAmoRjQCHY8d8b5LvE0BxaHy8OA59x9GzDPzOYA/cJ1c9x9bri/58K6SiqyMzM46peQkgbv3g2lxTDib5CqK8Ai8ZTo36j7geuBZrGFZvYPYDBBAvh5WHwH8JaZXQU0Acqee5kHfBKz+cKwDGBBufL+cYxd6qOB1wcTZP73jmDm5NMeDd6LSFwkrMfSzIYAy919Yvl17n4+waWsmcCZYfFZwOPunk+QcJ40sxSCvpeddlFFeUWxXGxmE8xswooVK3b/YKR+OfyncPzdMONV+Nd5wQzKIhIXiRwGMwAYambzgeeAo83sqbKV7l4CPA+cFhZdCLwQrvsYyALaELRA9orZbz7BpbHKynfi7g+7e6G7F+bk5NT8yKTuO+xKOOk++PI1eP5HULQ16ohE6oWEJRV3v8nd8929gKCj/R3gHDPrDDv6VE4Byob0fAscE67bnyCprABGAyPNLNPMOgJdgM+A8UAXM+toZhnhZ4xO1PFIPdT/YhjyR5g9Bp77IRRtiToikTov2b2UBjxhZtnh8lTgsnDdz4G/m9lPCS5jnefuDkw3sxcI+l+KgSvCVg5mdiUwhmBI8WPuPj2pRyN1X+EFkJIOo6+CZ86Es56DjMZRRyVSZ1nwd7vhKCws9AkTJkQdhtQ2U58LbrLc+zD44fOQ2TTqiERqDTOb6O6F1amrW4tFAA4cCaf+Hb79GJ46DbaujzoikTpJSUWkTM/T4fRHYdEEeOpUTWQpsgeUVERi9RgBZzwBi6fAk8Nh8+qoIxKpU5RURMrbfwic+RQsmw7/HAqbVkUdkUidoaQiUpGuJ8JZz8LK2fDEKbBRN82KVIcmPhKpTOdjg5Fgz4yEx08OJqFMy4o6qoareT40ahF1FLILSioiVek0CH70Ijz9A3jkmKijadhy9ofLPtLzcGo5JRWRXSk4HC77AJZ+EXUkDdeSKfD+72Huu9BZyb02U1IRqY5WnYKXRGO/E2DSk/Dp35RUajm1I0Wk9kvLDKbUmT0GVn0ddTRSBSUVEakbCs8P5mn77O9RRyJVUFIRkbqhWXvoMRymPA3bNkQdjVRCSUVE6o7+l8K29cEEoFIrKamISN2RXwh5BwUd9qWlUUcjFVBSEZG6pf+lsGo2zH0n6kikAkoqIlK3dB8OTdsFrRWpdZRURKRuScuAg86H2W9peHEtpKQiInWPhhfXWkoqIlL3NGsfPPtm8lMaXlzLKKmISN3U/1LYvgGmPBt1JBJDSUVE6qb8gyCvED7T8OLaRElFROqu/pfAqjkaXlyLKKmISN2l4cW1jpKKiNRdaRnh7MUaXlxbKKmISN12UNnw4oejjkRQUhGRuq5ZOzjgVJis2YtrAyUVEan7+l+i4cW1hJKKiNR9eRpeXFsoqYhI/dD/0mB48dcaXhwlJRURqR+6DwuHF/816kgaNCUVEakf0jKg8EKYMxZWzok6mgZLSUVE6o+y2YvHa/biqCipiEj90bTtd8OLt66POpoGSUlFROqXsuHFUzW8OApKKiJSv+QdBPkHB/OBaXhx0impiEj90/9SWP01fP121JE0OAlPKmaWamaTzey18P2jZjbVzKaZ2Ytm1jSm7g/MbIaZTTezZ2LKzzWz2eHr3Jjyg8zsczObY2YPmJkl+nhEpA7Yfyg0ba/ZiyOQjJbKNcDMmPc/dfcD3b0X8C1wJYCZdQFuAga4ew/g2rC8FXA70B/oB9xuZi3Dff0FuBjoEr5OTPzhiEitVzZ7sYYXJ11Ck4qZ5QMnA4+Ulbn7+nCdAY0AD1ddBDzo7mvCesvD8hOAse6+Olw3FjjRzDoA2e7+sbs78E9geCKPR0TqkELNXhyFRLdU7geuB77XW2Zm/wCWAt2AUWHxfsB+ZvahmX1iZmWtjjxgQczmC8OyvHC5fPlOzOxiM5tgZhNWrFhRw0MSkTqhaVs44DSYouHFyZSwpGJmQ4Dl7j6x/Dp3Px/IJbgsdmZYnEZwCWsQcBbwiJm1ACrqJ/EqyncudH/Y3QvdvTAnJ2d3D0VE6qr+F8P2jTDlmV3XlbhIZEtlADDUzOYDzwFHm9lTZSvdvQR4HjgtLFoIvOruRe4+D/iKIMksBPaK2W8+sDgsz6+gXEQkkHcQ5PcLLoFpeHFSJCypuPtN7p7v7gXASOAd4Bwz6ww7+lROAb4MN3kFOCpc14bgcthcYAxwvJm1DDvojwfGuPsSYIOZHRLu68fAq4k6HhGpo/pfouHFSZSW5M8z4Akzyw6XpwKXhevKkscMoAT4hbuvAjCzXwPjw3p3uvvqcPky4HGCDv83w5eIyHe6D4MxN8PLFwWzGDdUjVrBBYn/E2nBwKmGo7Cw0CdMmBB1GCKSTDNfg89fiDqKaGU1h6Gjdl2vAmY20d0Lq1M32S0VEZHk239I8JKE0zQtIiISN0oqIiISN0oqIiISN0oqIiISN0oqIiISN0oqIiISN0oqIiISN0oqIiISNw3ujnozWwF8s4ebtwFWxjGceFN8NaP4akbx1Uxtjm8fd6/WFO8NLqnUhJlNqO5UBVFQfDWj+GpG8dVMbY+vunT5S0RE4kZJRURE4kZJZffU9oddK76aUXw1o/hqprbHVy3qUxERkbhRS0VEROJGSUVEROJGSaUCZnaimX1lZnPM7MYK1mea2fPh+k/NrCCJse1lZu+a2Uwzm25m11RQZ5CZrTOzKeHrtmTFF37+fDP7PPzsnR6zaYEHwvM3zcz6JjG2rjHnZYqZrTeza8vVSer5M7PHzGy5mX0RU9bKzMaa2ezwa8tKtj03rDPbzM5NYnz3mdmX4ffv32bWopJtq/xZSGB8d5jZopjv4eBKtq3ydz2B8T0fE9t8M5tSybYJP39x5+56xbyAVOBroBOQAUwFupercznw13B5JPB8EuPrAPQNl5sBsyqIbxDwWoTncD7Qpor1g4E3AQMOAT6N8Hu9lODGrsjOH3Ak0Bf4Iqbsd8CN4fKNwL0VbNcKmBt+bRkut0xSfMcDaeHyvRXFV52fhQTGdwdwXTW+/1X+ricqvnLrfw/cFtX5i/dLLZWd9QPmuPtcd98OPAcMK1dnGPBEuPwicIyZWTKCc/cl7j4pXN4AzATykvHZcTQM+KcHPgFamFmHCOI4Bvja3fd0hoW4cPf3gNXlimN/xp4Ahlew6QnAWHdf7e5rgLHAicmIz93fcvfi8O0nQH68P7e6Kjl/1VGd3/Uaqyq+8O/GD4Bn4/25UVFS2VkesCDm/UJ2/qO9o074i7UOaJ2U6GKEl936AJ9WsPpQM5tqZm+aWY+kBgYOvGVmE83s4grWV+ccJ8NIKv9ljvL8AbRz9yUQ/CMBtK2gTm05jxcQtDwrsqufhUS6Mrw891gllw9rw/k7Aljm7rMrWR/l+dsjSio7q6jFUX7cdXXqJJSZNQVeAq519/XlVk8iuKRzIDAKeCWZsQED3L0vcBJwhZkdWW59bTh/GcBQ4F8VrI76/FVXbTiPNwPFwNOVVNnVz0Ki/AXYF+gNLCG4xFRe5OcPOIuqWylRnb89pqSys4XAXjHv84HFldUxszSgOXvW/N4jZpZOkFCedveXy6939/XuvjFcfgNIN7M2yYrP3ReHX5cD/ya4zBCrOuc40U4CJrn7svIroj5/oWVllwTDr8srqBPpeQwHBgwBzvawA6C8avwsJIS7L3P3EncvBf5eyedGff7SgFOB5yurE9X5qwkllZ2NB7qYWcfwv9mRwOhydUYDZSNtTgfeqeyXKt7Ca7CPAjPd/Q+V1Glf1sdjZv0Ivs+rkhRfEzNrVrZM0KH7Rblqo4Efh6PADgHWlV3qSaJK/0OM8vzFiP0ZOxd4tYI6Y4DjzaxleHnn+LAs4czsROAGYKi7b66kTnV+FhIVX2wf3YhKPrc6v+uJdCzwpbsvrGhllOevRqIeKVAbXwSjk2YRjAy5OSy7k+AXCCCL4LLJHOAzoFMSYzucoIk+DZgSvgYDlwKXhnWuBKYTjGb5BDgsifF1Cj93ahhD2fmLjc+AB8Pz+zlQmOTvb2OCJNE8piyy80eQ3JYARQT/PV9I0Ef3NjA7/NoqrFsIPBKz7QXhz+Ec4PwkxjeHoD+i7GewbDRkLvBGVT8LSYrvyfBnaxpBouhQPr7w/U6/68mILyx/vOxnLqZu0s9fvF+apkVEROJGl79ERCRulFRERCRulFRERCRulFRERCRulFRERCRulFREKmBm7czsGTObG06R8bGZjYgolkFmdljM+0vN7MdRxCKyK2lRByBS24Q3Pr4CPOHuPwzL9iGY1iVRn5nm303QWN4gYCPwEYC7/zVRcYjUlO5TESnHzI4hmIp8YAXrUoF7CP7QZwIPuvvfzGwQwXTrK4EDgInAj9zdzewg4A9A03D9ee6+xMzGESSKAQQ36M0CbiGYhn0VcDbQiOAGzBJgBXAVwezKG939/8ysN/BXghs6vwYucPc14b4/BY4CWhDccPd+ODnmP8LPSAFO88onMxTZbbr8JbKzHgSTSlbkQoJpZQ4GDgYuMrOO4bo+wLVAd4K7oQeE87SNAk5394OAx4C7Y/bXwt0HuvvvgQ+AQ9y9D8E07Ne7+3yCpPFHd+/t7u+Xi+efwA3u3ovgDvLbY9aluXu/MKay8kuBP7l7b4K78yucIkRkT+nyl8gumNmDBNPjbAe+AXqZ2enh6uZAl3DdZx7O4xQ+ya8AWEvQchkbTieWSjBlR5nYyQTzgefDeasygHm7iKs5QVL6X1j0BN+fdblsstGJYSwAHwM3m1k+8LJaKRJvaqmI7Gw6wZP6AHD3KwguOeUQzFt2Vdhq6O3uHd39rbDqtph9lBD802bA9Jj6Pd39+Jh6m2KWRwF/dveewCUEc8zVRFk8ZbHg7s8Q9A1tAcaY2dE1/AyR71FSEdnZO0CWmV0WU9Y4/DoGuCy8rIWZ7RfOIFuZr4AcMzs0rJ9exUO/mgOLwuXY581vIHh09Pe4+zpgjZkdERadA/yvfL1YZtYJmOvuDxD04/Sqqr7I7lJSESnHg9Erw4GBZjbPzD4juLR0A/AIMAOYZGZfAH+jisvIHjym9nTgXjObSjCj72GVVL8D+JeZvU/QoV/mP8AIM5sSk0DKnAvcZ2bTCB5IdecuDu9M4Ivw8lw3gj4ZkbjR6C8REYkbtVRERCRulFRERCRulFRERCRulFRERCRulFRERCRulFRERCRulFRERCRu/h/xOMorRrNdDgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.ylabel('Total Costs')\n",
    "plt.xlabel('Generations')\n",
    "plt.plot(mean_total_costs,label='Mean')\n",
    "plt.plot(min_total_costs,label='Min')\n",
    "plt.title('Neighborhood ' + str(num_houses) + ' houses')\n",
    "plt.legend()\n",
    "plt.savefig('plots/' + str(num_houses) + '_household_onlyPV.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we already stated above we can also optimize the cost for each house seperately. To check the results of the genetic algorithm above, we do the seperate optimization, which gives us the true result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEYCAYAAADmugmLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4VVXWwOHfIgFCCzUkoffeBBTUkUFRRGXEgjOIIzgWPuvo6NjGsYzojL33wthQ7IIIKKIOFlBCJ7SEEgglCQRCCunr++Oc6DUmJKSdW9b7PHnI3ffcc9eN27vP3nudvUVVMcYYY/xRPa8DMMYYY8pjjZQxxhi/ZY2UMcYYv2WNlDHGGL9ljZQxxhi/ZY2UMcYYv2WNlDHGGL9ljZRLRLaLyKmlyi4Rke+8iqksIjJZRJJEJFtEPhGRVl7HFIoCob6ISKyIzBGR3SKiItLF65hCWYDUmbNE5DsROSgie0XkZRFp5mVM1kgFEBHpD7wIXAxEAznAc54GZfxZMbAAON/rQEzAaA7cB7QD+gIdgIe9DMgaqaMgIn1F5Bv3KiNeRM72ee4bEbnc5/HPV0jieFxEUkUkQ0TWiMgA97mGIvKIiOwQkRQReUFEGpUTwkXAp6q6WFWzgDuB87y+0jFl87q+qGqKqj4HLKvlj2pqiB/UmbdVdYGq5qjqAeBl4MTa/dRHZo1UJYlIfeBT4AugLXAdMFNEelfi5WOBUUAvoAXwJ2C/+9yDbvkQoAfQHrirnPP0B1aXPFDVLUC++3rjR/ykvpgA4qd1ZhQQX8lja4U1Ur/2iXsFc1BEDvLrobSRQFPgAVXNV9WvgLnAhZU4bwHQDOgDiKpuUNU9IiLAFcDfVDVdVTOBfwOTyjlPUyCjVFmGe25T9/y9vhj/EzB1RkROA6bi8UWQNVK/do6qtij5Aa72ea4dsFNVi33KknCuSo7IrWzPAM8CKSLykohEAlFAY2C5T6Vd4JaXJQuILFUWCWRW4rOZmufv9cX4n4CoMyIyEngbmKiqm4/i89U4a6QqbzfQUUR8/2adgF3u79k4laFEjO+LVfUpVR2GM2TXC7gZ2AccBvr7VNzmqtq0nBjigcElD0SkG9AQ8LQSmTL5Q30xgcUv6oyIHAPMAS5V1UXV/VDVZY1U5f2IU0luEZH6IjIa+AMwy31+FU4SQ2MR6QFcVvJCETlWREa4Y87ZQC5Q5F4xvQw8LiJt3WPbi8jp5cQwE/iDiJwkIk2Ae4GP3C688S/+UF8QkQicCxmAhu5j4588rzNussUC4DpV/bRWPuVRskaqklQ1HzgbOAPn6uQ5YIqqbnQPeRwniSEFeB2nQSkRiVNRDuB03/cDj7jP3QokAktF5BDwJVDmRKmqxgNXuudOxRmDvrqsY423/KG+uA7jDBMDbHQfGz/kJ3XmJpyhwFdFJMv98TRxQmzTQ2OMMf7KelLGGGP8ljVSxhhj/JY1UsYYY/yWNVLGGGP8VrjXAVRVmzZttEuXLl6HEZSWL1++T1WD6gZRqy+1K9jqjNWX2nU09SVgG6kuXboQFxfndRhBSUSSvI6hpll9qV3BVmesvtSuo6kvNtxnjDHGb1kjFeLyCotIy8zzOow6VVysZBwu8DoMY0JCZm4BB3Pyq/x6a6RC3KyfdjLqoa9J2p/tdSh1Ztqby5n00lKvwzAm6Kkqt364hgnPfk9uQVGVzmGNVAjLyS/k6a8SGdShOZ1aNa74BUEiOrIhKYdyvQ7DmKD32g/bmbd2L5OP60RE/bAqncMaqRD22g/b2ZeVx82n98bZdiY0xERGkJ6dX+UrO2NMxVbuOMC/523g1L7RTBvVrcrnsUYqRGUcLuCFb7ZwSp+2DO/Syutw6lRMc2chcOtNGVM7DmTnc83MFURHRvDoBYOrdRFsjVSIennxVg7lFnLT2NDbeT62eSMA9mZYI2VMTSsuVv723ir2ZeXz3EVDad64frXOZ41UCNqXlceM77cxflAs/ds19zqcOhfT3Nleaa/1pIypcc//bwvfbErjzvF9GdShRbXPZ41UCHru6y3kFRZz42mh14sCiHF7UnusJ2VMjfphyz4e/WITfxjcjj+P7Fwj57RGKsTsOniYt5YmMXFoB7pFheau400bhtOsYbgN9xlTg1IO5fLXd1bStU0THjhvYI0lYwXsskimap5elADAX0/t6XEk3opuHmGNlDE1pLComOveXkl2XhFvXzGSJg1rrmmpsCclIjNEJFVE1vmUTReRNSKySkS+EJF2pV5zrIgUichEn7KpIpLg/kz1KR8mImtFJFFEnpJQyoWuY4mpWbwXt5OLRnaifYtGXofjqdjmEeyxOSljasTDX2zip+3p/Oe8gfSKblaj567McN9rwLjSManqIFUdAswF7ip5QkTCgAeBz33KWgF3AyOA44C7RaSl+/TzwDSgp/tT+r1MDXn0i000qh/GtSf38DoUz8VERpBiPSljqu2L+L28+L+tXDSiE+cc077Gz19hI6Wqi4H0UmWHfB42AdTn8XXAh0CqT9npwEJVTVfVA8BCYJyIxAKRqrpEVRV4AzinSp/EHNGqnQeZv24vV4zqRuumDb0Ox3MxzSNIzcylsKjY61ACzprkgxQVa8UHmqCXtD+bm95fzcD2zblzfL9aeY8qJ06IyP0ishO4CLcnJSLtgXOBF0od3h7Y6fM42S1r7/5eury895wmInEiEpeWllbV0EOOqvLg/I20btKAy0+q+p3fwSSmeQTFCmlZobW4bnWlZuYy8fklPLZwk9ehGI/lFhRx5VsrqCfCcxcNrfKyRxWpciOlqneoakdgJnCtW/wEcKuqll5vpqx5Jj1CeXnv+ZKqDlfV4VFRQbO/Wq37NmEfS7bu59pTetC0Bic0A1msu+qEJU8cnTd+SKKguJgLhnX0OhTjsbtmr2PDnkM8/qfBdKzFtT9r4hvrbeAznDmn4cAsN/ehDXCmiBTi9JBG+7ymA/CNW96hVPnuGojJuIqLlQcXbKRDy0ZMHtHJ63D8RnSkNVJHKye/kLd+TGJsv2i6tGnidTjGQ+8u28F7cclcd0oPTukTXavvVaWelIj45i+fDWwEUNWuqtpFVbsAHwBXq+onOEkUY0WkpZswMRb4XFX3AJkiMtLN6psCzK76xzGlfbZ2D/G7D3Hjab1oGF473fHSyskIbSUiC93szoUliTMicpGbKbpGRH4QkcE+rxknIpvczM/bfMq7isiP7rneFZEGRxtjrN3Qe9Q+WJ7MwZwCrrAh45C2blcGd86O53c92nDDqbW/IEBlUtDfAZYAvUUkWUQuAx4QkXUisganwbn+SOdQ1XRgOrDM/bnXLQO4CngFSAS2APOr+mHMrxUUFfPoF5voHd2MCUNqPuvmCF7jt1matwGLVLUnsMh9DLAN+L2qDsKpIy/Bz1mizwJnAP2AC0WkZGb2QeBx91wHgMuONsCWjevTILyeLTJbSUXFyqvfbWNIxxYM69yy4heYoJSRU8BVM5fTukkDnpw0hLB6tX/HUIXDfap6YRnFr1bidZeUejwDmFHGcXHAgIrOZ47eOz/tYPv+HGZcMrxOKlMJVV0sIl1KFU/glyHf13GGe29V1R98jlnKL8O/xwGJqroVQERmARNEZANwCjDZ51z34NzKUGki4twrZT2pSlm4fi9J+3O4dVyfkNrWxfyiuFi58b1V7M3IZda04+ssS9iWRQpSmbkFPPllAiO6tuLk3m29Dgcg2h3exf23rKAu45eedHkZoa2Bg6paWKr8NyrKBo2OtFUnKuvlb7fRsVUjTu8f43UoxiPPfZPIoo2p/POsfnXam7ZGKki9vHgr+7Pzuf3MvgFx5SsiJ+M0UreWFJVx2FFlhFaUDRrbPMJWQq+E5UkHWJ50gMtO7FqnPXLjP75NSOPRhZuZMKQdU46vmYVjK8saqSCUeiiXl7/dxlkDYxnSsfpL5deQFPfmbdx/f77ZW0QG4cxLTlDV/W5xMuCb51yS+bkPaCEi4aXKj1qMu36fcx+5Kc8r324lMiKcC4Z7k3YuImEislJE5rqPXxORbe6ybKtEZIhbLu7SaoluIs5Qn3OUuSybqdjug4e5ftYqerZtyn9qcOHYyrJGKgg9sSiBgqJibj69t9eh+JoDlHw5TMXN4hSRTsBHwMWqutnn+GVATzeTrwEwCZjjrkzyNTCx9LmOVkxkBPlFxaRn51fl5SEhaX82n8fv5c8jO9fooqFH6XpgQ6mym1V1iPuzyi07g1+WV5uGO09ZwbJs5gjyCou4auYK8guLef7Pw2jcoO7rgDVSQWZLWhbvLtvJ5BGdPLuXpbyMUOA0EUkATnMfg7NaSWvgOfeqOA7AnXO6Fuf2hQ3Ae6oa777mVuBGEUl0X1thIk9Zfr6h14b8yvXS4q2E16vHJSd08eT9RaQDcBZOT7siE4A31LEUp8cdSznLstVa0EHkvrkbWL3zIA9PHER3j7b2seUHgsxDCzYSEV6Pv47xbiuOcjJCAcaUcezlwOXlnGceMK+M8q04V8TVEuOzjXwo7lBckdTMXN5fnsz5wzrQ1r352QNPALcApZfWvl9E7sK9nUFV8zjy8mtllf+KiEzD6YHRqZPd+P7h8mTeXJrE/43qxhkDYz2Lw3pSQSRuezqfx6cwbVR32tgishWKcb94LQ29bDO+205hUTH/N8qbm3dFZDyQqqrLSz11O9AHOBZoRQ0l29iya7+I353BPz5ey8hurTyfNrBGKkioKvfP20DbZg25YlRXr8MJCFHNGhJWT+yG3jIcyi1g5tIkzhwY6+USSCcCZ4vIdmAWcIqIvKWqe9whvTzgv/zSqy4v2aa8clOGjJwCrnprBS0bN+DpC4cSHuZtM2GNVJCYt3YvK3cc5O9je3syuRmIwuoJbZs1tJ5UGd5ckkRmXiFX/r67ZzGo6u2q2sFdZm0S8JWq/tknS1RwtvYpWX5rDjDFzfIbCWS49+SVuSxbXX+eQFBcrPztvVXsyTjMsxcNJaqZ9yMy9m0WBPIKi3hwwUb6xDTj/GEdKn6B+Vl0ZIT1pErJLSjiv99vY1SvKAa098u5upkiEoUzjLcKuNItnwecibPEWg7wF3CWZRORkmXZ4NfLshkfT3+VyFcbU5k+ob/fLH9ljVQQeHNJEjvSc3jj0uPsZsujFNs8goTULK/D8Cvvx+1kX1Y+V4/2rhdVmqp+g7OUFqp6SjnHKHBNOc+VuSyb+cXXm1J5YtFmzjumPX8eWbc37B6JDfcFuIM5+Tz9VSKjekUxqldoT/ZWRckNvcZRWFTMi4u3ckynFozo2srrcEwdSdqfzfXvrKRPTCT3n1v3N+weiTVSAe6ZrxLJzC3gH2f28TqUgBQTGUFWXiGZuQVeh+IX5q7ZQ/KBw1w9uodffVGZ2nM439lhV0R48c/DaNSgbrb0qSxrpALYjv05vL5kOxcM60ifmEivwwlIMe4NvTYv5WzH8dRXCfSObsaYPn6xKLGpZarKPz5ey8a9h3hi0hA6ta69HXaryhqpAPaf+RsIr1ePG8fW/sZjwco2P/zF3DW72ZqWzV/H9KSezW2GhDeWJPHxyl3cMKaXv+yW8BvWSAWon7alM3/dXq4a3f3nrdDN0YuxbeQBpxf19FeJ9IpuyhkDbDuOULBsezrT565nTJ+2XHdKD6/DKZc1UgGouFiZPnc9sc0jbCvvamob6dwHEuqN1Ly1e0hMzeK6U6wXFQpSDuVy9cwVdGjZiMf+NMSv/5tXZvv4GSKSKiLrfMqmu0vhrxKRL0SknVt+kVu+RkR+EJHBPq8ZJyKb3GX0b/Mp7yoiP7pL6L/rrnhtjuDjlbtYuyuDW8b19rtJzkATUT+M1k0asCeE56SKi5WnFiXQo21TzvRwjTZTN/ILi7l65gqy8wp58eLhNG9U3+uQjqgyPanX+O2KwQ+r6iBVHQLMxVnJGmAb8HtVHQRMB14CZz8Y4FmcpfT7AReKSD/3NQ8Cj6tqT+AAzsZ3phw5+YU89PlGBndozoTBZW5Ia45SdGQEKSHck5q3bg8JqVlcd0oPu88uBEyfu57lSQd4aOIgeseUXrfX/1TYSKnqYiC9VNkhn4dNcBdrVNUf3KXwAZbirJEFztpaiaq6VVXzcdbhmuAua3IK8IF73Os4y5yYcry0eCsph/K4c3w/v+6iB5LY5hEhmzhR0ovqHtWE8YPaeR2OqWXvx+3kzaVJTBvVLWD+e1d5TkpE7heRncBF/NKT8nUZMN/9vbyl8lsDB929g3zLy3vPaSISJyJxaWlpVQ09YO3NyOXF/23lrIGxDO9iN1rWlJgQ3kZ+QfxeNqdk8dcxPa0XFeTWJB/kjk/WcXy31tziXxuiHlGVGylVvUNVOwIzcTan+5mInIzTSNXIEvo+7xnSS+k/uGAjRarcdobduFuT2rVoRHp2Pjn5hRUfHESKipXHFm6mm/Wigt6+rDz+783lRDVtyDOTj/F8ZfOjURORvg2cX/JARAbh7KI5QVX3u8XlLZW/D2f3zPBS5aaU5UkH+HjlLq44qSsdW/nfDXeBrJu7FcWW1GyPI6lbH61IJjE1i7+P7W29qCBWUFTMNTNXkJ6dz4sXD6N1gO01V6VGSkR8t309G9jolncCPgIuVtXNPscsA3q6mXwNcJbdn+MuCPk1MNE9biowuyoxBbPiYuVfn8YTHdmQq0f77/0MgapntDN5nJCa6XEkdSevsIgnvkxgYPvmdl9UkLv/sw38uC2dB84f6K+r2h9Rhaugi8g7wGigjYgkA3cDZ4pIb6AYSOKXpfLvwplnes5d96vQHZ4rFJFrcfZwCQNmqGq8+5pbgVkich+wEni1pj5csPhgeTJrkjN44k9DaNLQFq6vaZ1bN6Z+mLA5JXRWQ5+5dAe7Dh7mwfMH2Rp9QezD5cm89sN2Lj2xK+ceE5jb+FT4jaeqF5ZRXGZDoqqXA5eX89w8nP1eSpdv5ZedNU0ph3ILeOjzjQzt1IIJQ2zeoDbUD6tHtzZNSQyRnlRWXiHPfJ3IiT1a87uebbwOx9SSNckHuf3jtRzfrXVAL0AdOLNn5cjOK+T9uJ0kpATnF8wzXyWyPzufe87ub1e8tahndNOQ6Um98u1W0rPzueX0wP3iMkeWlhm4iRKlBW7krsJi5eYP1rBwQ4rXodS4rWlZ/Pf7bVwwrAODOrTwOpyg1rNtM3YeyOFwfpHXodSq/Vl5vLx4K2cMiGFwR6tTwSi/sJir3lrOgZx8XpoSeIkSpQV8I9W8UX06tmpE/O5DFR8cQFSVez5dT8PwMP4eQPc0BKpe0U1RhS1pwd2bevbrLRwuKOKmsVangtXdc+KJSzrAwxMH079d4CVKlBbwjRRAv9hINgRZI7Vg3V4Wb07jb6f1om0zW+W8tvWMbgrA5iAdNgbYti+bN5c6+4/1aNvU63BMLXhraRLv/LSDq0d35w+Dg2MOOygaqf7tmrNtfzbZecFxM2Z2XiH3zl1Pn5hmTD2+s9fhhITOrZtQP0xISA3entT9n22gYXgYN51u+48Fo5+2pXPPnHhO7h0VVD3loGik+sVGogob9wZHb+qprxLYk5HLfecMCOgJz0BSkuEXrAk43yXs48sNKVxzcg/rmQehnek5XPnWcjq1bsyTFx4TVDdnB8U3YP/2ztbpwTAvlZCSyavfbmPisA62Pl8d6xHdNCh7UoVFxdw7N56OrRrxlxO7eB2OqWHZeYVc8UYchUXFvDJlOJER/r31xtEKikYqJjKClo3rE78rsBspVeXO2eto3CDM1ufzQK+2zdiRHnwZfu8s28nmlCzuOLMvEfVt/7FgUlys3PjeKjanZPLM5KF0iwq+ucagaKREhP7tmrN+T2A3UnNW72bp1nRuHteHNgGeNhqIegZhhl9GTgGPfbGJkd1acXp/W/4o2Dzx5WY+j0/hjrP6MapXcC66HRSNFED/dpFs2ptJQVGx16FUSUZOAdPnbmBQh+ZMPq6T1+GEpF5uhl8wreH35KIEDh4u4M7x/exm8CAzd81unvoqkT8O78ClQTyMGzSNVL92keQXFQfsVfC/523gQE4+/z53YFBNegaSkgy/YFl5YuPeQ7yxZDuTju0YFPfLmF+sST7ITe+tZljnlkw/Z0BQX4AETSPVv52bPBGA81I/bNnHu3E7ufykrgG5SnFpIjJDRFJFZJ1PWSsRWSgiCe6/Ld1yEZGnRCRRRNaIyFCf10x1j08Qkak+5cNEZK37mqekhv4PrR9Wj65tmpAQBI1UcbHyj4/WEtmovi1/FGT2ZuRyxRtxtGnakBcvHkbD8OCeZwyaRqprm6ZE1K8XcPNSuQVF/OOjtXRq1ZgbxgTN/SuvAeNKld0GLFLVnsAi9zHAGUBP92ca8Dw4jRrOivsjcBYgvrukYXOPmebzutLvVWU9o5sFxXDfO8t2sGLHQe44sy8tmzTwOhxTQw7nF3HFG3Fk5Rby6iXDQ2LuOmgaqbB6Qp+YSOJ3Z3gdylF5clEC2/fn8J/zBtKoQXBcEanqYiC9VPEE4HX399eBc3zK31DHUpxNMGOB04GFqpquqgeAhcA497lIVV3i7kf2hs+5qq1n26YBn+GXmpnLA/M3cny31pw3tL3X4VSbiISJyEoRmes+7ioiP7o97HfdPeoQkYbu40T3+S4+57jdLd8kIqd780mqp7hY+fv7q1m3O4OnLjyGPjGRXodUJ4KmkQJnXmr97kM4313+b/3uQ7y0eCsTh3XgxB5Bv2VCtKruAXD/beuWtwd2+hyX7JYdqTy5jPIa0Su6WcBn+E2fu4G8gmLuPzdo5iquBzb4PH4QeNztlR8ALnPLLwMOqGoP4HH3OESkH85Gq/1xet3PiUjAXRE+uSiBz9bu4fYz+jCmb7TX4dSZoGqk+reL5FBuIckHDnsdSoUKi4q57aM1tGxcnzvO7Ot1OF4q61tUq1D+2xOLTBOROBGJS0tLq1QwPdsGdobfN5tS+XT1bq45uUdQ3DMjIh2As4BX3McCnAJ84B5Sulde0lv/ABjjHj8BmKWqeaq6DUgkwPaw+2TlLp5clMAFwzpwxUndvA6nTlXYSJUzCT7dneReJSJfiEg7t9zTSfB+sU73NxDmpZ7/ZgtrkjO45+z+oTJnkOIO1eH+m+qWJwMdfY7rAOyuoLxDGeW/oaovuTtDD4+Kqtw9JF3aNCG8ngRk8kROfiF3zl5H96gmXDk6aL7IngBuwdkFHJydvw+qaslCnb496Z973+7zGe7x5fXKA8JP29K55YM1HN+tNfefOzBYeseVVpme1Gv8dmL6YVUdpKpDgLk428aDx5PgfWIiqSf+vzxS/O4MnlyUwPhBsYwfFBwrFVfCHKDk4mQqMNunfIp7gTMSyHCHAz8HxopIS7eujAU+d5/LFJGR7gXNFJ9zVVtJhl8gpqH/e94Gkg8c5t/nDgyKjC8RGQ+kqupy3+IyDtUKnqtU77sqPe/atn1fNv/3ZhwdWjXihT8Po0F4UA1+VUqFn7isSXBV9W0FmvDLf3BPJ8EbNQijW1RT1vtxI5VXWMSN766mZZMGTJ8wwOtwaoWIvAMsAXqLSLKIXAY8AJwmIgnAae5jgHnAVpwhmJeBqwFUNR2YDixzf+51ywCuwhn+SQS2APNrMv5eAZjh99XGFN5auoMrTurGiG6tvQ6nppwInC0i24FZOMN8T+B8r4S7x/j2pH/ufbvPN8f57iqvV/4rVel516aDOflc+toyAGZMPZbmjYNrTb7KCq/4kLKJyP04V7EZwMluca1OgovINJxeF506lb0qQ/92kSzbVjqxzH88vjCBTSmZ/PeSY4N2mE9VLyznqTFlHKvANeWcZwYwo4zyOKDWWvie0U2Zt24PuQVFAbHW3b6sPG75YA19Yppx09iguY0BVb0duB1AREYDf1fVi0TkfWAiTsNVulc+FecCaSLwlaqqiMwB3haRx4B2OCM2P9XlZzlaeYVFXPnWcpIPHGbmFSPo0qaJ1yF5psp9R1W9Q1U7AjOBa93iWpsEd9+zwiudfrGR7M7I5UB2/hHj98LypHReWryFScd25OQ+bSt+gfFEz7bN3K1f/L83parc9uFaDuUW8sSkIUExzFcJtwI3ikgizpzTq275q0Brt/xG3HvxVDUeeA9YDywArlFVv73HoLhYueWDNSzdms5DEwdxbIjvhlATA5xvA+e7v9faJHhllSz/4m/JE9l5hdz03mpimzfijrNCOpvP743s1oqwesIX8Xu9DqVCs5bt5MsNKdw6rk9Q3zejqt+o6nj3962qepyq9lDVC1Q1zy3PdR/3cJ/f6vP6+1W1u6r2VtUaHR6uaQ9/sYnZq3Zz8+m9OeeYgMnvqDVVaqREpKfPw7OBje7vnk+C9ytZHsmPbupVVe78ZB1J6Tk8+sfBNAuy/V6CTeumDTmhe2s+W7vHr++525qWxb2frud3PdrwlxO6eB2OqQFvLU3i+W+2MHlEJ64e3d3rcPxCZVLQy5wEF5F1IrIGp8G53j3c80nwVk0a0L5FI37adqA6p6lR7y9P5qOVu7h+TE9GBs+kdlAbPyiWpP05fpspmpNfyFVvrSCifj0euWAw9WxR4oD35foU7pq9jjF92nLv2f1DLtW8PBUmTpQzCf5qGWV+Mwl+xoAYXl+ynQPZ+Z4nJ2xOyeSu2es4oXtrrjulZ8UvMH5hbL8Y7vh4HZ+u2e13i/6qKrd/tJbNqZm8celxxDS37eAD3fKkA1z7zgoGtm/O05OPITws9FLNyxOUf4nzhnagoEiZu6Za01vVlpNfyDUzV9C0YThPTBpiW3AEkJZNGnBijzZ8tsb/hvxe/2E7s1ft5u9je3NST+9TpU31bNqbyaWvLSMmMoJXph5L4wZVTroOSkHZSPVrF0mfmGZ8uGKXp3HcPTuexLQsnvjTMbRtZle7gWb8oFiSDxxmTbL/zG/GbU/nvs82cGrfaK76vc1ZBLqd6Tlc/OqPRNSvx5uXjSCqWfCvan60grKRAjhvaHtW7TzIVo8WCn132Q7eX57MdSf34Hc9g37x2KA0tl8M9cOEz9bu8ToUwFnd/OqZK+jQshGP/tHmoQJdWmYeF7/6I3mFxbxx6Qg6tmrsdUh+KWgbqQlD2lNP4OOVdd+b+mlbOv/8ZB0n9WwpuHe1AAAgAElEQVTDX8fYPFSgat64Pif1jPKLIb+c/EKueGM5h3ILeOHiYTRvZBmigSzjcAFTZvxEyqE8ZlxyLL1jmnkdkt8K2kYqOjKCE3u04eOVuygurrsvmJ3pOVz51nI6tmrMM5OH2gRogDtrYCy7Dh5m5c6DnsVQWFTMtW+vZG3yQZ6aFDr7CAWrQ7lOA5WYmskLFw9jWOeWFb8ohAX1N+j5QzuQfOAwy7bXzTJJWXmFXP56HEXFyqtTj7Wr3SBwWv9oGoTV47M13gz5qSp3fLyOrzamMv2cAYztH+NJHKZmZOYWMHXGT8TvyuDZyUP5fS9LfKlIUDdSY/tH07hBGB/VQQJFUbFyw6yVJKZl8ezkoXQN4bW2gklkRH1G9Ypi3to9ddojL/HElwm8G7eT607pwUUjOtf5+5uak5VXyCX/Xcba5AyemTzULjgqKagbqcYNwjljQCzz1jqLhdYWVWX63PV8uSGVu//QzxIlgsz4QbHsycglLqlubxCf+WMSTy5KYOKwDtx4WvAsHBuKsvMK+ct/f2LVzoM8feExjBtgDVRlBXUjBU6WX2ZeIQvXp9TK+VWVB+Zv5LUftnP577oy5fgutfI+xjun9oumVZMG3PfZegqKiit+QQ2Y8d027vh4HaN7R/Gf80Jvo7tgciA7n4te+ZEVO5w5xTMGxnodUkAJ+kZqZLfWxDaP4J2fdtRKhtZjCzfz4uKtTDm+sy0cG6SaNgzn/nMGsCY5g2e/TqzV91JVnvwygXvnruf0/tG8ePEw6lvyTcDak3GYC15cwvo9h3h28lDOGmQN1NEK+tofVk+4/KRu/LBlP2//tKNGz/30ogSe/iqRScd25J4/2FpbweyMgbGce0x7nvkqkTXJtZPpp6rc/9kGHv9yM+cNbc+zk4eGytYbQWlLWhYTn1/C3oxcXv/LcTbEV0VB30gB/OWELpzUsw3T564nIaX6ewSpKk8tSuDRhc6Xyb/PHWg3VoaAe87uT5umDbnxvdU1PseZW1DErR+u4ZXvtjH1+M48MnGw3b4QwNYkH+SCF5aQW1DErGkjOb67LSxdVSHxf0G9esKjfxxMkwbhXPfOymp9wWTmFnDVWyt4bOFmzjumPQ9PtDv/Q0XzRvV5+IJBJKZm8fDnm2rsvImpmZzz7Pe8F5fMX0/pwT1n97c6FcA+XJ7MBS8soXGDMD646gS/W6A40IREIwXQtlkEj1wwmI17M3lg/saKX1CGhJRMJjzzPQs3pPDPs/ry6B8H26KxIeaknlFMOb4zr363ja83plbrXKrKrJ92MP7p70jNzGPGJcO5cWxvGzYOUPmFxdw9ex03vb+aoZ1a8sk1J9qtKDUgpJbbPblPWy49sSszvt/GST3bMKZvdKVep6p8smoXd3y8jsYNwph5+QjbFyqE3XZGH5Zu3c+lry/j0hO78vexvWnU4Ojmjnbsz+GBBRuYt3YvJ/ZozeN/HELbSFuEOFClHsrlmrdXsGz7AaaN6sYtp/e24doaElKNFMCtZ/Rm6db9XD9rFVeN7s4lJ3ShScOy/wyqysL1KTzxZQLr9xxiWOeWPDt5qO3fE+IaNwjno6tP5MH5G3n1u218tTGVhyYO4tgurSp8bWJqFs99k8jsVbsJE+GWcb25clR3G97zY/mFxWzdl0Xv6Ga/6eXmFhTxxpLtPPfNFvIKinn6wmP4w+B23gQapEKukWoYHsZLU4Zxz5x4Hv58E69+t42rR3fnzyM70zC8HhmHC0jNzGNzSiYv/G8L63YdonPrxjxywWDOGdLOro4M4KSlTz9nAGcMjOGWD9bwxxeXcOaAWIZ1bsngjs3pF9ucRg3CSM/OJzE1i4TUTL5P3Mf8dXuJCA/jkhO6cMVJ3eyCJwAs3bqfKTN+omOrRpzaN5rT+kYztHNLPlqxi6cWJbD3UC6jekXxz7P60ivaFoqtaVLRvUMiMgMYD6Sq6gC37GHgD0A+zpbvf1HVgyJSH2cr+KE4DeAbqvof9zXjgCeBMOAVVX3ALe8KzAJaASuAi1U1v6LAhw8frnFxcUf/iX2s2HGAx77YzHeJ+2jaMJz8omLyC3+5WbNTq8Zcd0oPzj2mfUg1TiKyXFWHex1HTaqJ+lKe7LxCHv1iM/PW7mHvoVwA6glENqrPwZyCn49r3qg+fx7ZiUtP7ErrpsG1b1Cw1Rnf+pKenc/n8XtZuD6F7xL3kV9YTFg9oahYGdqpBbeM62PD/0fpaOpLZRqpUUAWToNT0kiNBb5S1UIReRBAVW8VkcnA2ao6SUQaA+uB0cBOYDNwGpAMLAMuVNX1IvIe8JGqzhKRF4DVqvp8RYHX5JfOki37mbN6N5ER4bSNjKBts4bENI9gSMcWIXkjZbB94UDtNlK+Ug7lsjY5gzW7MkjLzKN7VBN6tG1Kz+hmxEZGBO2wXrDVmfLqS05+IYs372Pp1v2c2KMNp/Zta4kuVXA09aXC4T5VXSwiXUqVfeHzcCkwseQpoImIhAONcHpah4DjgERV3eoGOAuYICIbgFOAye7rXwfuASpspGrS8d1b230MpkZER0YQ3S+CU/tVLinHBJbGDcIZNyDGbsytQzXRTbgUmO/+/gGQDewBdgCPqGo60B6nN1Ui2S1rDRxU1cJS5WUSkWkiEicicWlpaTUQujHGGH9WrUZKRO4ACoGZbtFxQBHQDugK3CQi3YCy+sN6hPIyqepLqjpcVYdHRdk+LMYYE+yqnN0nIlNxEirG6C8TW5OBBapaAKSKyPfAcJxeVEefl3cAdgP7gBYiEu72pkrKK3rvacA+EUmqavx1pA3OZ/R3peMMuo2Lli9fniUiNbdMRO2xOuMHrL7UuCrXlyo1Um6m3q3A71U1x+epHcApIvIW0BgYCTyBk0DR083k2wVMAiarqorI1zhzWrOAqcDsSoQwTVX9vislInGBMJkcKHFW06ZA+IyB8t8iUOKsBqsvNag6cVY43Cci7wBLgN4ikiwilwHPAM2AhSKyys3KA3gWaAqsw8ng+6+qrnF7SdcCnwMbgPdUNd59za3AjSKSiDNH9WpVPogxxpjgU5nsvgvLKC6zIVHVLOCCcp6bB8wro3wrzlyWMcYY8yuBehPQS14HUEkWp/8IlM9ocfqHQPl8QR9nhTfzGmOMMV4J1J6UMcaYEGCNlDHGGL/l142UiIwTkU0ikigit5XxfEMRedd9/sfSyzfVUYwdReRrEdkgIvEicn0Zx4wWkQw3E3KViNxV13G6cWwXkbVuDL9ZmEwcT7l/zzUiMtSLOKvD6kyNxxrUdcbqS43HWvP1RVX98gdntfQtQDegAbAa6FfqmKuBF9zfJwHvehBnLDDU/b0ZzkK6peMcDcz1g7/pdqDNEZ4/E2eJK8G5x+1Hr2O2OuP53zRo64zVl8CoL/7ck/p5UVp1tu6YBUwodcwEnEVpwVk3cIzU8ZLEqrpHVVe4v2fi3AdW7vqDfm4Czmr3qqpLcVYDifU6qKNgdabuBXKdsfpS9466vvhzI1XeorRlHqPODcMZODcEe8IdCjgG+LGMp48XkdUiMl9E+tdpYL9Q4AsRWe4uLVVaZf7m/szqTM0L5jpj9aXm1Xh98eedeSuz+OxRLVBbm0SkKfAhcIOqHir19Aqgs6pmiciZwCdAz7qOEThRVXeLSFuc1UI2qupin+f95u9ZRVZnal4w1xmrLzWvxuuLP/ekkil7UdoyjxFnD6vmQHqdROdDnB2JPwRmqupHpZ9X1UPqrMaBOitv1BeRNnUcJqq62/03FfiY3670UZm/uT+zOlPDgrzOWH2pYbVRX/y5kVqGuyitiDTAmbScU+qYOTiL0oKzSO1X6s7O1RV3fPpVYIOqPlbOMTEl49gichzO331/3UUJItJERJqV/A6MxVlj0dccYIqbgTMSyFDVPXUZZzVZnalBIVBnrL7UoNqqL3473KfO1vQli9KGATNUNV5E7gXiVHUOzn+4N8VZnDYdp5LVtROBi4G1IrLKLfsH0AlAVV/AqdxXiUghcBiYVNcVHYgGPnbrcTjwtqouEJErfeKch5N9kwjkAH+p4xirxepMjQvqOmP1pcbVSn2xZZGMMcb4LX8e7jPGGBPirJEyxhjjt6yRMsYY47eskTLGGOO3rJEyxhjjt6yRMsYY47eskTLGGOO3rJEyxhjjt6yRMsYY47eskTLGGOO3rJEyxhjjt6yRMsYY47eskTLGGOO3rJFyich2ETm1VNklIvKdVzGVJiIni8haETkoIvtF5GMRCZStuoNKINQXXyLyXxFREenhdSzGHA1rpALLeuB0VW0BtAMSgOe9Dcn4OxH5HdDd6ziMqQprpI6CiPQVkW/cnky8iJzt89w3InK5z+Ofr6rdXSgfF5FUEckQkTUiMsB9rqGIPCIiO0QkRUReEJFGZb2/qqaUbM/sKgLsythPeV1f3OPDgaeBa2vvkxpTe6yRqiQRqQ98CnwBtAWuA2aKSO9KvHwsMAroBbQA/sQvWzs/6JYPwWlw2gN3HSGOTiJyEGf3zb8DD1Xl85ja5S/1BfgbsFhV11ThYxjjOWukfu0T96r3oNsQPOfz3EigKfCAquar6lfAXODCSpy3AGgG9MHZDXmDqu4RZ5/lK4C/qWq6qmYC/+YIW1Sr6g53uK8N8E9gYxU+p6kZfl1fRKQj8H8cuREzxq9ZI/Vr56hqi5If4Gqf59oBO1W12KcsCedK9ojcL6hngGeBFBF5SUQigSigMbDc54tugVte0TnTgdeB2e6Qjql7/l5fngDuVdWMo/5kxvgJa6QqbzfQUUR8/2adgF3u79k4XyAlYnxfrKpPqeowoD/OcM3NwD6cYbv+Pl92zVW1aSVjCscZSoo86k9japs/1JcxwMMisldE9rplS0RkcrU+mTF1yBqpyvsR54vlFhGpLyKjgT8As9znVwHniUhjN833spIXisixIjLCnafIBnKBIvcq+2XgcRFp6x7bXkROLysAETlPRHqLSD0RiQIeA1a6vSrjXzyvLziN22Cc+ashbtkfgI9r8HMaU6uskaokVc0HzgbOwLmifQ6Yoqolc0KPA/lACs4w3Eyfl0fifLkcwBny2Q884j53K5AILBWRQ8CXQHmT6+1xhncygbVAMXBuDXw8U8P8ob6oaqqq7i35cYv3qerhmvmUxtQ+UVWvYzDGGGPKZD0pY4wxfssaKWOMMX7LGiljjDF+yxopY4wxfitgbwJt06aNdunSxeswgtLy5cv3qWqFNxQHEqsvtSsY64zxD9VqpETkepxlWgR4WVWfEJFWwLtAF2A78EdVPeAu6fIkcCaQA1yiqivc80zFWeIH4D5Vfb2i9+7SpQtxcXHVCd+UQ0SSvI6hpll9qV3BWGeMf6jycJ+7KvMVwHE4NwyOF5GewG3AIlXtCSxyH4Nzv0hP92ca7hYTbqN2NzDCPdfdItKyqnEZY4wJHtWZk+oLLFXVHFUtBP6Hc2PpBJybE3H/Pcf9fQLwhjqWAi1EJBY4HVjoLph5AFgIjKtGXOYo7Dp4mJcWbyG/sLjig40BXvjfFpL2Z3sdhgkR1Wmk1gGjRKS1iDTGGcbrCESr6h4A99+27vHtgZ0+r092y8or/w0RmSYicSISl5aWVo3QTYn7P1vPYws3sy8rz+tQTAD4NiGNB+ZvZO6aPV6HYkJElRspVd2As7fNQpylelYDhUd4iZR1miOUl/WeL6nqcFUdHhVlc7TV9X3iPuat3cs1o3vQrkW5++YZA0BBUTH/+nQ9nVs35rLfdfU6HBMiqpWCrqqvqupQVR0FpONsZ57iDuPh/pvqHp6M09Mq0QFnpejyyk0tKigq5u458XRq1ZgrRnWrk/cUketFZJ27S+0Nbtm7IrLK/dkuIqvc8i4ictjnuRd8zjNMRNaKSKKIPOUm5Zha9saSJBJTs7jzrH5E1A/zOhwTIqrVSPmsxNwJOA94B5gDTHUPmQrMdn+fA0xxt8YeCWS4w4GfA2NFpKWbMDHWLTO16PUftpOYmsVd4+vmC6e8RBtV/ZOqDlHVIcCHwEc+L9tS8pyqXulT/jxO8k1JIo7NYdayfVl5PLFwM6N6RTGmb9uKX2BMDanufVIfikhrnJ1Er3FTzR8A3hORy4AdwAXusfNw5q0ScVLQ/wLO5n0iMh1Y5h53r209UbtSM3N54ssERveu0y+cnxNtAESkJNHmIfexAH8ETjnSSdzeeaSqLnEfv4GTnDO/9kI3Dy/YxOGCIu4a3w/ruJq6VK1GSlVPKqNsP85ma6XLFbimnPPMAGZUJxZTeQ/O30ReYZ1/4awD7ncvag7jXLD43rh0EpCiqgk+ZV1FZCVwCPinqn6Lk1ST7HPMERNtcHpcdOrUqaY+R8hZvfMg7y3fyeW/60qPtpXdj9OYmhGwK06YqlmedIAPVyRz1ejudIuquy8cVd0gIiWJNln8NtHmQpzh4hJ7gE6qul9EhgGfiEh/jjLRBngJYPjw4bYnTRUUFyv3fBpP6yYN+euYnl6HY0KQrd0XQoqKlbvnrCM6siHXntyjzt+/nEQbRCQcZ07zXZ9j89xeOaq6HNiCs9NsMk5yTQlLtKlFH63cxcodB7l1XG+aRdT3OhwTgqyRCiFv/7SDdbsO8c+z+tGkYd13ostJtAE4Fdioqsk+x0aJSJj7ezecBImtbrJNpoiMdOexpvBLco6pQYdyC3hg/kaO6dSC84d2qPgFxtQCG+4LEenZ+Tzy+SZO6N6a8YNivQrjN4k2bvkkfj3UBzAKuFdECoEi4EqfhJqrgNeARjgJE5Y0UQue/DKB/dl5zLhkOPXqWbKE8YY1UiHioQUbyc4r5F9n9/csO6usRBu3/JIyyj7ESUkv6/g4YECNBmd+ZXNKJq/9sJ1Jx3ZkUIcWXodjQpgN94WAVTsP8m7cTv5yYhd6RjfzOhzj51SVe+bE07RhODef3sfrcEyIs0YqyBUVK3fNXkdU04Zcf2ovr8MxAWD+ur38sGU/N43tRasmDbwOx4Q4a6SC3LvLdrImOYM7zupLUw+SJUxgyckv5L656+kbG8nk4+zeMuM9a6SCWHp2Pg99vpERXVtx9uB2XodjAsBzX29hd0Yu/zq7P+Fh9vVgvGe1MIg9OH8jWbmFTD9ngC1lYyq0bV82Ly3eyrnHtOe4rq28DscYwBqpoLU8KZ1343Zy2e+60suSJUwFVJV/fRpPg/B63H6GJUsY/2GNVBAqLCrmn5/EE9s8wpayMZXyxfoUvtmUxt9O60XbyAivwzHmZ9ZIBaE3liSxYc8h7hrvzcoSJrAczi/i3k/X0zu6GVOP7+x1OMb8in2DBZnUQ7k8tnAzv+8VxbgBMV6HYwLA898ksuvgYd6dNtKSJYzfsRoZZO77bAP5RcWerixhAsf2fdm8sHgr5wxpx4hurb0Ox5jfsEYqiHyXsI85q3dz1e+706VNE6/DMX7u52SJsHr848y+XodjTJmskQoSuQVF3Dl7HV3bNOGq0d29DscEgM/jU/jakiWMn7M5qSDx/Ddb2LYvm7cuG0FE/TCvwzF+LjuvkH99Gk/f2EhLljB+rVo9KRH5m4jEi8g6EXlHRCJE5FsRWeX+7BaRT9xjR4tIhs9zd/mcZ5yIbBKRRBG5rbofKtRsScvi+W+2cM6QdvyuZxuvwzEB4KlFCezJyOW+cwZYsoTxa1XuSYlIe+CvQD9VPSwi7wGTfLdjEJEP+fWGdN+q6vhS5wkDngVOw9l1dZmIzFHV9VWNLZSoKnd+so6I+vW446x+XodjAsCmvZm8+t02Jh3bkWGdW3odjjFHVN1LqHCgkbv9d2N8tvEWkWbAKcAnFZzjOCBRVbeqaj4wC5hQzbhCxierdvHDlv3cekYfopo19Doc4+dKLmqaRYRz6zhbWcL4vyo3Uqq6C3gE2AHsATJU9QufQ84FFqnqIZ+y40VktYjMF5H+bll7YKfPMclu2W+IyDQRiRORuLS0tKqGHjQO5uRz39wNHNOpBRceaytWm4p9tGIXP21P57Yz+tDStuEwAaDKjZSItMTp8XQF2gFNROTPPodcyK+3BF8BdFbVwcDT/NLDKutmHi3rPVX1JVUdrqrDo6Kiqhp60Hhg/kYOHi7g/nMG2vbepkIHc/L597wNDO3UgguGdfQ6HGMqpTrDfacC21Q1TVULgI+AEwBEpDXOMN5nJQer6iFVzXJ/nwfUF5E2OD0n3/9jOuAzbGjK9uPW/cxatpPLT+pKv3aRXodTKSJyvZtkEy8iN7hl94jILp+EmjN9jr/dTabZJCKn+5Rbok0V/HxRc65d1JjAUZ0U9B3ASBFpDBwGxgBx7nMXAHNVNbfkYBGJAVJUVUXkOJwGcj9wEOgpIl2BXcAkYHI14gp6eYVF3P7xWjq2asQNYwJjt10RGQBcgXPxkg8sEJGSi5jHVfWRUsf3w6kL/XF66l+KSMmHtUSbo7Rsezqzlu1k2qhu9I0NjIsaY6AajZSq/igiH+AM4xUCK4GX3KcnAQ+UeslE4CoRKcRp1CapqgKFInIt8DkQBsxQ1fiqxhUKnv9mC1vTsnn90uNo1CBg7onqCyxV1RwAEfkfzrxleSYAs1Q1D9gmIok4DRy4iTbueUoSbayRKkd+YTF3fLyW9i0accOptiq+CSzVuplXVe8G7i6jfHQZZc8Az5RznnnAvOrEEioSU7N47ustnD24Hb/vFVDzcuuA+92h4MPAmTg97/3AtSIyxX18k6oewEmeWerzet+EmtKJNiPKekMRmQZMA+jUKXQTS175biubU7J4ZcpwGjew+/dNYLG7+AJIcbHyj4/XElG/HneOD6x7olR1A/AgsBBYAKzG6YE/D3QHhuBkiT7qvqS8hBpLtDkKO9NzeGpRAqf3j+bUftFeh2PMUbNGKoC8F7eTn7al848z+wbkPVGq+qqqDlXVUUA6kKCqKapapKrFwMv8MqRXXkKNJdpUkqpy1+x1hIlwz9n9K36BMX7IGqkAkXool/vnbeC4rq344/DATB8Wkbbuv52A84B3RCTW55BzcYYFAeYAk0SkoZtU0xP4CViGm2gjIg1w5j/n1NVnCCRz1+zh601p3Di2N7HNG3kdjjFVYgPUAeKu2fHkFRbzwHkBnT78oTsnVQBco6oHRORNERmCM2S3Hfg/AFWNd5faWo8zLHiNqhYBWKJNxQ7m5POvT+MZ1KE5l5zQxetwjKkya6QCwIJ1e1gQv5dbxvWmW1RTr8OpMt91HX3KLj7C8fcD95dRbok2FfjPvI0cyCngjUtHEBa4FzXG2HCfv8s4XMBds+PpFxvJFSd18zocEwCWbNnPu3E7ueKkbgFzo7cx5bGelJ97YP4G9mXl8erUY6lvWyqYCuQWFPGPj9fSqVVjrh9j90SZwGffen5syZb9vPOTc0U8sENzr8MxAeCZrxLZti+bf587MJBu9DamXNZI+amc/EJu+2gNnVs35oZTA2PpI+OtDXsO8cL/tnDe0Pa2+aUJGjbc56ce/nwTSftzmDVtpF0RmwoVFhVzywdraN6oPv+0zS9NELFGyg/FbU/ntR+2M+X4zozs1trrcEwAeOW7bazdlcGzk4fSyvaJMkHEhvv8TG5BEbd8sIb2LRrZzqmmUramZfH4ws2c3j+aMwfGeB2OMTXKelJ+5rGFm9m6L5uZl4+gSUP7z2OOrLhYufXDNTQMr8f0CQMQsXuiTHCxnpQfWbHjAK98u5XJIzpxYg+b+DYVe3NpEsu2H+DO8f1oGxnhdTjG1DhrpPxEbkERN7+/mpjICG4/w4b5TMV2pufw4IKNjOoVxcRhHbwOx5haYeNJfuKhBZvYkpbNW5eNoFlEfa/DMX6uZJhPgH+fa8N8JnhZT8oP/LBlHzO+38aU4zvb/S2mUmb+tIMftuznjrP60aFlY6/DMabWWCPlsczcAm5+fw1d2zThNhvmM5WwMz2H/8zbwEk923DhcYG5bYsxlVWtRkpE/iYi8SKyTkTeEZEIEXlNRLaJyCr3Z4h7rIjIUyKSKCJrRGSoz3mmikiC+zO1uh8qkEyfu549GYd59I+DbWtvU6HiYuXmD1YTJsKD5w+yYT4T9KrcSIlIe+CvwHBVHYCzt88k9+mbVXWI+7PKLTsDZ+O6nsA0nG3DEZFWwN3ACJxdWe8WkZaVjSM9O58b31vFdwn7qvpRPPPl+hTei0vmyt93Z2inSn9kE8LeXJrE0q3p/HN8X9q1sI0MTfCr7nBfONBIRMKBxhx5G+8JwBvqWAq0cHdlPR1YqKrpqnoAWAiMq2wAjRuEMXvVbpZu3V/1T+GBfVl53PbRWvrGRtrafKZSkvZn88D8jfy+V1TA7s5szNGqciOlqruAR4AdwB4gQ1W/cJ++3x3Se1xEGrpl7YGdPqdIdsvKK6+UiPphdGvThI17D1Xxk9Q9VeXWD9ZwKLeAx/80mAbhNjVojqyoWLnpvdWEhwkPnD/QhvlMyKjOcF9LnN5RV6Ad0ERE/gzcDvQBjgVaAbeWvKSM0+gRyst6z2kiEicicWlpaT+X94mNZMOezKp+lDr35tIkFm1M5fYz+tAnxjalMxV74X9biEs6wPQJA4htbsN8JnRU5xL+VGCbqqapagHwEXCCqu5xh/TygP/izDOB00PyHaPogDM8WF75b6jqS6o6XFWHR0VF/VzeN7YZuw4eJuNwQTU+Tt3YnJLJ/Z9tYHTvKC45oYvX4dQpEbneTbKJF5Eb3LKHRWSj2/P+WERauOVdROSwTwLOCz7nGSYia90knKckyLsVa5MzeHzhZs4aFMuEIe28DseYOlWdRmoHMFJEGrtfEmOADe48E27ZOcA69/g5wBQ3y28kzvDgHuBzYKyItHR7Z2Pdskrr6/ZGNu31795UbkERf31nJc0iwnl44uCQGrIRkQHAFTgXLYOB8SLSE2cOcoCqDgI24/TES2zxScC50qf8eZzkm5JEnErPYQaaw/lF3PDuSto0bcj959hNuyb0VGdO6kfgA/LgFoQAAA/xSURBVGAFsNY910vATBFZ65a1Ae5zXzIP2AokAi8DV7vnSQemA8vcn3vdskrrG+s0Uv4+L/XQgk1s3JvJwxMHE9WsYcUvCC59gaWqmqOqhcD/gHNV9Qv3McBSnJ50udyLoEhVXaKqCryBczFUaU8vSuD2j9Ye/SfwwIMLNrIlLZtHLhhMi8a2BYcJPdW6MUdV78ZJH/d1SjnHKnBNOc/NAGZUNY7oyIa0aFyfDXv8t5H6emMqM77fxtTjO3Nyn7Zeh+OFdTgJNa2Bw8CZQFypYy4F3vV53FVEVgKHgH+q6rc4STXJPseUm2gjItNwelx06tTp5/LUzDzmrNrF9An9CQ/z36SVxZvTeO2H7fzlxC62EokJWf77f+hREBH6xDTz2+SJ3QcP87f3VtE3NpLbz+zrdTieUNUNwIM4w3sLgNVASQ8KEbnDfTzTLdoDdFLVY4AbgbdFJJKjSLQpbw5zeJeWZOcXsdGPh4f3ZeVx0/ur6dm2qe0rZkJaUDRS4Az5bdqbSXFxmd9XnikoKua6d1ZSUFjMcxcNJaJ+6G4Fr6qvqupQVR0FpAMJ4Kw4AowHLnJ73Khqnqrud39fDmwBeuH0nHyHBMtNtCnP8C6tAGcHZH9U7KabZxwu4OnJx4R0nTEmeBqpmEgOFxSRlJ7jdSi/8sjnm1iedID/nD+Irm2aeB2Op0SkrftvJ+A84B0RGYdzm8LZqprjc2yUiIS5v3fDSZDY6ibbZIrISDc5Zwow+2jiaN+iEe2aRxCXdKBGPldNe/W7bfxvcxp3ju9ntyiYkBc0i8X9nDyx55DfNAaLNqTw4uKtXDSiE2cPttRh4EN3TqoAuEZVD4jIM0BD/r+9e4+uqroTOP795QFJIIRAYkKAhNQElEcRDC9R1GF4uVg6nSkVKY9aQEGZ2tJlx2VbXcPUaWuLdVqZKgVsRbBVXsUKlrS1SpcmGsDII4T3K0kBCQkJECDh1z/OCU3TRAL35t5zb36fte7K5Zx9z/0l7Ht+9+y9z96Q645cy3NH8o0CFohILVAHzGkwoGYu8EsgFtjoPq5JTq8u5B88hap6asRc4dEKfvj2bsb1S2HqsPSrv8CYMBc2SSo7pSMRAkV/rWLCgG7BDoeSivN8841Cbu7Wie9O7BvscDxBVe9oYltWM2VXA6ub2VcA9PcllpxeiawvLOXY6fP07OKNpS6qai7xn69tI6VTDM/+R9u6RcGY5oRNc19MdCSZSR08McKv5lIdc5ZvobZO23w/lFflZLj9Uoe90S+lqnx77Q5KKs7zf5NvISHOFr40BsIoSYHT5Bfse6VUle+s28H2kkqe+9JAzzQ9mn/UJzWe+PZRFBzyRr/UL98/xPrCUuaP6X1lYIcxJgyT1NHy81TVBG96pFfzj7BqyzG+9i9ZjO2XGrQ4zGeLjBAGZSR6Ikl9sP8U33uriLF9U5h7543BDscYTwmrJHVTajwQvOmRthwuZ8GbO7m7T7ItvxEChmQkUny8ispzwftSU1Jxnnkrt9KraxwLvzSQiAjrhzKmobBKUvUj/IqCkKROnKlh7qtbSescy/P3D7KTTQiob1bbeiQ4V1P1fZcXai+zeHoO8THWD2VMY2GVpLolxNApJirggydqLtUxe/kWqmpqeWnardbpHSJu6dmZqAjhoyDc1Fs/UGJ7SSU/uf8WbkzuGPAYjAkFYTMEHZzpkW7u1ondAUxS9bMDfHKsgpem3mo3X4aQ2HaR9OueEJR+qZ+/u5/VW4/x2OhsxvRNCfj7GxMqwupKCgI/PdLC3GLe2l7GkxNutoESIWhIRiKFxyq4UFsXsPdct62EZ98u5t6BaTw2Ojtg72tMKAq7JHVTajxnL9Zx9HTrT4/0RsFRFr2znweGpjPrjsxWfz/jfzm9unCh9jI7SgJz9f3+vk95fFUhwz/XhR9N+rz1XRpzFWGXpK4MnmjlGdE/2H+KJ9du5/asJBbc189mBwhRt2YkAoGZbHb3X8/w8PItZCZ14KVpObSPspu8jbmasEtSvVPiiZDWXQCx8GgFs18pIKNrBxZ9eTDRHl6TyHy25Pj2ZCZ1aPXJZksrzvPgyx8R1z6Slx8cSkKsDa4xpiXC7uwa2y6SvmmdeLOwlNq6y34//q7SM0xf9iGJHaJZPtNONuEgJyORgkPlrVJfwElQD/wij+qaWpZ9ZQjdO8e2yvsYE47CLkkBzLs7i/0nz7JmW4lfj7v3eBVTl+YT1y6SlbOG0y3BTjbhYMKAVE6fu8SK/CN+P3Z9giqvvsgrM4fSLy3B7+9hTDgLyyQ1rl8qA3sk8HzuHmou+WfU1sFPzzJlST6REcLK2cM9M3O28d3dfW7g9qwknsvdw+mzF/123MYJalB6ot+ObUxb4VOSEpFviMhOEdkhIq+JSIyIrBCRYnfbMhGJdsveJSKVIvKx+3iqwXHGu6/ZJyJP+PpLiQjfGn8TpZU1fvl2vKv0DA8szqPusrJy1jCbNDbMiAjfndiX6gu1/OQPe/xyzBJLUMb4xXUnKRHpDnwNyFHV/kAkMBlYAdwEDMBZlG5Wg5dtVtVb3McC9ziRwCJgAtAXeEBEfF6AaWRWErdnJbHonX1UX6i97uP8ufgEk158H4CVs4eRnRLva2jGg/qkxjN1WDqv5h32edDNhwfLue+Fv1iCMsYPfG3uiwJiRSQKiANKVXWDuoAPgR5XOcZQYJ+qHlDVi8Cvgft8jAuAx8f1ofzsRZZsPnBdr1+Zf4SZv3JG8a17dKTNJhHmvjGmN51io1nw5i6c6nttVJVXPjjElF/kER8TzZpHbrMEZYyPrjtJqWoJ8GPgCFAGVKrqpvr9bjPfNODtBi8bISKFIrJRRPq527oDRxuUOeZu+yci8pCIFIhIwcmTJ68a48CenZnQP5Ulmw9Sfg19DZfqLvP9jUU8uXY7d2Qn8fqcEaQmxLT49SY0dY5rx/wxvXl//yk27Tp+Ta+tuVTHt1Z9wlO/3cmo3smse3SkXXUb4we+NPcl4lzxZAJpQAcRmdqgyP8D76nqZvffW4EMVR0I/AxYV3+oJg7f5NdYVV2sqjmqmpOcnNyiOL85tjfnLtby/Q1FLZr6ZvPek4x//j1eevcAXx6WzpLpOXRsH1ZTHAaViDzm9lfuFJGvu9u6iEiuiOx1fya620VEfur2VX4iIoMbHGeGW36viMzwV3xThqbTO6Uj33trFyerLly1/OXLyobtZdzz08284a4jtmR6jt2aYIyf+NLc96/AQVU9qaqXgDXAbQAi8jSQDMyvL6yqZ1S12n2+AYgWkSScK6eeDY7bAyj1Ia5/kHVDPF8dmckbW44x5rn3+N0npU025RwtP8ec5VuYtvRDai8rS2fk8MwXBhBlN+r6jYj0B2bjNPEOBCaKSDbwBPBHVc0G/uj+G5x+ymz38RDwc/c4XYCngWHusZ6uT2y+ioqM4L/v7U9ZRQ13/ugdFm4q5kwTi2iqKu/uOcm9i/7CIyu2EiHCy18ZwvyxfWyqI2P8yJdLhCPAcBGJA84Do4ECEZkFjANGq+qVuyNFJBU4rqoqIkNxEuQpoALIFpFMoARn8MUUH+L6J9+Z2JdRvZP53w1FzFu5jSU9DzJ9RAZllTXsP1HNvpPV7C6rIjJCeHxcH2benklMtE1Z0wpuBvJU9RyAiLwLfAHnivwut8yvgD8D/+Vuf8Xt38wTkc4i0s0tm6uq5e5xcoHxwGv+CHLEjV3JnX8nCzcV87M/7ePVvMM8fOeNJMZFc/jUOQ6Xn2Pv8Sr2HK+mR2IsCycN5N8GdSfSkpMxfnfdSUpV80VkFU4zXi2wDVgMnAUOAx+489mtcUfyfRGYKyK1OEltsnvyqRWRecDvcUYILlPVnT78Tk0a1TuZkVlJrNl6jB9vKmb+64UApHaKIeuGjsy4LYMHR2aSZrMBtKYdwDMi0hWnDtwDFAApqloGoKplInKDW765/soW9WOKyEM4V2Ckp6dfU6CZSR14YcpgHh5VybO/380PNu4GICpC6JEYS3rXDkwbnsH9Q9JpF2VX28a0Fp86W1T1aZxml6seU1VfAF5oZt8GYIMvsbREZIQwKacnEz+fxoFPq0nvEmeroQaQqhaJyA+BXKAaKMT5gtOc5vorW9SPqaqLcb44kZOTc11rtwzokcDymcPYe7yK9lGRpHWOsSZgYwKoTX7aYttF0i8twRJUEKjqUlUdrKqjgHJgL3DcbcbD/XnCLd5cf2Wr9mM2JTslnvSucZagjAkw+8SZgKpvyhORdODfcfqR1gP1I/RmAL91n68Hpruj/Ibj3OZQhtM0PFZEEt0BE2PdbcaYMGNjq02grXb7pC4Bj6rqaRH5AfC6iMzEGZAzyS27Aaffah9wDngQQFXLReR/gI/ccgvqB1EYY8KLXM+d9cHmdog/gzNAw8uSgE+DHUQLNI4zQ1VbdiNaiBCRKqA42HG0gNUZYxoI1Suph0LhAyEiBaqaE+w4riZU4vRRcSj8jqHyfxEqcZrQZ31SxhhjPMuSlDHGGM8K1SS1ONgBtJDF6R2h8jtanMY0EJIDJ4wxxrQNoXolZYwxpg2wJGWMMcazPJ2kRGS8iBS76wk90cT+9iLyG3d/voj0CkKMPUXkHREpctdIeqyJMneJSKWIfOw+ngp0nG4ch0RkuxtDQRP7m12/KVRYnfF7rGFfZ4zHqaonHzgzou8HPge0w5mMtG+jMo8AL7rPJwO/CUKc3YDB7vN4YE8Tcd4F/M4Df9NDQNJn7L8H2IgzgetwID/YMVudCfrfNKzrjD28//DyldRQYJ+qHlDVi8CvcdYXaug+nPWHAFYBo8VdHyRQVLVMVbe6z6uAIppYNiJEXFm/SVXzgPr1m0KF1ZnAC/U6YzzOy0mqJWsGXSmjqrVAJdA1INE1wW06GgTkN7F7hIgUishGEekX0MD+ToFNIrLFnVqqsRat0+RhVmf8L9zrjPE4L0+L1JI1g1q0rlAgiEhHYDXwdVU902j3Vpy5zapF5B5gHc6S6IE2UlVL3ZnIc0Vkt6q+12C/Z/6e18nqjP+Fe50xHuflK6mWrBl0pYyIRAEJOGsUBZSIROOcbFao6prG+1X1jKpWu883ANEikhTgMFHVUvfnCWAtTvNYQwFfp8nPrM74WRuoM8bjvJykPgKyRSRTRNrhdHKvb1Sm4TpEXwT+pKoB/Rbn9mcsBYpU9blmyqTW93uIyFCcv/upwEUJItJBROLrn+OswbSjUbHm1m8KFVZn/KiN1BnjcZ5t7lPVWhGZh7OYXSSwTFV3isgCoEBV1+N80JeLyD6cb8OTgxDqSGAasF1EPna3PQmkA6jqizgnw7kiUgucByYH+sQIpABr3fNeFLBSVd8WkTkN4mxy/aZQYXXG78K+zhjvs2mRjDHGeJaXm/uMMca0cZakjDHGeJYlKWOMMZ5lScoYY4xnWZIyxhjjWZakjDHGeJYlKWOMMZ71N1+Ku/8x3uqkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum cost: 43534.36377500013\n",
      "Genetic algorithm minimum: 43579.274812500225\n"
     ]
    }
   ],
   "source": [
    "min_sum = 0\n",
    "\n",
    "for l in range(num_houses):\n",
    "    data_consumption = df_cons.iloc[:,l].values\n",
    "    g = Grid(num_houses=1,num_storages=1,max_capacity=0,num_pvtypes=len(pv_peakpowers),pv_peakpower=pv_peakpowers)\n",
    "    g.set_costs()\n",
    "    total_costs = []\n",
    "    # For all possible PV types\n",
    "    for i in range(data_production.shape[1]):\n",
    "        # Set pv type of house to this specific one\n",
    "        g.change_pvtype(num_house=0,pv_type=i)\n",
    "        # Simulate the grid for one year\n",
    "        res_dict = g.simulate(data_cons=data_consumption,data_prod=data_production)\n",
    "        # Objective function\n",
    "        total_costs.append(objective_function(res_dict))\n",
    "    plt.subplot(2,3,l+1)\n",
    "    plt.plot(pv_peakpowers,total_costs)\n",
    "    plt.title('House ' + str(l))\n",
    "    min_sum += np.min(total_costs)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "print('Minimum cost: ' + str(min_sum))\n",
    "print('Genetic algorithm minimum: ' + str(min_cost_genetic))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see the minimum cost calculated manually is close to the optimal cost, calculated with the genetic algorithm. But the run of the genetic algorithm takes a huge computational amount.\n",
    "\n",
    "A next step might be adding a grid system with storages. Hence we would have optimized the size of the pv systems first and then checked if a grid system is feasible. We can also think of doing both a the same time which takes even more computational power, but might result in a better/cheaper solution of the energy cost problem."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}