Testing_Battery_scheduling_algorithms.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import csv\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "import matplotlib.cm as cmx\n",
    "import matplotlib.mlab as mlab\n",
    "import matplotlib.ticker as mtick\n",
    "import matplotlib.patches as mpatches\n",
    "import matplotlib.ticker as ticker\n",
    "import itertools\n",
    "import math\n",
    "import random\n",
    "import networkx as nx\n",
    "import colorsys\n",
    "import pandas as pd\n",
    "from pyomo.opt import SolverFactory\n",
    "from pyomo.core import Var\n",
    "import pyomo.environ as en\n",
    "\n",
    "import seaborn as sns\n",
    "from scipy import stats \n",
    "from six.moves import cPickle as pickle\n",
    "from mpl_toolkits.axes_grid.inset_locator import inset_axes\n",
    "from __future__ import division\n",
    "\n",
    "%matplotlib inline\n",
    "sns.set_style({'axes.linewidth': 1, 'axes.edgecolor':'black', 'xtick.direction': \\\n",
    "               'out', 'xtick.major.size': 4.0, 'ytick.direction': 'out', 'ytick.major.size': 4.0, \\\n",
    "              'axes.facecolor': 'white','grid.color': '.8', 'grid.linestyle': u'-', 'grid.linewidth': 0.5})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# testProfiles are 1 year hourly data\n",
    "testProfiles = np.genfromtxt('TestData.csv', delimiter=',')\n",
    "# prices are 7 different sctructures of price --> see paper\n",
    "prices = np.genfromtxt('./IntermediateData/testPrices.csv', delimiter=',')\n",
    "# entries 0-10 are demand\n",
    "# i.e. first demand profile is testProfiles[0]\n",
    "# entries 10-20 are generation\n",
    "# i.e. first gen profile is testProfiles[10]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# define the cycle charging scheduling algorithm\n",
    "# standard logic - charges battery with RE until full then discharges\n",
    "# includes efficiency and limited by demand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cycleCharging(noUnits, storageProperties, load, price, billingOption):\n",
    "    # the input should have generation as +ve\n",
    "    \n",
    "    # make sure the load is in kW, then storage and price are adjusted using width of time period\n",
    "\n",
    "    seriesLength=len(load)\n",
    "    # load defined as -ve when importing from grid and +ve when producing --> EE convention\n",
    "    demand = -load\n",
    "    original_demand = np.zeros((seriesLength))\n",
    "    original_demand[ np.where(load<0)[0] ] = -load[ np.where(load<0)[0] ]\n",
    "    \n",
    "    ### ------------------------- ESS properties ----------------- ###\n",
    "    maxSOC = noUnits*storageProperties[0]# kWh\n",
    "    # scale maxSOC appropriately\n",
    "    maxChg = noUnits*storageProperties[1] # kW\n",
    "    maxDisChg = noUnits*storageProperties[2] # kW\n",
    "    etaChg = storageProperties[3]\n",
    "    etaDisChg = storageProperties[4]\n",
    "    ### ---------------------------------------------------------- ###\n",
    "\n",
    "    ### ------------------------- CHOOSE POLICY OPTION ----------- ###\n",
    "    policyOption = billingOption[0]\n",
    "    ### ---------------------------------------------------------- ###\n",
    "    \n",
    "    priceSolarExport = billingOption[1]\n",
    "    solarExports = np.zeros((seriesLength))\n",
    "    exports = np.zeros((seriesLength))\n",
    "    buy_price = np.zeros((seriesLength))\n",
    "    for i in range(seriesLength):\n",
    "        buy_price[i] = price[i]\n",
    "        if demand[i] < 0:\n",
    "            solarExports[i]=-demand[i]\n",
    "            exports[i] = solarExports[i]\n",
    "            demand[i] = 0\n",
    "    # now alter the price in the periods where there is solar available (if FIT policy)\n",
    "    if policyOption == 1:\n",
    "        buy_price[ np.where(solarExports>0)[0] ] = priceSolarExport\n",
    "    \n",
    "    # get the storage profiles\n",
    "    SOC = np.zeros((seriesLength))\n",
    "    deltaSOC = np.zeros((seriesLength))\n",
    "\n",
    "    # calculate initial costs\n",
    "    oldPrice = np.zeros((seriesLength))\n",
    "    for i in range(seriesLength):\n",
    "        oldPrice[i]=(price[i]*demand[i] - solarExports[i]*buy_price[i])\n",
    "    oldTotalPrice = np.sum( oldPrice )\n",
    "\n",
    "    # cycle through the series\n",
    "    for i in range(seriesLength):\n",
    "        bottleneck = np.zeros((3))\n",
    "        if exports[i]>0:\n",
    "            # establish the min bottleneck for charging\n",
    "            bottleneck[0] = maxSOC - SOC[i-1]\n",
    "            bottleneck[1] = maxChg\n",
    "            bottleneck[2] = exports[i]*etaChg\n",
    "            # act on the min bottleneck\n",
    "            actual_bottleneck = np.min(bottleneck) \n",
    "            exports[i] = exports[i] - actual_bottleneck/etaChg\n",
    "            SOC[i] = SOC[i-1] + actual_bottleneck\n",
    "            deltaSOC[i] = 0 + actual_bottleneck\n",
    "        else:\n",
    "            # establish the min bottleneck for discharging\n",
    "            bottleneck[0] = SOC[i-1]\n",
    "            bottleneck[1] = -maxDisChg\n",
    "            bottleneck[2] = demand[i]/etaDisChg\n",
    "            # operate on the actual bottleneck\n",
    "            actual_bottleneck = np.min(bottleneck)\n",
    "            demand[i] = demand[i] - actual_bottleneck*etaDisChg\n",
    "            SOC[i] = SOC[i-1] - actual_bottleneck\n",
    "            deltaSOC[i] = 0 - actual_bottleneck\n",
    "            \n",
    "    newPrice = np.zeros((seriesLength))\n",
    "    for i in range(seriesLength):\n",
    "        newPrice[i] = (price[i]*demand[i] - exports[i]*buy_price[i])\n",
    "    newTotalPrice = np.sum ( newPrice )\n",
    "    #print oldTotalPrice - newTotalPrice\n",
    "    \n",
    "    #now calculate the new demand\n",
    "    newDemand = demand-exports\n",
    "    SOC = SOC\n",
    "    \n",
    "    return newDemand, oldTotalPrice, newTotalPrice, SOC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def priceDriven2(noUnits, storageProperties, load, price, billingOption):\n",
    "    price1 = np.copy(price)\n",
    "    seriesLength=len(load)\n",
    "    demand = -load\n",
    "    # -ve demand is grid import, +ve is generation\n",
    "    original_demand = np.zeros((seriesLength))\n",
    "    original_demand[ np.where(load<0)[0] ] = -load[ np.where(load<0)[0] ]\n",
    "\n",
    "    ### ------------------------- ESS properties ----------------- ###\n",
    "    maxSOC = noUnits*storageProperties[0]# kWh # scale maxSOC appropriately\n",
    "    maxChg = noUnits*storageProperties[1] # kW\n",
    "    maxDisChg = noUnits*storageProperties[2] # kW\n",
    "    etaChg = storageProperties[3]\n",
    "    etaDisChg = storageProperties[4]\n",
    "    ### ---------------------------------------------------------- ###\n",
    "    ### ------------------------- CHOOSE POLICY OPTION ----------- ###\n",
    "    policyOption = billingOption[0] # solar export price, but must be used locally first\n",
    "    marketArbitrage = billingOption[2]\n",
    "    ### ---------------------------------------------------------- ###\n",
    "    priceSolarExport = billingOption[1]\n",
    "\n",
    "    solarExports = np.zeros((seriesLength))\n",
    "    exports = np.zeros((seriesLength))\n",
    "    buy_price = np.zeros((seriesLength))\n",
    "    for i in range(seriesLength):\n",
    "        buy_price[i] = price1[i]\n",
    "        if policyOption == 1:\n",
    "            if demand[i] < 0:\n",
    "                solarExports[i]=-demand[i]\n",
    "                exports[i] = solarExports[i]\n",
    "                demand[i] = 0\n",
    "    # now alter the price in the periods where there is solar available (if req.)\n",
    "    buy_price[ np.where(solarExports>0)[0] ] = priceSolarExport\n",
    "    # changed later if not policy 1\n",
    "    # copy the buy price to create speed up\n",
    "    buy_price1 = np.copy(buy_price)\n",
    "    \n",
    "    # get the storage profiles\n",
    "    SOC = np.zeros((seriesLength))\n",
    "    deltaSOC = np.zeros((seriesLength))\n",
    "\n",
    "    # calculate initial costs\n",
    "    oldPrice = np.zeros((seriesLength))\n",
    "    for i in range(seriesLength):\n",
    "        oldPrice[i]=price1[i]*demand[i] - solarExports[i]*buy_price[i]\n",
    "    oldTotalPrice = np.sum( oldPrice )\n",
    "    \n",
    "    # boolean variables for charging/discharging availability\n",
    "    removeMINH = np.ones((seriesLength))\n",
    "    removeMAXH = np.ones((seriesLength))\n",
    "    # make faster --> we know dont want battery to discharge when there is export\n",
    "    removeMAXH[ np.where(solarExports>0)[0] ] = False\n",
    "\n",
    "    if policyOption==1:\n",
    "        # alter the price slightly to discharge as soon as load>PV (reality)\n",
    "        for j in range(len(price1)):\n",
    "            price1[j] = price1[j]+(8760-j)/100000000.\n",
    "            buy_price1[j] = buy_price1[j]+(j)/100000000.\n",
    "            \n",
    "    else:\n",
    "        # alter the price to make the algorithm faster if there are lots of similar prices (buy early sell late encouraged)\n",
    "        for j in range(len(price)):\n",
    "            #if demand[i] >= 0:\n",
    "            price1[j] = price1[j]+(8760-j)/100000000. #price1[j] = price1[j]+np.random.uniform(low=-0.5, high=0.5)/10000000\n",
    "            buy_price1[j] = price[j]+(j)/100000000. \n",
    "\n",
    "    while np.any(removeMAXH) == True:\n",
    "        #for iiii in range(1):\n",
    "        matrix = np.zeros(( 2, len ( np.where(removeMAXH==True)[0] ) ))\n",
    "        matrix[0,:] = np.where(removeMAXH==True)[0]\n",
    "        matrix[1,:] = price1[ np.where(removeMAXH==True)[0] ]\n",
    "        #print removeMAXH\n",
    "        #print removeMINH\n",
    "        \n",
    "        # find the maximum available price \n",
    "        indici = np.where( matrix[1,:]==np.max( matrix[1,:] ) )[0]\n",
    "        #print indici\n",
    "        maxh = np.int(matrix[0,indici[0]])\n",
    "        #print 'Maxh = ', maxh\n",
    "\n",
    "        # find the last hour before maxh when storage was full\n",
    "        r1 = np.where( SOC[0:maxh+1] == maxSOC )[0]\n",
    "        if r1.size == 0:\n",
    "            r1 = 0\n",
    "        else:\n",
    "            r1 = np.where( SOC[0:maxh+1] == maxSOC )[0][-1]+1\n",
    "        # find the first hour after maxh when storage is empty\n",
    "        r2 = np.where( SOC[maxh:] == 0 )[0]\n",
    "        if r2.size == 0:\n",
    "            r2 = len(SOC)\n",
    "        else:\n",
    "            r2 = np.where( SOC[maxh:] == 0 )[0][0]+maxh-1\n",
    "\n",
    "        #print 'r1 = ', r1\n",
    "        #print 'r2 = ', r2\n",
    "\n",
    "        # find the minh in the time range\n",
    "        range_price = buy_price1[r1:r2+1]\n",
    "        range_remove = removeMINH[r1:r2+1]\n",
    "\n",
    "        # if there is no hour in the range then remove maxh and skip to the end\n",
    "        if np.any(range_remove) == False:\n",
    "            removeMAXH[maxh]=False\n",
    "        else:\n",
    "            matrix = np.zeros(( 2, len ( np.where(range_remove==True)[0] ) ))\n",
    "            matrix[0,:] = np.where(range_remove==True)[0]\n",
    "            matrix[1,:] = range_price[ np.where(range_remove==True)[0] ]\n",
    "\n",
    "            indici = np.where( matrix[1,:]==np.min( matrix[1,:] ) )[0]\n",
    "            minh = np.int( matrix[0,indici[0]] + r1 )\n",
    "            # calculate the marginal cost of operation\n",
    "            MoC = buy_price[minh]/(etaChg*etaDisChg)\n",
    "            #print 'Minh = ', minh\n",
    "\n",
    "            if ((MoC<price[maxh]) and (minh!=maxh)):\n",
    "\n",
    "                bottleneck = np.zeros((5))\n",
    "                bottleneck[0] = deltaSOC[maxh]-maxDisChg\n",
    "                bottleneck[1] = maxChg - deltaSOC[minh]\n",
    "                if maxh > minh:\n",
    "                    bottleneck[2] = maxSOC - np.max(SOC[minh:maxh+1])\n",
    "                else:\n",
    "                    bottleneck[2] = np.min(SOC[maxh:minh+1])\n",
    "\n",
    "                ############### ARBITRAGE ###########################   cannot output more than local consumption\n",
    "                #if marketArbitrage == False:\n",
    "                if demand[maxh]>0:\n",
    "                    bottleneck[3] = demand[maxh]/etaDisChg\n",
    "                else:\n",
    "                    bottleneck[3] = 0\n",
    "                #else: # the market arbitrage allowed case is yet to be tested\n",
    "                #    bottleneck[3] = maxSOC\n",
    "\n",
    "                # if there is export initally cannot store more than export\n",
    "                if exports[minh] > 0:\n",
    "                    bottleneck[4] = exports[minh]*etaChg\n",
    "                else:\n",
    "                    bottleneck[4] = maxSOC\n",
    "                #print bottleneck\n",
    "                actual_bottleneck = np.min(bottleneck)  \n",
    "\n",
    "                # Update the aggregated demand\n",
    "                demand[maxh] = demand[maxh] - actual_bottleneck*etaDisChg\n",
    "                # if not buying from solar increase the demand at minh\n",
    "                if exports[minh] <= 0:\n",
    "                    demand[minh] = demand[minh] + actual_bottleneck/etaChg\n",
    "\n",
    "                # update the storage charging schedule\n",
    "                exports[minh] = exports[minh] - actual_bottleneck/etaChg\n",
    "                deltaSOC[minh] = deltaSOC[minh] + actual_bottleneck\n",
    "                deltaSOC[maxh] = deltaSOC[maxh] - actual_bottleneck\n",
    "                if maxh>minh:\n",
    "                    SOC[minh:maxh] = SOC[minh:maxh] + actual_bottleneck\n",
    "                else:\n",
    "                    SOC[maxh:minh] = SOC[maxh:minh] - actual_bottleneck\n",
    "\n",
    "                # account for small rounding errors\n",
    "                if ((demand[maxh] < 0 + 0.0001) and (demand[maxh] > 0 - 0.0001)):\n",
    "                    demand[maxh] = 0\n",
    "                if exports[minh] < 0 + 0.0001:\n",
    "                    exports[minh] = 0\n",
    "                    # update the price\n",
    "                    if policyOption==1:\n",
    "                        buy_price1[minh] = price[minh]+(minh)/100000000.\n",
    "                        # remember buy_price1 is the version used for scheduling and buy_price for calculating cost\n",
    "                        buy_price[minh] = price[minh]\n",
    "\n",
    "                # check if at the charge or discharge operation is at capacity at either\n",
    "                # maxh or minh and remove that hour from the price distribution.\n",
    "                if deltaSOC[maxh] <= maxDisChg:\n",
    "                    removeMAXH[maxh] = False\n",
    "                    removeMINH[maxh] = False\n",
    "                if deltaSOC[minh] >= maxChg:\n",
    "                    removeMINH[minh] = False\n",
    "                    removeMAXH[minh] = False\n",
    "\n",
    "                ######################## ARBITRAGE #####################\n",
    "                # check whether exports are allowed or not\n",
    "                #if marketArbitrage == False:\n",
    "                #    # if not remove periods where all demand is met\n",
    "                if demand[maxh] <= 0:\n",
    "                    removeMAXH[maxh] = False\n",
    "\n",
    "            else:\n",
    "                # if there is no price incentive remove the hours\n",
    "                removeMAXH[maxh] = False\n",
    "                removeMAXH[minh] = False\n",
    "\n",
    "    newPrice = np.zeros((seriesLength))\n",
    "    for i in range(seriesLength):\n",
    "        newPrice[i] = price[i]*demand[i] - exports[i]*buy_price[i]\n",
    "    newTotalPrice = np.sum ( newPrice )\n",
    "    newDemand = demand-exports\n",
    "    \n",
    "    return newDemand, oldTotalPrice, newTotalPrice, SOC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# These are the storage device properties\n",
    "cycles = 3000\n",
    "maxLi = 30\n",
    "installedCap = 1 # kWh\n",
    "maxDOD = 0.85 # % equivalent full cycle (EFC)\n",
    "EoL = 0.8 # end of life % equivalent full charge \n",
    "maxSOC = np.mean([EoL*installedCap, maxDOD*installedCap])# kWh\n",
    "minSOC = 0 # kWh\n",
    "maxChg = 0.5 # kW\n",
    "maxDisChg = -0.5 # kW\n",
    "etaChg = 0.9\n",
    "etaDisChg = 0.9\n",
    "\n",
    "# Make a list of the storage properties to pass to the functions\n",
    "storageProperties = []\n",
    "storageProperties.append(maxSOC-minSOC) \n",
    "storageProperties.append(maxChg) \n",
    "storageProperties.append(maxDisChg) \n",
    "storageProperties.append(etaChg) \n",
    "storageProperties.append(etaDisChg) \n",
    "\n",
    "# billingOption is a list that describes what pricing scheme in in effect \n",
    "# 1=FIT and 2=NEM\n",
    "# If a FIT is specified then an export price is required.\n",
    "billingOption = []\n",
    "billingChoice = 1 # 1 means feed in tariff\n",
    "priceExport = 0 # cents/kWh\n",
    "# choose whether device is allowed to buy and sell back to grid\n",
    "marketArbitrage = False\n",
    "billingOption.append(billingChoice) \n",
    "billingOption.append(priceExport)\n",
    "billingOption.append(marketArbitrage)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# most challenging case is highly variable prices and a FIT (prices[6] and billingOption[0]=1)\n",
    "# negative prices havent been tested yet, but no reason logic approach cannot work,\n",
    "# though with negative prices bugs may need ironed out\n",
    "load = testProfiles[8]-(testProfiles[13]*0.25)\n",
    "price = prices[6,0:len(load)]\n",
    "# 0 austin FIT\n",
    "# 1 austin TOU FIT\n",
    "# 2 San Diego TOU NEM\n",
    "# 3 PGE flat rate NEM\n",
    "# 4 PGE FIT\n",
    "# 5 Xcel Energy NEM\n",
    "# 6 Austin DA NEM\n",
    "billingOption[0] = 1 # 1 means feed in tariff, 2 means NEM\n",
    "billingOption[1] = 1.5 # solar export price if specifying policy option 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA60AAACUCAYAAAB4OM8nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8W9X5+PGPJG87duJsh0zGSQKEHUbYECibFuiALy2F\nFn5QSmmBQpkByl6FMgo0QFktOwRIAiQBAknIInuckD3kFe+p/ftDsiPZ2rqWZPt5v168iK+uzj06\nupLuc885zzF5PB6EEEIIIYQQQoh0ZE51BYQQQgghhBBCiFAkaBVCCCGEEEIIkbYkaBVCCCGEEEII\nkbYkaBVCCCGEEEIIkbYkaBVCCCGEEEIIkbYkaBVCCCGEEEIIkbYkaBVCCCGEEEIIkbYkaBVCCCGE\nEEIIkbYkaBVCCCGEEEIIkbYkaBVCCCGEEEIIkbYyUl2BUEwmkyfVdRBCCNGzeDweU6rr0J3Jb7MQ\nQgijRfPbnLZBK4DHk/hvo9VqpaSkxIDaiDbSpsaS9jSWtKexelJ7mkwSrxpBfpvTj7Sn8aRNjSXt\naaye1J7R/jZHFbQqpQ4CPgae0lo/q5QaDrwBWIBS4HKtta3Dc54CjgE8wJ+01kuir74QQgghhBBC\nCBHFnFalVD7wT2CO3+b7gOe01icAm4ArOzznJGB/rfWxwFXAM4bVWAghhBBCCCFErxFNIiYbcDZg\n9dt2MjDd9+9PgNM7POc0YBqA1no90E8pVZhQTXspa2Ujf3v+O3aWN6S6KkIIIYQQohtbvK6Mu19c\ngN3hSnVVhIhJxKBVa+3UWrd02JzvNxy4Ahja4fEhQKXf35W+bSJGL3ywijWbq3ju/ZWprooQQggh\nhOjG7p+6iOUbK1myvjzVVREiJkYkYopm9mzEfZRSU4B7/LdZrdbgO8egoaHBkHJSpcXWCoDNZov7\ndXy9ooKSAbkcsE8fQ+rU3ds03Uh7Gkva01jSnkIIIYRItXiD1kalVK6vB3YYgUOH8f3t37Nagjdh\nU0ha6ynAlLa/TSaTx4isWN09u1Z21laggaysrLheR6vdyX/neHNgffLEBYbUqbu3abqR9jSWtKex\npD2FEEIIkWrRzGkNZjZwke/fFwGzOjz+BXAxgFLqcMCqtZZJmSngdsuSekIIIYQQQojuK2JPq1Lq\nCOAJYBTgUEpdDFwGvKaUugbYDvzHt+//gN9qrRcopZYppRYAbuAPXVR/Q7TYnOyuaGS/4X1TXZWQ\nDFgWTwghhBBCCCG6nYhBq9Z6Gd5swR1NDrLvL/3+fVtCNUuiu19cwIbtNTz9l5MZM6wo1dURQggh\nhBBCCOET7/DgHmXD9hoAdlc0prgmoZmiSXclhBBCCCGEED2MEdmDewwPPW8MrkmiXSGE6NGUUo8C\nJ+D9TX8IWAK8AVjwJkG8XGttU0pdBtyId9rOS1rrqSmqshBCCBET6Wn1k87zRtO5bkIIIVJDKXUK\ncJDW+ljgJ8A/gPuA57TWJwCbgCuVUvnA3cDpeKf8/FkpVZyaWgshhBCxkaA1zUlHqRBCiDDmAZf4\n/l0L5OMNSqf7tn2CN1A9Gliita7zLVc3H5iU3KoKIYQQ8ZGg1Y90ZgqRPlwuN4vXltFqd6a6KkKk\nLa21S2vd5PvzKmAGkK+1tvm2VQBD8a6dXun31LbtQgghRNqTOa3+wozBfXPmer76YRcv3XYaFkvy\nYn0ZFiy6m5r6Vsqqmhk3OrGRh5/N38rLH6/h1COH8+dfHW5Q7YTomZRSF+ANWs8AfvR7KNR4najG\n8SilpgD3+G+zWq1x1DBQQ0ODIeUIL2lP4/X0Nq2priaZL6+nt2ey9cb2lKA1Su/M3ghAdb2Ngf1y\nU1wbIdLX7x74ErvTzZv3/oSiguy4y9m0qxaANVuqjKqaED2SUupM4A7gJ1rrOqVUo1Iq1zcMeBhg\n9f03xO9pw4DvI5WttZ4CTGn722QyeUpKShKus9VqxYhyhJe0p/F6epv2Ky5O6uvr6e2ZbL2xPWV4\ncJqTOa2iu7E73QA0tjgSKqdtkIF8BIQITSlVBDwGnKu1rvZtng1c5Pv3RcAsYBFwlFKqr1KqAO98\n1m+TXV8hhBAiHhK0dmOtNidvzFxPZU1LwLYn3l7GZl8vlRDdndy4ESKsXwADgHeVUl8rpb4GHgB+\no5T6FigG/uPrdb0N+BxvUHuv1rouRXUWQgghYiLDg/10t+mjH32zmXdnb2S5ruDJG08C4PNF2/l6\n2S4WrCrlg4fPlV4qEZfq+lZK9zRx4Jj+qa6KECIMrfVLwEtBHpocZN/3gfe7vFJCCCGEwaSn1U93\nS3pU29AKQHl1c/s2u8MV8H8h4vG7B77ktue+o8Z3jqVEN/s8CiGEEEKIriFBawLmLt3Bs++tSLic\nTbtqOe+mj1m8tiym58k1vegqjrZ5qc3xz0v1GHQXyCTjBYQQQgghejUJWgPEdpH91H+X8/n322lu\nTSzhzMfzNgPw8ser43q+LUyvqgS2IhHuFA4/6G4jH4QQQgghRNeQoNVPd71IttldfDZ/a6qrIWLU\n0GzH5U7zky6F1fO0HVw6WoUQQggherW4EjEppa4CLvfbdKTWusDvcQcw3+/x07TWMsmyC308bzPn\nTBrdabtc76enhmY7l941kwPH9OeGn45JdXVCCteLnyxyDgshhBBC9G5xBa1a66nAVACl1EnAzzvs\nUqe1PjmxqgmQ+Xw9VYUvedbaLVVA+gatXy3byQEj+sX13FAjF9xuD3WNNvoV5iRQs8TZHS4yM8yY\nZE0dIYQQQoi0ZsTw4LuB+w0oJ+VSPTw42PE9Miu1R+ou72pbQiYjPf3Ocn597+dstaZuicjGFgcX\n3fYpj725LGV1EEIIIYQQ0UkoaFVKHQXs1Fp3THubo5R6Wyk1Xyn1l0SO0R1E21PjcrnZWd7QxbUR\nIr3NXboTgA3ba8Lv2IWRvbWyEYBvV+zuuoMIIYQQQghDxDU82M/vgNeCbL8ZeBPvZec8pdQ8rfXS\ncAUppaYA9/hvs1qtCVYPGhoaoi6nprYGqzV8k5SXl+FoyQ7YVlpaSk6WJWL5r87cwvfrqrjxYsW4\nkYXt21uaWwBwuVyd6mqz2QGw2+2dHmtqamr/t8vpxGq10tCwNyi2Wq0B67Ua0Z4QW5t2Fzsrmvlw\n3k5+85PR9C3I6vLj7anc+97F0p53Tl2FGl7I5WeMiv2YdTacTjdD+udG/Zympua43+vKigosrtA3\naWpra8OW3dwS+nMRTjTtWVW1t/172rlstJ74eRdCiF6vuwz5EsIn0aD1ZOCPHTdqrf/V9m+l1Bzg\nYCBs0Kq1ngJMafvbZDJ5SkpKEqye94I02nL69e0bcd/Bg4cwsF/gRf/QoUPJzY7clN+vWwLAnkZT\nwHFy80oBaGp1djp+TvY2ALKyshg6dChvzFzPxAOHMHZkMfl5lUAlAJaMDEpKSigs3HsxXlJSQqvN\nGfC3EWJp0+7inte+pKyqmbkra7n+kkO7/HjN7tr2f/fp0yfq9qysXUJlbSW3XnFczMe85omPAfjk\niQuifk5eXl7c7/XAQYMoGdwn5ONFRUVhy87N9QZKGb5zO1rRnJ/+7d/TzmWj9cTPuxBC9HqSzkF0\nM3EPD1ZKlQCNWmt7h+3KNzTYpJTKACYBaxOsZ1Ika05rqFVOWmwuVmysCPm8H3fW8t6cH7nlmW+7\nqGa9V9vcTXe6L0FjAE8MJ3pmRvwzCGI5Tlc8Pxz5rRZCCCGE6D4SmdM6FGiPsJRStymljtVaa2An\nsBjvsjcztNaLE6tmciQrXOl4Me5/Ab10feigNdzyIy02J3f9awGrN+1JtHoiTRkVxMVSzOFjB8V9\nnESz8rZV078Yp8vN3KU7aGpxJFR2MtgdLppb07+eQgghhBDpLu7hwVrrZcBZfn8/7PfvWxOsV48W\nLmiIJTDx37O2wcaKhsrOO0mXUrcWalmWJ99exuFjB3Py4fsEbN9WWo/H42F0SVHIMmMJfc0GLAez\ndksVL01bzV1XHs2AvtHPpw3mk2+38MonaznmoFLu+O3RcZeTjGVu/u+eWbTYnDENxxZCCCGEEJ0l\nOqe1W/p43mbyczI5feIIQ8qL9vLXZPIGrOECU3eq193ppdKx2R1ONxfd9ikH7zuAB6+bFFDHr5bt\n4qtluzoFrX98/CsgwrxVj4dwZ61RS9HY7N6RAfe/soimFgcffr2Jqy88OLAaMdrty/q7cUeEzMMR\nJGNp1ha/+eRCdDWl1EHAx8BTWutnlVKvAUcAVb5dHtNaf6aUugy4EXADL/nWXRdC9DZpeN0jRDi9\nMmj998drADoFrV0duJjwfkeEPYx8iQiftqGlqzcbO+Tb/xRrbnXQandRXJjTvu2GJ7425Dizvt/G\n9cPDJLWK9IHzPRwswEz0s5qMnlYhkkUplQ/8E5jT4aG/aa0/7bDf3cBEwA4sUUp9pLWuTlplhRBC\niDgktE5rb7Rldy1VdS3xPdl3oRyup1ViVhGKUeeG/+n3qztn8Jt7Pzeo5EBtPa3RRJgd52vX1Lcy\nf1XbMit7A8y2YDPRtpCQVfQwNuBsINLaREcDS7TWdVrrFrx5JyZ1deWEEEKIREnQGiDypfDfX13M\nFfd9EVfp0Vwohxoe7PF4JKLthlptTv705Nd888OuVFfFz94TKZnJkkOd/98u383Ft33Ktyt2t297\ncdrqsGUknJRKolbRg2itnb4gtKPrlVJzlVL/U0oNAIbQtk6aVwXepIpCCCFEWuuVw4PTWsdr8TAX\n11FduEugm1LLdAVbdtfx+FvLOMl//mk870sKsgd35bIzbT5bsBWAmQu2ccKhwwBobN67klbASF6D\ngs10ilm3l9Wzo7SBEw4bluqqiJ7lDaBKa71CKXUb3nXQF3TYJ+JHQSk1BbjHf5vVGqlDN7KGhgZD\nyhFe0p7G6+ltWl1TTTJfXk9vz2Trje0pQaufeK/Pu/Sy3r/wdLrS7nG8DW3UXMequhb+8d/lqFH9\nDCnPSEm7j2EKf7zympaoguJgMWtXxtIrNlawYFUppxwxnHGji7vuQD7XP+ZNnjVh/wEUFWR3+fFE\n76C19p/fOh14AXgfb29rm2HA9xHKmYI34AXAZDJ5SkpKEq6f1WrFiHKEl7Sn8Xp6mxb3K07q6+vp\n7ZlsvbE9JWjtJkwmk/SadiNvzdrAih8rWfFjkGWIgLomW8QyuipZULzLKhnto683kZ8T/VdQVV0L\nMxZsA7o2aL3rxYUAzFy4LanL1dgd7qQdS/R8SqkPgFu01luAk4E1wCLg30qpvoAT73zWG1NWSSGE\nECJKErSmgM3u4sm3l3HBifuy7z59Ax4LdS2ejGGaPYnD6aKqrpUh/fNTcvxwSxftKKtnysthOzeC\nMuwMSNKpFE3I/dWynfTtszdz8Z7aFqZ9s5mm1s7Lxfz91cV+f8X+IuoabRTmZ2EymeK+IbCrooGs\nDAuDivPier4QXUEpdQTwBDAKcCilLsabTfgdpVQz0Aj8Vmvd4hsq/DneD9G9Wmtj1rgSQnQrHukJ\nEd2MBK1+4v34RhtQtl0nz/p+O3aHiwWrS3n/oXPDlyVDguNyxwsLWL+tmn/fMZnBMQQYRnVuhjsl\nftxZa8xB4pToz9T2snrWbanirONGh93P7fauq9reFlG07RNvL2PN5qqAbW0B5jZr/d6yO7wIl9vD\n7MU7mDh+cNByd5Y3cN2jczn1yOH8+VeHR65ICNc+MheIsA6uEEmmtV6Gtze1ow+C7Ps+3mHCQggh\nRLch2YP9Jak30+5b3qNtSRBZM9J467d5lx3cXdmY4pp0ZjZH9353PC3SpbP9+se+4vkPVrG9rD7s\nft8s38X/e3gOLbbOvaZ7Bb7ImvrQw6adLr/hsx0a49vlu3j2vRUhe7D19hoA5i7dGbbOiVi3tYqV\nG4MPBw+mudUR//JZaaKsqglXMlNQCyGEEKJXkqA13cn1YFIkMyA0x3GTYntp+AAxFsFGBsxevD3E\nzqHLaQtGW+1O1m6pimsIe1VdC2u37O1ZjTag73ikylpv8LfFWofL1XluaLTlJuLWZ7/jzhc7JmcN\n7df3fh738lnpYM3mPfz+wdk8887yVFdFCCGEED2cBK0GarU7saZhz55IL9EGUB9+tan939c//hXx\n3MGoqG5mZ3lDxP2efmdFzGW3efzNZdz23HcsWVcedj9TkPHBrb7RBuDtWQ4Wz0cT4/vHy6/O3Nrp\n8Y5tHm2AvWbznqj2i4fN77V3R22jGeYu3YnH4+Hz77d3+55jIYQQQqSnXhG0/rChgunfbu60feOO\nmoC/E+1su+WZb7nm4TkxX7iFvYBOsIOot3XU2hwuVvpn7PXAGzPXc+0jc6Iaxlhd38prn66lscVh\neN3Wbqmi1eaMqqfV5XLz/twfEz7mVQ98yXWPzg3YZlSvcturWLS2DICt1vD5XKJJ+hDt6R7uNSzR\n1Z22WcK0uc0ROnj82/Pzo6xR7+P/HixaW8az763g1me/S12FhBBCCNFj9Yqg9Z6XF/LytDWdtt/0\n9DxDj7PNN4SzojpU0Br5kjxd5i12V8++u4I7/xU4RPPd2RvZVdFIY7M94vOXrCvng6828dbM9TEf\n2+328Mm3W6isCf7+3/bcdzz0nyVR9bR25WkQS9nhbqh0xVzsYGUG66GNdShyY0vge+//7LtjGNIb\nr8+/387/vtRR7790fTl/f2URDmf3WAanpr4VgPLq5hTXRAghhBA9Ua8IWtska9mYWNOISyKm2Kzc\nWMlrn64N+n4uWVdmyDFqGiOvo9rR/FVWXpq2mr89/13Ic+0HXYElzvmVhp2+Hg8bd9SwbEN5h80e\n3DEk1alttAX2akewXFfyuwe+jHr/9noF+TzF0hQ19a08/8GqkAWs29q5ZzaS1Zv3sGV39CuFPPve\nCt6atSHsPm1fA3WNNu799/csWlvGDxvCD7lOJf+vLSO/w57/YCUfz+s8MkYIIYQQvVdcS94opU4G\n3gPW+jat1lr/0e/x04EHARcwQ2t9f4L1NITbA5YuiA+NjIVDBTvBtqZrr+zrM9YxbGABpx01okvK\nb0t2c8bRIykZWBDwWPzLFiVYKaC2wRvollc3M350ccj9kpEUKBwPwUcZPPDq4k5Bf7hmuX/qok77\ndhxy7y+aubXRxj7+71djs503wvSMW/c0dX5+dIcJ6XbfsOFolr6J9WbZjU99s/e5sVUrZYy87zZz\nwTYALjhxX+MKFUIIIUS3lsg6rd9orS8O8dgzwJnAbuAbpdQHWut1CRzLGB4P4YboGhUEOhxdP6Rv\nxcYK7M70TOTy3hzvXMxIQavL5cbt8ZCZYYnrOI4gWWI7vofR9nrXduxZjeNc8I9Fwz091sRC0ZQZ\ni1Dnedvc1Hi9NWtDxN7ESIL12AVtC4+HxhYH9/37+7g+B6ECSf9MxkaJ9XtlT+3e4eVG3pjaWd7A\n1Olr+MPFhzKwX27Mz2+1Ofls/lZOnziCooLsgMdktIgQQnQv6drxIUQohg8PVkqNAaq11ju11m5g\nBnCa0ceJR7DP5/PvrzT8OHe+uKB9OZCVGyu544X5NLcmltjHRGAAdteLC/l62a6Eyky1K//+JT+7\n9dP4C4jiCzepX8p+F+6NzaHf72BzNI228sdKlq4PPrQ0WcPk4xFt7ONwurn35YWs31bN5l3RD9MN\np6nFwW3PGZ9IKF1a+8m3l7FsQwWvfro28s5BvDtnI699to5//M+7xI3/aeQIk8xKCCGEECJRifS0\njldKTQeKgXu11m2T1YYA/hPdKoCI47yUUlOAe/y3Wa3WBKrn1dCwd0ii1WolwxIYp89cuK3933V1\ndVEf03+/0tJS8nICm3Ldxu0M7Z/LnS8uAeCjOWsIdvlqtVppbt6bvKS5uTmgbJvN2wNot9up2hNb\nL5DVag1YUsRqtVJdb6MwP7NTO8SioaEBq9XKw2+voyA3g+t/ekDI44dT7UveEu/7XFFRQaYncIkh\ntzuw97W6em+blZWV0ZSXGVXZLS0tMderrm5v8LR+a+j3yr9OAPX19bwyYwsHjS5i4rj+AEET8JQG\nqU+oOnZMRuW/X1lZ9D2q1VXVWK3JCbtsNlvQ4MfhcAR9nRu2hx6K7L9/VVVDp8cqKpo7batpCJ6o\nK9x5EOyxjtv837dwZZWXl2NvzgrYVlNtXPs3t3pf37crdnPmEcVMX7Cbc44pobSqlQn79g0517rt\n877D6p37u81ai9Vqpb5h79rB//podfu/jfjeNrIcIYQQQnR/8QatPwL3Au8CY4CvlFL7aa2DXfVF\n1XeitZ4CTGl/ksnkKSkpibN6e/lf+AwdOjTsUNSiwkJKSkpYs3kPuyubOPOYkSH39a/bkKFDKcgN\nDIYGDhxIyZDC9r/7FhUBu+gYuJaUlJCXVwZ4A5n5a/Zw46VHk5PtfWuys7cBkJmVxYABA4DoM5CW\nlJS09/gCZOf3429PfMHYkf147IYTA/b1eDx8PG8zhx4wiFFDCzsWFcBqtVJSUsLW0iXtxwl1/Gjr\nGY8BAwdSUlIUsM1k+iHg7+Li/nhPVxgyZEinYY2hZOfkRFUvm8PFtK83cfrEEfTtu3eIsdkc+qZA\ncf/+AX/bPdksWl/FovVVXHjawd5tQYK3IUOHdtqWX+gtK9Lr8n8tgwcPCbuvv37F/eJ+f2KVnZ2N\ny9O5hzozMzPmOvjvX2erBvYOWy5vyCAzp7DT/tkhlqoKd+xgj3XcNtTvfQtX1uDBgxnQN3DYbr/i\nYkpKOr/v8cjK1ID3Nd451RtkLl7vDUSvOv8gLjwp+L3Fts97fr73e8pkNlNSUkJhn0a8M0ACGXW+\nJOu8E0IIIUT6iyto1VrvBt7x/blZKVUGDAO2Ala8va1thvm2pVy0oyLb1macPHFEVElzvlm2k+nf\nbgk8VsedTKaohz6+PnM9V1/oDV78nxJrVmIIHApaVuXtXQrWQ7V5Vx1Tp68F1kaVXCZdhWuhm56e\nx/N/PZWszPjm0Abz8TebeXPWBpauL+fUI4cbVm4wrbbOgez/3TMLiC4hUJtoEiKlTNA5rbGf9/OW\n76IoP5tDDhjY6bG7X1oYV9Xi5T/aIV1FWmMXaO+J9cSQYVoIIYQQwghxjRFVSl2mlLrZ9+8hwGB8\nt9y11tuAQqXUKKVUBnAu8IUx1U1MpEuteC/F/vXR6qAZSuO1vbQ+8k4Ga0pwzm0qBE/QE3r/8upm\n1hicaKdtiPPuysYOAVdiF/bBnh1qjmqsbn9hfvT1SHJ8Euy+zvayhpjnhD/25rL2LNPR3izqqmRC\n1z06t0vKTba29nGn8ZxoIYQQQvRM8Q4Png68rZS6AMgCrgUuVUrVaa0/8v39X9++72itNyZe1cQl\nu4egrCrxQNa/xslI4NOdBO2Bi3BB7QwyV7SiujnIntFpe0camh2s2FgR03PCCfbaekOC1lCv8bVP\n40s+XlXXElWDT/tmM7nZxvTAu9weXH6ZrdtubBhp8doyPp63mbOOG8XxhwwzpEyH08XS9RUcMXYQ\nmRnmTkF826iTSF+jNfWt/ON/y/nNOeMZM6wo/M7CMEqpg4CPgae01s8qpYYDbwAWoBS4XGttU0pd\nBtwIuIGXtNZTU1ZpIYQQIkrxDg9uAM4L8/g84Nh4K9VVYg1ZEwpxPfD7B2cnUkKQIqWHw1/QtWsj\nPSdIMDjl38YMF12wqjSq/eINPrdEMYQzlHgzBif7jAt1Y8a6pzHo9kj+/soirvnZhIj7TZ2+Jq7y\n2/jfoLr+sbnsqghfX2eQ5ZqCqawJfkPl/le8a+Su2rSH4utzGD+6f9D9/EU67/735Ubenb2R/NxM\nmlocvHz76Qzpn9/p+e4IUes7szfyg65gR3kDr951RsR6icQppfKBfwJz/DbfBzyntX5PKfUgcKVS\n6nXgbmAiYAeWKKU+0lpXJ73SQojUkktK0c0YvuRNOkvmUh8djxRvoNLdOteSupxKFB2tHevznxnr\nO2Xm3R0hwAilorqZT+dvjeu5EQV5bdPnbem8MUoVNcGTDEWsRjqMDwZW/rgnruJ2lDXg6fplkwNu\nUEUKWAF+e190MyZe/nhNwLqtwVTVRdeTG26khsfjYaNvrntTi3co9rcrApMsWXxfYpHOCZcvqA02\nqkF0GRtwNoH5I07GOyoK4BPgdOBoYInWuk5r3QLMByYlsZ5CiDQRzxrnQqRSrwpaY5bABXtdo63T\ntu4SgJZVNQVdciUayY1ZYz/YzvIG5izZEbgxQvKfUBfps77fFrpuCbaD0c34uwe+jLxTGjD6M+L2\nwJuz1idURrBMzomqDfL9EMrcpTt58cNVIXs4Dbux0KHxOxZrahseHKGnta0YGRmSPFprpy8I9Zev\ntW470SqAoQRfks6Y9NRCiG7lH/9bntD0KCGSLZF1WrudSNd2sVz7RZq/2HGdzFikOs/J7x+czUH7\n9ueh646P+blGVN3hdNHQ7KC4MCf8sYIezBPmL6+G5uDrcQYre/bi7Tz9zgqeuelkbnjia+668mgm\nHuhNjv3+3B8jlhPMW7M2RN4pDl0RXCW9o9XgibtOl5tVm+LrpW1z+wvzefyGEyP2eCYq1Et/Y6Y3\n6D7m4KEcsn/nbMiuaOfqh2naaErYG4xGt2Oi5878VVYmTZBlbwwS6t2P+IHryjXUZS1e40h7Gq8n\ntanb7aGx1Ulhh7XqF67YzFFjI08vMUJPas900Bvbs3cFrTHuP+2bzVx06v5BH7vrxeQsm9FVMUNd\noy1wbc8OB1qzOc4suwZEOX/5xzy2ldbz37+f3Wn92644fLirtqffWQHADU98DXjnEk5//HxKq5rC\nHquxJXS2245LDlmCLKsUT+9Zx+Gchkh60Jrc40VD+96v396f2iTooUY/TJ+3mVOOSHC5pSDvc6f3\nIspgtD24TfDcefg/S7r18ltpoFEplevrgW1bei7YknTfhyukK9dQl7V4jSPtabye1KZ3v7iA5Rsr\nee3uwDwDxf2Kk/Yae1J7poPe2J69anhwpECg43C21z6LL2NpMLFci3fVhbvDb/6C0Umi2hgR42zz\nLflT2xB+rl7b+/ndyt386YmvaW51xHX8WNt75sJtXPPQnIj7JZvTZXyEmewhnl217Eyinn9/pSHl\ndMUc4U3mZIm1AAAgAElEQVS7AhN0xXuMjuv3dhoeHGVf6973MPZ6fPhV4OiFpM+p7llmAxf5/n0R\nMAtYBByllOqrlCrAO5/12xTVTwiRJMs3emcF7K6ML4eHEOmgVwWt/4kjCJ29eEfknZJg085a7ngh\n/iHHENg73GJzJlqloIy8xowUwLQd6pHXl7LFWufN3puEi9xFa8oMLc+oC/NgPbbdickEaw1eR9co\nMxduM6ScRN7qaOL5Lxdt5/ybp8e11nOkhE6maHta4xwevKOsnlc7LG10/s3T+fz77bEV1AsppY5Q\nSn0NXAH8yffve4HfKKW+BYqB//h6XW8DPscb1N6rtY4/LbkQQgiRJL1qePDn32/n+ksOjek5T7+z\nnNMnjjCmAlH2IgVc7Pn+HfW8tRRwuz2YTG1BZvz13LCtmqpY1rTscCibw9X56HFWZ9HaspABlLkL\ng8OK6mYGFefRao99fqrFYny9pKMrsk+/iz6rs/cGRWLvU8ceUX//+mg1AHOW7uTK8w4MeCzRTuz6\nJu9c8EhfRW03m2I9dUKd8y9/vJozjxkZY2m9i9Z6Gd5swR1NDrLv+8D7XV0nIYQQwki9qqcVvENJ\nQ6mqDR4wffT1pq6qTtpIZBjor++dxc3PzAM6X3g2tThoDJL4aNOu2k7bbvnntzz8nyVBj7FwtZUn\n314WkLnUHSSiirTkTSxue+67oNuN7tF84l3d/u+rHvgSl9vDb+79POZywgUz8ZKYNbIXfYFiNNZs\nrmL15vgTQ7lcbq57dG7IxzN8Ny5cQdaBDXfWRvqYrNtaxdylO6Pauf04MX72yquCZ7GUGydCCCGE\n6HVB6yOvLw352IchgtNXPlmb+IFj6OZYvXlPt5rLVddoZ+MObxDacTj1L++cwa/umgnA9rJ6lq4v\nB+DPT30TsVz/FnvwtSV8tWwXG7ZXt2/r2ET/+nBVpzLe/kJ32hb+SJEZ3dO6py5w+ZM7Xpgfcxn/\n/ULz3pz4shmHlcTzsNXW89eMu/PFBdz+fPD3N5qmnvJy2Jw5WMzer/Tp326huTV0MrBOx45we2LF\nxr2rpESdPTjMLq12J/NXWnE43bTavVMVvlvVBYnEhBBCCNEj9LqgNZKp09ck/ZjB4tmVP1Z23tiV\ndTBohUxnkB6eNtc/9hX3/vv7oL1AISrVya3PBu/9DGXL7s7TtTweKN3TxMOvL6GqriXmYZNdOTwY\n4pvX+fbnXbOMTjJHpesdNZF3SjNNYbJExypyVl4TKyJ8L/gPEZ+1sMNc0DAn+o6y4L30bdMS/Ovm\n8XhobHEwf1XwVPsLV5cC0NzqDHiOv6nT1/Lw60t48u1lXPK3z/jPZ+tCvv40zc0lhBBCiCTqVXNa\nozHtm81JP2awi7WGZu/FcLKytxp1nGiyvzo7RELLNpRjrWwy5PjReGPmemYv3kFpVROWOK6I43lO\nd/XGjHWcemSCy6n0YL+8c0aqqxDA/35KsOHzoWwLkrjprc838MbM9Tx89SEB3w9Ol4dfhXndlTWB\n69n+oCu45yVvErjbrziKYw8uYdbCbQB8t9Ib+L4/90eOPXho1PUVQggRO6M6KIRIhR7f07p5t/Hz\n/OIVLNaZuXBb0IvLeDNw+nt39kZe+KDzkNmuFOrr0L/3cEuHJTqmvPw9L02Lfl5gG4/Hw4Zt1ZF3\nDKK0yhsk252umL/Cu7qnNZ3siZBRVhgn4o2jMKfdL++cEXaUQ4SnB9U2f1zvrE9ocvPrM/ZmBH7w\ntSXYHcGHgfeie0FCCJESyV7GTggj9fie1qr6zkmAUiHU9ZhR6z8G88bM9V1WdiihLjz9kxr99dno\nlgWMdEfQgzd5UyLiCUDl4lp0iQSuJZpaHFTXt+L/TdNxSG68563LldhlTse5ym/NCj6UPdTnvRtN\n7xdCCCFEF+nxPa3dNcBI1yEcpXuamPLqatZvDdXDmbx6G5GsasGqUuzOKOfY+vSmnlYIXG5pT21L\nmD1FImIZzhuMCVO3+L4LOXe5G9RdCCGEEKkRd0+rUupR4ARfGQ9prT/0e2wbsBNou8V+mdZaUkPG\nIskXcNFeL//vS01pdSuPv70s6ONGXjS3lbU5yPI4kLoeGHN3iAwMdN/U77nzt0czZ8kOnuvCkQEi\nPdmdbn7clUiSrMAPaqhEY73rUyWEEKknw4VFdxJX0KqUOgU4SGt9rFKqP7Ac+LDDbmdprRsTrWCi\nokkMlAy7KmJrimTU+tvlu1Ej+zGoOM+Q8jweD5kZezvvv/lhV4Llef9/Y6jlcVL0XZsu51Sy/LCh\ngp/d+kmqq9HjRc4eHN7UT7om8/n/5u6IvFMIP+gKJBwVQoj00HEUX7qO6hMimHiHB88DLvH9uxbI\nV0pZjKlSz/TxvM3Y7NGvQ5mMuOjRN5dyzcNzDCvv/JunBywx8/hbwXtjoxXpDmCq7hD2sphVJEmk\nRErT5oXPbD5/pbU9eVJHbrcHR4zD4I3QljU4GnK/Xwghupb0rIruLK6eVq21C2hbo+QqYIZvm79/\nKaVGAd8Bf9Nap+ST0l3ji1a7i53lDV0+BDbShXKsvlwcf69MR5Fe+87y1HTkhwoMhEjEdY/O5ZMn\nLgj5+A8bKiKWUdNgC7r9glumx12vZDFijroQQojQ5GtWdGcJZQ9WSl2AN2g9o8NDdwOzgGpgGnAR\n8H6EsqYA9/hvs1qDL14fi5bW7pk45sm3f0jasaxWK1VVdUG3d9Tc3AyAy+ns8nqVl5cze+HGkI/H\ns0yOEYwMzIXwt3OXcVP/6+vrDfkOTZTD6YhqvwWrSoNu93g8afE6hBCiu5ORYqI7SyQR05nAHcBP\ntNYBEY/W+nW//WYABxMhaNVaTwGmtP1tMpk8JSUl8VavXd7G+Nbx7E1KSkr4z5ed558Ga/+8vDKg\nCktGBtC1ywkNHDiI91/tmnl6QqSjwYOHGFZWnz6FQT/DyVZRE7z3N1oOpzstXocQQvQ0m3bVsnrz\nHq7+6cFkWHr8giKim4s3EVMR8Bhwuta6Oshj7wLnaa3twElECFhF6oXq5Qilorq5i2qylwwXFL1N\nawzz3iPpSZ8fl9uDpZctNZUopdTJwHvAWt+m1cCjwBuABSgFLtdaJ3ZXQQjRbXRMvPTh15sAGD+m\nPycfvk8qqiRE1OLtaf0FMAB4VynVtm0usFpr/ZGvd/V7pVQL3szCqQta5Tqn2+o5l9xCRCfU8k7x\n6ElTr10uNxaz5PqLwzda64vb/lBKvQo8p7V+Tyn1IHAl8ELKaieESKpQiZi++H47Jx02rNetjiC6\nl3gTMb0EvBTm8aeBp+OtlEhftSESvXSJHnTRLUQ0DD3lfT2t20vrjSxVdG8nA//P9+9PgJuRoFWI\nXm/15j0sWV/OxPHGTVERwmgJJWLqDmQNKmNV1iYvsVV5EoYgC5FWDIxaXb6g9dn3VhhXaIrI/au4\njVdKTQeKgXuBfL/hwBXA0JTVTAiRViprumfiUtF79Pig9aVPw69tKNLX/a8sSnUVhEgqI9fQa5vS\napHkGr3Vj3gD1XeBMcBXBP7mR3VHt6sy+zc0NEhWaANJexqvJ7ZpVVVVyMfq6uq69PX2xPZMpd7Y\nnj0+aJU1NY0m7SlEVzEyd1JPSsTUk15LsmitdwPv+P7crJQqA45SSuVqrVuAYUDEK56uyuxvtVol\nK7SBpD2N1xPbdED/AYAO+lhRUVGXvt6e2J6p1BvbU27BixjJcGshugOJ83o3pdRlSqmbff8eAgwG\nXsW7bjq+/89KUfWEEClg5GgeIZKtR/e0uqSX1XCSWE6IrmNkj2JbWT3iMytf5fGYDrytlLoAyAKu\nxZvN/3Wl1DXAduA/KayfECKdyJ1OkeZ6dND67Yrdqa6CEEJEzchLBrln17tprRuA84I8NDnZdRFC\nCCES1aOHB1fXSSY0I+2pbaEqidmDhRDx60nzQHvOKxHCy+5wpboKohcKu6JGjxiWI3qyHh20Sk+D\nsX57/xc0tTpTXQ0hei5DEzEZV5YQwjg1Da1cdNunPPPO8lRXRfQCi9aUtv877JxW+dEQaa5HB609\nqadBiO7ilv87ItVV6LbWbQ29HEGs3G1zWntA8jT5Lhc9yZbddQB8uXhHimsieoO/v7o41VUQwhA9\nOmiV5W4S8/wHK1myrizV1RDdzImH7ZPqKnRbRs7D98j3nxBCCD+7K5tSXQUh4tajg9YZC7alugrd\nQnOrI+j2mQu2cd/URRK4djODivPifu4BI/rGtP+9Vx/Lf+45M+zjHzx8LtMfPz/qMseO7MdV5x8Y\nUz0Abvj5oXzyxAWYu3HHYkWNcXPGV26qNKwsIYQQ3d/z769MdRWEiFuPDlqr61tTXYVu4YtF28M+\nft/URazevCdJtRGJ6pOXGdV+g4vz+PCR89oD1f5FOdx91TFRPfe2Xx/Fw384nsPVIIoLc/jVGQqA\nyRNHAPDBw+cy7bHzOVwNIivTginKBA998rJ47IYTufCk/aLa/8+/Orz930P653v/IckkANhZ3pjq\nKhhGRgcLIUTXkq9Zke56dNAqojN1+tqI++yq6DkXwD3d5Ikjo9qvuDCHzAwzI4cUAmAymdrnQYYz\npqSISYeUcOCY/u3bLj1zLJ88cQE3/OIwALIyLVg6dHk+e8sp/Pz0A7j398e2bzv5iMChxHddeXTY\nYz//11MD/vZ4PFxxzngKcjPZb7g3+O7OPa1Gs1Y2SgwvRJqRmzDCaE6Xm8Zme6qrIUSXkqBVRKVs\nj8yD6A5GDOnDcROGRtwvPzeTP/3SG2BeeuZYDhzTnzuumEhBbvhe2tennMmTfz4prrqNHFLI5WeN\n4/Cxg5g8cQT7D+/LTZcGJm3yz2yYm20JeOy/95/F8MF9eO6WU9q3ZWVauOjU/fnv388mN7tHLzsd\nl2sensOqTdGNkrj41P057ICBXVyj+Mg1vhBChPbnp77hV3fNxOGUpZREzxX3VZ5S6ingGLzXE3/S\nWi/xe+x04EHABczQWt+faEXjMXJIH7aXNaTi0D3Oh19vSnUVutSHj5zL7som/vj4V4aXPeX3x/D1\nsl0sWGVl2KACtlrrEyqvX59sahpsQR8zm0z065PDh4+cy+zFO3j+g1Xtj1170QTGjSpmxJDCgF7Q\nAX1zefgPx7f/Pe3R88BkYs6SHUwcP4S+fbI576aPfcfOSajubdp6ZDvKztwbqD5786ms3VrFk2//\nAEBBXhYAI3w9wwA5WYGBLUS/1JUa0Q+X201mhoVHrj8et9vD/FVWnnz7B1y9MInRUeMHc/lZ47jg\nlumprooQ3VJzq4Nbn/2OEw8bxiWnHRDVcxqb7e3fbULEa1up97qixeYiM6Pz76IQPUFcPa1KqZOA\n/bXWxwJXAc902OUZ4CJgEnCGUmp8QrWM033XHMe5x5YAcOYxI9vn3XX05I0n8uB1kxg3qtjQ4w8q\nzuNwNcjQMts8eN0khg8uAOCBa4/j9xcexOv3nMnPTo5uLmA6euDa4xg1tDDiflPvnBx1mT8/3Xvh\n8HqYZEEAZrOZUUMLueS0/Rk+uIBD9h8Qdv+rzj+Qv19zHGceE3oo7jM3ncy0x87niLGDuemyI3j/\n4XN56saTGDawIORzBvTNZfrj53PwvoHHLy7M9h33IF6f8pNOw2Tb/NT3/mdmWBjrdz6/MeUnnH3c\naEaXFHUattuRxWLGYjZxxtEj6dtn73Fv+81RYZ+XqFFDCxkzrKj970HFeZxyxHCe/+up/PuOwPf8\nT784lOGD+3RqJ4DrLtiPQf1yedavR7ajMSVF3Pabo3jqzyfz6B9PwGQyYbGYOfGwffjgkfMYMaRP\n+77THz+f1+85k0euPz5gaHMknzxxQcCc22Duuupo7rs6+jK70uDiPMxmE2/dd1ZcibDa3H1V+CHe\noTx+wwm8dd9ZXHrm2M4PynjKXqW7Zv7/8OtNbCut5/UZ66N+TqgbkELEQ6aDiJ4s3p7W04BpAFrr\n9UqpfkqpQq11vVJqDFCttd4JoJSa4dt/nSE1jkFxYQ7nHTeMay7ee7HddkH0zDvL29dI2394PwAe\n/eMJ7T1K4bw+5Ux+PeXzsPtccOK+XHX+gZhMpvYy33vwHC65/bOI5e8/vC8mEzz8hxOob7Jx3aNz\naW51AvDbcw/kZ6d4A5Pn/3pa+3Mm7Ocd1nf52eOYdEgJ+w4rorHFQVFBNjX1rfQrzOH/7plJXWP6\nznmYsN9A/nnzKRHfg4F9c6Mu8/KzxnH5WeMAuPf3x3LPywuD7tcWx/367PH8+uzxuNwebHYni9eV\ns3C1lWt/dgiXT5nVvv/pR42gIC+LQw4YyPWXHMrvHviS8upmAP5y6eGccsTwTsfwBkYmCsIkStpT\n24LJZOIvlx7O428t49dnj2P91mpOOXI4/fpktyc0Gj64D5eeOZa3P98AeAMkhzPwDuvokiL+cPEh\nHDimf3vwGa8LT9o3oedH4+xJo4MmbBo+uE+nbadPHMnpIebuHjymL2ceH/o+2fN/PTVomW0sZhP3\nXHUMb8xazxXnjMdkMtGvMId+hTm43R7OO2EMsxdv57wT9uXd2RvDvqZMy977gsFGfkwcPyTs87tS\nVoaZf958CnWNdvrkZ9K/yPu5KszP4sKT9mtPhrV7926GDRvGwtVWHnytfUAN1140gRf8evLBm6Dr\nqPFD+OSJC2hotlOQm4nHA2azCZfLjcVixuPxMHPhNqZ9s5nJE0ewatMe/nDxIe2JtH51huLnpx9A\neXUTL364mh90RXIaRKSFTTtr+fM/vuHGXx7GaUeNSHV1YtJqi25opv+NQ1mDWMSiqq6FxhYH/Yty\naW5xdFotoLbBRh+/nvsWmzPqsuVUFOku3qB1CLDM7+9K37Z63//911qoALr+ijdGvzpjbNCFvd9/\n+Fwuvu3TkM8bO7If/frkkJ+bSVOLg4f/cDyrN+/hrVkbAvb7zTnjOl2Am6PMEHPBifty0uHeBDX9\ni3J5/q+ncsV9XwDQrzB88JFhMXPACG8QXlSQ7XuOd0jneceP4c0O9ewKpx45nLlLd3ZZ+eEy0T57\n8ync+NQ3OF3uTo8dPjZ0r3fHMi1mE3k5mZx8+D6cfHjndUc7dgRM+f0xXPvIXCD4kNWAY/n9+6Hr\nJnH7C/M7/Vj4D9kdP7o/wZx02DDe/nwD119yKEDQIUE/OXZU2Lqkg5suO4LXPl3LcQdHnoubiBf/\ndholA0L3cvsbVJzXab4teD/DV194MFdfeDAejweTCb5ctJ3qem9vye8vPIipH6/hjz8/rH3/Ntf8\ndAK3vzDfgFdijLt/dwwlAwsoiTCNte2zcezBJeRmW2ixuTj/xDGcfdxo9h/el3e+3Mi+w4q4+LQD\nyMzYG6S3XTi1fbQsvgDeZDJx9nGjOfu40QBBh1FazCZKBhSQHeGzJHqeOUu9v8svT1vd7YLWDMve\nz3urzYnFYsJmd3Ua/pvhdzNL4gQRi7ZrwTY3/vIwPp63uf3v6x6dy3sPncPcpTuZu3Rne4eH6Oy/\nX2hK9zRyyWkHhL2RLdKHUZlLwkVjUUVqSqkpwD3+26xWawJV8mpoaAhaTq1fj2Oo45hNnYOTP1+8\nH1arlUeunoDD6SYvx8aJB/bhrb2dcEw+YjCVFeWdyisvK+Wso4cyc1Fp2DrnmFs71alfQSY1jQ5a\nGuvjbpdjx+bz5qzw+2Rnmrls8ihembElrmMA5GbElwig7XXdf+XB3PXK6rD7vXjTUVzzxJJOj2V6\nGnn82kNZvL6KooLMqNrqmvP2jblNG2oraazbe2r7j7Ovrq4hXHHDBmSxYTucc0wJxbl2Hvv9WG5+\nae/NhFjq8q+/HInJZDLks5IqBwwx8+DvDqapvoqmxKb7AqE/89jrsSY4n7ijUycUsW5zTnvQeuS+\nuRz5F+/IDqvVSk1tdfu+mXROZhbqfbv6vH156RPvhcioIflsK4s+EVpbYNnRpaeN5O05e5e3srfU\nYbUGX6PZn397XnzicN74chsTRuZitVrJt8CVP/GOKqisMH4954kH9KFPjof62j3U1xpevEgTZVVN\nPPTaEv5wySFk+W6+OZydbzy28d4wSr9xkBl+N20uuf0zBvXLpaKmhemPnx9Q39lL9t4wdzhCv85k\nStc2TbamFge3Pvst28sauPS0kfyqpKTLjzl36Q4sZnN7Z4W/Pz3xNaVVjbz74Lls2Fbd6fF//G95\np22X/C3yiL5gXpq2muMmDG0fdWMkt9vD/DWVnNVvYMiEjy63h2Xry5m5cBu/mHwAY0caO2Wvo7aR\nal8t28UnT1zQpccSxog3aLXi7VFtUwKUhnhsmG9bWFrrKcCUtr9NJpOnxIAvC6vVSrBySoBLz7Rx\n0Jj+lJQEzos767hRDCnO54ITx3DhXz9p356fkxG0LIDxo4tZt9X7hXLDpYFrXR594BAWrS1j2LBh\nXPfzfbjsbBvfLN/F5l11nXokH7puEgcFmad351W5fLtiN5MnjQvozYiVGrEFvaMm5OMWi5kLTjmI\nmYvKKa2KPWPwkeMGUzK4mCje8k7a2rZvsQMIHbSGOy/aHtt3dOfhuUD7RYS/U49R5OVEXtv0olP2\n44OvNvHMTSczrKQo5H7FxcWUlITuNbz2ksFMPKico8YPITPD3ClwMeK8781Cfua7qF3POMbF6i0/\nBD3G9ioT4A0+x+43ktfuHsQX32/n/a82cd1FE9r3P2fSaD6bv7X9eeedfFB70PrPW07npqe/ITPD\nwtotVQHlv3z76VTWtnD783t7cN954Bwee3MZ367YHbDvz8+cwEFqH775YRfjRxdz1IToerH82/Pn\nJSVcdMaEiPOijVJSUsLpxyXlUL1KuESKqfDmzA1ssdZx09PzmDTBe67ZnW6WbSjniLGDAXA4Xaz8\ncQ/NrQ4ee3MZv7/gIM4/Mf5BXAtWWSkuymHfYX1ZuNpKY4uDI8cO7jTcMphbnpnHhu01HDV+cMDa\n1h1/m9t+a6rqWsnPzaTF5qSh2R7wu//Mu8t57IYTA5LQJdsv75zB6JJCHrru+IDt9U12zKa9SfAc\nTpfvZr3399LlcuN0e9heWk9lTQuTDon8Hdtqc1Ld0EphXhatdhcDIkz3qWu0sWZLFRPHD2bL7jpU\nkECmsdlOWXUz++3TN2gZ9U12Zi/ezobtNYwfXcy+w/piNpv4atlORgzpw+SJI9uzzy9cXdo+jePt\nOduZdPgYZi7cxqlHDm+fSuavpqGVvJxMsjMtNDbbabW7sJhN2BwuCvOzIl5bNDTbeeq/3sDz8beW\ncdaxozho3/44XR4WrS1li7UOgItu+xS7o+szA19x3xcce/BQbvi5dwRXTnYGm3bVkp1p4an//sCe\n2hZ+MVkxYb8BmEwmKqqbOWLc4IDfhJqGVjIzLCxcZeWEQ4exo7yBZevLefuLbaza2sR9Vx/L7spG\nhg0swGQy4fF4ePtzzf++1O1lLF1fztQ7JvOXp7+h1e5i0oQSJh44hEP3H0h+h6B3q7WOphYHo4YW\nBoxqqG+yY3e4KMjNJMf3/rbanLwxaz16W+B1sN3hIiuFn8FgPB5P+/SaUOwOF+u3VdPY4uCYA4e0\nj2iKVovNyfbS+oAcKOnMFM98CqXUccC9WuvJSqnDgWe01sf7Pb4WOAfYBSwELtNah5/81bFiJpPH\niLkeoS5go7VoTSmDivPIz/V+KbUNue3I5XIzZ+lOGpsd7XNO/bndnqAn3mufruWDr7yZeafeOZlB\n/SL/YCbij49/1Z5lLpjjJgzlb7+ZyAsfrGTGgm0xlf3y7afTrzCHuUt38vz7K2OuW9udLo/Hwz/f\nXRF0+Lb/Hetgc18j3S1ranFQVdfCiCGF3Pvv71m6vpxpj50f9UV4uLvR59/8MR4PPHjtJA7eL3wi\nJ39WqzWg11ju+CXG/zP/xaLt/PPdFdx06eGcHGSesVHem7ORI8cNZnSHmxl1jTb+7x7v8IZI7+u2\n0nrenLkeNbIfl5x2QPv5/ckTF7Sfd2/MXI/L5UaNLMbucLXfmV+yrowvF+/gcDWofUh4bYON3ZWN\ntNicFBfmBCS5ikWi36HpxHeB1Ku7k3yJFG/RWp+rlBoHvOJLqhgVo36bb/3nV6zbVs8D1x7He7N/\nZMWPlUH3m3RICfNXBr8J+uLfTmNnWQNL1pdzxtEjaW51cNeLC7nuogkcf+gwMjPM7K5oZGC/PCxm\nE4vWlrHfPkUMHVDAz279JGiZ7z98LovXlLHf8L58sWg7FTXNzFu+m/NOGMPVFx5MeXUzv3vgy/b9\nJ08cwWU/GcubMzcE9KDGa799ijhuQgkjhxRyyAEDWbelCqfLzRFjB+Nye3B7PKzQFcxZupOFq719\nBbnZFg7drx/ZOTkcOLo/Jxw6jPLqZtZs3oPFYuYwNZBVP+4hPyeTI8cPJjc7A2tlI+/P/ZE+eVnt\nqwNMvWMy28vquW/qIkaXFLZnuj/liH3YtKuWneXe9dr/evmRVNW1dFrjvbgwm5svO5KSgflY9zQx\n7evNrNpUicVipl+fbHZVNJJhMQedvnPG0SO56JT9yMvJZNHaMsYMK+T5D1axaWfgEIvzTxzDxafu\nT3amBYfTzXPvr2xvB/AmlLvktP2xOVy88+VG6puMy+ORlWHGHmIEQKjXlZlh9ga0Ld5RLRPHD2Ho\ngHyysyxUVDfz9Q+7DKtfdzRpQgnzVxk3Uuzmy45g3Khinv9gJcs2BOZDMJnCz9t97I8nkJ+byfKN\nFVTWtNC/KBe9vZrvVlp55c4zcLndbNhWzdylO1m+sZKJ44dwxbnjcbrcrN1SRYvNyUWn7M+aLXt4\nedoafn7aAfTJz+SuFxcyoG8uxx08lLVbqxg3qphD9h/IqKGFvs9pFWXVTcz7YRejSoq47MyxVNW1\n8PwHq8jJsvDeQ+cC3tEnH3z1Ixu2VbNsQwVF+ZnUNQWOlvrwkfOYPm8zP+gKbvvNUdQ22HC63GRm\nmOmTl0XpniYqapo5cEx/Zi7cxjtfekOzR68/gXGji9lV0UBRQTZut4esTAvZmRbMZhMOp5ubn57H\nFudk89cAAAyuSURBVGsdD103if1H9DP0Rlu0v81xBa0ASqmHgRMBN/AH4DCgTmv9kVLqROAR364f\naK0fj7X8dAlak2H24u0sXV/Brb8+ssuH57w3Z2N7ZsP7rj6WgrxM/vKPeYwZVsTZx43mhENLyMvJ\nDBm05mZnhJzY33ZR3mp3csPjX8fUU5ubncG7D54TsC1SUNr2+LhRxfxyssLudHHMQdHPi/R4PLjd\nnpjvTIWys7yBpevLufCkfWN6H61WK3e/urY9kZMErYnp+JlvbnVE1ZPeVZZtKGdI//ywWaODWbKu\nDJPJxJHjBndRzaLTHb5DoyVBKyil7gN2aK3/7ft7AzBRax3V2Hkjfpunf7+GV+bNwWSSGZ1CCNEd\nDSscwrPX/NKQsqL9bY57TqvW+rYOm1b6PTYPSI91HLqBcJlQjfazU/bn+EOG0b8op30oRLAg6dIz\nx2KtbGJbaT21jd75euceP5pfnTGWHTt3U9RvAKVVTRwwvB8PvrY4oHc5JyuDf912Gq/PWMdB+w7g\nsAMG8v7cHzly3GAqa1sY2j+fgf1y+cUdMzjmoCGcf+K+7BPkgr5vn2xqG2w8cO1xOF0e+hcGrhF6\n6ZljsdmdXHFufMtztGXzNcrwwX3insz/wq2n8fK01fxicnRr+4nopTJgBdqHN8bqqBRmFhY9WrhE\nikmxuuYHskZ2fVJAIYQQXaPcuRmX+xIs5uQNq467p7Wr9aae1nRW12hj6vQ1/HKyosQXWBrVpk6X\nG4vZFLJX0u32Zmft6ckh5Bw1lrSnsXpSe0pPKyilXgI+01p/7Pv7O+DKUFN4OiZJ3LhxI7t37w62\na9RsLjtryjaQk5ODB2hqdVKQk4HD5cFiNtFxtoYH7++B/zQOt9sDZhMO31y02kYHKzbVMG5kIQOL\nsmmxu3A4PWRmmNha2sTAomz6FmTi9kCr3UWL3UVBTgbVDXayM80U5mVS12SnrslJYX4GdY1OzCZv\ncsK+fbLIz7Gw2drEjvImxo8qYtiAXOxON1tLmyirbmXk4DxGDsn3/qYBNocLs9nUvuxVW1LHVrsL\nE94hvQBOtweH040J79w1j8eD0w0ZZrA53Hg80NzqpMXuIjvLgt3hpm9BJnnZFpptbnaUN1FVb2Ng\nYQYtDuhbkMWuimbcwIhB3mHRFrMJD1CQk4HN4cLp8o4yysvNwGI2kWE2kZlhpqnVhQcPedkW3G6w\nWKC+2UmfvAzysjOoqbeTnWWmstZGfk4GBbkZtNhd2Bwu3G7va9pa2sjIwQU4XC4sJhMNLU42Wxsx\nm0yMGprPwKJsMjJMWMxm3G637/8eb8pO3/w9m8MFHu/Q2pxsC07fkFzvtYCHmkaH9/oAb/Iri9lE\nbYMdu8tNn9wMMJmw2104XB6abU4sZrP3Pc7PJCfTgsvtITfbggdwOt2YzHuXJ3M43WRnev/d0tJK\ndnYOrQ4XjS1OMswm+vbJaqsqHs/eeYeYwGzytnOr3RUwDrWp1UVuloUMi4mcLAuNLU5yss3YHW6c\nLsjKMGH2vU8ZGSaaWl24XB5ys83kZmfgcLjJyDDhdnuvm5pbXdidLgrzMzHhfZ7Z7M1J4nC6cbrc\nmDBhMoHb48Hp9OB0ebA5XJRWtWKxmCjKz6RvQSZ98jLxeDy43N7PlMeXz9piNmMxQ6vdTYvdO5c2\nN8tCXo4FjweaWpzUNztwe7z1d7qgINdCn7xMbA7vZ897fYf3fPN48LgcFBbk4fZ437uaBjubrY3k\nZWcwfHAe+TkZOJxu72sxm7HZXWT5znm320OGxYTb433PXB4PTpeb/JwMzCYTTpeHxhYndqebDDMU\n5Gbi8nhosbnIyTSTlWWhrsHOrj0tmEywb0kB9c0OqursDO2fg9Pl4cddDeRkWciwmBlQlE3fPpnk\nZWdgNkGzzYXL7abV7qaqzobbAzlZZrIyLJhMYHe6KK+2Udtop39hFiOHFJCdaaa8upVhA3PZUdFM\neXUrY0oKKO6TSYvNTUVtK3vqbPQtyKIwP4OivExqGhxsr2hizNAC9hmUx+7KZsqqW9uHno8pKWDY\ngFxcLmhsbiE/L7f981HX5KCq3kZRvnepuRa793PU2OrNZN7S6sTtgfzcDBqaHBQXZuFyQ262mV0V\nLdQ22ikuzMLtxtfW3nOhpsGBCVAjChnl+44zm2FwXjFj+huz6sOwYcO6tqdV9A5FBdn8JcjSH0bI\niDAsN9olgoQQQkQtXCLFTroqSWK2JcvwmyE/7+rxXYd3cfkJ6Ek3l9KFtKmxpD2N1Rvb05jJfEII\nIYToDr4ALgbwJVK0aq0bUlslIYQQIjwJWoUQQoheQmu9AFimlFoAPIM3kaIQQgiR1mR4sBBCCNGL\nBEmkKIQQQqS1tA5ae3oCHiGEEKK7kd9mIYQQyZa22YONopTyaK3lF9ZA0qbGkvY0lrSnsaQ9RVeQ\n88pY0p7GkzY1lrSnsXpje8qcViGEEEIIIYQQaUuCViGEEEIIIYQQaas3BK33proCPZC0qbGkPY0l\n7WksaU/RFeS8Mpa0p/GkTY0l7WmsXteePX5OqxBCCCGEEEKI7qs39LQKIYQQQgghhOimJGgVQggh\nhBBCCJG2JGgVQgghhBBCCJG2JGgVQgghhBBCCJG2JGgVQgghhBBCCJG2MlJdga6klHoKOAbwAH/S\nWi9JcZXSmlLqUeAEvOfFQ8AS4A3AApQCl2utbUqpy4AbATfwktZ6qlIqE3gNGAm4gN9qrbck/1Wk\nF6VULrAGuB+Yg7Rn3Hzt9FfACdwNrELaMy5KqQLgdaAfkI03dX4Z8ALe78tVWutrffveAlzi236v\n1nqGUqoIeBsoAhqBS7XW1Ul/IaJbkt/m2Mhvs/Hkt9k48ttsHPltDq/H9rQqpU4C9tdaHwtcBTyT\n4iqlNaXUKcBBvvb6CfAP4D7gOa31CcAm4EqlVD7eL6XTgZOBPyulioFLgVqt9fHAA3h/WAXcCbR9\nYUh7xkkp1R+4BzgeOBe4AGnPRFwBaK31KcDFwNN4P/N/0lpPAoqUUmcppUYDv2Rvuz+plLLgvfD4\n2teeHwK3puA1iG5IfptjI7/NXUZ+mw0gv82GuwL5bQ6pxwatwGnANACt9Xqgn1KqMLVVSmvz8N6x\nAagF8vF+sUz3bfsE75fN0cASrXWd1roFmA9MwtveH/n2ne3b1qsppcYC44HPfJtORtozXqcDs7XW\nDVrrUq311Uh7JmIP0N/37354L95G+/V4tbXnKcBMrbVd///27h9EijMO4/j3LMSYFCkEo4KxSPj1\ngiAYLLQQURFUSBFJIoopTBrRSgULC0GsrARREAPWFkHFJIUQtUyT8BPFfyDEK5JwSaFRz+Kdg7k7\nb7O7N3c3t34/1e07czD3sNwz7+zMu5nDwCPKe7qe59i+Ujfs5t7YzQ2zmxtlNzfLbu5gkCetHwHD\ntdfD1ZjeIjNfZea/1cu9wA/A+5n5vBp7Bixjcq6TxjPzNTAaEQtn49hb7DRwsPbaPPu3ClgcEVci\n4mZEbMQ8+5aZl4GVEXGPclJ8CPiztkvXedbGpG7YzT2wm2eE3dycVdjNjbGbOxvkSetEQ3N9APNB\nRGynFOO3EzZNlV+v4++EiPgSuJWZD6bYxTx7M0S5+riDcvvMBcZnYp49iIjdwOPM/ATYAFyasEsv\nub3TWWrafP90wW5uht3cOLu5QXZzZ4M8aX3K+Ku3yykPhGsKEbEJOAJszsy/gX+qxQoAVlAynZjr\npPHqwfqhzHwxW8feQluA7RFxG9gHHMM8p+MP4JfMfJmZ94ERYMQ8+7YOuAaQmb8C7wFLatu7zrM2\nJnXDbu6R3dwou7lZdnOz7OYOBnnSep3yEDMRsRp4mpkjc3tI7VWtOHYK2FpbaewGsLP6eSdwFbgD\nrImID6tVztYBNyl5jz13sw34ebaOvY0y8/PMXJOZa4FzlBUKzbN/14ENEbGgWvjhA8xzOu5RnjEi\nIj6mnGj8HhGfVdt3UPL8CdgSEQsjYjmlBH9jfJ5j2UvdsJt7YDc3y25unN3cLLu5g6HR0dG5PoYZ\nExEngfWU5bUPVFct9BYRsR84DtytDX9F+ae+iPKQ957M/C8idgGHKctsn8nM76tVy84BnwLPga8z\n88ks/gmtFRHHgYeUq2cXMc++RMQ3lNvjAE5QvvbBPPtQnTScB5ZSvkbjGGVZ/bOUi5l3MvNgte93\nwBeUPI9m5o/V71+i3Bb2F7C7+gRI+l92c/fs5pljNzfDbm6O3dzZQE9aJUmSJEnz2yDfHixJkiRJ\nmuectEqSJEmSWstJqyRJkiSptZy0SpIkSZJay0mrJEmSJKm1nLRKkiRJklrLSaskSZIkqbWctEqS\nJEmSWusNGHH2sicgELoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f21dde024d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# simulate a 10 kWh battery\n",
    "fig = plt.figure(figsize=(16,2))\n",
    "ax1=fig.add_subplot(121)\n",
    "plt.plot(np.arange(0,len(price)), load)\n",
    "ax2 = fig.add_subplot(122)\n",
    "plt.plot(np.arange(0,len(price)), price)\n",
    "if billingOption[0] == 1:\n",
    "    plt.plot(np.arange(0,len(price)), np.ones(len(price))*billingOption[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# how big is the battery\n",
    "noUnitsSim = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if billingOption[0]==1:\n",
    "    \n",
    "    modelPrice1 = np.zeros(len(price))\n",
    "    for j in range(len(price)):\n",
    "        modelPrice1[j] = price[j]+(8760-j)/1000000.\n",
    "    modelPrice2 = np.zeros(len(price))\n",
    "    for j in range(len(price)):\n",
    "        modelPrice2[j] = price[j]+(j)/1000000.\n",
    "    priceDict1 = dict(enumerate(modelPrice1))\n",
    "    priceDict2 = dict(enumerate(modelPrice2))\n",
    "\n",
    "    # split load into +ve and -ve\n",
    "    posLoad = np.copy(load)\n",
    "    negLoad = np.copy(load)\n",
    "    for j,e in enumerate(load):\n",
    "        if e>=0:\n",
    "            negLoad[j]=0\n",
    "        else:\n",
    "            posLoad[j]=0\n",
    "    # print posLoad\n",
    "    # print negLoad\n",
    "    posLoadDict = dict(enumerate(posLoad))\n",
    "    negLoadDict = dict(enumerate(negLoad))\n",
    "\n",
    "    # want the MILP model to charge early when PV>load (this is likely reality), so provide a very small incentive with the FIT\n",
    "    FITseries = np.ones(len(price))*billingOption[1]\n",
    "    for j in range(len(price)):\n",
    "        FITseries[j] = FITseries[j] + j/1000000.\n",
    "    FITDict = dict(enumerate(FITseries))\n",
    "    \n",
    "    # model\n",
    "    m = en.ConcreteModel()\n",
    "\n",
    "    # set\n",
    "    m.Time = en.RangeSet(0, len(price)-1)\n",
    "\n",
    "    # variables\n",
    "    m.SOC = en.Var(m.Time, bounds=(0,noUnitsSim*storageProperties[0]), initialize=0)\n",
    "    m.posDeltaSOCGrid = en.Var(m.Time, bounds=(0,noUnitsSim*storageProperties[1]), initialize=0)\n",
    "    m.posDeltaSOCPV = en.Var(m.Time, bounds=(0,noUnitsSim*storageProperties[1]), initialize=0)\n",
    "    m.negDeltaSOC = en.Var(m.Time, bounds=(noUnitsSim*storageProperties[2],0), initialize=0)\n",
    "    m.posEInGrid = en.Var(m.Time, initialize=0.)\n",
    "    m.posEInPV = en.Var(m.Time, initialize=0)\n",
    "    m.negEOut = en.Var(m.Time, initialize=0)\n",
    "    m.posNetLoad = en.Var(m.Time, initialize=posLoadDict)\n",
    "    m.negNetLoad = en.Var(m.Time, initialize=negLoadDict)\n",
    "\n",
    "    # m.helperVar = en.Var(m.Time, initialize=10.)\n",
    "    m.Bool_char=en.Var(m.Time,within=en.Boolean)\n",
    "    m.Bool_dis=en.Var(m.Time,within=en.Boolean,initialize=0)\n",
    "\n",
    "    # parameters\n",
    "    m.price1 = en.Param(m.Time, initialize=priceDict1)\n",
    "    m.price2 = en.Param(m.Time, initialize=priceDict2)\n",
    "    m.posLoad = en.Param(m.Time, initialize=posLoadDict)\n",
    "    m.negLoad = en.Param(m.Time, initialize=negLoadDict)\n",
    "    m.FIT = en.Param(m.Time, initialize=FITDict)\n",
    "    m.etaChg = en.Param(initialize = etaChg)\n",
    "    m.etaDisChg = en.Param(initialize = etaDisChg)\n",
    "    m.ChargingLimit = en.Param(initialize=5.)\n",
    "\n",
    "    # objective function\n",
    "    def Obj_fn(m):\n",
    "        #return sum((m.price1[i]*m.posNetLoad[i]) + (m.FIT[i]*m.negNetLoad[i]) for i in m.Time)\n",
    "        return sum((m.posEInGrid[i]*m.price2[i]+m.negEOut[i]*m.price1[i]) for i in m.Time) + \\\n",
    "                sum((m.FIT[i]*m.posEInPV[i]) for i in m.Time)\n",
    "    # ...\n",
    "    m.total_cost = en.Objective(rule=Obj_fn,sense=en.minimize)\n",
    "\n",
    "    # constraints\n",
    "    def SOC_rule(m,t):\n",
    "        if t==0:\n",
    "            return (m.SOC[t] == m.posDeltaSOCGrid[t]+m.posDeltaSOCPV[t]+m.negDeltaSOC[t])\n",
    "        else:\n",
    "            return (m.SOC[t] == m.SOC[t-1]+m.posDeltaSOCGrid[t]+m.posDeltaSOCPV[t]+m.negDeltaSOC[t])   \n",
    "    m.Batt_SOC = en.Constraint(m.Time,rule=SOC_rule)\n",
    "    \n",
    "    # boolean constraints\n",
    "    # not actually needed - no need to formulate as MILP, LP is sufficient\n",
    "    def Bool_char_rule_1(m,i):\n",
    "        bigM=500\n",
    "        return((m.posDeltaSOCGrid[i]+m.posDeltaSOCPV[i])>=-bigM*(m.Bool_char[i]))\n",
    "    m.Batt_ch1=en.Constraint(m.Time,rule=Bool_char_rule_1)\n",
    "    def Bool_char_rule_2(m,i):\n",
    "        bigM=500\n",
    "        return((m.posDeltaSOCGrid[i]+m.posDeltaSOCPV[i])<=0+bigM*(1-m.Bool_dis[i]))\n",
    "    m.Batt_ch2=en.Constraint(m.Time,rule=Bool_char_rule_2)\n",
    "    def Bool_char_rule_3(m,i):\n",
    "        bigM=500\n",
    "        return((m.negDeltaSOC[i])<=bigM*(m.Bool_dis[i]))\n",
    "    m.Batt_cd3=en.Constraint(m.Time,rule=Bool_char_rule_3)\n",
    "    def Bool_char_rule_4(m,i):\n",
    "        bigM=500\n",
    "        return((m.negDeltaSOC[i])>=0-bigM*(1-m.Bool_char[i]))\n",
    "    m.Batt_cd4=en.Constraint(m.Time,rule=Bool_char_rule_4)\n",
    "    def Batt_char_dis(m,i):\n",
    "        return (m.Bool_char[i]+m.Bool_dis[i],1)\n",
    "    m.Batt_char_dis=en.Constraint(m.Time,rule=Batt_char_dis)\n",
    "\n",
    "    #ensure charging efficiency is divided\n",
    "    def pos_E_in_Grid_rule(m,i):\n",
    "        return m.posEInGrid[i] == m.posDeltaSOCGrid[i]/m.etaChg\n",
    "    m.posEInGrid_cons = en.Constraint(m.Time, rule=pos_E_in_Grid_rule)\n",
    "    def pos_E_in_PV_rule(m,i):\n",
    "        return m.posEInPV[i] == m.posDeltaSOCPV[i]/m.etaChg\n",
    "    m.posEInPV_cons = en.Constraint(m.Time, rule=pos_E_in_PV_rule)\n",
    "    # ensure discharging eff multiplied\n",
    "    def neg_E_out_rule(m,i):\n",
    "        return m.negEOut[i] == m.negDeltaSOC[i]*m.etaDisChg\n",
    "    m.negEOut_cons = en.Constraint(m.Time, rule=neg_E_out_rule)\n",
    "\n",
    "    # ensure battery can only discharge to demand level\n",
    "    def E_netLoad_rule(m,i):\n",
    "        return m.negEOut[i]+m.posLoad[i]>=0\n",
    "    m.demand_cons = en.Constraint(m.Time, rule=E_netLoad_rule)\n",
    "\n",
    "    # ensure PV charging can only come from PV\n",
    "    def E_PV_rule(m,i):\n",
    "        return m.posEInPV[i]+m.negLoad[i]<=0\n",
    "    m.EPVchargingLimit_cons = en.Constraint(m.Time, rule=E_PV_rule)\n",
    "\n",
    "    # ensure charging rate obeyed\n",
    "    def E_charging_rate_rule(m,i):\n",
    "        return m.posDeltaSOCPV[i]+m.posDeltaSOCGrid[i]<=m.ChargingLimit\n",
    "    m.chargingLimit_cons = en.Constraint(m.Time, rule=E_charging_rate_rule)\n",
    "\n",
    "    # calculate the net positive demand\n",
    "    def E_pos_net_rule(m,i):\n",
    "        return m.posNetLoad[i] == m.posLoad[i]+m.posEInGrid[i]+m.negEOut[i]\n",
    "    m.E_posNet_cons = en.Constraint(m.Time,rule=E_pos_net_rule)\n",
    "\n",
    "    # calculate export\n",
    "    def E_neg_net_rule(m,i):\n",
    "        return m.negNetLoad[i] == m.negLoad[i]+m.posEInPV[i] # EInPV only (PV rule) and battery cannot discharge (netLoad rule)\n",
    "    m.E_negNet_cons = en.Constraint(m.Time,rule=E_neg_net_rule)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if billingOption[0]==2:\n",
    "    modelPrice1 = np.zeros(len(price))\n",
    "    for j in range(len(price)):\n",
    "        modelPrice1[j] = price[j]+(j)/1000000.\n",
    "    modelPrice2 = np.zeros(len(price))\n",
    "    for j in range(len(price)):\n",
    "        modelPrice2[j] = price[j]+(8760-j)/1000000.\n",
    "    priceDict1 = dict(enumerate(modelPrice1))\n",
    "    priceDict2 = dict(enumerate(modelPrice2))\n",
    "    \n",
    "    #modelPrice1 = np.zeros(len(price))\n",
    "    #for j in range(len(price)):\n",
    "    #    modelPrice1[j] = price[j]+(j)/1000000.\n",
    "    #priceDict1 = dict(enumerate(modelPrice1))\n",
    "\n",
    "    # split load into +ve and -ve\n",
    "    posLoad = np.copy(load)\n",
    "    negLoad = np.copy(load)\n",
    "    for j,e in enumerate(load):\n",
    "        if e>=0:\n",
    "            negLoad[j]=0\n",
    "        else:\n",
    "            posLoad[j]=0\n",
    "    # print posLoad\n",
    "    # print negLoad\n",
    "    posLoadDict = dict(enumerate(posLoad))\n",
    "    negLoadDict = dict(enumerate(negLoad))\n",
    "\n",
    "    # model\n",
    "    m = en.ConcreteModel()\n",
    "\n",
    "    # set\n",
    "    m.Time = en.RangeSet(0, len(price)-1)\n",
    "\n",
    "    # variables\n",
    "    m.SOC = en.Var(m.Time, bounds=(0,noUnitsSim*storageProperties[0]), initialize=0)\n",
    "    m.posDeltaSOCGrid = en.Var(m.Time, bounds=(0,noUnitsSim*storageProperties[1]), initialize=0)\n",
    "    m.posDeltaSOCPV = en.Var(m.Time, bounds=(0,noUnitsSim*storageProperties[1]), initialize=0)\n",
    "    m.negDeltaSOC = en.Var(m.Time, bounds=(noUnitsSim*storageProperties[2],0), initialize=0)\n",
    "    m.posEInGrid = en.Var(m.Time, initialize=0.)\n",
    "    m.posEInPV = en.Var(m.Time, initialize=0)\n",
    "    m.negEOut = en.Var(m.Time, initialize=0)\n",
    "    m.posNetLoad = en.Var(m.Time, initialize=posLoadDict)\n",
    "    m.negNetLoad = en.Var(m.Time, initialize=negLoadDict)\n",
    "\n",
    "    m.Bool_char=en.Var(m.Time,within=en.Boolean)\n",
    "    m.Bool_dis=en.Var(m.Time,within=en.Boolean,initialize=0)\n",
    "\n",
    "    # parameters\n",
    "    m.price1 = en.Param(m.Time, initialize=priceDict1)\n",
    "    m.price2 = en.Param(m.Time, initialize=priceDict2)\n",
    "    m.posLoad = en.Param(m.Time, initialize=posLoadDict)\n",
    "    m.negLoad = en.Param(m.Time, initialize=negLoadDict)\n",
    "    m.etaChg = en.Param(initialize = etaChg)\n",
    "    m.etaDisChg = en.Param(initialize = etaDisChg)\n",
    "    m.ChargingLimit = en.Param(initialize=5.)\n",
    "\n",
    "    # objective function\n",
    "    def Obj_fn(m):\n",
    "        return sum(( (m.posEInGrid[i]+m.posEInPV[i])*m.price1[i]+m.negEOut[i]*m.price2[i]) for i in m.Time)\n",
    "    # ...\n",
    "    m.total_cost = en.Objective(rule=Obj_fn,sense=en.minimize)\n",
    "    \n",
    "    # boolean constraints\n",
    "    # not actually needed - no need to formulate as MILP, LP is sufficient\n",
    "    def Bool_char_rule_1(m,i):\n",
    "        bigM=500\n",
    "        return((m.posDeltaSOCGrid[i]+m.posDeltaSOCPV[i])>=-bigM*(m.Bool_char[i]))\n",
    "    m.Batt_ch1=en.Constraint(m.Time,rule=Bool_char_rule_1)\n",
    "    def Bool_char_rule_2(m,i):\n",
    "        bigM=500\n",
    "        return((m.posDeltaSOCGrid[i]+m.posDeltaSOCPV[i])<=0+bigM*(1-m.Bool_dis[i]))\n",
    "    m.Batt_ch2=en.Constraint(m.Time,rule=Bool_char_rule_2)\n",
    "    def Bool_char_rule_3(m,i):\n",
    "        bigM=500\n",
    "        return((m.negDeltaSOC[i])<=bigM*(m.Bool_dis[i]))\n",
    "    m.Batt_cd3=en.Constraint(m.Time,rule=Bool_char_rule_3)\n",
    "    def Bool_char_rule_4(m,i):\n",
    "        bigM=500\n",
    "        return((m.negDeltaSOC[i])>=0-bigM*(1-m.Bool_char[i]))\n",
    "    m.Batt_cd4=en.Constraint(m.Time,rule=Bool_char_rule_4)\n",
    "    def Batt_char_dis(m,i):\n",
    "        return (m.Bool_char[i]+m.Bool_dis[i],1)\n",
    "    m.Batt_char_dis=en.Constraint(m.Time,rule=Batt_char_dis)\n",
    "\n",
    "    # constraints\n",
    "    def SOC_rule(m,t):\n",
    "        if t==0:\n",
    "            return (m.SOC[t] == m.posDeltaSOCGrid[t]+m.posDeltaSOCPV[t]+m.negDeltaSOC[t])\n",
    "        else:\n",
    "            return (m.SOC[t] == m.SOC[t-1]+m.posDeltaSOCGrid[t]+m.posDeltaSOCPV[t]+m.negDeltaSOC[t])   \n",
    "    m.Batt_SOC = en.Constraint(m.Time,rule=SOC_rule)\n",
    "\n",
    "    #ensure charging efficiency is divided\n",
    "    def pos_E_in_Grid_rule(m,i):\n",
    "        return m.posEInGrid[i] == m.posDeltaSOCGrid[i]/m.etaChg\n",
    "    m.posEInGrid_cons = en.Constraint(m.Time, rule=pos_E_in_Grid_rule)\n",
    "    def pos_E_in_PV_rule(m,i):\n",
    "        return m.posEInPV[i] == m.posDeltaSOCPV[i]/m.etaChg\n",
    "    m.posEInPV_cons = en.Constraint(m.Time, rule=pos_E_in_PV_rule)\n",
    "    # ensure discharging eff multiplied\n",
    "    def neg_E_out_rule(m,i):\n",
    "        return m.negEOut[i] == m.negDeltaSOC[i]*m.etaDisChg\n",
    "    m.negEOut_cons = en.Constraint(m.Time, rule=neg_E_out_rule)\n",
    "\n",
    "    # ensure battery can only discharge to demand level\n",
    "    def E_netLoad_rule(m,i):\n",
    "        return m.negEOut[i]+m.posLoad[i]>=0\n",
    "    m.demand_cons = en.Constraint(m.Time, rule=E_netLoad_rule)\n",
    "\n",
    "    # ensure PV charging can only come from PV\n",
    "    def E_PV_rule(m,i):\n",
    "        return m.posEInPV[i]+m.negLoad[i]<=0\n",
    "    m.EPVchargingLimit_cons = en.Constraint(m.Time, rule=E_PV_rule)\n",
    "\n",
    "    # ensure charging rate obeyed\n",
    "    def E_charging_rate_rule(m,i):\n",
    "        return m.posDeltaSOCPV[i]+m.posDeltaSOCGrid[i]<=m.ChargingLimit\n",
    "    m.chargingLimit_cons = en.Constraint(m.Time, rule=E_charging_rate_rule)\n",
    "\n",
    "    # calculate the net positive demand\n",
    "    def E_pos_net_rule(m,i):\n",
    "        return m.posNetLoad[i] == m.posLoad[i]+m.posEInGrid[i]+m.negEOut[i]\n",
    "    m.E_posNet_cons = en.Constraint(m.Time,rule=E_pos_net_rule)\n",
    "\n",
    "    # calculate export\n",
    "    def E_neg_net_rule(m,i):\n",
    "        return m.negNetLoad[i] == m.negLoad[i]+m.posEInPV[i] # EInPV only (PV rule) and battery cannot discharge (netLoad rule)\n",
    "    m.E_negNet_cons = en.Constraint(m.Time,rule=E_neg_net_rule)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#print m.total_cost.expr\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13.2353301048\n"
     ]
    }
   ],
   "source": [
    "# set the path to the solver\n",
    "opt = SolverFactory(\"cplex\", executable=\"/opt/ibm/ILOG/CPLEX_Studio1271/cplex/bin/x86-64_linux/cplex\")\n",
    "\n",
    "t = time.time()\n",
    "results = opt.solve(m)\n",
    "elapsed = time.time() - t\n",
    "print elapsed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "SOC = np.zeros((len(price)))\n",
    "posDeltaSOCGrid = np.zeros((len(price)))\n",
    "posDeltaSOCPV = np.zeros((len(price)))\n",
    "negDeltaSOC = np.zeros((len(price)))\n",
    "posEInGrid = np.zeros((len(price)))\n",
    "posEInPV = np.zeros((len(price)))\n",
    "negEOut = np.zeros((len(price)))\n",
    "posNetLoad = np.zeros((len(price)))\n",
    "negNetLoad = np.zeros((len(price)))\n",
    "j = 0\n",
    "for v in m.component_objects(Var, active=True):\n",
    "    #print (\"Variable component object\",v)\n",
    "    #print (\"Type of component object: \", str(type(v))[1:-1])\n",
    "    varobject = getattr(m, str(v))\n",
    "    for index in varobject:\n",
    "        if j == 0:\n",
    "            #print j, index, varobject[index].value\n",
    "            SOC[index] = varobject[index].value\n",
    "        elif j == 1:\n",
    "            posDeltaSOCGrid[index] = varobject[index].value\n",
    "        elif j == 2:\n",
    "            posDeltaSOCPV[index] = varobject[index].value\n",
    "        elif j == 3:\n",
    "            negDeltaSOC[index] = varobject[index].value\n",
    "        elif j == 4:\n",
    "            posEInGrid[index] = varobject[index].value\n",
    "        elif j == 5:\n",
    "            posEInPV[index] = varobject[index].value\n",
    "        elif j == 6:\n",
    "            negEOut[index] = varobject[index].value\n",
    "        elif j == 7:\n",
    "            posNetLoad[index] = varobject[index].value\n",
    "        elif j == 8:\n",
    "            negNetLoad[index] = varobject[index].value\n",
    "    j+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "netLoad = posNetLoad+negNetLoad\n",
    "# print posEInGrid\n",
    "# print posNetLoad\n",
    "# print negNetLoad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(SOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.130748987198\n",
      "2.87181401253\n"
     ]
    }
   ],
   "source": [
    "t = time.time()\n",
    "[newDemand1, oldTotalPrice1, newTotalPrice1, SOC1] = cycleCharging(10, storageProperties, -load, price, billingOption)\n",
    "print time.time() - t\n",
    "t = time.time()\n",
    "[newDemand2, oldTotalPrice2, newTotalPrice2, SOC2] = priceDriven2(10, storageProperties, -load, price, billingOption)\n",
    "print time.time() - t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f21ddb21110>]"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAHaCAYAAADsYJkHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8FOX9wPHP7JWTQDgEAgIV8VHReuJ94FW1av3Vu56t\nra3Wqq21tv3VKt6t1p9atVarrRU8QOqBqFUBxQNQQC4BH+5zIYGEJJtjs9f8/thNssnuJrubPZLd\n77svS3Z25plnn539zjPPPPM8hmmaCCGEyA+WbGdACCFE5kjQF0KIPCJBXwgh8ogEfSGEyCMS9IUQ\nIo9I0BdCiDwiQV8IIfKIBH0hhMgjEvSFECKPSNAXQog8Yst2BsIZhiFjQgghRBJM0zTiWa9XBX2A\nZMcCcjqdVFRUpDg3fZ+USyQpk+ikXCL1lTIxjLjiPSDNO0IIkVd6VNNXSh0EvAU8qrV+Uim1NzAZ\nsAI7gKu01i09z6YQQohUSLqmr5QqAZ4AZoctvgd4Smt9IrAOuLZn2csdLnczHq8329kQQuS5njTv\ntADfBZxhyyYCM0J/vw2c3oP0c8qP37qVq6b+LtvZEELkuaSbd7TWPsCnlApfXBLWnFMFDO9B3nKO\naW/KdhaEEHkunb13ur2drJSaBNwVvszpdEZfuRsulyvpbTMpkTzurKvjv2uWcPmhx+Ow25PaX18p\nl0ySMolOyiVSLpZJqoN+g1KqSGvdDIygY9NPBK31JGBS62vDMMxku0f1la5VieTxd3OewW2vomC1\ng1+fdVFS++sr5ZJJUibRSblEysUySXWXzVnAhaG/LwT+m+L080qzpRqAysaaLOdECJErkq7pK6WO\nAB4BxgBepdRFwBXAC0qpnwGbgX+nIpNCCCFSoyc3chcT7K3T2RlJ5yZNql31PPzhVK47/lzGDpV7\ny0KI/JUXT+Q+9MGrbPB+xd0f/j3bWRFCiKzKi6Df4G0AoIWGLOdECCGyKy+Cft8ng48KIVJDgr4Q\nQuSRPAv6UmMWQuS3PAn68Y81LYQQuSxPgn7X6hobZQRMkXOce/Zw62t/Y33ljmxnRfQieRH0u6vn\nXzfzNhkBU/RKO2r38OzH7+H1+xPe9k8fvMi2wAru/fAfaciZ6Kt63XSJ2SIjYIre6LczH8Ntr8L7\nsZ8bTzs3oW2bfI1gA48px7Zolx81/QTmjxSiN2m2VQGwo76qB6nI8S/a5UXQ783WV+7EHXY/Yd7a\nb1i1fVsWcyR6I+l3JlJFgn4Wra/cwe8/vpufTr2vbdljXz3OpM/uz2KuhBC5TIJ+Fq3eEazRu+09\nuXQXfcmD777KvTNfynY2RB5L6Y1cpVQp8CJQDhQAd2ut30/lPnom/RfJPr+fRreb/iUlad+X6HuW\nuOaG/roizi0MwMRM4tiVJiERTapr+j8EtNb6FOAi4PEUp58RN099nJunJpf1H79yL9fNvA2Xu7nb\ndc04f5XJ/OCFECKaVHfZ3A18O/R3eeh1rxFv6NzJN0nvo9leCUBVXS39CouSTgekz4UQIvVSWtPX\nWr8KjFJKrQM+AW5LZfrJMvpA+Jy5dGG2syB6ox5d5MkVooiU6jb9K4EtWuuzlFKHAM8DR3ax/iTg\nrvBlyc4839Ws9T6fL/RJzS7T7/xesnkBqKraTZE/8kcXnmZdXV3b31NWTOPwvUZ0WM8kWNv3erxp\nKZd8s3rnDl78eiaX7Xt6trOS8Hfi8XgS3sYfCABgml0f963kWImUi2WS6uad44H3AbTWy5RSFUop\nq9Y66jPkWutJwKTW14ZhmMnOPN/VrPU2W/vH7Cr9zu8lmxeAvfYaTEWUqRnD0+xfuR0qo79XUVHR\ndn1id9iTzktX5ZJvfjfnGdwFVUxdN5uzj7szq3lJ9DtxOBwJb2OxBC/kLYYlrm3lWImUi2WS6hu5\n64CjAZRSo4GGWAE/K3p/K49Io4AZPBQDBLKck9SY+81Kpn85L2P7e2PRfKrq67pfUfRqqQ76zwBj\nlFJzgZeB61OcfpKyEe3lDNPb9IV7O4l4atmTTNs4OSP7mrl0Ia+sf5Fb33ooI/sT6ZPS5h2tdQNw\nSSrT7M3qGhtp9HioKC+P8m73N9HMHKlxikxI5qZs6m7kbtgVHJ7Z46hJWZoiO+SJ3B74ydu/4Zcf\n/G+2s5FTXvh0FvfNfDnb2WDKvI+5esqd1DXFP0JlXWMjvjiHQA4EMnnC7/kVjoxZmDsk6PeAYXRV\nk4r+KwkEAry/fAlNLe70ZKqPe9f5H5Y3fprtbDBj61Tc9l1MX/RZXOvXNTZy3czb+Mkr93W/MhCI\n98k8IVJMgn6avLJwTtTl//jkfZ5f/Sy3v/F0hnMk0mlHXS0ATfadWc5JekhFP3fkRdDPxnj6sWqr\netcGAKr8Gzu+YZhsq66OnphUCrPGjLNGXlpYmFi6CecjwQ2S2UkXZE6K3JEXQb8vCNgbuXXWHZ2W\nyg8tPVJ/FrUYif2UzLjb9DseAx6vN+77AWan5scGt5ufv/IIc1Yuj3Pf6TX9y3ls2d2rRmrJC3kR\n9Htr6JRm3d6lrqkp4l7LBzvf5Kapj3VY9u6yxbyzdFGHZQEzsRuzid7GNTHx+f1c+frN/OSVBxLc\nOuhfn33Abss6/v51Mk2Lqf0Vfb5mNdM2TuY3/03us4jk5ewcuY9/+CZ2q42fnxo+r2j2o2z2cyBi\nue7tX2P6rbx2+ZMdlleiO7x+4ZvnADjn0JgjjES1aMO6tr/jbTYK1+z1ANBg257wtgAenyep7SD1\nvXe21uwCwLR3PxqtSK2crel/XvM+H+96J9vZEL1KMHJ12efKmr4HyB9a+EjS2yZzkmiVingtbfq5\nI2eDfm+V2E9Hrgv6jChf1Y9eupefvvzn6Ksn2ByUzJHQ2Hlk8wTvO4jclBdHQbyP3/9g8m/Z09iY\n5tyITGhqcTP9y3kdJp3PtEabk1rrpqjvBeKN4mb7sRuIsVEgEKCytpbfv/4cm3cHm01mLl0Ittbm\nnBQ8nJXiNn25csievAj68fI76pm59Iu07iOxGlv+/jB6+sTq3TNfZNrGyTz4btjTvb2oOBOfDS32\n+gHT5IEPJrPeu5h73n8WgNU7Nvcgd5EkSOeOvA/6ndtKe9J2KlKnp0+sbmsKBr0tDVtSkZ2sMzG7\n7O7V6HUB0BxoaF+/k57E7WyE/Kdmv83dM6ZkYc+5TYJ+hoN8+I9HTjCpV+2q59Ipv8YbGhjMMDMT\nrhI9ScVqquksntxv35O9QdAeeOdVrn3pfjbtqkp52nN3v8vK5s9Tnm6+S3nQV0pdoZRappRarJQ6\nJ9XpJ0UuTfucZE+HM5Z+gWnvbpC0XnCyTfSEb8bu2/+b2Xd23SMJg0AgwLKq1D6UNX+dZmnDXBps\n27j3g+dSmrZIn5QGfaXUIILTH54AnAucn8r0e6sZS75gzqoV2c5GTon/idVO28WIfvPWfkOLfVcP\ncpRa8X+61m6mZtz3Oeav0+xo6FjzfmfZYprtlTG2iD8f4WoaXG1/twQ6doDI7CiiIhGprumfDszS\nWru01ju01j9Ncfo91H3tKplDdcqaF/j7ir+lKAfRtukFNdMe+u1/nuX6Vx6Oe/1Uh4zXls5OcYod\n9fQ7csc1vELsfYSH5EcXP8b2wNcd3t9eG/uEF2u/by5ewI9fuh+XuznqvmPd3P1o9Qoue+1GJn/+\nUcx9iuxJddAfAxQrpWYopT5VSp2W4vSTIo072bfRt4Qay4as7b+3zZoVHmgb3G6ufv1mfvpK9D79\nANvc65n07j9jvt/dKSfWSWntzh1c9tqNPPTetIj3Xl73b1y2bby84OOo28Yq0WlLPwDgnY2xH47s\nbd9HPkn1MAwGMAj4PjAa+EgpNVprHfWIU0pNItgc1CbZmedjzVrvdDrxer1gjZ6+r1Mtp7GhIem8\ndF4/2va+sLzU1dd3mZZJsEC9Xm/KyyVb4s2L0+mk0G5POP3GxoYOrwOBAE6nE5/P11bFMU2zy3zE\nOo66WrartjbmutG23VlZSYsr2DyyLnQTtN62JXLd0C8nYG/sMBxE5/UCfj/YoncOCJgmO2s61vRb\nt5/2ZbA2vqj+I853HRI1r/UuFzaLNWLb+rDjNxBWpn5/sKw7l/PnG9YyY8NsfjXhB7hc9RHptdpc\nXc0H65bEfD+Tov1+ttXu4T29mKsOPzmpYzTbUh30K4F5WmsfsF4p5QKGAFFv7WutJwGTWl8bhmEm\nO/N8rFnrKyoqsNsdbe0FndfpPNNRSWkp1EVftzud14+2vd1ub8tLWb9+MUomuK0Rtk2qyyVb4s3L\nsGFDKS5IbLhigJJ1pRB2LrVYrVRUVFBgd0DoqzYMo8t8LKmKDDKxjq1WXoc95rqzNqzm6hM6XvTu\ntddeDCkrA6A27OnciooKVmzZzKLNa/jRiWcQqz7deR+NBcE8R2tysRgGRUVF0N4E37Z9cXERtASX\n2QoKon7OoqJCHDY7NHTcdsDuSnC27sPSttxqtbXlJTy91z5+AAr8TFv9GYePOACqo3+WX89+HL+j\n/UvM5vEb7ffz29lP0+LYzfRVhdx+dt+bHTbVzTsfAKcqpSyhm7ql0PlZ8F6uU03J5/fzjTO5Aa6i\nJp+ylHJb3E+sdhKrGaO7h4umfdE+Q9a/9fNx78/lbmbu6pVddr+duX16xLKu7gHcO/9PvOd8na+3\npu4Zg3gervrbgujNMcHvInL78DQ9jmoa3F3PBtd5qOdYwgN+b+S27gFgd1Mtry74hLnfrMxyjhKT\n0qCvtd4OTAcWAO8BN2mtM3Ybf/LnH/HVpvUdlnWY4zSJZsRbpz/JnZ/eFzGUbib19fbPXWHNAPPX\naa6echffbN/W5TYvL/iIldt6HvRMAjS43d2Odz9900tJpX/L9Ed4avmTfPTNsqS2h9jj8Ozp1FSV\nPCPmSSl8aaVnR4x1ouev83nk0Q87ntxindhy6fmU1ze/wlPLnux+xV4k5f30tdbPaK2PCv03I9Xp\nx1LlquftbdP40xd/6bD8ujd/26N0d/INAAu3rOpROq3CfyfxPtrel3vvvPDpLG58r/07+Ov8f+K2\nV/Ho3K4nP59VOYO7P38w7v0EAgG21UTOPNZoc3LtW79K24mzwR68Cly/e2tC25lhlzKxYmCg7a5O\ndnW+6npj0XwgsjKyu3lP3GmGH/rfbN/G9C/ndbvN11u38OSsjIWU2OK8YumtcmY8/caWluhvWH3t\nf0f5rrL59eVSjQdgY1UlowcPwWJpr0t8tuXLTkdZej7zTdMeY5exFlVwVNT30z12TKKfKtDDcvj7\nR+8mtH58n799nY1V7X36O9f0P1j7Od8/8lgs3aQZ74n2zs/uB2DCPuMYPXhIzPXu/uzPGJYA6utR\nnHHQoXGl3Zt9uX4tJQUFjB85KqP7zfthGDqLeUmasvR7j1SOQPnal5/x248m8ce3XuiwPFOfd5ex\nFoBtDdGbjSxhASgdV07haT4xawaXvHRLh/ef/+SDjuvHyMJVk++MSDmaOVVvJ5S/WAHY74/edPOX\nOe3NXZ0flGurrHQT9AP2Rq6cfAdef3xzFNR00Zz1mV6NYQnmY3dDltv8UzS0x18W/V9CV7Opkl9B\nv0ffVaoDRXYv26d/OY+rX7854RpjLPM3B9u01zYtTUl6yUvuRm4q9/tp9XthwxoHvb/jjdhbhp0B\nWhztXSvnrl0CtvQNDV3X2MjC+ugPrfkC7VfInZt3Wk9wnYu09XOEL/Y4qlmzo2NHiGR+SU8ujP/m\netr18eadvAj6qfi9u70xmo+SZfOweOvq1KaZgPfXfQLA3G1zY65T7XLR7Il3ir1470+kV6z0Pf7s\njasfTTxz6n7dlKrBxkw6l8xNUx/jupm3JZlaMK3ZunPnBjPs/6NtmPwPMVVXZ3WNjVzy8s3c+3Zi\no3c+/dE7PPrBmynJQ7blRdDf6o9/XJxYl90bfUto9nhYX7mDa6ZMSrqbVvhhv8azMKk0MuWGd2/n\nmtduT2wjo8uXaRe7eS71p5vwq6SEb89kvLLY8ZvoPO9vVzqXnYlJg9vNmpYvk8xLrF496e/o9+na\nVWD1siLBE+pHVTOZv+f94IsMjdyaLnkR9NtFu5Mb/69vT2Mjj388jWZ7JX9blNzlZjYuDHfV1/Oj\nKffy3+VfJZ4TW3xXOL1lIFO/GV+N/h9z3+dfn37Yo32Ft6snelIJH4q5pzd145Fo76Wu1jdNkxZv\n/JOsG71omsbwz/XVpvU89uGbVLtccd93CCYS/L42+ZZ0s2Lv1Hu+jQzaVl3d9iUn0oMmYAbiuixP\nNX/Azw8m3879MyO7OW6rjuym2Nkzn8yk0e7kn6v/kY7sdSnd4WxH7Z44f7Adg9iHO9/kPefr6clU\notLci6vnVzmdt+86vc6ni2g1eCNK6Il7YM4UldefvvgL82reD17RvvyHlKTZF+Rd0F+1fRu3zrqD\nG16NPbhVq7e+6tnUiTdPfbxH27fa46nG73CxrPHTDsunzPuYW2fdwb0zu36wKLMnqvh+kHvMbT0e\nfte5Zw+3vP+//OSVexPcMvGg8bOXH+LtJalrjuuNz150qN2H/bmo/iNmVb7V9ro35j1esW7o+xx1\nGc5J9uRd0F+2JfjEbr0t+sM04Qf0S2tf6PxmQlof7EqX+VuCzTVf137VzZrpl/DDTzYvryz4pEf7\nXLszOPBLfOPE9yxQ7bFuZPKaLka5TLD2aXZo3km/eJrfWgp2ccnUG1i2eWO363Y1TEbnt1qbdwxr\n8IrM5XHRp0mbvkhEKg+Xvv5w1/rqrodi6FaUwox18ulttdNMfnceRzW+QPxt1lO/mtPl+z29WV7j\njz7cQybueKXk9yddNvu2RL4+X8BPJWsytr92sQJZ6N9upgfsKsCkf1yf+NuD31+exhtjvex32tOJ\n3xO1uib+7sEbG7ru2eOy7OCOmX+P+X7nIyoyyJh4w54D6KlAIMAPp9zD5ZN/z5odyQ/DXFVfx3Nz\n30/spm7IMx+/x2WTf0O1q3cPFgcS9BOKBa8snJOlXirRc1nj7z3j5Mcr4Ij91OXzq59NKC1LAies\ndNf0e9k5JUKsQdOi8TtcXVYGDEuAPdbIJqBYlYvOezbtzbyx+ZW48xORv0772VVfT5N9Bz5HLXfM\n6fpeXVcP6d0+43E+2PkmV0z/Bb949dGE8jS7cgYBRwNvfDU/oe2yIf+CfjdRu6ta8a7G7nvKZJJp\nJF4jSZdkrhi21PdsFM1Xl/43qe3SE6CT76ifif7p7gTnB07mJLnLF7u5rvMwFNEETJMPVnT/RPeM\nrVNjv2mLvytpZ43W9ntDVUZyV/S9/eQPaQr6SqkipdR6pdQP05F+tqTiC+0tt4B6w+2AettWFm1Y\n3/2KMbSOtxMuvvbm1H/4Tb7Ehp8IvxHaG+cQT2ZqS9MS/RkJCzBva3xDk09ZFnu4inDLtmzqdvz+\naPr6MOWpkK6a/h1ATZrSzprOk00nI12x9vb/PNPjLpCp0c0n9BZ0eLlpd/c9bz5avYI5q1bgcjfj\nSXKQuN3uxGq66ZauG7lZDWrW2N9NPGMfJVIm98//M9e/dl/c68en599JX+hckfKgr5TaHzgQiD0r\nci/2wc6+Ob7GJt9SvunBTax4vbl4QfT+6jF+05E/gcSD0tPL/8bfV/yNH791K1e+fjMA7y5bHHVd\nH9Ev75vtOxPebzRrd8bqedITqQsU2e6l5PP72W1ZF7E83m89kZOWxxFsbjUs8W+TyYH34nXT1MfS\nkI/Y0jGe/iPAL4BrulsxlROjN4bNkBUrDZPgRMet3l4wj4OGj0hqf9HEM6F2+MTo8aYVCGsL+HLV\nSkYOKA++6HR8VVVVMqDTI+8ulwt32FwDrekGAu1PJHdX5uHvv7zu3wAcMbRjuXk93uBk2BF573zf\noWOm6131CU9SvlR/w7srP41ajgFH933Aw2tjiR5vc1ak5pmIXbt2URrKRnVN6i6Ks13T/NdHkfdZ\ndlXvptaytdtwXl1dHfVqNfgdRf9cTqeT6oaGKOtHV1tXG/O9aGl3nhi9u+OlsbEp4WOqEp3Ryd9T\nGvSVUlcD87XWG5VS3a6fyonRN+5uv3yPlYZBaDLy0G9s5obPOP3wXya1v2i6mzwbwBY2MXq8aVnC\najKPLn6RV64K9VDo9CsaNGhwxP6cTicFBQXQ3DFdiyUYMbubJDzaZ4i2zOFwgK/9vWumTMJuFLTt\np02nPJf1K+ty/9Hee3Dhs5j25i7z3JXw2l6ix9uAsrKUNFwOGjSobd/bmht7nmBI+muyXfuwembE\nsiGDBmPE0be9fNDAyOOF1u8o9gTxtvq6KOvH2EdNFcQ55XVFRUXExOjdHS/FJcVxHVP//KTjuE+Z\nnPw91c075wDnK6UWAD8B/qiUOj3F+4gqfNAq555Y07Z1PPCidTtLh/DaS2L3BYIHelNY00SmJo3u\naY2x2V5JvW1Lt7W7RLoStm3Tg4DfU93NFhWv8E+d7SaZtIt3WlDTTPKORPStrpkyiV9N653z1/53\nR8dxn34w+bdJPR+QjJTW9LXWl7b+HWq62aS1npXKfcQSHqPmfrM87u1c7uwFkK4kc1M2Zi0vakzp\nOtB0jvlev7/DBOdR9t5lesnqep+ZZ41SE+2xFMZ8n6V3Hs/xCHQ1tkMX3DHmfGi2V7Ld7NhRIN03\nupOtLPkd9Tj31HQ5XWSq5Ew//fDa0k5X/P3pF6xP7/g4PZP6GuDdM6YQCAQIdFPD7rznn73yIL/8\n4H8j1ttWU82S8LFaOl3GR36CxD5T+KTqqeK1trcB3/DKI8xdHf/cCKkKGeF98/0p7aff+64aZq9O\n79hQt83+Y1rTf31R+wNX3VXGenLVlqmeV2mbGD3UXp8x4WfY+XtijJNujbx8mrLmhST3GP8XlMwj\n9wHTxCTZy93YVjZ/zrvLFQ22YMOmz1HH9C/nsXnPTn595gVt69351r86bNdgj94QeuuHdwAwxnZY\n1PcbbD0cXycdwh7gqbas46nl8TcBvLvm47huxHcn/JDI9s3XdNu4J75j4K9L/wr2xNL+8Ovun48I\nBAJYLMH6bTLNc6+uf7E9rW6+q558l6lqOux2PxnZSybEWdi+DPdlT7aLn0n7qIRdrpQEV3PHsXqm\nbZzMF7Uf0tTi5oZX/sIHK5ayNsFZvTa546wt9+CJyd6gztqzp4hbhQeHlMb8Pl6+ifrHyme6XefG\nqf+Xsv1191VF+y7rmpq4ZdoTfKZX8+6yxfhitd1n6B582mr6vVWXj3CnwWMfv8pjF9+c8HZmGk9O\nnccuaTV5/kdUW9bz3KoknpINBRvDkli+A6bJVVPuRA3YP/F99mHhHQ9yvaZf76lLydVRsqotYcdz\nF7Vpw9Lxe4g2P/SqbVu6PMZXVK4GzsXn92OzBj/03z6ewQ5zFX9dugqANVVnRt3WkqEZxnIm6Ldk\n6M534rL/g46Y4zRGkPH4otcS73zrhVRnqc1/1v4H7C0sb+xdT8ym21y9nENGjQEyM11iNtVaN2c7\nC0m55o1bOL7fGR2W3b/goS63qTG38tZXX7TNxfHLw29hSc0XHZqt1lZvjD4seIZq+jnTvDNtdUY6\nCXVp0YbIJxEhuVr7Xz74T8LbxHvQxBrga0td9Kaob9yJzSD2mY5/GF/s8c3Bm2s+r/yo/UVux/yU\nSOQm57Uv3R+1lg6JB9Yvdic4W5oBr696r+3lowufBLu78ypRJTJqbE/kTNCvt6bj8fjEPLTwkShL\nk/silzbMjWOt5NKOVbPc7FuWVHqd/XXpX1OSTi4Lv/rqrieVSEyDbRszlvRsqtNWiffGMfGY7d1m\nu70vFyaR4SR6ImeCfm+WrZ90IBDgQ72CZm/HmkY2JncXnYRd4ZhJ9k8XibMkHPIS+27K2TvFKaZe\nzrTpZ78oo0vzoyBdvvv6ogW8UxU5gJwEmd6htc/35KVv5dQvMdV6x+ix8RnZr4JVtd09JxRj+O8M\n/S7lUEuz9I6F0nXaG2LMQRur906ipn85LyXp5Ks/vPlP1nkWY8ivsEvJPOcS8x5Agj9Hn8OV0CaZ\nngYzGdK8k6Tu5qXtsG6KD4QnZ83o0fapGutl2sbJKUknX633Ls7S9Jt9SzJH62sbX07JvuMZKC7c\nyuZ53Y4nFXtaSanpixg+qX6P8au/FbHciLOf7+zKnp00hMgkMxAAM7GzY6xgnYkeMl3NAw3Zb4iW\nmn4GpGMO1F2dhpMVIlc9Nis1Extt3r2LJ5Y9kZK0eiZGm36GzgYS9NOsEs31rz2Y8nT9ZiCizpKJ\nCbZbuZry63F/kT2L6uekJJ0X5r2fknR6qsMTwmEy1asud4J+gm1vmeS2V6U8TX8ggNnpoY9M3kP6\ny5To0xUKIZJT25i6yXS6ko45ch9SSs1XSi1USl3Q/RYiGZHTEKZ6iN6ufVO1IWP7EiJgiT3pevx6\nb8UQ4M8fP5eR/aQ06CulTgEO0lofC5wFZHbG3zwSrWtY5LL03bQyRy1JW9pCRLD1fLiO3j5DWYsj\nM+NPpbqm/wlwcejvWqBEKZXF8fVyV7QD2O/vXNPv3Qe5EJkkT6IHpXq6RD/Q2jD1Y+Dd0LKoQlMq\n3hW+LJOzwvdlDQ2R3cJ21+zGWVjUvk5j3506T4hUcjqdfONamvAkLZmWifiXln76SqnzCQb973S1\nXmh2rUmtrw3DMDM5K3xfVlxSAp2mjy0fUE54+VltDuitI04LkUEVFRV9YkTXTMS/lAd9pdSZwB+A\ns7TWGetMnujkHX1dtKf6Mjoue0FmehoIIVIrpUFfKdUfeBg4XWtdk8q0RUfRwnvngankCX8hRGep\nrulfCgwGpimlWpddrbVOzcSiok20B7EyPf+vEKLvSfWN3GeBZ1OZpoguak2/84lAqvpCiE5y4olc\ntzcVD270fZ27bAbkLq4QopOcCPpPzX4r21nIuKg3cjst85kyPo4QoqOcCPrOhtSPbdPbRWve8UcM\nzSDtO0KIjnIi6Kd6kpK+Krymv0BvRDfJoGhCiI4k6PdVUT6z32yv6f9l3tMYVmnTF0J0lBNBf7sn\n+vjU+SYQNrGypdiVxZyIVLF4SrOdBZFjciLo02lc+XwQbcC1N+au5a1PgifAfLz4yU1yX0akVm4E\n/Ty0ZOdqTdwvAAAgAElEQVSKiGUux2aee3spK7dtSXhOUZFew40D41rv0NKTOLb8zLbXBmB4izAD\nFv58yiSmXfp0XOmMsR0ad96snrK40xV9nwT9PqrOGvmQs23IdoqOnMXdnz+IYZGqfra9fNGTbX9f\nfsSZXawZdFi/k7n1Oxdx/cTvgj/43OT5487lpR88zNRLn+Bbew2Na79DzHE8dOHPul2vyBtMz8iD\nMPD1VhkUoFVaRtkUQoDN2j6VhNXoeloJm6c/v//uZW2vp13+BIFAAIslMiDfdMhNPLnoeQopo9m+\nM+L9e865Lq78jR80nkX1lRxYfnBc6/dln69fme0s9BoS9IVIQpF3WNSA28ZXAISaZuzNjBg4MOaq\n4xwTuOH08yOWRwv4ACfufyAn7v8Ic1ev5KnlT2LxlBJwRM6v0NmDE+9k6qKPWdrwCQA3nXY+63Ye\nxfiRe3e7bV83Z9uHvX4s/UzJ/es60aeEt2f3ZveedSMHFR8fsXwoivtPvpMXLvoTAH///r384djf\nMnxAOWdXXEBFqG3//L0vbdvm/u9fy8hBgxLOw8kHjOeysVfzl3N+1+mdYNPewZ3yN3bocP73nB+0\nvS5yODh41OiYJ5eunDfykoS3ySbT3v2EQjZP/wzkJDa7J3bFIJXSMTH6o6GJ0ecppSakOv1o4r1J\nJnq/Cw8/gYGBfXqcTplvFC9e8Ff2sR+e9ImkyDsMw1sU9b1Rgwdz53lXdlh2+4Rf88Slv2TcsOEU\nFxQCUF5SwiGjxgDwoxPP4LFLbuLVi5/iiuMmJpWnzi448lhGDux4wihyOAD443lX8uDEO4Hg/YKe\nsnn68+sjf8WrFz/FVcefwrRLn2Z8UeSJL14nDz6nx3lKJQtWfnXEL7O2/9Gl38rIflI9MfrJwLjQ\nxOg/Bv6ayvRjeeTCn3PV6KsYFBjLeSMvZqTlYKyefm3v27xlvHrxU4wvOq7H+/r+qMu6XykOR5ad\nwlNn/xmrp4zbJ/yae074A/93+n0dV/I5UrKvVBhq7Mco67d7nM7gwL78cP+fxHw/YAb47WlXMygw\nlqGotmaSzu4/+U4enHgnIy0HdVg+xnYY1x90A89d8XsK7Xb+dMF1/Oo7/wOe6MEbiOi5YgYslPpG\n8PTFv+WBM37TvtwfbJc/tvyMiDRevfgpjtxn35j7CNdas759wm1cf9ANcW0Tr/7+UW0nHAjW7l+9\n+KkO9wtO2etcDi2NPAko6wRO3eu8qOkOY38eOOtWjh67X4crg7u+137ie+miJ5l26dO8evFTceX1\nxtPOjfneuIIjOGZAZDmn00l7n8Cx+ypevfgpvrf3pRR7h8e97VX7XdtxgS/YljTEHBex7hBzP6Zd\n+nTECSZT3axT3aZ/GvAmgNZ6tVKqXClVprWu72a7HrFZrRwxagznHRMM6ldxKgCXTG39QRlYLBZ+\ndNzZ3DZ7XtL7OW/kJfzg2JN5Y8urAIy2HsJm/7Iutyn1jeSZy37Hsx+/y9zd7wJQYd2f288OXh6/\nctWfO6z/9Hf/THlJadsPKxAIcP2rD1Nr3ZR0vlPBarFwzvgTeHr58h6l87cf/DqYntXC8ysjR+H2\n+v2MGzacp39wW9uy26Y/zX6DxzCrckbbMrvVypghe/F/F98Y9j3DQxf+NKl8jbYdwpDigZQX9+d/\nDjuWIWVlQDBomgEDw2IyruhQHvh+xxPW7RNuo8pVm1QTyZH7jE0qr9GYAQuGJcDgwiER73XO2w2n\nRK9h33DcWVRUVDBn6tsdln+75ETuOPfymPs+c/j3aWxxYw/duLZYLNx0yE08segfWEw7AUdyDwpm\nYiLzUv9w7jnrRvoVFtG/uBgI5v/K4yZy5XETcTqdVFRUcOnk2zAdwdniCjxDaHHsaktjb+vBnHfY\nBM47LHrDxry13/DXBc9TbhsKGDxy0S8AOHZfxdH7PMV/VyzhhW+eS+8HDZPqoD8MCB/wZVdoWVqD\nfiz9PKNxOTZjMYMfc9TgwfTz7o3LvjXq+scNPJN5Ne8DwfbQFU2fd3j/imM71o5KHMXQRVOhGTD4\n5xV/AOC6k85i7uvvtr4Tc5tB/co6vLZYLIwu25vaxk2xd5QBpmkyvKy865V8drDFN8z1mQcdxvNx\ndqj4y0XBoD5ranvQt1nae8NcNvZqXl3/YjepBJ9bMHwOTFvk6KMPX3h9zC3HFR7BOs8iTtznsIj3\nUhm4e2LCgIksqp/DWQcem9J0/+/0+7q93/DjkyKnwg7ebH6UZo+H29/4G5XoDu8PDnS8KjJNMMIe\nLbFarFEfQEy1Kw49L6J5LJrLDriEV9b/C4BnL72Du95+gU0tK/jR+B9y9iFHdLntceP257hxD0d9\nz2KxdDiWMyHdvXe6fEJIKTUJuCt8WbKzwbtcroht+7kUdaWVlNaNb3tvSPP4qEHfu3Ucl5x4FEdX\n70NVQz0TRu/Dz59349gveA4LuIvYuXNn27q2YZs5YJBiVfP8tjT89QOxlrXPEmnffniHPPnWHol1\nnyUMC+yf0Oe8ZPyxLPl4CZbi6D009ms+m9UtC7EO2B13mt4twZnN7KM6/hhNnw3D5otYf7B3X8oM\nC+a28Rgj26O16XVg2D14Nh/Ak5dfwM2vv4htaHv52rYewY8mfpvnvnmVI8qO7vC5/TV7YR1YRaC5\nBEtRI77qYZQEzJhl06KPoEAtxvTZ8Dc14WwJTnR9TMVoXlo8BIi9bUHlt3GP+JJvm9/h3ANH8ceX\nVlFwwJf4tu/b7Xfx8wlnUt14AkP69Uv6+Ey3Kw85nkv8x+CwWnv8G/LuGIN9+CY8Gw7C0tLS48/8\n+xMv4eb3H8VS3MApRRdz9iFjcdjtwX0598E6yMng2qOpGTq3bZsRjMXq6RgMzYABpoFhDeCvL8dw\nuLEUBmtdgZYizKZ+WMs7jrjbemy16vwbLWgp7PLztZbJEGs/TK+Dgrp92bN7N7889lwg2DzV0/IZ\nbCvA8BQzcsCQjBxfqQ76ToI1+1YVwI5YK2utJwGTWl8bhmEmOxt862VYuEtPPIYHXrBwxZVHtr13\n8kGHsHI6nH3cGN6bt6nD+hUVFR3SGFO6L+u/DF4uP3zziVRUBO+u+3aMxbdjLMfdfiCHuu/ntr9+\nGjVPT/3mFCqGtdfcrz/jdJ6YNpjzf3FUW1pxe/Zkmj3RB1AbdORIKnb3Y9OaxC+ofDvju3k06tR9\nqKio4LBBRzH/y8gufoUOKyNHjsS7eTzezePblg8sK+CMI4/ijCOPitjGs+7wDq8vOnVcxHcYLlA3\nhOYvz2JgWQGjr+iYB8+aYG0r1vaTLjufye/tx00/OJzG+mrwFtKy/KQutwk3sts1+r7W35Bv6/74\ntu4PxFc28Wj5+gQAhn9vNGNGj25b7tu2H75t+/HH353G9X9qv+9yyi0H4/ebvPXkWW3LLj19P6bO\nWgNAeb8C9rhaIvaz99BStla2V45Kiuw0Nne8+vzRuQfyr5mrKC2yc8xFCps1dtNca5lUVFTQr/QO\nvlXRn7KS1N5rq6ioYOoBj6Q0za6kOuh/ANwNPKOUOhxwaq2zNvLXsQcP5/U/n4vd1l5jOPPo0Ywd\n0Z99RvRn/vId1Da0cOh+Q/j+xMibcLdfeSQ/+9NsAAb3bz8gn7ztFFZurGbvocGbxTP+8j0qa5q4\n7oFZbetMvf+7FBd27Bh8xlGjOOnQERQWJF7sZiD2pa7NauGH5x7IpH8siCutw9VenHPCt7j3+S/i\n3n+BI1iGFYNLor7/11+fEnV5VzenvnfSPujNe9Cb9wDgsHXdLj7tgXN4f8FmzjxmdJfrRTN25AAm\nXRds+mjMSmNj31FUYKW5xc+oYf26XzlOZx87hvfmb2L8PtErO1aLwduPnM+eejeNbi8j9wru+zdX\nHsHDU4JX28MGtR97j/7qZBatruLJ15a2LTvliJHcevkRnPfr4KRKg/oX8sgtJ/HwlMWs3FANwNXf\nPYALThnH+Sfvi9WS2FAlh4yLvF/SF6V6jtx5SqnFSql5QAC4MZXpJyM84ANYLAb7jQq2Tb9w15kY\noWXRVAwp5ZX7vktldSODB7QH/dHDyxg9vL0GbxgGwwaV8IuLD207CDsH/Nb1kgn4AB5f+02tkw8b\nydwl29peH3vw8ITu/E+67hgMI7EDvsAeLMfLzlD856N1Ee8Pj3EyOOeE2FcS150ffBJ0o7OO12av\n5Xsndd0+XlRg439O7h1t6LnssV9N5Knpy/j5RYekLM0bLvw2l5+5PwP6Re+N1Vp7Li8rpLysvffR\nSYeNZP/RA5n/9Q4mHjGSx6cuAWBQ/yLOPGY0I/cKjkI6bFAxA/oFt3vlvu9SXdvM0EHFFDps/OnG\nEyL2l2jAzyUpb9PXWnd+UqTXiueLLy2yUzpyQFzpnXnMaBqaPB1OCKly0anjmD5nLQA/v+jbNLq9\nFBfauPa88QzqX8TCVV08HdpJa8B/5JaT+PXjn8S1zb57B8ugsMBGcaGNJndku3+4a88bzwFjBrLP\niO4fePlWRX9uv+rIuPIRyzO/Oy2h9U85YiQfLd7W/Yp5qGJIKfffkHz/+2gMw4ga8F//83k0NHui\nVpJa7TWwmPNDFYInbzuFstL25pXx+0TehC0tslNaJI/fxmL0pglIDMMwk81PtDb9XGKaJo1uHyWF\ntqi1dK/Pzz3Pf8HSNe1dye76yTEMLvbSv3wwtQ0tfL7cya49zfzqB+1t6fNX7KC+0cOJh1bw9qcb\nOObg4VTVNDF25ACuuft9TjpsBN+fuC/7hp343pu/ib9NX8Z3jh7N6RNGMah/IXsNDHZ3e3PuenZW\nN3L9BT3v058uTqeTYcOGM/PzDZx4yIgONct8luu/oWT0lTIxDAPTjG9oXQn6OWbZml1M/u9q7vrJ\nMfQrdvSoXNwtPgoc1qgnGb8/gLWLG2C9mRwr0Um5ROorZZJI0JcB13LMIfsN4ZD9UnPDqav7D301\n4AuR7+SXK4QQeUSCvhBC5BEJ+kIIkUck6AshRB7pdTdyE31oSAghRPx6VZfNnlBKmVprOWN0IuUS\nScokOimXSLlYJtK8I4QQeUSCvhBC5JFcCvp3ZzsDvZSUSyQpk+ikXCLlXJnkTJu+EEKI7uVSTV8I\nIUQ3JOgLIUQekaAvhBB5RIK+EELkEQn6QgiRR3rdMAzJUEo9ChwDmMAtWuuFWc5S2imlHgJOJPgd\nPggsBCYDVmAHcJXWukUpdQXwS4JzFj+rtX5eKWUHXgBGA37gR1rrDZn/FKmnlCoCvgbuBWYjZULo\n894O+IA7geXkcbkopUqBF4FyoIBgt8ydwNMEY8hyrfUNoXV/A1wcWn631vpdpVR/4GWgP9AAXK61\nrsn4B0lSn6/pK6VOBsZprY8Ffgz8NctZSjul1CnAQaHPfBbwGHAP8JTW+kRgHXCtUqqE4I/8dGAi\n8Cul1EDgcqBWa30CcD/Bk0auuANo/QHmfZkopQYBdwEnAOcC5yPl8kNAa61PAS4CHif4G7pFa308\n0F8pdbZS6lvAZbSX3f8ppawET4wfh8rkdeC3WfgMSevzQR84DXgTQGu9GihXSqV+ZvLe5ROCtQ+A\nWqCE4A91RmjZ2wR/vEcDC7XWdVrrZuBz4HiCZfZGaN1ZoWV9nlJqf+BA4J3QoonkeZkQ/MyztNYu\nrfUOrfVPkXLZDbTOqF5OsJLwrbAWgtYyOQV4T2vt0VrvAjYTPL7Cy6R13T4jF4L+MGBX2OtdoWU5\nS2vt11o3hl7+GHgXKNFat4SWVQHDiSybiOVa6wBgKqUcmch7mj0C3Br2WsoExgDFSqkZSqlPlVKn\nkeflorV+FRillFpHsAJ1G7AnbJW4yyRsWZ+RC0G/s5waEa8rSqnzCQb9X3R6K1YZJLq8z1BKXQ3M\n11pvjLFK3pVJiEGwVnsBwWaNf9Hxs+VduSilrgS2aK33BU4FpnRaJZHP3ufKIxeCvpOONfsKgjen\ncppS6kzgD8DZWus6oCF0ExNgBMFy6Vw2EctDN+oMrbUnU3lPk3OA85VSC4CfAH9EygSgEpintfZp\nrdcDLsCV5+VyPPA+gNZ6GVAEDA57P+4yCVvWZ+RC0P+A4M0YlFKHA06ttSu7WUqvUO+Bh4Fzw3oN\nzAIuDP19IfBf4AtgglJqQKjHwvHApwTLrPWewHnAR5nKe7porS/VWk/QWh8DPEew905el0nIB8Cp\nSilL6KZuKVIu6wjew0ApNZrgiXC1UuqE0PsXECyTOcA5SimHUqqCYIBfRccyaS2/PiMnBlxTSv0J\nOIlgV7MbQ2fvnKWU+ikwCVgTtvgagsGukOANpx9prb1KqYuA3xDscvaE1vqlUA+E54BxQAvwQ631\n1gx+hLRSSk0CNhGszb1InpeJUupnBJsBAe4j2L03b8sldFL7JzCUYJfnPxLssvkMwYrwF1rrW0Pr\n3gRcQbBM7tBazw5tP4Vgs1ktcGXoartPyImgL4QQIj650LwjhBAiThL0hRAij0jQF0KIPCJBXwgh\n8ogEfSGEyCO9apRNwzCkK5EQQiTBNM24ng7uVUEfINkupE6nk4qKihTnpu+TcokkZRKdlEukvlIm\nhhH/aBDSvCOEEHlEgr4QQuQRCfpCCJFHJOgL0QfsWfwVzrdnZjsbIgf0uhu5QohIq+65H4CK887N\nck5EXyc1fSGEyCMS9IUQIo9I0BdCiDwiQV8IIfKIBH0hhMgjEvSFECKPSNAXQog8IkFfCCHySFwP\nZymlHgJODK3/ILAQmAxYgR3AVVrrFqXUFcAvgQDwrNb6+bTkWgghRFK6rekrpU4BDtJaHwucBTwG\n3AM8pbU+EVgHXKuUKgHuBE4HJgK/UkoNTFfGhRBCJC6e5p1PgItDf9cCJQSD+ozQsrcJBvqjgYVa\n6zqtdTPwOXB8SnMrhBCiR7pt3tFa+4HG0MsfA+8CZ2qtW0LLqoDhwDBgV9imrcuFEEL0EnEPuKaU\nOp9g0P8OsDbsrVhTtnQ7lYtSahJwV/gyp9MZb5Y6cLlcSW+by6RcIvXlMklnvvtyuaRLLpZJvDdy\nzwT+AJylta5TSjUopYpCzTgjAGfov2Fhm40AFnSVrtZ6EjCp9bVhGGayU5P1lWnNMk3KJVJfLJON\noX/Tme++WC7plotlEs+N3P7Aw8C5Wuua0OJZwIWhvy8E/gt8AUxQSg1QSpUSbM//NPVZFkIIkax4\navqXAoOBaUqp1mXXAM8ppX4GbAb+rbX2KqV+B7wPmMDdWuu6NORZiLxlmmZCk2AL0Vk8N3KfBZ6N\n8tYZUdadDkxPQb6EEEKkgTyRK4QQeUSCvhBC5BEJ+kIIkUck6AshRB6RoC+EEHlEgr4QQuQRCfpC\nCJFHJOgL0ZeYZrZzIPo4CfpCCJFHJOgLIUQekaAvhBB5RIK+EELkEQn6QgiRRyToCyFEHpGgL4QQ\neUSCvhB9ifTTFz0kQV8IIfJIvBOjHwS8BTyqtX5SKfUCcARQHVrlYa31O0qpK4BfAgHgWa3182nI\nsxBCiCR1G/SVUiXAE8DsTm/9Xms9s9N6dwJHAR5goVLqjbDJ1IUQQmRZPM07LcB3AWc36x0NLNRa\n12mtm4HPgeN7mD8hhBApFM/E6D7Ap5Tq/NYvlFK3AlXAL4BhwK6w96uA4SnKpxBCiBSIq00/islA\ntdZ6qVLqd8AkYF6ndYzuElFKTQLuCl/mdHZ3QRGdy+VKettcJuUSqS+XiXPHDgxLevpf9OVySZdc\nLJOkgr7WOrx9fwbwNDCdYG2/1QhgQTfpTCJ4wgDAMAyzoqIimSzhdDpJdttcJuUSqS+WycbQvxXD\nh2NYrWnZR18sl3TLxTJJqsqglPqPUmqf0MuJwNfAF8AEpdQApVQpwfb8T1OSSyGEECkRT++dI4BH\ngDGAVyl1EcHePFOVUk1AA/AjrXVzqKnnfcAE7tZa16Ut50IIIRIWz43cxQRr8539J8q60wk28wgh\nhOiF5IlcIYTIIxL0hRAij0jQF0KIPCJBXwgh8ogEfSGEyCMS9IUQIo9I0BeiDzFlEhXRQxL0hRAi\nj0jQF0KIPCJBXwgh8ogEfSGEyCMS9IUQIo9I0BdCiDwiQV8IIfKIBH0h+hLppy96SIK+EELkEQn6\nQgiRR+KaGF0pdRDwFvCo1vpJpdTewGTACuwArtJatyilrgB+CQSAZ7XWz6cp30IIIZLQbU1fKVVC\ncE7c2WGL7wGe0lqfCKwDrg2tdydwOsHpFX+llBqY8hwLIYRIWjzNOy3AdwFn2LKJwIzQ328TDPRH\nAwu11nVa62bgc+D41GVVCFG3fEW2syD6uG6DvtbaFwri4Uq01i2hv6uA4cAwYFfYOq3LhRApsuqe\n+7OdBdHHxdWm3w0jweVtlFKTgLvClzmdzugrd8PlciW9bS6TconUV8rENE08W7biqOhYd0pX3vtK\nuWRSLpZJskG/QSlVFLoCGEGw6cdJsLbfagSwoKtEtNaTgEmtrw3DMCsqKpLKkNPpJNltc5mUS6S+\nUiY1Xy5k058eZvBJJ3ZYnq6895VyyaRcLJNku2zOAi4M/X0h8F/gC2CCUmqAUqqUYHv+pz3PohD5\nybV2HQC7P/0syzkRuaTbmr5S6gjgEWAM4FVKXQRcAbyglPoZsBn4t9baq5T6HfA+YAJ3a63r0pZz\nIYQQCes26GutFxPsrdPZGVHWnQ5M73m2hBBCpIM8kSuEEHlEgr4QQuQRCfpCCJFHJOgLIUQekaAv\nRG8nY+iLFJKgL4QQeUSCvhBC5BEJ+kL0UobR7fBVQiRMgr4QQuQRCfpCCJFHJOgLIUQekaAvhBB5\nRIK+EELkEQn6QgiRRyToCyFEHpGgL0RvFaOf/ufnX0jVnI8ynBmRKyToC9EHrX38yWxnQfRRSU2M\nrpSaCLwGrAwtWgE8BEwGrMAO4CqtdUsK8iiEiML0+zGs1mxnQ/QxPanpz9VaTwz9dxNwD/CU1vpE\nYB1wbUpyKISIyvn2O9nOguiDUtm8MxGYEfr7beD0FKYtRP7pZuwd15o1GcqIyCVJNe+EHKiUmgEM\nBO4GSsKac6qA4d0loJSaBNwVvszpdCaVGZfLlfS2uUzKJVJfKRNXfX2X77ub3Sn9HH2lXDIpF8sk\n2aC/lmCgnwbsA3zUKa24hgfUWk8CJrVtZBhmRUVFUhlyOp0ku20uk3KJ1FfKxNevH7VdvF9YVJjS\nz9FXyiWTcrFMkgr6WuvtwNTQy/VKqZ3ABKVUkda6GRgB5NbpUYhM625oZZlQSyQhqTZ9pdQVSqnb\nQn8PA4YC/wIuDK1yIfDflORQCBGdTKMokpBs884M4GWl1PmAA7gBWAK8qJT6GbAZ+HdqsiiEECJV\nkm3ecQHnRXnrjJ5lRwgRP6npi8TJE7lCCJFHJOj3Ul6Xi63TpuNraMx2VkSWeHZXd/m+NOmLZEjQ\n76U2/uOfbHnpFTb+S26N5CNvfT2VH87KdjZEDpKg30u5d1YC0LJrV5ZzIrKhpZtaPiBVfZEUCfq9\nnNFdX22Ro+IJ6BL0ReIk6AvRGwUkoIv0kKAvhBB5RIK+EL2QGU97vVwMiCRI0O+15BedLyo/nMXu\nefOznQ2RJ3oytLLIANeatdnOgkizdU8+DcDgt/7TvjCOmn5cVwNCdCI1/SzZ8I/n+fz8C3Hv3Bn1\n/dYftL+pKZPZEr2FBHSRJhL0s2THzHcBWPPYE1nOiei75MQgEpfTQb9y1hzmXXQZLdU12c5KTD6X\nK+pyw8jpryZv+VtaCPh8mH5/l+t1935wpRRlSuSVnG7TX/fEUwBUz5tHxXnnZjk30Zm+6D9uwypB\nPxctuOTytr9HX30l1uLittdbXn6VgUcfxe7PPmf7629mI3sZ17K7mj2LFjP0jNMwrNZsZycv5GRk\nqVm0mObt7RN3Nbe08MA7r7CtOvqj7a4mDx5vHDWrHrh56uPcN/PliOXunTvZ8vKrbJ4SfC8QCESs\ns6ubuVJFfPwtLWx5ZSp1cz5OWZqLti/nppl/ZE9zHS3V1QS83g7v+5qaWfvE32jaui1i280vTmHD\n359te7116mssu/U3cQf8PYsWd3kzNxAIsHLblqjHVDIaN26K7wokASv+9w7WP/0M8y64hIDPl9K0\nRXQ5V9P3Nzez+t4HOiybq5exdN9drH1vHf+68o8ANKxbj2ma9Bu3L5f/8T3KShy8cNd3sFut/O71\nf2Bg8OAFP+lyX3NWraCmwcVFRx3Xtuy1T5fw9op53HnBBew7bDimaVL54WzKNyzh69GFeP2XMu2L\nTxkTls7Wqa8BEDj1HH435x76+/pzzarVbe8/8OHz/Oj4s3D7WtinfBSDSwa2vefxebBZbVhS0Bzk\na2qmas5H+CccwJBBFRTaCnqcZjQvLpnOtvod/Pr4n1Fgc3R4b97ab9hRW8OFE46L2M40TfwBPzZr\n9MPWbA1uhtEWDC2WYLks2PoVzieeY6gOjmnUcOwxlI7dB9M024a62Fm7B5vVSmmhg6raJjZsaWT8\n2IGUFDooLrTj9XuxWqxUu1yYJgwpK+Ohz4I9b+Yun8WQB4In7sOf/ztFg4ewp6meyn+/TPWs2VTN\nmk3Jgfv3sOQiffzq2xx93smYDjsGBsWOIgKBAM61K3lz7dd80vgJB5ccxx3nXEmjtwnt3Mn+w0fw\nzKLJnDVuImPKR1JkK2Tn+x+w+cUpDPrbk9j79+twPO1ZspRVk+4FwHHcMRz0w2vw1dbST+1HwAwk\nfewt+vILWiqr2l77m5qxlPXrWYEAPr+PBduWcETFwRTZC+PaZsHWr6h113PWuIk93j9AQ0sjJY7i\nqMOoLHauYMyAkQwqLk/JvhJl9KZuX4ZhmMnkZ2udk0nv/okm04PNZ3L9f3ZHrDPjpP5sHFnA/hua\nOXNBezv6jgd+zLSv3457X/tvbOaor5v4aEI/Tl7kwu43mXZGOe7iAgKGvy3wOLwmHoeFEZUeLpod\nnN562bgiPp4QPKhvebkqIu03ThnAoFofg2t9HLjR3bb867GFfHxkP/zW4AFUXtSf+077DV9sW8qL\nS3cpPsAAACAASURBVKcDYDUsqMFjuf2EG7DUNzFl9dusdG2mtrmOACbNPjcFVgctfk9bukXWYjyB\nFvxmsPY2cXEDh+gmVowtZM7RZcF0/SaHbPLxzb4leC0mA6ua8FoNGgcVY3c106/RzxGrm9jTz8a8\nQ0r4/vhzWFn1DetqNlHc4KXFbjCw3o8lYDKiKlgL3r6XHecQe8QcsDafyYSVjfisBn4LbK4ooHpA\ne4Af7WxhjNPD/G+X4HEEA40RMLEGYOIiF+M3BMusqtzG4FofFhO2DrXz3vH9sftMrnynGnuooqpH\nF6A2t6BHFzDr6DIKPAFsfhOf1WD0Dg8bRhZQXu9nZGWwvOw+kyX7F9NcaKGkyY/fYuAubA92R61o\n5NgV2RsGe86EflQOsnHQumYOXtd+7GwZaidgMdhVbmPDyAJsfpP6Eitqk5uh1V6aiiz4LQaHrmkG\nYMW+hdSXWNFjCmkqtHDM8kaOXB27B1l9iYXXTx2A32rQXGChwBPguGWNWAIms48uw281sPlMDtzU\nQlGzH+cQO7X9rOy32c0JSzuW1+unDmDrsI6VAEyTnxx2CY0+Hws2rGZ77SoACkr74fJ0X97F9iK+\nt/8Z7GqsYd7WRTR73VgNK2VeK94iG/0KyjhgwLdYsfwTdg20R+wbYOy2Fk5e1MCqsYUsOWQAHr8X\ni0nb7zGc3RqsHIyo9FCx28vCA4vBMLBbbJiAL9B+NWMETAwTAqF0XrzgUQrjPEl1ZgQrOnEN1JXy\noK+UehQ4huBtplu01gvj3TbZoP/Tt37LNf9cB8DME8s499P4m0PeOaGMo1c0Uugx+eKgEk5b6KK5\nwOClswcywOXH5jdxFVs5cUkDXxxUwqUf7olIo6rcxsYRDgKG0eGH/8qZ5fzPR7UUedo/05fjizlq\nZc+6YVaW29h41nj4eh0BAxYfWEx/l59xW1qiBp7VYwpYqoo5VDdR1hhg5dhCvrXdw65yG4sOLGZQ\nnY/L39tD+BHz0ZGlnLKoIWYe1o90MHabJ2L52yf1p2KXlyO6CBStWmwGb546gJGVHorcAQ7XzRHr\nvHFKfxqLrIxf18xhazq+3xq4s8ljM3D4ek/Fqbd577gyzp4X/+/xvePKGFnlYb/NLdT2szK0Jhgk\nnz9/EEd93cjB64MntJoyK1/tX8y+W1toLjD47LBSrnsj2Hy7Ymxh23qt1owqYO2oAvbd0oLa0n7M\nrNu7gH23Bl+/dXJ/No0oYN8tbs75LHqeF+9fxBHfNLetbxpQ4DFpKLYwfr2byoE2DlrvZkhtMN+T\nvzuQgAUKW0y89mBlprafFZsfbpwWHEH300NLWDW2iP858n+44MCz4y6rcFkL+kqpk4HfaK3PVUod\nAPxTa31svNsnG/Tvu/eaLgOUEEL0dnPPP5D/vfbepLZNJOin+kbuacCbAFrr1UC5UqosxfuI4LHL\n8MNCiL6tdM2mjOwn1TdyhwGLw17vCi2Leq2klJoE3BW+zOl0Rlu1S9bxx/L4t5aCaVLsDuCzGZiA\n12Z0bDc2TSwmBCxG22uAAq9JwIDjlzZSV2pl+X5FlNcHL88qdnkJWAw2VgTbGofv9uIqtjK41kd9\nqRWv1WBoTXCd0iY/lgDsLrexY7AdtcmN2uRmzlH92DnIzrgtLUxc7OLL8SUs26+IQKi90zBNvHYL\nmME2vhJ3ABNoLLaCaVLSHAh9FrD5AEyKWoLbDXAFG6kLPCbV/W3U9A/m7bhlwc9iCZg0FFvZZ3sL\nlQNtVA+w4bcY+K0GdSUWHD6TsgY/7gILDcVWmgsMrIFgG3aLw8BvMShsCeC3GlT3t7JXjQ+7z6S0\nOUB1fyumYVDW4CdggXKXn6YCC65iC3Y/HLKmiUG1PpbtV8zm4Q6aCi1ggCVg4rca2L0mBmCYUOwO\n4Cq2UFdqpcQdIGCAxYS9anw4vAGsAdjTL/jZHN5gu7TNH/y+S5oD+KwGTYUWqgbaKPCY7BoYbEMN\n3lsx8NgM7D4T0wi2xZrAAJefvfb4qCmz4ioJ6y4YdrHpt0JZQyC0HZQ2BShqCbCr3IbNBwYmDm9w\ngxa7BY8j2LZd1BKgyB2grtSKu8CCETrWCj0m5XU+MAzqSi147BYsARO3w0K/Jj/WQPC4dXhNHN4A\nxe4Ao3Z6KGoxcQ62s2Z0AR67gdUPhZ4AhS0mpgUaiyzYvSZDa7z4rAbOIXb8FiP0Gb34LcF8BSxQ\n099GgScQ/C0YwTxZ/SY1/W24HQYFXpOh1T4O/6YJ5xA7a0cV0GI3sASC34k1YGLxB8sjYAmWZbE7\nQIHXxO0w8NgNyuv9lDX68VkN9pRZsfrBNMBdYDCo1s8RqxtxDnGwcYQDtbmFhiILFhP23epm1wA7\nphFsvttdbqewJYDFNClymxy4oZkCr8nOQfa242K/UHPNmlEFbB3qYGC9D7fDwsHrmqkrtbJUFVHu\n8rP/Rjd2r0lzoYVNFQ4G1vkZssfL/ENKKWv08+01zRS2BO/JvX9cGdX9bfRr9AfvS/hNyuv9DK71\n0VxgobglgN+AoTU+mgotlNf78FkNLCa02A0aioPH55gdwSZQt93AMKG5MPg9NxZa2DHYzqqxRcGy\nNU2aB41OKv4lKtXNO88C72it3wq9/gy4Vmu9Jq7MJNm8A7Bt+zaGDhuKaZr4TD8Ww4LFsFDf5Kas\nqICAGcButdPi8+DzWQgYLTS5/dgsNjChX0kB9lBPD5/pxzRNrIaFnXsaMIwA/UuLKLQ7ME0Tt9eD\nz2fSv7gIl7sZDAOP14/DasVZU0eBzcGAfgWUFhUQCASwhE48VsPKnqZmSgvtmKH/OawOAgE/GAYW\nDHymHwODAquDek8DFsNCsa2QFp+HNVtrKCiwUlpQSG1TE36/Sf+SQoYPKsXvN1m7fTd7Gty0eD24\nqWffoRU01zZR6XHjN9z4TR/7DNqbYf0HMLC0mAa3m4Wb1lBgFDO8/yAa3V72GV7O5p112KxWahpc\nDOpXisNuZ9TQflgNCxgGW3fVsm7nTkYPGoLDZsc0TPYePADDMKmqbWDrrjr2Ki+l1l3HARV7Y7da\n8Pp9OKx23L4WrBYrFgxaPAGaWrwUF1mpqW9mrwElWK1WAgE/hmFgGBb8AT8Oqx3//7d353FylHXi\nxz9V3XPP5JjckwPCkS8SFlRAcAMYFFZQ+KEGBAkoxyryAxR5yS6/VTGgu7r4wptFEETlWFFQNxEE\nDeAKhCPcV/hCDghJ55hMMvdMX1W/P6pmMleYTqcnM931fb9eeaX7qaeqn/pO17eqn3q6H98L/x6t\n1NaUUV1WQdb3cIC4G6etM4Xn+4yvqegd4ZPJZoi5MTw/OImmUh5VFXE2JhJMmzaV/31+AwsOm0VZ\n3MX3PeJuvOdj8qBRF+lsBgeIuTEcxyGVTRN3YzjhnRDHcchkBww5dBxcx+kdMhlzY6Szwc1s13GJ\nucGJxsfHddze9X384Mt54bEQj8VJZ9M4jks8XKenruM4ZH2PMjdO2gva6Pk+ZbF47+u6jovn7xy2\nGY/Fe0c39exnJpth8+YtNDTMwPd8dEMT82bV44VnwJ59dRyHrJft3abjODiEZa6LFx43Pj7pbKb3\nPYMPG7Zvo6G+HnyfslgZjuOwo62L5vZuZk+toyOVwvFjeL7PuOpy8KEzmWZ7Rzvd6RRzJk9m/bZG\nKsvLaJhQz6aWJuoq6vB9n47ubqaMr6W5s5MyN05TeweT6mqoKi9jU8sOulMpXNenYeJk4m6c9u5u\nyuMxJlTXkvGypDIpXMehPF6O53mkveC9szGxkdkNs/Dwg2HdbpaYE/wNytw49MYjeL86QHcmGfx9\nHZeUl6bMLaM8tvMmcWe6CwcnGIXm+2R9j1Qmy/iqmt3IeAPfaqPXp78E2KSqN4XP1wKHqerQXzsd\n2Jg9SPqJRIKGhoa81i1lFpfBLCZDs7gMViwxGc0+/b8ApwOIyPuBRK4J3xhjzMgraNJX1RXAsyKy\nAvgxcEkht2+MMWbPFPwbuap6VaG3aYwxpjDG3M8wDPW1ZWOMMYUxpn6GYU+IiK+qdsYYwOIymMVk\naBaXwUoxJiX5K5vGGGOGZknfGGMipJSS/jWj3YAxyuIymMVkaBaXwUouJiXTp2+MMWZ4pXSlb4wx\nZhiW9I0xJkIs6RtjTIRY0jfGmAixpG+MMRFiSd8YYyJkzP32Tj72ZDL2YiUi1wHHEvwNvwOsBG4H\nYsAm4FxVTYrIYuBywANuVtVbRaQM+CWwD5AFzlfVtXt/LwpPRKqAV4BvAQ9hMSHc338BMsDVwEtE\nOC4iUgv8GpgIVBCMxd8M3EiQQ15S1YvDulcCZ4Tl16jq/SIyHrgLGA+0A2er6va9viN5Kvor/XAy\n9gPDCdgvJPhJ55ImIscDh4T7fBLwQ+Ba4AZVPRZYDVwgIjUEB/kJwELgKyJSD5wNNKvqMcC/E5w0\nSsXXgZ4DMPIxEZFJBFOSHgOcApyGxeU8QFX1eIL5P35EcAx9WVUXAONF5GQRmQucxc7YfV9EYgQn\nxr+FMfk98K+jsA95K/qkzyhNxj7K/k5w9QHQDNQQHKhLw7JlBAfvUcBKVW1R1S7gcWABQcz+ENZd\nHpYVPRE5CDgYuC8sWkjEY0Kwz8tVtU1VN6nqF7C4bAMmhY8nElwkzO3TQ9ATk+OBP6tqSlUbgbcJ\n3l99Y9JTt2iUQtKfTjABe4+eydhLlqpmVbUjfHohcD9Qo6rJsGwrMIPBsRlUrqoe4ItI+d5o+wi7\nHriiz3OLCewLVIvIUhF5VEQ+QsTjoqq/AeaIyGqCC6ivAjv6VMk5Jn3KikYpJP2BSupnUN+NiJxG\nkPQvHbBoVzHY3fKiISKfBZ5Q1XW7qBK5mIQcgqvaTxF0a9xG/32LXFxE5BxgvaoeAHwYuGNAld3Z\n96KLRykk/QT9r+wbCG5OlTQR+SjwNeBkVW0B2sObmAAzCeIyMDaDysMbdY6qpvZW20fIx4HTRORJ\n4J+Bb2AxAdgCrFDVjKquAdqAtojHZQHwIICqvghUAZP7LM85Jn3KikYpJP3ITcYejh74HnBKn1ED\ny4FF4eNFwAPAU8CRIjIhHLGwAHiUIGY99wROBR7ZW20fKap6pqoeqapHA7cQjN6JdExCfwE+LCJu\neFO3FovLaoJ7GIjIPgQnwlUicky4/FMEMXkY+LiIlItIA0GCf43+MemJX9EoiV/ZFJHvAscRDDW7\nJDx7lywR+QKwBHijT/HnCJJdJcENp/NVNS0ipwNXEgw5+4mq3hmOQLgFOBBIAuep6jt7cRdGlIgs\nAd4iuJr7NRGPiYhcRNANCPBtguG9kY1LeFL7BTCNYMjzNwiGbN5EcCH8lKpeEda9DFhMEJOvq+pD\n4fp3EHSbNQPnhJ+2i0JJJH1jjDG5KYXuHWOMMTmypG+MMRFiSd8YYyLEkr4xxkSIJX1jjImQMfUr\nm47j2FAiY4zJg+/7OX07eEwlfYB8h5AmEgkaGhoK3JriZ3EZzGIyNIvLYMUSE8fJ/dcgrHvHGGMi\nxJK+McZEiCV9Y4qAl8mQ6ewc7WaYEmBJ35gi8PwlX+Kpz5w72s0wJcCSvjFFoHvzltFugikRlvSN\nMSZCLOkbY0yEWNI3xpgIsaRvjDERYknfGGMixJK+McZEiCV9Y4yJEEv6xhgTITn9yqaIXAccG9b/\nDrASuB2IAZuAc1U1KSKLgcsBD7hZVW8dkVYbY4zJy7BX+iJyPHCIqn4QOAn4IXAtcIOqHgusBi4Q\nkRrgauAEYCHwFRGpH6mGG2OM2X25dO/8HTgjfNwM1BAk9aVh2TKCRH8UsFJVW1S1C3gcWFDQ1hpj\njNkjw3bvqGoW6AifXgjcD3xUVZNh2VZgBjAdaOyzak+5McaYMSLnmbNE5DSCpP9PwJt9Fu1qypZh\np3IRkSXAN/uWJRKJXJvUT1tbW97rljKLy2DFHJORbHcxx2WklGJMcr2R+1Hga8BJqtoiIu0iUhV2\n48wEEuG/6X1Wmwk8+W7bVdUlwJKe547j+PlOTVYs05rtbRaXwYoxJuvC/0ey3cUYl5FWijHJ5Ubu\neOB7wCmquj0sXg4sCh8vAh4AngKOFJEJIlJL0J//aOGbbIwxJl+5XOmfCUwGfisiPWWfA24RkYuA\nt4FfqWpaRK4CHgR84BpVbRmBNhsTWb7v79Yk2MYMlMuN3JuBm4dYdOIQde8B7ilAu4wxxowA+0au\nMcZEiCV9Y4yJEEv6xhgTIZb0jTEmQizpG2NMhFjSN8aYCLGkb0wx8f3RboEpcpb0jTEmQizpG2NM\nhFjSN8aYCLGkb4wxEWJJ3xhjIsSSvjHGRIglfWOMiRBL+sYYEyGW9I0pJvblLLOHLOkbY0yE5Dox\n+iHA/wA/UNWfisgvgcOBprDK91T1PhFZDFwOeMDNqnrrCLTZGGNMnoZN+iJSA/wEeGjAov+nqn8a\nUO9q4ANAClgpIn/oM5m6McaYUZZL904S+BiQGKbeUcBKVW1R1S7gcWDBHrbPGGNMAeUyMXoGyIjI\nwEWXisgVwFbgUmA60Nhn+VZgxrttW0SWAN/sW5ZIDHduGVpbW1ve65Yyi8tgxRyTxKZNOO7I3Ior\n5riMlFKMSU59+kO4HWhS1RdE5CpgCbBiQB1nuI2o6pJw3WAFx/EbGhryalAikSDfdUuZxWWwYozJ\nuvD/hhkzcGKxEXmNYozLSCvFmOSV9FW1b//+UuBG4B6Cq/0eM4En82+aMcaYQsvrc6KI3Csi+4VP\nFwKvAE8BR4rIBBGpJejPf7QgrTTGGFMQuYzeORy4HtgXSIvI6QSjee4WkU6gHThfVbvCrp4HAR+4\nRlVbRqzlxhhjdlsuN3KfJbiaH+jeIereQ9DNY4wxZgyyb+QaY0yEWNI3xpgIsaRvjDERYknfGGMi\nxJK+McZEiCV9Y4qIb7+nb/aQJX1jjIkQS/rGGBMhlvSNMSZCLOkbY0yEWNI3xpgIsaRvjDERYknf\nGGMixJK+McZEiCV9Y4qJfTnL7CFL+sYYEyE5zZErIocA/wP8QFV/KiKzCSZHjwGbgHNVNSkii4HL\nAQ+4WVVvHaF2G2OMycOwV/oiUkMwPWLfydCvBW5Q1WOB1cAFYb2rgRMIZtr6iojUF7zFxhhj8pZL\n904S+BiQ6FO2EFgaPl5GkOiPAlaqaouqdgGPE0yObowxZowYNumraiZM4n3VqGoyfLwVmAFMBxr7\n1OkpN8YUyPNfumK0m2CKXE59+sNwdrO8l4gsAb7ZtyyRSAxdeRhtbW15r1vKLC6DFXNMuhOJEWt7\nMcdlpJRiTPJN+u0iUhV+AphJ0PWTILja7zETePLdNqKqS4AlPc8dx/EbGhryalAikSDfdUuZxWWw\nYolJpr0dvf6HzDnr0/3KR6rtxRKXvakUY5Jv0l8OLALuCP9/AHgKuEVEJgAZgv78ywvRSGOiKPGn\n+2l+7nlaX3l1tJtiSsiwSV9EDgeuB/YF0iJyOrAY+KWIXAS8DfxKVdMichXwIOAD16hqy4i13JgS\n52ezAHiZzCi3xJSSYZO+qj5LMFpnoBOHqHsPcM+eN8sY08u+hWsKyL6Ra4wxEWJJ3xhjIsSSvjHG\nRIglfWPGOuvTNwVkSd8YYyLEkr4xY50z7JfbjcmZJX1jxjrr3jEFZEnfGGMixJK+McZEiCV9Y8Yo\nx/ryzQiwpG+MMRFiSd8YYyLEkr4xxkSIJX1jjIkQS/rGjFV2I9eMAEv6xhSZTX9+kFRz82g3wxQp\nS/rGFJm1P7uZVd/6j9FuhilSec2RKyILgd8BPZN3vgxcB9wOxIBNwLmqmixAG40xA7SvXjPaTTBF\nak+u9P9XVReG/y4DrgVuUNVjgdXABQVpoTFmSF2bNo12E0wRKmT3zkJgafh4GXBCAbdtjBng7dvv\nHO0mmCKUV/dO6GARWQrUA9cANX26c7YCM/a0ccYYYwor36T/JkGi/y2wH/DIgG3lNNZMRJYA3+xb\nlkgk8mpQW1tb3uuWMovLYMUSk7a2tndd3tXVVdD9KJa47E2lGJO8kr6qbgTuDp+uEZHNwJEiUqWq\nXcBMYNhIqeoSYEnPc8dx/IaGhnyaRCKRIN91S5nFZbBiiUlm3DjebWBmVWVVQfejWOKyN5ViTPLq\n0xeRxSLy1fDxdGAacBuwKKyyCHigIC00xhhTMPl27ywF7hKR04By4GLgeeDXInIR8Dbwq8I00Rgz\nJJtRy+Qh3+6dNuDUIRaduGfNMcYYM5LsG7nGFC270je7z5L+GOal06PdBDOGWe+OyYcl/TFq22OP\n88TpZ9H0xFOj3RQzCnzf553/vnv4isbsJkv6Y1Ri2X0AbPqzDYKKos717wxfyS71TR4s6RszBvnZ\nzGg3wZQoS/rGjEV2EW9GiCX9scoOejMse5OY3WdJf4xzbMq8aLL+ejNCLOkbM8qy3d2Dhuf6uSR9\nOy+YPFjSH+NyOvhNUXvyzMU8/bkBcw7l8He394bJhyX9UdK1aRNblj+E73m7qBEc0C0vvrT3GmVG\nTbajs3+BJXQzQvZkEhWzB5774qUApJq2M/vMMwYtt6s4Mzx7j5jdV9JX+qntO9iy/OExnUAbH31s\nyHK7gWuMGQklfaX/ytVL6HpnA/GaaiZ98OjRbs6Q/Owuuncs6ZekzQ/8BbeykmRjI1OPX4gT23nd\n1bUxQeWM6WTaO2j8+6PDb2zsXsvsFt/37SJnLyrppN/1zgYAktu2kc5mKYvFRq0tj+kqJtbWMn/m\n7H7l3YkEyW1N+Jk0ldOn95Y77s5k4HkerlvSH8r2Gt/z8DOF+7Zre6qDl7e8zlGz3ofrDP036t66\nlYrJk3FclzU33tRbvv6Ou/rVe+7/XsbkYxaw7bHHc3rtQu7HaEks/RPrbr2N2Wd9mjmfOXO0mxMJ\nJZdJ/GyWZz5/Metu2zmHy+sbN7D4nku58ZH7hlzn7uXKoy9s7H3+zNo1PLN2TU6v5w24EdvU0sWq\nddv7t8n3+fELP+aax/5jyG08c+EXePaiSwC4/fFHeOCl52h9bdXO9j09dBfQSNgbXWHrmzfy1o4c\nflumgDrTXax79glWfPIM3rrsK6Rb333+2Vx977Gb+MGKW3ji7WdZ87Of0/zCi/2WN7/wIs9+/mLe\n+tXttK56fdjt5Zrwe7bd1bhtl8vfatzKxf99Pa9uWD942Y4NpDKp3ufJpu20P71y1+16/Ame/9IV\ntL75Zs7ty8W6W28D4J3f/JZMZ+cwtUdOKpumM901aq+/Nzljqb/bcRw/n/Z4vsdzq1+kvK6SqpTP\n1su+0W/5W7MnsGxBGZ7r8NH9j2Na2Xga7nsOP51BrrqSxd+5G6eyk0Pml/PF407ligeuBeC6f/o3\nPN+nO5NiXEUNnekuUi3NNCc2MPOQ9/HzX15PU00HnzjmfGbUTWXahHFcevuNxCZtZmrZPlyx8DNk\nW9vZ+uWrAbjz5HpOPuZcbn/5br78m82D9uP3nzyNdypXMHNrmtMf2jk76p0nT8SfOZW2ZDtHzDyM\nRQefzKxxM2jubmWZLmfWuOkcNOUAytw4U2omAdDc3UplvII176xl0pTJvdvyfZ90Nt37ySHmxPDx\nyWQzJFe9ydbv/oSXFx3GjCOP5rDpBxOPxXpPBDHHJeNlcR2XrJ/FwSHmxuh4TUnXVFA3ezaTq+tp\n7m7Bw6cz1UVlvALXcdnR3kUqk2HK+Bq++uC3Afjse09n/tR5bO9qJp1NE3Nj3PDX/8Lxfc754Nlk\nHJg7aQ6d6W4832Nj62bauls5eJrg4+OG7Xl5y+vUZBzKn36dqvIqNs7el1ZvC97Klzji1DNJjo9z\n8zN3cdHvGqlMB/vS+vlTGdeRoWNWPQfNP4rOdBdXP3w9AIdNO5jmlgzr3ulixpRq5tbPZO7MOja0\nbmJqzWTufe1+AD489x95eN0KAKY2pfnMgzsA6KyvZvOpR/LG9gQnLStskhyoudbliUNr2XHwLCZ6\nc/jAvvuQeHoF+y9/lRfmVbFyfjXxLHx63knctfZBAA4YfwCrW1ZT15Fl3tQDaHaynHjXy1S0J3nh\n9EPJHjCbI2YeysSKel55/QUavn9n7+t1lzkkP3I44/76HFuOPoBHpyVprY7z2fcvYtaEabSnOmhs\nTNCShPftN494LIaDg+u4rG/ZyDMbXqVh/GTqY3Ws/+nPOPCdZO+2Z93wHeLjxoEPLySUra0txPwK\nzjj8Q2S9LMlsimQmScbLEo/Fae5qxXGgpqwagG2dO/B8j+5Mkp8/cxeLD/sk/zDtIGrKq8l4GTzP\nI+NlyPoejhO0Kea4uI7Lxcv+DYCz/uH/UF1WhUzen5buNhqbGtl3+hw63lxN/X77Q0UFDg5N7e3M\nmjiZTDaLh8fm9kZcx2FG3VQ83+e2Z3/DoRMP5ODZ86ksqyDmuGR9j+auFv7zsRtZMPsIjtv3KLoy\n3ZTHyilz48ybvF/e7wPHcfB9P6c+soInfRH5AXA0QY/jl1V115cPAxuTZ9L/xbN3s2XZ/fgOrJpb\nyRfvHfrq56n51Rz1av+riR+dPRV8n/KMT9ZxOO65dp4/qIrmcf17vmJZn6wLF927jcqUj+eAGzb1\nkSNqeWleNfg+h6zpJlnmMGdzir8dUceRr3Zw1Cs7X/O/zpjMvLeTnPD07l1prmsoJzGljNf2q6Kz\nyuXsQz/BXS/9cVC9r9YdT/cf/8q2dCuPv7cW1/PpqnDZOqmMWNanLO3TXekyqTlDS22MTDx4n8ze\nnOLEJ1qp6/LYNCnO/ceO5/DXOmmvdjnmhQ7+elQdWdfhkDVdNE6M89j7alm4so2p2zNM2xF0M9x9\n4kQ2TymjpjPLzK1p5q/t4u3p5eyzOUVtp0d9a7a3natnlfPivGo2TC8nlvHJxmC/jSlO/XtLb52X\nD6jk4Q+MAyCe8bnkt40AdFY4rHhvLa/uV8nhqzo5eG13v20P1DQ+hg9Mbhm6zpuzK5jQlqW626Om\nO/jktnZmOfttTPWrlyxzeHX/SuZsSrF1Uhmvza2ktSZGOu5w0e93fcU91nRUur37uSut1S4rxEe+\nAwAACMRJREFU3lvLSStac9rm5klxOqpc9t+wM2YrDqth5fwajn2ujfe/vvMqel1DOXMTqUHbWD2r\nnPuOm0As4zO5JXh/zt6cojLl8fKB1bhZn0Pf7OLgtd08+I/jaJoQpyztUdvpsWN8nIktGdqrXVwf\n3rO2m631cSa0Zemocnl7Rjk4DlXdHrWdWQ5f1cmOcXFWzq/moHXdvO/1Tn534kRS5S5laY/pTRkO\nWtdNW7XbL2fc8bF6JrRlqEj5vLZ/Vf8dCHPXuA6P85c2AfCzRZNJVgzdofL+VZ0csL6be0+YSDbm\n8PUPfYlDp78np3gPNGpJX0Q+BFypqqeIyHuAX6jqB3NdP9+k//lbPs95y7YPX9EYY3aho9Ll2fdU\nc9zz7QXb5htzKpjYmmVCW4at9WXMbExz94kTmb0lxT++1NFbr63aJfGlT3DBBxfn9Tq7k/QL3af/\nEeCPAKq6CpgoIuMK/BqDHPZGNPrijDEjp6bbK2jCB5i3PsmU5gxlWZjZGPzUxpl/3dEv4QPUdXq0\nLNs7c2cUOulPBxr7PG8My0bUyvk16D4VtNYM3p3m2hirZ5Wzoy4YuZNxoaXGpb1q9O9hr55V0fu4\nubb/yKIt9XEe+kAdmyaV1gCr4T7HpWPQXuWSjDsky/ZsGF9nRbB+xg22mXWC/zsqXdqq+//9M7Gg\nXtYN2rAnOip3bjsdg41TymgaF6NpXIzsgF3qKndoqXFpGeK9O5RU3KFxQpzWAe1PxvtvOONCV8Xu\nxS8zeoPbDLBuZsXwlQpgpDPKu77rRGQJ8M2+ZYlEYrdf5OvHXcXXKq/b7fUGquvIkok5dPU5aHu/\nDp/LOOIBdR3PpyLlkyx38F0HfB/HJ3ico1cOqBq+0hDqmzO017g4HqTKHCpSPpm409uHPxLiGR/f\ngWwseI2KpEdV0ht0f2RvihFnRtUMWlPttGZ35LTOxLJ6dqS3U5GaxLTqyazP6M6FvsO8qkOpq/N5\nZdsbJP3u3kX1ZZPZnt5GmVPOzOoZ1JZX0ZHqZkeqmYmV49jSuY26sjp8D7amtvSuV+3WcvK+x7Jy\n+zq6vS58PLJZSHpJplXXM7eygRd2KIdNm8+4slqaupp4ZsvLzKmdRWcqzZqON6l0q2jPtuESo6Fm\nKnXltTR17WBr9zZq3HHEXEh6KaZVTWFTRyPlsTgHjt+H2tgEDp/+Hp7YspKXt71Bl7fzU3NFyiMT\nc3r/ngA1fj0dznYmuFNo9hpxcPHpf3+gilq6/PZdH/2+T3W3RybukCpzqUh6pMocyrI+dR0eO8YF\n92CGOk4qkh7Jcqff8Th1exofaJoQx+tZx/ep6wza1VYTnM1iWZ94xidVtvN4nNiapaUuhuc6lKc8\nKtI+HZUuXmyY48T3wXFwvOCYr0z5pOMOrueTKt+ZPxzPxwE8B8rTPlnX6Y1LdsBrnLf/OXnlv91V\n6D79JcAmVb0pfL4WOExVc7prmW+fPgQni4aGhrzWLWUWl8EsJkOzuAxWLDEZzT79vwCnA4jI+4FE\nrgnfGGPMyCto0lfVFcCzIrIC+DFwSSG3b4wxZs8UvLNVVa8q9DaNMcYUxugPYTHGGLPXjLnxgPZr\ne8YYM3LG1G/v7AkR8VXVzhgDWFwGs5gMzeIyWCnGxLp3jDEmQizpG2NMhJRS0r9mtBswRllcBrOY\nDM3iMljJxaRk+vSNMcYMr5Su9I0xxgzDkr4xxkSIJX1jjIkQS/rGGBMhlvSNMSZCxtzPMORjTyZj\nL1Yich1wLMHf8DvASuB2IAZsAs5V1aSILAYuBzzgZlW9VUTKgF8C+wBZ4HxVXbv396LwRKQKeAX4\nFvAQFhPC/f0XIANcDbxEhOMiIrXAr4GJQAXBsMzNwI0EOeQlVb04rHslcEZYfo2q3i8i44G7gPFA\nO3C2qhbNJN1Ff6UfTsZ+YDgB+4UEP+lc0kTkeOCQcJ9PAn4IXAvcoKrHAquBC0SkhuAgPwFYCHxF\nROqBs4FmVT0G+HeCk0ap+DrQcwBGPiYiMolgdrpjgFOA07C4nAeoqh5PMP/HjwiOoS+r6gJgvIic\nLCJzgbPYGbvvi0iM4MT4tzAmvwf+dRT2IW9Fn/QZpcnYR9nfCa4+AJqBGoIDdWlYtozg4D0KWKmq\nLaraBTwOLCCI2R/CusvDsqInIgcBBwP3hUULiXhMCPZ5uaq2qeomVf0CFpdtwKTw8USCi4S5fXoI\nemJyPPBnVU2paiPwNsH7q29MeuoWjVJI+qMyGftoUtWsqnaETy8E7gdqVDUZlm0FZjA4NoPKVdUD\nfBEp3xttH2HXA1f0eW4xgX2BahFZKiKPishHiHhcVPU3wBwRWU1wAfVVoO8EyjnHpE9Z0SiFpD9Q\nSf0i3rsRkdMIkv6lAxbtKga7W140ROSzwBOqum4XVSIXk5BDcFX7KYJujdvov2+Ri4uInAOsV9UD\ngA8Ddwyosjv7XnTxKIWkn6D/lX0Dwc2pkiYiHwW+Bpysqi1Ae3gTE2AmQVwGxmZQeXijzlHV1N5q\n+wj5OHCaiDwJ/DPwDSwmAFuAFaqaUdU1QBvQFvG4LAAeBFDVF4EqYHKf5TnHpE9Z0SiFpB+5ydjD\n0QPfA07pM2pgObAofLwIeAB4CjhSRCaEIxYWAI8SxKznnsCpwCN7q+0jRVXPVNUjVfVo4BaC0TuR\njknoL8CHRcQNb+rWYnFZTXAPAxHZh+BEuEpEjgmXf4ogJg8DHxeRchFpIEjwr9E/Jj3xKxol8YNr\nIvJd4DiCoWaXhGfvkiUiXwCWAG/0Kf4cQbKrJLjhdL6qpkXkdOBKgiFnP1HVO8MRCLcABwJJ4DxV\nfWcv7sKIEpElwFsEV3O/JuIxEZGLCLoBAb5NMLw3snEJT2q/AKYRDHn+BsGQzZsILoSfUtUrwrqX\nAYsJYvJ1VX0oXP8Ogm6zZuCc8NN2USiJpG+MMSY3pdC9Y4wxJkeW9I0xJkIs6RtjTIRY0jfGmAix\npG+MMRFiSd8YYyLEkr4xxkSIJX1jjImQ/w+6Bw0CIkBtFwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f21ff025910>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize = (6,8))\n",
    "ax1 = fig.add_subplot(311)\n",
    "ax1.plot(np.arange(0,len(load)), load)\n",
    "ax1.plot(np.arange(0,len(load)), newDemand1)\n",
    "ax2 = fig.add_subplot(312)\n",
    "ax2.plot(np.arange(0,len(load)), load)\n",
    "ax2.plot(np.arange(0,len(load)), newDemand2)\n",
    "ax2.plot(np.arange(0,len(load)), price)\n",
    "ax3 = fig.add_subplot(313)\n",
    "ax3.plot(np.arange(0,len(load)), load)\n",
    "ax3.plot(np.arange(0,len(load)), netLoad)\n",
    "ax3.plot(np.arange(0,len(load)), price)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Remember that bills can be negative if lots of solar produced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Min bill calculated cycle charge 40287.7584251\n",
      "Min bill calculated paper algorithm 30212.2645252\n",
      "Min bill calculated MILP 30212.2841169\n"
     ]
    }
   ],
   "source": [
    "if billingOption[0] == 1:\n",
    "    print 'Min bill calculated cycle charge', sum([newDemand1[x]*price[x] for x in range(len(price)) if newDemand1[x]>0]) + \\\n",
    "            sum([newDemand1[x]*billingOption[1] for x in range(len(price)) if newDemand1[x]<0])\n",
    "    print 'Min bill calculated paper algorithm', sum([newDemand2[x]*price[x] for x in range(len(price)) if newDemand2[x]>0]) + \\\n",
    "            sum([newDemand2[x]*billingOption[1] for x in range(len(price)) if newDemand2[x]<0])\n",
    "    print 'Min bill calculated MILP', sum([netLoad[x]*price[x] for x in range(len(price)) if netLoad[x]>0]) + \\\n",
    "            sum([netLoad[x]*billingOption[1] for x in range(len(price)) if netLoad[x]<0])\n",
    "else:\n",
    "    print 'Min bill calculated cycle charge', sum([newDemand1[x]*price[x] for x in range(len(price))])\n",
    "    print 'Min bill calculated paper algorithm', sum([newDemand2[x]*price[x] for x in range(len(price))])\n",
    "    print 'Min bill calculated MILP', sum([netLoad[x]*price[x] for x in range(len(price))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  65   66  151  213  243  344  390  410  598  741 2333 4443 6438 7108 7562\n",
      " 8162 8595]\n",
      "[  64  159  218  255  351  399  579  723 2354 4446 6439 7112 7560 8181 8182\n",
      " 8614]\n"
     ]
    }
   ],
   "source": [
    "A = netLoad - newDemand2\n",
    "print np.where(A<-0.1)[0]\n",
    "print np.where(A>0.1)[0]\n",
    "# times with different schedules algorithm vs MILP - more than one solution though when lots of periods with same price"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}