Model Building.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Gurobipy commands for getting the results:\n",
    "\n",
    "    * `printAttr('X')`:     Prints the values of all decision variables having non-zero values.\n",
    "    * `objVal` or `ObjVal`: Gives the Objective values of the model.\n",
    "    * `getVars()`:          Gives the list of values of all variables.\n",
    "    * `model.getAttr(GRB.Attr.X, Name_of_varibles or list of variables)`:  Gets you the value of the variables requested."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import random as rn\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as st\n",
    "from gurobipy import *\n",
    "import time\n",
    "import numpy as np\n",
    "import seaborn as sn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (10,6)\n",
    "plt.style.use(\"fivethirtyeight\")\n",
    "plt.rcParams.update({'figure.autolayout': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Defining the Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Minimum Hourly Mean Demand=  539.0 \n",
      "Maximum Hourly Mean Demand=  109517.0 \n",
      "Median Hourly Mean Demand=  8946.0 \n",
      "Mode Hourly Mean Demand=  ModeResult(mode=array([539.]), count=array([1])) \n",
      "Average Hourly Mean Demand=  18606.0 \n",
      "Standard Hourly Mean Demand=  25879.0\n"
     ]
    }
   ],
   "source": [
    "path = \"C:/Users/garg1/Desktop/IAC Jupyter Files/demand_profile.csv\";\n",
    "\n",
    "#importing demand profile\n",
    "demand_profile = pd.read_csv(path,header=None);\n",
    "demand_profile = demand_profile.iloc[:50,:].copy()\n",
    "nr_cust,nr_t = demand_profile.shape;\n",
    "\n",
    "#demand_profile = demand_profile*2\n",
    "\n",
    "#getting demand of the customers\n",
    "demand = demand_profile.sum(axis=1).values;\n",
    "\n",
    "#getting mean hourly demand of each customer\n",
    "mean_hourly = np.rint(demand_profile.mean(axis=1).values)\n",
    "\n",
    "#printing mean dats\n",
    "print(\"\\nMinimum Hourly Mean Demand= \",round(np.min(mean_hourly),ndigits=0),\n",
    "        \"\\nMaximum Hourly Mean Demand= \",round(np.max(mean_hourly),ndigits=0),\n",
    "        \"\\nMedian Hourly Mean Demand= \",round(np.median(mean_hourly),ndigits=0),\n",
    "        \"\\nMode Hourly Mean Demand= \",st.mode(mean_hourly),\n",
    "        \"\\nAverage Hourly Mean Demand= \",round(np.mean(mean_hourly),ndigits=0),\n",
    "        \"\\nStandard Hourly Mean Demand= \",round(np.std(mean_hourly),ndigits=0))\n",
    "\n",
    "# Loading Data Related to fuel\n",
    "path = \"C:/Users/garg1/Desktop/IAC Jupyter Files/fuel_properties.csv\";\n",
    "data2 = pd.read_csv(path);\n",
    "data2['Eff- kWh/$/CO2'] = data2['kWh/Unit Fuel']*100/data2['Unit Price']/data2['CO2 lbs/Fuel Unit']\n",
    "data2['Eff- kWh/$'] = data2['kWh/Unit Fuel']*100/data2['Unit Price']\n",
    "\n",
    "#getting length of the data\n",
    "nr_fuel = len(data2.index);\n",
    "\n",
    "alpha_rate = list(data2[\"kWh/Unit Fuel\"].values) #capturing kWh/fuel source\n",
    "beta_rate = list(data2[\"CO2 lbs/Fuel Unit\"].values) #capturing CO2 lbs/fuel source\n",
    "\n",
    "fuel_supply=[]\n",
    "capacities = [.50, .60, .20]\n",
    "for i in range(nr_fuel):\n",
    "    fuel_supply.append(round(demand.sum()*(capacities[i]/alpha_rate[i]),ndigits=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* __Defining all the cost Parameters__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#other fuel properties\n",
    "cost_fuel = list(data2[\"Unit Price\"].values) #capturing cost of unit fuel source\n",
    "\n",
    "cost_ramp = 0.50; #ramping up or down cost\n",
    "\n",
    "caligraphy_n = 24; #maximum number of price change allowed\n",
    "\n",
    "#epsilon_rate = round(np.median(mean_hourly),ndigits=0); #maximum ramping allowed at no cost\n",
    "\n",
    "epsilon_rate = round(np.median(mean_hourly),ndigits=0); #maximum ramping allowed at no cost\n",
    "\n",
    "cost_gamma = [10,10,10,10,15,15,15,15,20,20,20,20,20,10,10,10,10,10,10,15,15,15,10,10]\n",
    "#cost_gamma = [5,5,5,5,10,10,10,10,15,15,15,15,15,5,5,5,5,5,5,15,15,15,10,10]\n",
    "\n",
    "\n",
    "cost_emission = 10;  #Penalty cost of CO2 emisisons\n",
    "\n",
    "#inconvenience cost to the customer\n",
    "cost_xi=pd.DataFrame(0,index=np.arange(nr_cust),columns=np.arange(nr_t));\n",
    "\n",
    "#Maximum overconsumption allowed more than psi\n",
    "omega_rate=pd.DataFrame(0,index=np.arange(nr_cust),columns=np.arange(nr_t));\n",
    "\n",
    "for i in range(nr_cust):\n",
    "    n = demand_profile.iloc[i,:].max()\n",
    "    m = demand_profile.iloc[i,:].mean()\n",
    "    for t in range(nr_t):\n",
    "        r = demand_profile.iloc[i,t]\n",
    "        cost_xi.iloc[i,t] = abs(round((1-(r/n))/5,ndigits=2))\n",
    "        omega_rate.iloc[i,t] = abs(round(r*(1-r/n)*1,ndigits=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1AAAAIwCAYAAACImIrfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde1jUZf7/8ScBoVkygjCKeMJyRDLNAzqFoNKWhy00D9RqbpSaVliJCnjMdDU13TVys00tSzL8oQVsrodtMURJ21VKU1nLNKHEPKBiiRzm94dfPusIKurYIL0e18V1wX2/577vz4eh5u19+LgUFBTYEBERERERkSu6xdkDEBERERERuVkogRIREREREakiJVAiIiIiIiJVpARKRERERESkipRAiYiIiIiIVJESKBERERERkSpSAiUiUkNs2rQJk8nErFmznD0UqcTBgwcxmUz06dPnV+lv1qxZmEwmNm3a9Kv0d7FLvR/79OmDyWTi4MGDThkX/O93MWrUKKeNQURuXkqgRET+j8lkwmQyXTam/MOfsz6UOlNiYqJxj8q//Pz8CAwM5OGHH2batGns2bPH2cOsUS6+3z4+PgQEBBASEsKoUaP4+9//TnFx8Q3pe9SoUTfte/3XTlZF5LfFzdkDEBGRm8vdd99tfDA9d+4cP/30E9nZ2fz5z3/mz3/+M4MHD2bu3LncdtttTh5pzREbGwtAaWkpp06dYt++fXz88cesWLECi8XCW2+9Rbt27exeM2LECPr374+/v78zhkyHDh3Ytm0b3t7eTun/cvz8/Ni2bRt169Z19lBE5CakBEpERK5KmzZtiI+Pr1CenZ3NqFGjSExM5NixY3z44YdOGF3NVNn9LigoYPr06SxZsoR+/frxz3/+kxYtWhj13t7eTk1ebrvtNlq2bOm0/i/H3d292o5NRKo/LeETEXGQL7/8kj/+8Y/cdddd+Pj4EBQUxHPPPceBAwcqxJYvj6psH8il9mdcuKQqMTGRsLAw/Pz8CAkJueSYunfvjpeXV6VjAHjvvfcwmUxMmzbtqq61Mu3atePjjz/G29ubtWvX8o9//KNCTH5+PnFxcbRv3x6z2UzTpk3p168fn332WYXY8iWDs2bNYseOHfTv358mTZrQpEkTnnjiCXJzcwHYv38/Tz75JC1atKBBgwb06dOHnTt3Vmjvm2++4eWXX6Zbt260aNECX19f7r77bkaPHs2hQ4cqxJfv4Rk1ahQHDx7kqaeeIiAgALPZTFhYGGvWrKn0Ppw+fZoJEybQunVrzGYznTp1IiEhAZvNdrW39LJMJhPz5s1j0KBBnDhxgpdfftmu/lJ7oDZt2kRkZCRBQUH4+vpy55130q1bNyZOnGiMsU2bNqxYsQKAhx9+2G4ZYbkrvR+vtCfPZrPxxhtv0KlTJ8xmM0FBQUycOJHTp09Xeq2XWo538XUmJibStm1bADZv3mw39vKxXG4PVH5+PuPGjaNt27b4+vrSvHlzBg0aRGZmZoXYa32PiMjNTQmUiIgDrF27lt/97nd88sknhISEEB0dzd133218sPzyyy8d1ldCQgIxMTEEBAQwfPjwyyZQw4YNo6ysjGXLllVav3TpUm655Rb++Mc/OmRsZrOZqKgoAFauXGlX9/XXX9O1a1feeustmjdvzvDhw+nTpw/bt2+nb9++vP/++5W2uWPHDnr37o2bmxtDhw4lKCiItLQ0IiIi2Lt3Lz169ODYsWM8/vjjdO3alc2bN9OvXz8KCwvt2klLS2Pp0qU0atSI/v37M2LECCwWC++//z49evQgLy+v0v4PHTpEeHg4hw4dIjIykn79+rF3714GDx5cIfErKioiIiKCv/71r5hMJkaOHElISAjz5s0jLi7uWm/rZZW3+49//KPS5ONC69ev55FHHiErK4uQkBCef/55HnnkETw9PXnrrbcoLS0FzidHd999NwCPP/44sbGxxtfFrub9eKH4+Hhee+01QkJCGDlyJJ6enixcuJC+fftSVFR0NbfATps2bRg5ciQAjRs3thv7lcZ28OBBunfvzttvv42/vz/PPfccPXv2ZNOmTTzyyCO89957lb7uat4jInLz0xI+EZGLXO4Uu++//75CWWFhIc8++yzFxcWkpKQQGhpq1L333nuMHj2akSNHsmXLFlxcXK57fJmZmaxfv5577rnnirH9+/dn8uTJLF++nPj4eG699VajLjs7m+zsbH73u9/RrFmz6x5Xua5du/Laa6/x73//2ygrLS3lj3/8IydPniQtLc3ug+zhw4cJDw9n3LhxPPTQQ/j6+tq1t379epYtW0ZERARwfuZiwIABfPrppzz44IPExcXx7LPPGvEvvPACy5Yt4/3337ebYYiMjOTZZ5/Fw8PDrv0NGzYQGRnJa6+9xp///OcK15OZmcnEiRMZN26cUTZw4ED69+9PQkICYWFhRvkbb7zB9u3b6d27N8uXL+eWW87/O+VLL71Et27druY2VllAQACNGjUiLy+P7OxsunbtesnY9957D5vNRlpamjFLU+748eO4uZ3/WPDss8+yc+dOdu3axR/+8IfLtnk178cLbd26lU2bNtG4cWMApkyZwhNPPMGaNWtYuHAhY8aMuar2yt1zzz14enqyaNEimjRpUunyx0t56aWX+OGHH4iLi7NLeJ9//nkeeOABxo0bR48ePSrsK7ua94iI3Pw0AyUicpHZs2df8quypV5r1qzh+PHjRERE2CVPAEOHDqVdu3bs2bOHbdu2OWR8Q4cOrfKH1Vq1ajFkyBB++ukn/v73v9vVLV26FMCYMXKUBg0aAHD06FGjbP369XzzzTc8/fTTFWYBGjRoQHR0NGfPniUlJaVCeyEhIUbyBODi4sLAgQOB8/t8Ll6GFRkZCVBhGZ+fn1+F5Angd7/7Ha1ateJf//pXpdfTpEmTCh/mw8PDady4Mdu3b7crT0xMxMXFhWnTphnJU3kbzzzzTKXtO0Jl97wy5WOq7IAPLy+va+r7at6PFxo5cqSRPAG4uroybdo0XFxcWL58+TWN5Xrk5eXxr3/9Cz8/vwq/76CgIJ566imKiopISkqq8NqreY+IyM1PM1AiIhcpKCi4ZF2fPn3YvHmzXVn58ryLk6dyYWFhZGdn8+WXX9K5c+frHl/Hjh2vKv6pp57ijTfeYOnSpTz66KPA+X06q1atwt/fn4ceeui6x1SZC2fbtm7dCkBubm6lM3z79+8H4L///W+Fuso+nJcnDEFBQRVm9crrfvjhB7tym83GypUr+eCDD9i1axcFBQXGkjXAbnbuQm3atMHV1bVCeaNGjeyS4tOnT7N//34aNGjAXXfdVSH+/vvvr7R9R7rSDOegQYNITU0lPDycfv360bVrVzp16kTTpk2vuc+rfT+Wq+x+3HXXXfj6+rJ//35Onz7NHXfccc3julpfffUVAF26dKn0vdCtWzcWLlxY6XLcqr5HRKRmUAIlInKdTp06BVBh6Vk5s9lsF3e9LtXPpTRr1ozf/e53rFu3jv/+97+0bNmSpKQkzpw5wwsvvFDpB7/rcfjwYQC7E+COHz8OQGpqKqmpqZd87ZkzZyqUVfYhunzMl6u7+PlIEyZM4M0336RBgwaEh4fTsGFDatWqBcAHH3xQ6ewicMmjrl1dXSkrKzN+Lv/9+vj4VBp/tb+3q1HZPa/M73//e1atWkVCQgIrVqww9sa1bt2a2NhYu5m+qrrW67rU63x8fMjPz//VE6jr+Tuu6ntERGoGJVAiItep/MPTkSNHKq3Pz8+3i4P/LaW6cAak3MmTJy/b37Xsoxo2bBjr1q3jnXfeYdasWbzzzjvGoQyOVn4aWqdOnYyy8mt/7733eOSRRxze55X89NNPvPXWW7Ru3Zp169ZV+GC+atWq6+6j/Bp/+umnSusv9f64Xvv37ycvLw83N7cKz4KqTHh4OOHh4fzyyy/85z//4Z///CdLlizhySefrLA/rSqudV/fkSNHKp2pK79/F/6OXFxcKv1bgSv/vVTVtfwdi8hvk/ZAiYhcp/LN+BcfF10uIyMDwO7Dbflx0OVHcV9ox44djh4iDzzwAM2bN2fFihV89tlnfP311/Tu3dtY7uYo+fn5xt6qQYMGGeXlyVRWVpZD+6uqAwcOUFZWRvfu3SskT3l5eZc85v1q3HHHHQQEBJCfn88333xTof7ipZ+O8uqrrwLQu3fvq5qxqV27NiEhIbz88stMnz4dm81md+x2+UzejZpBqex+7Nu3jyNHjhAQEGB3LSaTqdK/Faj87+Vaxl6+VHTr1q2cO3euQn35aXpVSVJFpGZTAiUicp369OmDl5cXKSkpFT4UJiYmsmPHDgIDA+1mZMq/f/fdd+2eD/T9998ze/Zsh4/RxcWFp556ioKCAuMwg6eeesqhfWRnZ9OvXz+OHz9Oz5496dmzp1HXu3dvAgICeOeddy75bJwvv/zSWOrnaE2aNAHg888/t5vJKCws5IUXXqCkpMQh/QwePBibzcaUKVPsPrx///33vPXWWw7po1xBQQExMTGsXLkSk8lU4TlQldm4cSM///xzhfLy2ZXyJY3wv+WAl1raeL0WLVpk13ZpaSlTp07FZrMxePBgu9hOnTqRm5vL+vXr7cqXLVtm7K+7UL169XBxcbmqsTdq1Ijw8HDy8vJYsGCBXd2ePXtYunQpHh4edv8wICK/TVrCJyJynerUqcNf//pXhg4dSt++fXnkkUdo1qwZu3btYv369Xh6evLmm2/aLXXq1asXFouF1atXk5eXR3BwMIcPH+Yf//gHDz30kEOWlF1syJAhzJw5k8OHD9OiRYtrPlp5586dxkEQxcXFHD16lB07dhin3j3++OPMmzfP7jXu7u4sX76cRx99lD/84Q907NiRtm3bUqdOHfLy8vjqq6/Yt28fGRkZ13wa3OWYzWb69+/PqlWr6Nq1K927d+fUqVOkp6dTq1Yt2rRpU+nDd6/W888/zyeffMKaNWvo2rUrDzzwAKdOneKjjz7CarVW+nDhqii/32VlZZw6dYp9+/aRlZXFL7/8gsVi4a233iIgIOCK7UyaNInvv/+e+++/nyZNmlCrVi2+/vprPv30U7y8vOyeB9ajRw8WLFjAK6+8wp49e4xZ0wuP6r4eXbp0oWvXrvTr14+6deuyYcMGdu/eTfv27Xn++eftYkePHs0///lPhgwZQt++ffHx8TGO4X/ooYdYt26dXXydOnXo0qULWVlZREZG0q5dO9zc3Ljvvvsue5jH/Pnz6dmzJ3/605/IyMigU6dOHD58mI8//phffvmFBQsWVDjCXER+e5RAiYg4QM+ePVm/fj3z58/ns88+IyUlBR8fHx5//HHGjx9f4TlLHh4epKSkMGXKFDZs2EB2djYtWrRg5syZhIWF3ZAEql69evTu3ZtVq1bx5JNPXvPelV27drFr1y7g/DIwT09PWrRowYsvvsjAgQMJCgqq9HWtW7dm8+bNvPnmm6xZs4YVK1Zgs9kwm820atWK6OjoSvfEOEpCQgLNmjVj9erVLF68mPr169OrVy8mTJjAE0884ZA+PDw8+Pjjj3n11Vf56KOPjGcRxcTE8PDDD19zAlU+K+nu7s7tt9+On58fERER9OnTh549e+Lu7l6ldmJiYvjkk0/YsWOHseTUz8+PUaNG8eyzz9olB2FhYcyZM4d33nmHxYsXGw+3dVQCNXPmTNLS0li2bBnff/899evX59lnnyU+Pr7CcfMhISEkJSXx6quvkpqayq233sp9993Hhg0bSElJqZBAwfkZrokTJ7JlyxY2bNhAWVkZsbGxl02gmjZtysaNG3nttddYu3Ytn3/+OXXq1OH+++9n9OjRl30eloj8drgUFBTYrhwmIiI3O5vNZiyF2r179w2Z6REREanptAdKROQ3IjU1lW+++Yb+/fsreRIREblGTkug5s+fT/fu3WncuDEtWrQgMjKS3bt3XzL+hRdewGQykZCQYFdeVFTEuHHjCAgIwM/Pj8cee4y8vDy7mEOHDhEZGYmfnx8BAQGMHz++wgk7mZmZhIWFYTabadu2rXGK1IUWL17MPffcg9lsJiwsjC1btlzHHRAR+XXMnTuXCRMmEB0dTe3atR22BEtEROS3yGkJVGZmJk8//TTr1q0jNTUVNzc3+vbty4kTJyrEpqSksH37dho2bFihLj4+nrS0NJYsWcKaNWs4ffo0kZGRxilLpaWlREZGUlhYyJo1a1iyZAmpqalMnDjRaOPAgQMMGjSI4OBgMjIyGDNmDOPHjyclJcWIWb16NXFxccTExJCRkUFwcDADBw68YacTiYg4yp/+9Cf+9re/0bhxY5YvX15hP5aIiIhUXbXZA1VYWEiTJk1ITEykV69eRvn333/PQw89xMcff8yAAQMYMWIE0dHRwPmH5915550sXLjQOFY0NzeXNm3akJycTHh4OBs2bGDQoEHs3LnT2ByblJTE6NGj2bdvH3Xr1mXq1KmkpaWxfft2o9/o6Gj27t3Lhg0bgPMPHgwKCuL11183Ytq3b09ERARTp0694fdHREREREScr9rsgSosLKSsrMw4JhWgpKSEYcOGMXbsWCwWS4XXZGdnU1xcTI8ePYwyf39/LBaL8VyIbdu2YbFY7E4WCg8Pp6ioiOzsbCPmwjbKY3bs2EFxcTHnzp0jOzu7QkyPHj0qff6EiIiIiIjUTNUmgYqLi6NNmzYEBwcbZbNmzaJevXo8/fTTlb7myJEjuLq6Gg/7K+fj48ORI0eMGB8fH7t6b29vXF1dLxvj4+NDSUkJx44d49ixY5SWllYaU96GiIiIiIjUfNXiOVATJkzg888/Z+3atbi6ugLn90h98MEHxnMqrobNZrN7vsmlnnVyuRibzWaUX/j95foREREREZGazekzUPHx8axatYrU1FS7jc2bNm3i8OHDWCwWvL298fb25tChQ0ydOpXWrVsD4OvrS2lpKceOHbNr8+jRo8Zska+vb4VZootnlCqLOXr0KG5ubnh5eVWYsaqsHxERERERqfmcmkDFxsaSnJxMamoqLVu2tKsbNmwYmzdvZtOmTcZXw4YNefbZZ43T8dq1a4e7uzvp6enG6/Ly8sjJyaFz584ABAcHk5OTY3e0eXp6Oh4eHrRr186I2bhxo13/6enp3Hvvvbi7u3PrrbfSrl07u37KY8r7ERERERGRms9pS/jGjh1LUlISy5cvx2QykZ+fD0CdOnW4/fbb8fHxqTC74+bmhtls5q677gLA09OTJ554gilTpuDj40O9evWYOHEiQUFBdOvWDTh/0ENgYCAjR45kxowZnDhxgilTpjB06FDq1q0LQFRUFG+//TZxcXFERUWxdetWPvjgAxYvXmz0/dxzz/HMM8/QoUMHOnfuzNKlSzl8+DBRUVG/wt0SEREREZHqwGkJVHlyEhERYVceGxtLfHx8lduZOXMmrq6uREVFcfbsWUJDQ1m0aJGxl8rV1ZWkpCTGjh1Lz549qVWrFgMGDGDGjBlGG82aNWPlypVMmDCBpUuX0qBBA2bPnm03tkcffZTjx48zd+5c8vPzCQwMZOXKlTRp0uR6boOIiIiIiNxEqs1zoERERERERKq7anEKn4iIiIhIdXTmzBlKSkqcPQxxsDp16uDmdm2pkBIoEREREZFKFBUVAef33UvNYbPZKCgo4I477rimJMrpx5iLiIiIiFRHZ8+e5bbbbnP2MMTBXFxcMJlMnDlz5pperwRKREREROQSXFxcnD0EuQGu5/eqBEpERERERKSKlECJiIiIiIhUkRIoERERERGRKlICJSIiIiJSQ4waNYrIyEhnD6NGUwIlIiIiIiJSRUqgRERERER+Iw4dOsTgwYPx9/fH39+fIUOGkJeXB0BhYSH169fn3//+txHfunVrgoODjZ/T09Px8/OjuLi40vZNJhMpKSkAWK1WZs2aBdSsmTE9SFdEREREpIouTgKSkpJ+1frrYbPZGDx4MLVq1SI1NRUXFxfGjRvH4MGDSU9P5/bbb6dt27Zs2rSJjh078u2333Lq1CmOHz/O4cOHadCgAZmZmQQHB+Pu7u6wcd1sNAMlIiIiIvIbsHHjRnbt2sXbb79N+/btuffee1m8eDFffvkln332GQAhISFs2rQJgMzMTLp06UL79u3JzMw0ykJCQpx2DdWBEigRERERkd+AnJwcGjZsSNOmTY2yZs2a0bBhQ/bu3QucT6C2bt1KcXExmZmZdO3alZCQEDIzM/n555/Zvn27EihnD0BERERERG48m82Gi4tLpXXl5VarlaKiIrZv387mzZuNBGrTpk1s3boVd3d3OnTocMk+brnlFsrKyiqUl5aWXrLvm432QImIiIiIVNGV9iTd6Prr0apVK3744QcOHjxozEIdOHCAH3/8kVatWgEY+6DeffddCgsLadu2LcXFxeTm5rJy5cor7n/y9/fn+++/r1D+3XffYbFYbsyF/co0AyUiIiIiUoOcOnWKr776yu7r4MGDdOvWjbvvvpsRI0aQnZ3Njh07GD58OG3btiU0NNR4fUhICCtXrsRqteLq6kqtWrXo0KEDK1euvOLyvb59+7Jo0SI+//xzbDYbZ86c4a233uKLL77g97///Y2+9F+FZqBERERERGqQrKwsu4QI4JFHHuG9994jMTGR2NhYI5kJCwtjzpw5dsvrunbtyoIFC+ySpa5du7Jly5YrJlDx8fEUFRUxaNAgTp06xd69e2nevDmLFi2iZ8+eDrxK53EpKCiwOXsQIiIiIiLVzcmTJ/H09HT2MG5KJSUldOnShT59+jBt2jRnD6dS1/r71RI+ERERERFxKDc3N9zd3alVq5azh+JwSqBERERERESqSHugRERERETE4bKyspw9hBtCM1AiIiIiIiJVpARKRERERESkipRAiYiIiIiIVJESKBERERERkSpSAiUiIiIiIlJFSqBERERERESqyGkJ1Pz58+nevTuNGzemRYsWREZGsnv3bqO+uLiYqVOnct999+Hn54fFYmHYsGEcOnTIrp2ioiLGjRtHQEAAfn5+PPbYY+Tl5dnFHDp0iMjISPz8/AgICGD8+PGcO3fOLiYzM5OwsDDMZjNt27Zl6dKlFca8ePFi7rnnHsxmM2FhYWzZssWBd0RERERE5ObQqFEjEhMTjZ9NJhMpKSlOHNGvx2kJVGZmJk8//TTr1q0jNTUVNzc3+vbty4kTJwD4+eef+fLLLxk7diyfffYZH3zwAXl5eQwYMICSkhKjnfj4eNLS0liyZAlr1qzh9OnTREZGUlpaCkBpaSmRkZEUFhayZs0alixZQmpqKhMnTjTaOHDgAIMGDSI4OJiMjAzGjBnD+PHj7d4Eq1evJi4ujpiYGDIyMggODmbgwIEVEjoREREREWcZNWoUJpOJ6OjoCnVTpkzBZDIRGRnp8H5zcnLo2bOnw9utjlwKCgpszh4EQGFhIU2aNCExMZFevXpVGrN37166dOnC5s2bCQoK4uTJk9x5550sXLiQQYMGAZCbm0ubNm1ITk4mPDycDRs2MGjQIHbu3Im/vz8ASUlJjB49mn379lG3bl2mTp1KWloa27dvN/qKjo5m7969bNiwAYDw8HCCgoJ4/fXXjZj27dsTERHB1KlTb9RtEREREREnOXnyJJ6ens4exlUZNWoUmzZtoqCggJycHOrUqQNASUkJQUFBuLu7ExQURFJS0nX106hRI+bMmcPgwYMdMWynuNbfb7XZA1VYWEhZWRkmk+mSMadPnwYwYrKzsykuLqZHjx5GjL+/PxaLha1btwKwbds2LBaLkTzB+WSoqKiI7OxsI+bCNspjduzYQXFxMefOnSM7O7tCTI8ePYx+RERERESqg6CgIAICAvjoo4+MsnXr1uHh4UFISEiF+OXLl9O5c2fMZjMdOnRg4cKFlJWVGfX79++nT58+mM1mOnbsyNq1ayu0ceESvoMHD2IymdixY8cVY1atWkXv3r1p0KABXbt2ZdeuXezevZsHH3wQPz8/evbsyYEDBy55rbNmzcJqtQKQkpJi5AmXGoMjuDm8xWsUFxdHmzZtCA4OrrT+3LlzTJo0iZ49e9KoUSMAjhw5gqurK97e3naxPj4+HDlyxIjx8fGxq/f29sbV1dUuplu3bhXaKCkp4dixY9hsNkpLSyu0c2E/IiIiIlLzfRv57TW9zqO5B/6v+ldalxuXS9F3RZXWtUhqcU39PfHEEyQmJjJkyBDgfJI0ePDgCsnIsmXLmDlzJnPmzKFt27bs2bOHF154AXd3d0aMGEFZWRlDhgzB09OT9evX88svvxAXF0dRUeXjvVqzZs1i5syZNGvWjDFjxjB8+HC8vb2ZNGkSPj4+jBo1itjY2OueMXOkajEDNWHCBD7//HPef/99XF1dK9SXlJQwYsQITp48yV//+tcrtmez2XBxcTF+vvD7C10uxmazVSnmUm2LiIiIiDjLwIED2bFjB99++y35+fl8+umn/OEPf6gQN3fuXKZNm0ZERATNmjWjV69evPjiiyxZsgSAjRs3snfvXv72t7/Rtm1bunTpwqxZs+zOJLgezz33HA8++CAtW7bk+eefZ8+ePYwYMYLQ0FACAwMZPnw4mZmZDunLUZw+AxUfH8/q1atJS0ujWbNmFepLSkp4+umn2b17N3//+9/x8vIy6nx9fSktLeXYsWPUr1/fKD969Cj33XefEXPxMrtjx47ZzSj5+vpWmEk6evQobm5ueHl5YbPZ7GasLoy5eFZKRERERMTZTCYTv//971m+fDmenp6EhITQuHFju5ijR4+Sm5vLSy+9RExMjFFeUlJiTCbk5OTg5+dn99qOHTtyyy2OmYcJCgoyvvf19a207MyZM/z888/cdtttDunzejl1Bio2Npbk5GRSU1Np2bJlhfri4mKioqL4+uuvSUtLw2w229W3a9cOd3d30tPTjbK8vDxycnLo3LkzAMHBweTk5NgdbZ6eno6Hhwft2rUzYjZu3GjXdnp6Ovfeey/u7u7ceuuttGvXzq6f8pjyfkREREREqpMhQ4bw4Ycfsnz5cmMp34XK9znNnz+fTZs2GV9ZWVl8/vnnwP9WZV2N8uTqwtcWFxdXGuvu7m58X76yy83NrULZhXuyLuTq6s5gOnIAACAASURBVFppXXnZjVgt5rQZqLFjx5KUlMTy5csxmUzk5+cDUKdOHW6//XZKSkr44x//yI4dO1ixYgUuLi5GTN26dalduzaenp488cQTTJkyBR8fH+rVq8fEiRMJCgoy9jT16NGDwMBARo4cyYwZMzhx4gRTpkxh6NCh1K1bF4CoqCjefvtt4uLiiIqKYuvWrXzwwQcsXrzYGO9zzz3HM888Q4cOHejcuTNLly7l8OHDREVF/bo3TkRERESc5lr3JF3OpfZGXa+wsDDc3d05duwYffr0qVDv6+uLn58f3333HY8//nilbbRq1YoffviB3Nxc41C2//znP5dMaABjZdjhw4eNsp07d17PpVySv78/P/zwQ4Ulhfv377cbiyM5LYEqT04iIiLsymNjY4mPjycvL481a9YAVDjgYeHChcaRiTNnzsTV1ZWoqCjOnj1LaGgoixYtMvZSubq6kpSUxNixY+nZsye1atViwIABzJgxw2ivWbNmrFy5kgkTJrB06VIaNGjA7Nmz7cb26KOPcvz4cebOnUt+fj6BgYGsXLmSJk2aOPzeiIiIiIhcLxcXFzZv3ozNZsPDw6PSmLi4OMaPH4+npycPPvggxcXFfPnll/z444+MGTOGbt260bJlS0aOHMnMmTM5e/YsEyZMsJsluljt2rXp1KkTCxYsoHnz5pw6dYpp06bdkGsMDw/HxcWFSZMmERgYCMDu3bt55ZVXuPfee+1O4nYUpyVQBQUFl61v2rTpFWMAatWqxdy5c5k7d+4lYxo3bnzFkztCQkLIyMi4bMywYcMYNmzYFcckIiIiIlId3HHHHZetHzp0KLfddhuvv/46r7zyCrVq1TIOb4Dzy/GWL1/O6NGjeeCBB/D392fGjBlG/aW88cYbjB49mh49etC8eXNee+01evfu7bDrKmc2m1mxYgXx8fEsWrQIOL8CrXfv3vzpT39yeH9QjR6kKyIiIiJSndyMD9L9Lfvwww8ZOXIkP/74I7Vr175i/E3/IF0REREREZFrVZ40VSV5uh5KoERERERERKpICZTIb1xycjJWqxUvLy+sVivJycnOHpKIiIjIVYuIiKjSGQrXy+kP0hUR50lOTmb69OkkJCRgtVrJysoiOjoagAEDBjh5dCIiIiLVjw6REPkNs1qtzJ49m9DQUKMsIyOD2NhYsrKynDgyERER59MhEjXbtf5+lUCJ/IZ5eXmRn59v9xTw4uJizGYzx48fd+LIREREnE8JVM2mU/hE5KpZLJYKM01ZWVlYLBYnjUhERESkelMCJfIbFhMTQ3R0NBkZGRQXF5ORkUF0dDQxMTHOHpqIiIhItaRDJER+w8oPioiNjSUnJweLxcLkyZN1gISIiIjIJWgPlIiIiIhIJbQH6tIaNWrEnDlzGDx4MAAmk4lly5YRERHh5JFVnfZAiYiIiIj8xo0aNQqTyWQ8luRCU6ZMwWQyERkZ6fB+c3Jy6Nmzp8PbrY6UQImIiIiI1CD+/v589NFHnDlzxigrKSkhKSkJf3//G9Kn2WzGw8PjhrRd3SiBEkNycjJWqxUvLy+sVivJycnOHpKIiIiIXKWgoCACAgL46KOPjLJ169bh4eFBSEhIhfjly5fTuXNnzGYzHTp0YOHChZSVlRn1+/fvp0+fPpjNZjp27MjatWsrtGEymUhJSQHg4MGDmEwmduzYccWYVatW0bt3bxo0aEDXrl3ZtWsXu3fv5sEHH8TPz4+ePXty4MCBS17rrFmzsFqtAKSkpGAymS47BkfQIRICnE+epk+fTkJCAlarlaysLGPqVwcKiIiIiJwX+e23dj8ntWhRadw/T53i7Z9+Mn4Or1uXET4+lcbG5ebyXVGR8fMsf38CrnM254knniAxMZEhQ4YA55OkwYMHV0hGli1bxsyZM5kzZw5t27Zlz549vPDCC7i7uzNixAjKysoYMmQInp6erF+/nl9++YW4uDiKLhjv9Zg1axYzZ86kWbNmjBkzhuHDh+Pt7c2kSZPw8fFh1KhRxMbGkpSU5JD+HEEzUALAvHnzSEhIIDQ0FHd3d0JDQ0lISGDevHnOHpqIiIiIXKWBAweyY8cOvv32W/Lz8/n000/5wx/+UCFu7ty5TJs2jYiICJo1a0avXr148cUXWbJkCQAbN25k7969/O1vf6Nt27Z06dKFWbNmUVJS4pBxPvfcczz44IO0bNmS559/nj179jBixAhCQ0MJDAxk+PDhZGZmOqQvR9EMlADnN/6VT3+Ws1qt5OTkOGlEIiIiInKtTCYTv//971m+fDmenp6EhITQuHFju5ijR4+Sm5vLSy+9ZPcMyJKSEmy28wd15+Tk4OfnZ/fajh07csstjpmHCQoKMr739fWttOzMmTP8/PPP3HbbbQ7p83opgRIALBYLWVlZhIaGGmVZWVlYLBYnjkpERERErtWQIUMYNWoUderUYcKECRXqy/c5zZ8/n86dO1faRnkidTXKk6sLX1tcXFxprLu7u/G9i4sLAG5ubhXKLtyTdSFXV9dK68rLyl/vSEqgBICYmBiio6Mr7IGaPHmys4cmIiIiUm1cas/TxR6oW5cH6tatUuyrN+hkvLCwMNzd3Tl27Bh9+vSpUO/r64ufnx/fffcdjz/+eKVttGrVih9++IHc3FzjBL///Oc/l0xoAOrXrw/A4cOHjbKdO3dez6Vckr+/Pz/88EOFJYX79++3G4sjKYES4H8HRcTGxpKTk4PFYmHy5Mk6QEJERETkJuXi4sLmzZux2WyXPGI8Li6O8ePH4+npyYMPPkhxcTFffvklP/74I2PGjKFbt260bNmSkSNHMnPmTM6ePcuECRPsZokuVrt2bTp16sSCBQto3rw5p06dYtq0aTfkGsPDw3FxcWHSpEkEBgYCsHv3bl555RXuvffeG3Jsuw6REMOAAQPIysri+PHjZGVlKXkSERERucndcccd1L3MTNjQoUN54403SEpKIiQkhF69erFs2TKaNm0KnF+Ot3z5csrKynjggQcYOXIkY8eOveIzn9544w0AevTowUsvvcSkSZMcd1EXMJvNrFixgi1btvDCCy8YfbZo0YIPPvjghvTpUlBQcPULG0VEREREariTJ0/i6enp7GFIFX344YeMHDmSH3/8kdq1a18x/lp/v5qBEhERERGRm1550lSV5Ol6KIESERERERGpIi3hExERERGphJbw1WxawiciIiIiInKDKYESERERERGpIiVQIiIiIiKXYLNpt0tNdD2/V6clUPPnz6d79+40btyYFi1aEBkZye7du+1ibDYbs2bNolWrVjRo0IA+ffqwZ88eu5iCggJGjBhBkyZNaNKkCSNGjKCgoMAu5uuvv6Z37940aNCAwMBAZs+eXeGmpaSk0LlzZ3x9fencuTNpaWlXPRYRERERqTlq1arFzz//7OxhiIPZbDYKCgqoU6fONb3+0o8QvsEyMzN5+umnad++PTabjZkzZ9K3b1+2bt1KvXr1AFiwYAELFy5k4cKF3HXXXcyZM4d+/frxxRdfcMcddwAwbNgwcnNz+X//7//h4uLC6NGjeeaZZ0hKSgLg1KlT9OvXj/vuu49//etf7Nu3j+eee47bbruN6OhoALZt28ZTTz1FfHw8Dz/8MGlpaTz55JOsW7eOjh07VnksIiIiIlJzeHh4UFJSwsmTJ509FHGwO+64Aze3a0uFqs0pfIWFhTRp0oTExER69eqFzWajVatWDB8+nLFjxwLwyy+/cNdddzF9+nSioqLIycmhc+fOrF27li5dugCQlZVFr169+OKLL7jrrrtYsmQJL7/8Mv/973+NM+Hnzp3L0qVL2b17Ny4uLkRFRXHixAk+/vhjYzwRERHUr1+fJUuWVGksIiIiIiJS81WbPVCFhYWUlZVhMpkAOHjwIPn5+fTo0cOIqV27Nvfddx9bt24Fzs8c3X777XTu3NmI6dKlC3Xq1LGLsVqtdg/UCg8P58cff+TgwYMAfPHFF3b9lMeUt1GVsYiIiIiISM1XbRKouLg42rRpQ3BwMAD5+fkA+Pj42MX5+Phw5MgRAI4cOYK3tzcuLi5GvYuLC/Xr17eLqayN8rryvi7XT1XGIiIiIiIiNZ/T9kBdaMKECXz++eesXbsWV1dXu7oLkyM4v+nr4oTpYleKKT9A4koxF5dVJUZERERERGoup89AxcfHs2rVKlJTU2nWrJlRbjabASrM8Bw9etSYCfL19eXo0aN2J+rZbDaOHTtmF1NZG/C/GSWz2XzZfqoyFhERERERqfmcmkDFxsaSnJxMamoqLVu2tKtr2rQpZrOZ9PR0o+zs2bNkZWUZe56Cg4MpLCxk27ZtRsy2bds4c+aMXUxWVhZnz541YtLT02nYsCFNmzYFoFOnTnb9lMeUt1GVsYiIiIiISM3nGhcX97IzOh47diwffvgh7777Lv7+/pw5c4YzZ84AcOutt+Li4kJpaSl//vOfufPOOyktLWXixInk5+fzl7/8BQ8PD+rXr8+///1vkpOTueeee8jLy+Oll16iffv2PPPMMwC0aNGCd955h507d3LXXXeRlZXFlClTePHFF43kp2HDhsycORN3d3e8vb1ZtmwZiYmJLFiwAD8/vyqNRUREREREaj6nHWNeftrexWJjY4mPjwfOL8d79dVXeffddykoKKBDhw689tprtG7d2og/ceIEsbGx/OMf/wCgV69ezJkzx679r7/+mrFjx7J9+3ZMJhNRUVHExsba7V9KSUlhxowZHDhwgObNmzNp0iQeeeQRo74qYxERERERkZqt2jwHSkREREREpLpz+iESIiIiIiLiOMnJyVitVry8vLBarSQnJzt7SDVKtTjGXERERERErl9ycjLTp08nISEBq9VKVlYW0dHRAAwYMMDJo6sZNAMlIiIiIjeEZkJ+ffPmzSMhIYHQ0FDc3d0JDQ0lISGBefPmOXtoNYZmoERERETE4TQT4hw5OTlYrVa7MqvVSk5OjpNGVPNoBkpEREREHE4zIc5hsVjIysqyK8vKysJisThpRDWPEigRERERcTjNhDhHTEwM0dHRZGRkUFxcTEZGBtHR0cTExDh7aDWGlvCJiIiIiMOVz4SEhoYaZZoJufHKl0fGxsaSk5ODxWJh8uTJWjbpQHoOlIiIiIg43KX2QOnDvNzsNAMlIiIiIg6nmRCpqTQDJSIiIiIiUkU6REJERERERKSKlECJiIiIiIhUkRIoERERERGRKlICJSIiIiIiUkVKoERERERERKpICZSIiIiIiEgVKYESERERERGpIiVQIiIiIiIiVaQESkREREREaoTk5GSsViteXl5YrVaSk5Md3oebw1sUERERERH5lSUnJzN9+nQSEhKwWq1kZWURHR0NwIABAxzWj0tBQYHNYa2JiIiIiIg4gdVqZfbs2YSGhhplGRkZxMbGkpWV5bB+lECJiIiIiMhNz8vLi/z8fNzd3Y2y4uJizGYzx48fd1g/2gMlIiIiIiI3PYvFUmGmKSsrC4vF4tB+lECJiIiIiMhNLyYmhujoaDIyMiguLiYjI4Po6GhiYmIc2o8OkRARERERkZte+UERsbGx5OTkYLFYmDx5skMPkADtgRIREREREakyLeETERERERGpIiVQIiIiIiIiVeTUBGrz5s089thjBAYGYjKZSExMtKsvLCxk3LhxtG7dmgYNGtCxY0cWLlxoF1NUVMS4ceMICAjAz8+Pxx57jLy8PLuYQ4cOERkZiZ+fHwEBAYwfP55z587ZxWRmZhIWFobZbKZt27YsXbq0wngXL17MPffcg9lsJiwsjC1btjjoToiIiIiIyM3AqQnUmTNnaN26Na+++iq1a9euUD9x4kTWr1/PokWL2Lp1KzExMUybNo0PP/zQiImPjyctLY0lS5awZs0aTp8+TWRkJKWlpQCUlpYSGRlJYWEha9asYcmSJaSmpjJx4kSjjQMHDjBo0CCCg4PJyMhgzJgxjB8/npSUFCNm9erVxMXFERMTQ0ZGBsHBwQwcOJBDhw7dwDskIiIiIiLVSbU5RKJRo0bMmTOHwYMHG2VWq5WHH36YCRMmGGW9e/cmKCiIuXPncvLkSe68804WLlzIoEGDAMjNzaVNmzYkJycTHh7Ohg0bGDRoEDt37sTf3x+ApKQkRo8ezb59+6hbty5Tp04lLS2N7du3G/1ER0ezd+9eNmzYAEB4eDhBQUG8/vrrRkz79u2JiIhg6tSplV5T5Lffnu+vRYtK6/956hRv//ST8XN43bqM8PGpNDYuN5fvioqMn2f5+xPg4VEhbn9REfG5ucbPzT08ePX/rvtif/vpJz49dcr4ebiPDw/UrXvZaymna9I1XUzXpGu6kK5J13QxXZOu6UK6Jl3TxarTNV0qply13gPVpUsX1q5dS+7/3cytW7eya9cuwsPDAcjOzqa4uJgePXoYr/H398disbB161YAtm3bhsViMZInOJ8MFRUVkZ2dbcRc2EZ5zI4dOyguLubcuXNkZ2dXiOnRo4fRj4iIiIiI1HzV+jlQs2fP5qWXXuLuu+/Gze38UOfMmUPPnj0BOHLkCK6urnh7e9u9zsfHhyNHjhgxPhdlxN7e3ri6utrFdOvWrUIbJSUlHDt2DJvNRmlpaYV2LuxHRERERERqvmqdQL311lts3bqVFStW0LhxY7Zs2cLkyZNp0qQJDzzwwCVfZ7PZcHFxMX6+8PsLXS7GZrMZ5Rd+f7l+RERERESkZqu2e6B++eUXmjRpwrvvvkufPn2MuOjoaL7//ntSUlL47LPPiIiI4JtvvqF+/fpGTJcuXXjkkUeYMGECf/rTn0hLS+Pzzz836o8ePcqdd95JamoqoaGh9OrVi6CgIF577TUj5uOPP2bYsGH8+OOP2Gw2GjZsyJIlS+jbt68RM3bsWHbv3s2aNWtu5K0REREREZFqotrugSouLqa4uBhXV1e7cldXV8rKygBo164d7u7upKenG/V5eXnk5OTQuXNnAIKDg8nJybE72jw9PR0PDw/atWtnxGzcuNGun/T0dO69917c3d259dZbadeunV0/5THl/YiIiIiISM3nGhcX97KzOi8sLGTv3r3k5+fz/vvv07p1a+rWrcu5c+fw9fVl06ZNrF27FovFgs1m45NPPmH+/Pk888wzdOzYkVq1anH48GHefvtt7r77bk6ePMlLL71E3bp1mTZtGrfccgvNmjUjLS2Nf/3rXwQFBbF3717Gjh3LwIEDefjhhwFo3rw5f/nLX/jpp59o3Lgxa9asYd68ecyYMYNWrVoBcMcddzBr1iwaNGhArVq1mDt3Llu2bOGNN97A09PTWbdQRERERER+RU5dwrdp0yYjibnQ448/zptvvkl+fj7Tpk0jPT2dEydO0LhxY4YOHcrzzz9v7D06e/YskydPJjk5mbNnzxIaGsq8efPsTt07dOgQY8eOJSMjg1q1ajFgwABmzJiBxwVHKGZmZjJhwgT27t1LgwYNePHFF3nqqafsxrV48WIWLFhAfn4+gYGBzJw5k/vvv/8G3R0REREREaluqs0eKBERERERkequ2u6BEhERERERqW6UQImIiIiIiFSREigREREREZEqUgIlIiIiIiJSRUqgREREREREqkgJlIiIiIiISBUpgRIREREREakiJVAiIiIiIiJVpARKRERERESkipRAiYiIiIiIVJESKBERERERkSpSAiUiIiIiIlJFSqBERERERESqSAmUiIiIiIhIFSmBEhFxguTkZKxWK15eXlitVpKTk509JBEREakCN2cPQETktyY5OZnp06eTkJCA1WolKyuL6OhoAAYMGODk0YmIiMjluBQUFNicPQgRkd8Sq9XK7NmzCQ0NNcoyMjKIjY0lKyvLiSMTERGRK1ECJSLyK/Py8iI/Px93d3ejrLi4GLPZzPHjx504MhEREbkS7YESEfmVWSyWCjNNWVlZWCwWJ41IREREqkoJlIjIrywmJobo6GgyMjIoLi4mIyOD6OhoYmJinD00ERERuQIdIiEi8isrPygiNjaWnJwcLBYLkydP1gESIiIiNwHtgRIREREREakiLeETERERERGpIiVQIiIiIjeQHpwtUrNoD5SIiIjIDaIHZ4vUPNoDJSIiInKD6MHZIjWPEigRERGRG0QPzhapebQHSkREROQG0YOzRWoepyZQmzdv5rHHHiMwMBCTyURiYmKFmG+++YYhQ4bQpEkTGjZsSGhoKDk5OUZ9UVER48aNIyAgAD8/Px577DHy8vLs2jh06BCRkZH4+fkREBDA+PHjOXfunF1MZmYmYWFhmM1m2rZty9KlSyuMZfHixdxzzz2YzWbCwsLYsmWLg+6EiIiI1ER6cLZIzePUBOrMmTO0bt2aV199ldq1a1eoP3DgAA899BBNmzYlNTWVrKwsJk2aRJ06dYyY+Ph40tLSWLJkCWvWrOH06dNERkZSWloKQGlpKZGRkRQWFrJmzRqWLFlCamoqEydOtOtn0KBBBAcHk5GRwZgxYxg/fjwpKSlGzOrVq4mLiyMmJoaMjAyCg4MZOHAghw4duoF3SERERK6XM0/BGzBgAJMnTyY2Nhaz2UxsbKwenC1yk6s2e6AaNWrEnDlzGDx4sFE2bNgwXFxcePvttyt9zcmTJ7nzzjtZuHAhgwYNAiA3N5c2bdqQnJxMeHg4GzZsYNCgQezcuRN/f38AkpKSGD16NPv27aNu3bpMnTqVtLQ0tm/fbrQdHR3N3r172bBhAwDh4eEEBQXx+uuvGzHt27cnIiKCqVOnOvx+iIiIyPW71Cl4SmJE5FpV2z1QZWVlrF27FovFQv/+/WnRogXdu3dn9erVRkx2djbFxcX06NHDKPP398disbB161YAtm3bhsViMZInOJ8MFRUVkZ2dbcRc2EZ5zI4dOyguLubcuXNkZ2dXiOnRo4fRj4iIiFQ/8+bNIyEhgdDQUNzd3QkNDSUhIYF58+Y5e2gicpOqtgnUTz/9RGFhIfPnz6d79+589NFH9O/fn+HDh7N27VoAjhw5gqurK97e3nav9fHx4ciRI0aMj4+PXb23tzeurq6XjfHx8aGkpIRjx45x7NgxSktLK40pb0NERESqn5ycHKxWq12Z1Wq1208tInI1qu2DdMvKygDo3bs3zz//PAD33HMP2dnZLF68mJ49e17ytTabDRcXF+PnC7+/0OVibDabUX7h95frR0RERKqX8lPwLnwOk07BE5HrUW1noLy9vXFzc6vwH7iWLVuSm5sLgK+vL6WlpRw7dswu5ujRo8Zska+vb4VZootnlCqLOXr0KG5ubnh5eVWYsaqsHxEREal+dAqeiDhatU2gbr31Vtq3b8++ffvsyr/55hsaN24MQLt27XB3dyc9Pd2oz8vLIycnh86dOwMQHBxMTk6O3dHm6enpeHh40K5dOyNm48aNdv2kp6dz77334u7uzq233kq7du3s+imPKe9HREREqh+dgicijubUJXyFhYXs378fOL9kLzc3l6+++op69erRuHFjRo8eTVRUFPfddx+hoaFs2rSJ1atXG8+L8vT05IknnmDKlCn4+PhQr149Jk6cSFBQEN26dQPOH/QQGBjIyJEjmTFjBidOnGDKlCkMHTqUunXrAhAVFcXbb79NXFwcUVFRbN26lQ8++IDFixcbY33uued45pln6NChA507d2bp0qUcPnyYqKioX/emiYiIyFUZMGCAEiYRcRinHmO+adMmHn744Qrljz/+OG+++SYAiYmJzJ8/n7y8PAICAhgzZozdfwTPnj3L5MmTSU5O5uzZs4SGhjJv3jy7U/cOHTrE2LFjycjIoFatWgwYMIAZM2bg4eFhxGRmZjJhwgT27t1LgwYNePHFF3nqqafsxrV48WIWLFhAfn4+gYGBzJw5k/vvv9/Rt0VERERERKqpavMcKBERERERkequ2u6BEhERERERuRrJyclYrVa8vLywWq0kJyc7vI9qe4y5iIiIiIhIVSUnJzN9+nQSEhKwWq1kZWURHR0N4NB9kFrCJyIiIiIiNz2r1crs2bPtnvuWkZFBbGwsWVlZDutHCZSIiIiIiNz0vLy8yM/Px93d3SgrLi7GbDZz/Phxh/WjPVAiIiIiInLTs1gsFWaasrKysFgsDu1HCZSIiIiIiNz0YmJiiI6OJiMjg+LiYjIyMoiOjiYmJsah/egQCRERERERuemVHxQRGxtLTk4OFouFyZMnO/xB2toDJSIiIiLiYMnJycybN8/4IB8TE+PwD/LiHJqBEhERERFxoF/rOG1xDs1AiYiIiIg40K91nLY4hw6REBERkRsuOTkZq9WKl5cXVquV5ORkZw9J5IbJycnBarXalVmtVnJycpw0InEkJVAiIiJyQ5UvZ5o9ezb5+fnMnj2b6dOnK4mSGuvXOk5bnEMJlBj0r4MiInIjzJs3j4SEBEJDQ3F3dyc0NJSEhATmzZvn7KGJ3BC/1nHa4hw6REIAbXYUEZEbR8uZ5Lfm1zpOW5xDh0gIoM2OIiJy4+j/MSJSk2gJnwD610EREblxtJxJRGoSJVACaLOjiIjcOAMGDGDy5MnExsZiNpuJjY3Vciap8bS3vObSHigB/vevgxfvgZo8ebKzhyYiIjXAgAEDlDDJb4b2ltds2gMlhuTkZObNm2dsdoyJidEfuYiIiMhV0r6/mk0JlIiIiIiIA3l5eZGfn4+7u7tRVlxcjNls5vjx404cmTiC9kCJiIiIiDiQ9pbXbEqgRERE5IbThnr5LdHJkzWbDpEQERGRG0ob6uW3Rg/Srdm0B0pERERuKG2oF5GaRAmUiIiI3FDaUC8iNYn2QImgtfkiIjeSNtSLSE2iBEp+88rX5s+ePZv8/Hxmz57N9OnTlUSJiDiINtSLSE3i1ARq8+bNPPbYYwQGBmIymUhMTLxk7AsvvMD/Z+/e46Is8//xv2bxEGoFKAwaHpNFyMpSGcZcTeyTaQcWG7LytxauZyUPgCOZMGezGAAAIABJREFURuoGo2KyaGyp+Klda02yzDKtXSmwRrBPuWxqpOvqesAhDqOioojz+8Pv3MsIKuDMfc3M9Xo+Hj4ect/X3Pc199yn93X08/NDVlaWw/KLFy8iOTkZvXr1QpcuXfDss8/ixIkTDmmOHTuGMWPGoEuXLujVqxfmzp2LS5cuOaTZtWsXhg4dCq1Wi/vvvx85OTkN8rB27Vrcd9990Gq1GDp0KL799ttb+PbkLjIyMpCVlYUhQ4agdevWGDJkCLKyspCRkSE6a0REXsFgMGDBggUwGo3QarUwGo3sUE9EHktoAHXu3DlEREQgPT0dvr6+1023ZcsWfP/99+jcuXODdSkpKdi6dSvWrVuHbdu24ezZsxgzZgzq6uoAAHV1dRgzZgyqq6uxbds2rFu3Dp988gnmz5+vbOPIkSN45plnEBkZifz8fMyZMwdz587Fli1blDSbN2/GvHnzkJiYiPz8fERGRiIuLg7Hjh1z4hEhEUpKSqDX6x2W6fV6lJSUCMoREbkKm+uKYzAYYDabUVlZCbPZzOCJiDyW0ADq0UcfxcKFCxETE4Nf/arxrPznP//BvHnzsHbtWrRq5Tjq+unTp/HnP/8ZixYtwrBhw9CvXz+89dZb2LdvH7766isAwM6dO3HgwAG89dZb6NevH4YNG4bXXnsN7777Ls6cOQMAWL9+PYKDg7Fs2TKEhYXhhRdewHPPPYdVq1Yp+1q9ejWef/55vPDCCwgLC8OyZcug1Wobrakiz8K2+URyYHNdIiJyBrfuA3X58mVMmDABSUlJjb7M7t27F7W1tYiOjlaWhYSEICwsDIWFhQCAoqIihIWFISQkREkzfPhwXLx4EXv37lXS1N+GPc0PP/yA2tpaXLp0CXv37m2QJjo6WtkPeS62zSeSA5vrEhGRM7j1RLppaWnw9/fH73//+0bXl5WVwcfHBx07dnRYHhgYiLKyMiVNYGCgw/qOHTvCx8fHIc3DDz/cYBuXL19GRUUFbDYb6urqGmyn/n7Ic3GyOyI5sLkuERE5g9sGULt27cJ7772HgoKCZn/WZrNBo9Eof9f/f303SmOz2ZTl9f9/o/2Q5zIYDAyYiLycvblu/clc2VyXiIiay22b8BUUFODUqVMICwtDx44d0bFjRxw7dgyvvvoqIiIiAABBQUGoq6tDRUWFw2fLy8uV2qKgoKAGtUQVFRUONUqNpSkvL0erVq0QEBDQoMaqsf0QEZF7Y3NdIiJyBretgZowYQJiYmIclj399NN4+umn8cILLwAA+vXrh9atWyMvLw9xcXEAgBMnTqCkpAQ6nQ4AEBkZieXLl+PEiRO46667AAB5eXlo27Yt+vXrp6T57LPPHPaVl5eHBx54QJk1vV+/fsjLy8Nvf/tbhzRPPfWUC749ERE5G5vrEhGRMwgNoKqrq3H48GEAwJUrV3D8+HEUFxfD398fXbt2bVC706pVK2i1WoSGhgIA7rzzTvzud7/DwoULERgYCH9/f8yfPx/33HOP0qcpOjoa4eHhmDJlCpYsWYKqqiosXLgQ48aNwx133AEAiI+Px5o1azBv3jzEx8ejsLAQ7733HtauXavse/r06Zg8eTL69+8PnU6HnJwcnDp1CvHx8U47Hrm5ucjIyFAe7ImJiXywExE5EZvrEhHRrRIaQP3www948sknlb/T0tKQlpaG5557DtnZ2U3axuuvvw4fHx/Ex8ejpqYGQ4YMwZ/+9Cf4+PgAAHx8fLBx40YkJSXhsccew2233QaDwYAlS5Yo2+jRowc++OADvPzyy8jJyUFwcDBMJpNDDdjo0aNRWVmJZcuWwWKxIDw8HB988AG6devmlGNhH143KysLer0eZrMZCQkJAMCHPRERERGRm9BYrVab6EzQ1ZGgTCaTQ+fm/Px8GI3GBnMUeSPWvhERERGRJ3DbQSRkI/PwupzckojI9XJzc6HX6xEQEAC9Xs97LBFRCzGAchP24XXrk2V4XU5uSUTkWiyoIiJyHjbhcxPX6wMlwwhRAQEBsFgsyoiHAFBbWwutVovKykqBOSMi8g6yNxMnInImBlBuRNZ+QHq9HqNGjcK2bduU727/mw92IqJbx4IqIiLnYRM+N2IwGGA2m1FZWQmz2SxF8AQAgwcPRmZmJsaOHYvjx49j7NixyMzMxODBg0VnjYjIK8jcTJyIyNkYQJFwu3btwqxZs7BhwwaEhIRgw4YNmDVrFnbt2iU6a0REXiExMREJCQnIz89HbW0t8vPzkZCQgMTERNFZIyLyOGzCR8KxaQkRqUXWptKA3N+diMiZWANFwrFpCRGpgSPRERGRMzCAIuHYtISI1CDzlAkMHomInIdN+MgtsGkJEbmazM2FOYw5EZHzsAbKjcg8S7ysIxASkXpkbi5cUlICvV7vsEyv16OkpERQjoiIPBcDKDfB5hVERK4lc3NhmYNHIiJnYxM+N8HmFUREridrc2F7IV1WVhb0ej3MZjMSEhKwYMECKb4/EZEztRKdAbqKzSuIiMhV7EGS0WhUgkcGT0RELcMmfG6CzSvkJXPfNyI1yd5Umn1NiYico1kB1JUrV3Ds2DFUVVW5Kj/Skrltvsxkf6EjOYkqNJB5GHMiInKeZvWBqq2tRefOnZGamooZM2a4Ml9SkrVtvszY941kI7IvjszDmBMRkfM0qwaqdevWCA4OhkajcVV+pMbmFfJh3zeSjchaIDaVJiIiZ2h2H6hx48bhvffeQ01NjSvyQyQVvtCRbEQWGrCptFjJycnQarXw8/ODVqtFcnKy6CwREbVIs0fh69mzJ2w2GwYOHIjnnnsOPXr0gK+vb4N0sbGxTskgkTezv9A11pyJyBvZCw3qN1tVq9DAHUaik7WpdnJyMtavX4/U1FSMHz8eOTk5SE1NBQAsW7ZMbOaIiJqp2fNA+fv733yjGg3bkxM1kawvVCQnmecjkvm7a7VaLFiwwKH/9KpVq7B48WJYLBaBOSMiar5mB1C7du1qUrrBgwe3KENEROTdZC00kHnQGD8/P5w8eRLt2rVTlp0/fx5dunSB1WoVmDMiouZrdgBF5AqyvlARkTxkHgWQNVDixMbG4quvvoLNZoNGo8HDDz+Mjz76SHS2iDwaJ9Il4TgXEhHJQOZBY8aNG4fU1FSsWrUK58+fx6pVq5Camopx48aJzppXi42NRV5eHuLj43H06FHEx8cjLy9P1X7qnCyevFGLaqC++uorvPPOOzhy5AiqqqpgszluQqPRYO/evU7LpCxkrYWRuVkLEclD5j5QwNWBJN59911cvHgRbdu2xbhx4ziAhIv5+/sjPj4eK1asUJbNmTMH69evR1VVlcv3L/s5T96r2QFUdnY25s+fj06dOmHAgAHw8/NrNN2bb77plAzKQuabjMzNWohILrIWlJEYfn5+OHr0KO68805l2enTp9G9e3dV+p6xgJS8VbMDqL59+6Jnz5748MMP0aZNG1flSzoy32Rk/u5ERESuIroGigWk5K2a3QeqoqICo0ePZvDkZCInlxSNk1sSERE538MPP4ycnBzMmTMHp0+fxpw5c5CTk4OHH35Ylf3L3O+PvFuzA6h+/frhP//5j1N2/s033+DZZ59FeHg4/Pz8sGHDBmVdbW0tXn31VQwaNAhdunRBWFgYJkyYgGPHjjls4+LFi0hOTkavXr3QpUsXPPvsszhx4oRDmmPHjmHMmDHo0qULevXqhblz5+LSpUsOaXbt2oWhQ4dCq9Xi/vvvR05OToP8rl27Fvfddx+0Wi2GDh2Kb7/91inHAZD7JmMwGLBgwQIYjUZotVoYjUYhk1uykysREXmTjz76CMOGDcP69evRvXt3rF+/HsOGDVNtFD4WkJK3anYA9Yc//AHvvfce8vPzb3nn586dQ0REBNLT0+Hr6+uw7vz58/jHP/6BpKQkfP3113jvvfdw4sQJGAwGXL58WUmXkpKCrVu3Yt26ddi2bRvOnj2LMWPGoK6uDgBQV1eHMWPGoLq6Gtu2bcO6devwySefYP78+co2jhw5gmeeeQaRkZHIz8/HnDlzMHfuXGzZskVJs3nzZsybNw+JiYnIz89HZGQk4uLiGgR0LSX7TcZgMMBsNqOyshJms1n14ImjAIrBwJWIyLU++ugjVFVVwWq1oqqqStUhzN2hgJTIFW7aByouLq7BsiNHjuBf//oX7r77bnTt2hU+Pj6OG9Vo8MEHHzQrI3fddReWLl2KsWPHXjfNTz/9hKioKHzzzTe45557cPr0afTu3RurV6/GM888AwA4fvw47r33XuTm5mL48OH48ssv8cwzz+Cf//wnQkJCAAAbN27ESy+9hIMHD+KOO+7Aq6++iq1bt+L7779X9pWQkICffvoJX375JQBg+PDhuOeee/DHP/5RSfPggw8iJiYGr776arO+6/Wwc7EY7IMlhswDpxAREZHnanWzBD/99BM0Gk2D5SEhIbh48SIOHTrUYF1j6Z3h7NmzAKCM/Ld3717U1tYiOjraIV9hYWEoLCzE8OHDUVRUhLCwMCV4Aq4GQxcvXsTevXsxZMgQFBUVOWzDnub9999HbW0tbDYb9u7di4SEBIc00dHRKCwsdNr3MxgM0r44igweZe5/JlJGRgaysrKUwHXIkCHIysqC0WiU9jogIiIi93fTAOqf//ynGvm4qUuXLuGVV17BY489hrvuugsAUFZWBh8fH3Ts2NEhbWBgIMrKypQ0gYGBDus7duwIHx8fhzTXdqgMDAzE5cuXUVFRAZvNhrq6ugbbqb8farnr1UQAUOVF2t7/rH4NlCz9z0Ri4EpERESeqEl9oNLT04VOjHv58mVMmjQJp0+fbtL8UjabzaEW7Ho1YjdKY58c+GZpnFnbJmt/kPo1Ea1bt1ZqIjIyMlTZv+z9z0SReeAUIiIi8lw3rYECgL/85S9YunQpAgMD8cgjj+Cxxx7DsGHD0KFDB1fnD5cvX8bvf/977N+/H59++ikCAgKUdUFBQairq0NFRQU6deqkLC8vL8egQYOUNNc2s6uoqHCoUQoKCmpQk1ReXo5WrVohICAANpvNocaqfppra6VaSnQtjEiiayLsx9doNCpNCNkPx/XsgWtjfaCIiIiI3FWTaqB+/PFHFBQUYMqUKTh8+DDi4+Nx9913IyYmBtnZ2fjXv/7lkszV1tYiPj4e+/btw9atW6HVah3W9+vXD61bt0ZeXp6y7MSJEygpKYFOpwMAREZGoqSkxGFo87y8PLRt2xb9+vVT0nz11VcO287Ly8MDDzyA1q1bo02bNujXr5/Dfuxp7Pu5VaJrYURyh5oIkaMAyoqjMxEREZEnuukofI2xWq348ssv8cUXX+Dvf/87rFYrevbsiREjRmDEiBF46KGH0KrVzSu3qqurcfjwYQDAiBEjMGvWLIwcORL+/v7o3Lkzxo0bhx9++AHvv/8+OnfurHzujjvuUIY9nzNnDj7//HNkZ2fD398f8+fPh9Vqxddffw0fHx/U1dXhN7/5DTp27IglS5agqqoKU6dOxRNPPIFly5YBuDqq4KBBgzBu3DjEx8ejsLAQiYmJWLt2LWJiYgBcHcZ88uTJyMjIgE6nQ05ODv7yl7/AbDajW7duzT2EDcg8WzdHYyMiIiIiT9GiAKq+K1euoLCwEF988QV27NiBAwcO4Pbbb8fDDz+MGTNmIDIy8rqfLSgowJNPPtlg+XPPPYd58+bh/vvvb/Rzq1evVoY7r6mpwYIFC5Cbm4uamhoMGTIEGRkZDqPuHTt2DElJScjPz8dtt90Gg8GAJUuWoG3btkqaXbt24eWXX8ZPP/2E4OBgzJo1C+PHj3fY79q1a5GZmQmLxYLw8HC8/vrreOihh5p1vK5H9qG0OYQ7yYbnPBERkWe65QDqWseOHcOOHTvw5ZdfYvDgwQ2G/qbGsRaGSB683omIiDyX0wMoajmWSBPJQfYaZyIiIk/W4j5QH374IY4cOYKqqiplyG9loxoNVq1a5bRMEhF5E5n7PBIREXm6Jg1jXt/XX3+N3/3udzh79ixuv/12+Pn5NUjjzLmRiIi8DSdvJiIi8lzNDqBefvll+Pv747PPPsO9997rijwREXk1zoFFRETkuZodQB06dAivvfYagyciohbi5M1ERESeq0kT6dbXvXt31NTUuCIvRETS4OTNRETkKrm5udDr9QgICIBer0dubq7oLHmVZgdQc+bMwfr161FVVeWK/BARERERUQvZp8owmUywWCwwmUxYvHgxgygnuukofJmZmQ2Wffzxxzhy5AhiY2MREhICHx8fx41qNHjppZecm1MiIiIiIrohTpXhejcNoPz9/Zu/UY2GQ/ESEVGjOOcdEZHrcKoM17vpIBL/+Mc/1MgHERFJwN605NoRCAEwiCIicgLZp8pQo5CuRRPpEhERtQSblhARudb1CqpkGO1Vre/epABq7NixGDRoEKKiotCvX78GfZ6IiIiagk1LiIhcT9am0moV0jUpgAoODsbFixeh0Wjg6+uLAQMGICoqCoMGDcLAgQPRrl07p2WIiIi8F2ugiIjIVdQqpGvSMObHjh3Djh07kJqaiqFDh2Lfvn1YunQpYmNj0b17d0RHR2P+/PnYunUrysvLnZY5IiLyLomJiUhISEB+fj5qa2uRn5+PhIQEJCYmis6aKjg3CxGR69j7f9Xniv5fLe4D9fPPP8NsNmP37t0oLCzEv//9b2g0GgBA7969UVRU5NSMEhGRd5C1aYnM/RKIiNTgVn2gbqaurg5ffPEFMjMzUVhYyGHMiYiIrsHmi0RErue2o/CdO3cOe/bsUWqg/u///g/nz59Hx44dERkZiaioKGVYWiIiIuIAGkRE3uKm80ABgMViwe7du5WAad++fairq0NoaCh0Oh2eeeYZREVF4e6773Z1fomIvIKszdhkJvvcLERE3qJJAVR4eDhatWqFgQMHIjo6GikpKdDpdPDz83N1/oiIvA4nk5WTfQCNxtrmExGR52jWMOZdu3ZFVFSU8i8iIkKNPBIReRX2hZEXax6JiDxfkwKo2tpa/PDDD9i9ezd2796NoqIiVFRU4I477sDAgQOh0+kQFRWFAQMGwNfXV418ExF5LPaFISLyfiww8V5NasLXunVrREZGIjIyEi+99BIAoKSkRBnC/P3338frr7+OVq1a4d5770VUVBRef/11l2aciMhTsS8MEZF3Y1Nt78ZhzImIVMb5gIiIvBubans3pw5jbrPZ0K5dO/Tv3x+ffPKJK/JLROQV2LSDiMh7sam2d7ulYcxtNhs6duyIoUOHQq/XQ6/X4/7770erVk3aLBGRtAwGAwMmIiIXEllQJXtTbW8vJGxSpNOnTx9oNBrYbDb06NEDBoNBCZhCQ0NdnUciIiIioiYT3QdJ5mkLRB97NTSpCd/cuXMxaNAg6PV6XLp0CV27dr1h+p07dyI6OtppmSQiIufy9tJBIpKbO/RBkvU+6w7H3tWa3QfqwQcfxOeffw6tVtvo+k8//RQTJkzAqVOnnJJBIiJyLg5iQSQPWV/i2QdJHBmO/a+a+4G2bdsiJiYGFRUVDdbl5uYiPj6+ybVP33zzDZ599lmEh4fDz88PGzZscFhvs9mQlpaGPn36IDg4GI8//jgOHDjgkMZqtWLSpEno1q0bunXrhkmTJsFqtTqk2bdvH0aNGoXg4GCEh4fDZDLBZnOMG7ds2QKdToegoCDodDps3bq12XkhIvIEGRkZyMrKwpAhQ9C6dWsMGTIEWVlZyMjIEJ01InIie2GJyWSCxWKByWTC4sWLkZubKzprLmfvg1SfTH2QRJLh2Dc7gPr4449RW1uLmJgYh0Dl3XffxeTJk/HUU0/hz3/+c5O2de7cOURERCA9Pb3RCXgzMzOxevVqmEwm7Ny5E4GBgYiNjcXZs2eVNBMmTEBxcTE2bdqE3NxcFBcXY/Lkycr6M2fOIDY2FkFBQdi5cyfS09ORlZWFVatWKWmKioowfvx4xMXFoaCgAHFxcXjxxRfx3XffNSsvRORZcnNzodfrERAQAL1eL8VLBXB1Hj+9Xu+wTK/Xo6SkRFCOiMgVZC4ssfdBys/PR21tLfLz85GQkIDExETRWfN6Mhz7Fg1jfvz4cYwcORJBQUH4+OOP8ec//xnz58/H2LFjkZWVBY1G0+yM3HXXXVi6dCnGjh0L4GqNT58+fTBx4kQkJSUBAC5cuIDQ0FAsXrwY8fHxKCkpgU6nw/bt2xEVFQXgaoQ7cuRI7NmzB6GhoVi3bh1SU1Px888/K0HasmXLkJOTg/3790Oj0SA+Ph5VVVX4+OOPlfzExMSgU6dOWLduXZPyQkSeReZmbDK0TyciOZpS3YiszRfdgbcf+2bXQAFASEgIPvnkE5SWluKhhx7CK6+8gokTJ2LVqlUtCp4ac/ToUVgsFofmgL6+vhg0aBAKCwsBXK056tChA3Q6nZImKioK7du3d0ij1+sdariGDx+O0tJSHD16FACwZ8+eBs0Ohw8frmyjKXkhIs/CklnvLh0kIjmaUt2IwWCA2WxGZWUlzGazV73A34zoFhbefuxvOoz5L7/80ujyDh064K233sLzzz+P559/HsnJyQ5pAwMDbyljFoul0e0EBgaitLQUAFBWVoaOHTs6BG0ajQadOnVCWVmZkqZLly4NtmFf16NHD1gslkb3Y99GU/JCRJ5F5mZs9geZ0WhUSgdlqHkjko3MQ2nLTIZhxEW7aQD161//+oa1SjabDe+99x7ee+89h+XOqhq+dt82m61BwNRYnm6Uxj6AxM3SXLusKWmIPI23V7Nfj+yTHHIiXyLvx8ISOdVvYQFAaWFhNBr52zvJTQOouXPnCgkS7MOkl5WVISQkRFleXl6u1AQFBQWhvLzcIZCx2WyoqKhwSGOvSaq/DeC/NUparbbRNPXX3ywvRJ5I5lIqlswSkQxYWCIfmVtYqOWmAVRKSooa+Wige/fu0Gq1yMvLw4MPPggAqKmpgdlsxqJFiwAAkZGRqK6uRlFRkdIPqqioCOfOnVP+joyMRGpqKmpqanDbbbcBAPLy8tC5c2d0794dADBw4EDk5eXhpZdeUvafl5enbKMpeSHyRDKXUrFkloiIvJHsLSzU0KJBJJyluroaxcXFKC4uxpUrV3D8+HEUFxfj2LFj0Gg0mDp1KlauXIlPPvkE+/fvx7Rp09C+fXvlBScsLAyPPPIIZs+ejT179qCoqAizZ8/GiBEjEBoaCuDqS5Kvry+mTZuG/fv345NPPsHKlSsxbdo0pdZqypQpyM/Px4oVK/Dzzz9jxYoVKCgowNSpUwGgSXkh8kSyl1J5eydXIiKZiR5IQRR3GCjI24/9TWugXOmHH37Ak08+qfydlpaGtLQ0PPfcc8jOzsbMmTNx4cIFJCcnw2q1on///ti8eTNuv/125TNr1qyB0WjE6NGjAQAjR47E0qVLlfV33nknPvroIyQlJWHYsGHw8/PD9OnTMWPGDCWNTqdDTk4OlixZgrS0NPTs2RM5OTkYMGCAkqYpeSHyNCylIiIibyRzE3XRLSxkOPYtmgeKyNvIOpCCzHMhERGR95J9vjuR7zUyHHuhNVBE7kCGkpLrEV1KRSQbWQtriNQmcxN10e81Mhx7oX2giNyBzBOqAuwHRKQW+0uNyWSCxWKByWTC4sWLva5vAJE7kHkSYdHvNTIcewZQJD0ZSkqISDzRLzVEMnGHgRREEf1eI8OxZxM+kh4HUiAiNYh+qSGSiTs0URfVZFf0e407HHtXYwBF0uOEqkSkBtEvNUSyETmJsMh+SO7wXuPtEzhzFD4isGM3EbkeR70kkofokej4XuNaDKCIBONNjkgevN6J5BAQEACLxYLWrVsry2pra6HValFZWSkwZ+QMbMJHJJDooUaJSF3e3qyFiK5ik13vxlH4iATiqFxERETeR4aR6GTGJnxEArGKn4iIyDuxya73YhM+IoFYxU9EROSd2GTXe7EJH5FArOInIiIi8iysgSISSIbJ5oiIiIi8CftAERERERERNRGb8BERERERETURAygiIiIiIqImYgBFRERERETURAygiIiIiIiImogBFBERERERURMxgCIiIiIiImoiBlBERERERERNxACKiIiIiIioiRhAkVvIzc2FXq9HQEAA9Ho9cnNzRWeJJMDzjoiIiJqLARQJl5ubi8WLF8NkMsFiscBkMmHx4sV8mSWX4nknJwbNRER0qzRWq9UmOhMkN71eD5PJhCFDhijL8vPzYTQaYTabBeaMvBnPO/nYg+asrCzo9XqYzWYkJCRgwYIFMBgMorNHREQeggEUCRcQEACLxYLWrVsry2pra6HValFZWSkwZ+TNeN7Jh0EzERE5A5vwkXBhYWENXl7MZjPCwsIE5YjUIrI5lejzjk3J1FdSUgK9Xu+wTK/Xo6SkRFCOiIjIE7l1AFVXV4clS5bgvvvug1arxX333YclS5bg8uXLShqbzYa0tDT06dMHwcHBePzxx3HgwAGH7VitVkyaNAndunVDt27dMGnSJFitVoc0+/btw6hRoxAcHIzw8HCYTCbYbI6Vc1u2bIFOp0NQUBB0Oh22bt3qui8vkcTERCQkJCA/Px+1tbXIz89HQkICEhMTRWeNXEh0HySR553o7y4r0UEzERF5B7cOoFauXIm1a9fCZDKhqKgI6enpWLNmDVasWKGkyczMxOrVq2EymbBz504EBgYiNjYWZ8+eVdJMmDABxcXF2LRpE3Jzc1FcXIzJkycr68+cOYPY2FgEBQVh586dSE9PR1ZWFlatWqWkKSoqwvjx4xEXF4eCggLExcXhxRdfxHfffafOwVCBqBJxg8GABQsWwGg0QqvVwmg0sk+CikT97hkZGcjKysKQIUPQunVrDBkyBFlZWcjIyFBl/yLPO9HfXVYsrCEiImdw6z5QY8aMgb+/P/70pz8py6ZMmYKqqips3LgRNpsNffr0wcSJE5GUlAQAuHDhAkJDQ7FPU72iAAAgAElEQVR48WLEx8ejpKQEOp0O27dvR1RUFICrJY4jR47Enj17EBoainXr1iE1NRU///wzfH19AQDLli1DTk4O9u/fD41Gg/j4eFRVVeHjjz9W8hITE4NOnTph3bp1Kh4V18jNzcW8efPQrl07HDt2DF27dsX58+eRnp7OQMaLiexUHxAQgOzsbKxcuRIlJSUICwvDrFmzMHXqVK/vg8T+V+Lk5uYiIyNDOecSExN5jyMiomZx6xqoqKgo7Nq1Cz///DMA4KeffkJBQQH+53/+BwBw9OhRWCwWREdHK5/x9fXFoEGDUFhYCOBqzVGHDh2g0+kcttu+fXuHNHq9XgmeAGD48OEoLS3F0aNHAQB79uxx2I89jX0bnm7hwoXw8fHBqlWrUFZWhlWrVsHHxwcLFy4UnTVyIZE1IcHBwZg7dy7OnTsHm82Gc+fOYe7cuQgODnb5vkVjUzJxDAYDzGYzKisrYTabGTwREVGzuXUANWvWLIwZMwY6nQ6dOnVCVFQUnnvuOUyYMAEAYLFYAACBgYEOnwsMDERZWRkAoKysDB07doRGo1HWazQadOrUySFNY9uwr7Pv60b78XQnT55Edna2w4t0dnY2Tp48KTpr5EIiO9VfuHABZ86cQU1NDTQaDWpqanDmzBlcuHDB5fsWjU3JiIi8HwcL8l6tRGfgRjZv3oy//vWvWLt2Lfr06YN//vOfmDdvHrp164Zx48Yp6eoHR8DVgSWuDZiudbM09gEkbpamsW17qmu/izd9N2qcvSak/rDOatWEVFVV4Y477sBtt90Gm82G2267Dbfffjuqqqpcvm/R7LUeRqNRaUqmdr8/NmUjInKd6zWRB8B7rRdw6xqohQsXYsaMGXj66adxzz334Nlnn8X06dPxxhtvAAC0Wi0ANKgFKi8vV2qLgoKCUF5e7jCins1mQ0VFhUOaxrYB/LcmSqvV3nA/nq5Lly6YMmWKQ4n4lClT0KVLF1X2z1IaMUTXhCQlJaG4uBhVVVUoLi5W+jLKQGRTMo4CSETkWhwsyLu5dQB1/vx5+Pj4OCzz8fHBlStXAADdu3eHVqtFXl6esr6mpgZms1np8xQZGYnq6moUFRUpaYqKinDu3DmHNGazGTU1NUqavLw8dO7cGd27dwcADBw40GE/9jT1+1Z5skWLFqGurg4zZsyAVqvFjBkzUFdXh0WLFrl833yZE0f0CIirVq1yCN7qj3xJrsMHOxGRa3HeOe/m1qPwTZ06FV9//TXeeOMN9OnTB8XFxZg1axaeffZZ/OEPfwBwdajzjIwMrF69Gr1798by5cvx7bffYs+ePbj99tsBXH1JPHnyJDIzM2Gz2TBr1ix07doVGzduBACcPn0aAwcOxODBg5GUlIRDhw5h+vTpmDt3rlLdWlhYiFGjRmH+/Pl44okn8Omnn+L111/H9u3bMWDAADEHyMmSk5Px7rvv4uLFi2jbti3GjRuHZcuWuXy/er0eJpPJoRlZfn4+jEZjg4725D0iIiJQXV0NPz8/HD9+HCEhIbBarejQoQP2798vOntejaMAEhG5Ft9tvJtb10AtXboUTz31FBITE6HT6fDKK6/ghRdewIIFC5Q0M2fOxLRp05CcnIxhw4bh1KlT2Lx5sxI8AcCaNWvQt29fjB49Gk8//TT69u2Lt956S1l/55134qOPPkJpaSmGDRuG5ORkTJ8+HTNmzFDS6HQ65OTk4P3338dDDz2Ev/71r8jJyfGa4Ck3NxdffPEFNm3ahF9++QWbNm3CF198oUotEEtp5LRo0SK0adMGwH/7HLZp00aVWk/ZcRRAIiLXEt1EnlzLrWugSD0iS0pYSiMvDmQghsj5v4iIZMFnnPdiAOVGRF5oIpv08GWOSH18sBOph9cbkXdx62HMZSJ6uEuRw1m7w5DORLIxGAy8xohUIPr5TkTO59Z9oGQielQs0W11RQ7pTETy4JQJpDbRz3eZ8XonV2ENlJsQPZACa4GIyNuxJoBEEP18lxWvd3Il1kC5CXcYFYu1QETkzVgTQCKEhYXBZDI51ISYTCaOeulivN7lpUbNIwMoNyG6CR0RkbdjTQCJMHjwYKxcuRJjx47F8ePHMXbsWKxcuRKDBw8WnTWvxutdTvaaR5PJBIvFApPJhMWLFzs9iGIA5SYMBgMWLFgAo9EIrVYLo9HIJnRERE7kDjX9JJ9du3Zh5syZ2LBhA0JCQrBhwwbMnDkTu3btEp01r8brXU5q1TxyGHMiIpICp0wgEUROEyIzXu9yUut64yASREQkBQ6WQyKInCZEZrze5aTW9cYaKCIiIiIXYU2IvDiBsvrUut5YA+VGeKGJw2NPRESuwJoQOXEYdTHUut5YA+VC/xrzryanLS0txTf//gb6DfpGI+bj847j4r8vtigfd2+8u9Hlv7z9C878/UyLthmSFoK2vdo2WH7mb2fwy5pfWrTNwImBuOOROxosv3j4Io6nHG/RNu8YfgcCJwU2us7++5SWluLgwYPo27cv/P39UVVVhR9//BGhoaHo3Llzg8+17dkWIekhjW6zJb9TaWkpDh8+jFnnZzUavLn6d7Lvv7q6Gh06dECvXr0a/d52on6n5nL272TnLddTc393gL+TN933mou/E3+nxvB3uv7v9M033yA8PBwBAQHKssrKShw4cADR/180fyc3+Z2u53rH0o41UG7i8OHDeHjYw9AN0QGAMmqI0WhkSYWLHT58GH379lVucgEBAejbty8OHDhw0xfKW1U/eLNss6heQnW94BGAy787icPfXU65ubk4/835ZgXNRNQy1dXV8Pf3d1jm7++P6upqQTkiZ+Iw5m6iurq6wYOM8xWoQ+RNrn7wJmKiv/r712g0SvB4+PBhVfZPYvB3l4+9OVF4eDgeffRRhIeH4+DBgygtLRWdNfJypaWl+Oabb7Bjxw5888030pxzHTp0QFVVlcOyqqoqdOjQQVCOyJkYQLmJDh06NLipcJQedYi8yTUWvKkZOLOETE4y/+4/7P1Bypc5+9woDJrlk5ubC71ejz/96U+qn/P22m4ZA/devXrhxx9/RGVlJWw2GyorK/Hjjz+iV69eorPm9eznfEBAAPR6vUvON/aBchMcpUcckcder9fDZDI5DLeZn58Po9HYYAJAb9w/iSHr7y7zfTYgIABvvvkmMjMzlY7VM2fOxLRp0zgXkRcTfc7Leq+x4wBV6lPrnGcA5UZ4oYkj6tiLfriJ3r9osl5zon93Ucdd5pe5iIgI1NXVYc2aNcpvPnHiRPj4+GD//v2is0cuIvqc5yTCpDa9Xo9Ro0Zh27ZtyjPG/rczz3kOIuFGDAaDFC9v7kjUsRc9vK3o/Ysk8xCzIn93kce9pKQEer3eYRn7mpIaRBUaiD7nOYkwqe2nn37C+fPnGzxj/vOf/zh1P+wDRSSYwWCA2WxGZWUlzGaz6i/vovcvir1PyJAhQ4QM4CGaqN9d5HG3v8zVJ8vL3KlTp7Bo0SIYjUZotVoYjUYsWrQIp06dEp01r2cvNDCZTLBYLDCZTFi8eDFyc3Ndvm/R53xiYiISEhKQn5+P2tpa5OfnIyEhAYmJiars/9q+MGocc7pK1LFv06YNJk6c6PCMmThxItq0aePU/bAGioikJLpkVlYij7v9Za6xpoveLiwsDJ07d3Z4mc7Pz5cieBStfqEBoO40JaLPeVlru2Un8thfunQJb7/9Nu677z5l32+//TYuXbrk1P2wBoqIpCS6ZFZWIo+7wWDAggULHGphZGmyKromQGYlJSUoLS1tMCqYGoUGBoMBjz76KOLi4hAYGIi4uDg8+uijqp7zMtZ2y07kse/Tpw+eeeYZh/v8M888gz59+jh1PxxEgoikJHogBVnxuIsj66ApookcwEPm640DWIgj8tirdc6zCR8RSUnmATRE4nEXhwMVyUdk80HROICFOCKPvVrPGNZAEREJwNoAIjmInINL5loYmWvfRJPh2LMGiohIZe7QuZkBHJE6wsLC0KVLFyEDeMhcC8PabnFkOPYcRIIIHOqU1CW6c7PIYZVJHN7nxBA5gIfsg4fIOk2HO/D2Y88mfCQ9Gaqayb2Iblaj1+thMpkcSqXz8/NhNBqdOlM7NSSq5o/3ObFE1vjKXNss83cn12INFLkFkSWjomsDSD6ih1AXOayyzETW/LnDfU7mGjCRpfHeXhNwPaxpJ1dy+wDq1KlTmDJlCu6++25otVrodDrs2rVLWW+z2ZCWloY+ffogODgYjz/+OA4cOOCwDavVikmTJqFbt27o1q0bJk2aBKvV6pBm3759GDVqFIKDgxEeHg6TyQSbzbFybsuWLdDpdAgKCoJOp8PWrVtd98UlIvomxwlVSW2im9UEBwdj4cKFDtfcwoULERwcrMr+ZSUyiBF9nxN9nyf5uEOhAXkvtw6grFYrRowYAZvNhg8++ACFhYVYunQpAgMDlTSZmZlYvXo1TCYTdu7cicDAQMTGxuLs2bNKmgkTJqC4uBibNm1Cbm4uiouLMXnyZGX9mTNnEBsbi6CgIOzcuRPp6enIysrCqlWrlDRFRUUYP3484uLiUFBQgLi4OLz44ov47rvv1DkYXkz0TU50bQDJR+YJXWUmMogRfZ8TfZ8n+YguNCDv5tYB1B//+EcEBwfjrbfeQv/+/dGjRw8MHTpUueHbbDZkZ2dj1qxZiImJQUREBLKzs1FdXa2UapWUlOBvf/sbVq5cCZ1Oh8jISLzxxhvYsWMHDh48CADYtGkTLly4gOzsbERERCAmJgYzZ87Em2++qdRCZWdn4ze/+Q2SkpIQFhaGpKQkDB48GNnZ2WIOjhcRfZMTXRtA4ohsUiSyWc2pU6fw2muvOQRwr732Gk6dOqVaHmQkMogRfZ8TfZ8n+YguNADkbrbq7dw6gPrss8/Qv39/xMfHo3fv3hg8eDDefvttJag5evQoLBYLoqOjlc/4+vpi0KBBKCwsBHC15qhDhw7Q6XRKmqioKLRv394hjV6vh6+vr5Jm+PDhKC0txdGjRwEAe/bscdiPPY19G9Ryom9yrA2Qk8xNiuoPq2wP4Lp06cJaVxcTGcSIvs+Jvs+TOKKCCNGFBjI/Y2Tg1gHUkSNHsG7dOvTo0QMffvghpkyZgtdeew1r1qwBAFgsFgBwaNJn/7usrAwAUFZWho4dO0Kj0SjrNRoNOnXq5JCmsW3Y19n3daP9UMuJvskB8nayFY2Dh4jhDtecjEQHMSLvczzn5CQyiBB9vcn8jJGBW0+ke+XKFTzwwAN49dVXAQD3338/Dh8+jLVr12LSpElKuvrBEXC1ad+1AdO1bpbGXst1szSNbZuaR4YJ16gh0ZPJytykiNccqY3nnJzqBxEAlCDCaDSq8tsbDAZh55jMzxgZuHUApdVqG1Tv//rXv8bx48eV9cDVWqKQkBAlTXl5uVJbFBQUhPLycodgx2azoaKiwiHNtTVJ5eXlAP5bE6XVahtNc22tFLWMyJsciSH6wWpvUlR/LiSZmhTxmlOf6EID0XjOyUfmIEL2Z4y3c+smfFFRUTh06JDDskOHDqFr164AgO7du0Or1SIvL09ZX1NTA7PZrPR5ioyMRHV1NYqKipQ0RUVFOHfunEMas9mMmpoaJU1eXh46d+6M7t27AwAGDhzosB97mvp9q4io6UQ/WNmkiNTGJj0kG5n7vvEZ493cugZq2rRpePTRR7F8+XKMHj0axcXFePvtt7FgwQIAV5vUTZ06FRkZGQgNDUXv3r2xfPlytG/fXinlCgsLwyOPPILZs2cjMzMTNpsNs2fPxogRIxAaGgrgaqmYyWTCtGnTkJSUhEOHDmHlypWYO3euUms1ZcoUjBo1CitWrMATTzyBTz/9FAUFBdi+fbuYg0Pk4USXzrFJEalNdKEBkdrsQcS1ta729zhvxmeMd9NYrVbbzZOJs2PHDixatAiHDh1CSEgIJk6ciMmTJzs0x0tPT8f//u//wmq1on///li+fDkiIiKUbVRVVcFoNOLzzz8HAIwcORJLly6Fn5+fkmbfvn1ISkrC999/Dz8/P8THx8NoNDr0cdqyZQuWLFmCI0eOoGfPnnjllVfw1FNPqXQkiLzL9Zoz8QFD3kqv18NkMjkUGuTn58NoNDYopSfyFrm5ucjIyFCCiMTERN7jyeO5fQBFRN6LD1aSCQsNxOG9hoiciQEUkWB8sBPJg9e7+hi4EpGzMYAiEogPdiIi12LTSSJyNrcehY/I23FULhJF5CTGRGri4B1E6vP2ZwwDKCKBZH+we/sN1l3Zaz5NJhMsFgtMJhMWL17M4+/lZL3eZB5Km8SS9ZqT4RnDAIpIIJkf7DLcYN0Vaz7lI/P1xvl4SASZrzkZnjHsA0UkkDv0gRLVqZ39EsQJCAiAxWJB69atlWW1tbXQarWorKwUmDNyFdmvNw7eQWqT+ZoT/YxR43pnDRSRQAaDAQsWLIDRaIRWq4XRaFQ9eBJVQuYOzRdlbV4hc82naKLOOXe43kQyGAwwm82orKyE2Wxm8EQuJ/M1J/IZo9Z7DQMoIsFEPthFVrOLfomXuXkFmzSJIfKcE329EclG5mtO5DNGrfcaBlBEgomsBRFZQib6JV6GNtrX4w41nzLW/Ik850Rfb0SykfmaE/mMUeu9ppVTt0ZEzXK9PlAAVLnR2EvI6rfRVquEzP79jEaj0k5ZzZd4mZtXAFePv4hmTCLP+bS0NJhMpkbXGY1GpKSkuHT/Is850dcbkWxkv+ZEPWPUeq9hDRSRQKJrQUSXkIlsviiqeUVaWhr8/Pyu+y8tLc2l+xdN5DmfkpICq9UKq9UKvLFP+b/VanV58ASIb9Ijcz8gWWs9SQz7fX7ChAk4cOAArly5ggMHDmDChAlS3OdFUuu9hqPwEQkkeqQaQN7RqdxhBEQA8Ft/Atb4u1Tbn2jucM4DYo67u5xzssnNzcW8efPQrl07HDt2DF27dsX58+eRnp7O406qkO0+L5oa7zVswkckkMgmdHaiqtlFk715hSjucM6LwnNOjIULF8LHxwerVq1SAteJEydi4cKFXj9dBJGM1HivYRM+IoFEN6EjUpvs57zMzehEOXnyJLKzsx2ajWZnZ+PkyZOq7F/mET+JvBVroIgEYom0OKIH8JAVz3kSQaPR3PBvV6rf7w+A0u/PaDTyvCfyUOwDRURScpdZ4tk2Xgwed3lERESgrq4Oa9ascWjC5+Pjg/3797t8/+7S74/E4f3G+7AJHxFJSfZhzIlksWjRItTV1WHGjBnQarWYMWMG6urqsGjRIlX2L3r0RSJyPgZQRIJxeF0x+FJDJAeDwYD09HS0b98eANC+fXtVR+CTvd8fkTdiAOVG+CItH3YuFocvNUTyEDl4h8FgwIIFC2A0GqHVamE0Gtnvj7yet7/TchAJN8EO7XJi52JxOJgBEalF1ukiSE4yvNOyBspN1H+Rtg+zmpWVhYyMDNFZIxdiPxwiIiLyJjK807IGyk3wRVpOYWFhMJlM+Oyzz5RakMcff5z9cFQgQwkZERGR2mR4p2UNlJsQ2aE9LS0Nfn5+1/2Xlpbm8jzIavDgwVi5ciXGjh2L48ePY+zYsVi5ciUGDx4sOmteT4YSMnfl7W3jiYhEE3mflWGQJgZQbkJkh/aUlBRYrVblH97Y5/B3SkqKy/Mgq127dmHmzJnYsGEDQkJCsGHDBsycORO7du0SnTWvJ0MJmTviwClERK4l+j4rwyBNnEjXjeTm5iIjI0NpypWYmCikKZGaE76lpaXBZDJdd73RaPTqAI4TLIrDiXTF4HEnItmofb9xh/usu7zTugoDKGpA5IuFbC817nCTk9X1+kCpPRKfbOe8uxQayHbciUgcte837nKf9WZswkckkAzV3O6Kc7OIIUPbeCJ3wP7N8uJ91vUYQBEJxJd4sUROrikrFhoQqcNd+jdz0Bj18T7reh4VQGVkZMDPzw/JycnKMpvNhrS0NPTp0wfBwcF4/PHHceDAAYfPWa1WTJo0Cd26dUO3bt0wadKkqzeTevbt24dRo0YhODgY4eHhMJlMsNkcWzdu2bIFOp0OQUFB0Ol02Lp1q+u+LEmDL/EkExYaEMlD9GAGsjIYDHj00UcRFxeHwMBAxMXF4dFHH+V91ok8JoDas2cP3nnnHdxzzz0OyzMzM7F69WqYTCbs3LkTgYGBiI2NxdmzZ5U0EyZMQHFxMTZt2oTc3FwUFxdj8uTJyvozZ84gNjYWQUFB2LlzJ9LT05GVlYVVq1YpaYqKijB+/HjExcWhoKAAcXFxePHFF/Hdd9+5/ssTkVPdqGkLm7W4HgsNiOTA6SLEyM3NxRdffIFNmzbhl19+waZNm/DFF18wcHUijxhE4vTp0xg6dCgyMzOxdOlSREREYNmyZbDZbOjTpw8mTpyIpKQkAMCFCxcQGhqKxYsXIz4+HiUlJdDpdNi+fTuioqIAXG0HOnLkSOzZswehoaFYt24dUlNT8fPPP8PX1xcAsGzZMuTk5GD//v3QaDSIj49HVVUVPv74YyVfMTEx6NSpE9atW6f+QXEhDiJBsuE5Lx8ed5KNiHOegxlcJeMofN7OI2qgZs2ahZiYGAwdOtRh+dGjR2GxWBAdHa0s8/X1xaBBg1BYWAjgas1Rhw4doNPplDRRUVFo3769Qxq9Xq8ETwAwfPhwlJaW4ujRowCu1oDV3489jX0bRERERPRfHMxADM5z6HpuH0C98847OHz4MObPn99gncViAQAEBgY6LA8MDERZWRkAoKysDB07doRGo1HWazQadOrUySFNY9uwr7Pv60b7Ic/EUYqIiIhcg4MZiCFr4KrmO51bB1AHDx7EokWLsGbNGrRp0+a66eoHR8DVgSWuDZiudbM09gEkbpamsW2T53CXUYqIZMJRudTFgiIShYPGiCFr4KrmO10rp23JBYqKilBRUeFQDVlXV4dvv/0WOTk52L17N4CrtUQhISFKmvLycqW2KCgoCOXl5Q7Bjs1mQ0VFhUOaa2uSysvLAfy3Jkqr1Taa5tpaKSIiur7rTWAMgC9VLpKSkuLw4sD+X6Qmg8HAa1tl9uNtNBpRUlKCsLAwBq5O5tY1UI8//ji+/fZbFBQUKP8eeOABPP300ygoKEDv3r2h1WqRl5enfKampgZms1np8xQZGYnq6moUFRUpaYqKinDu3DmHNGazGTU1NUqavLw8dO7cGd27dwcADBw40GE/9jT1+1YREdGNyToqF2uB5JacnAytVgs/Pz9otVqH6VjI+Xi9cbRTV3PrGij7iV5fu3bt4O/vj4iICADA1KlTkZGRgdDQUPTu3RvLly9H+/btlRMlLCwMjzzyCGbPno3MzEzYbDbMnj0bI0aMQGhoKICrJ5nJZMK0adOQlJSEQ4cOYeXKlZg7d65SazVlyhSMGjUKK1aswBNPPIFPP/0UBQUF2L59u4pHhIjIs8nauZm1QPJKTk7G+vXrkZqaivHjxyMnJwepqakAro74S87H641cza1roJpi5syZmDZtGpKTkzFs2DCcOnUKmzdvxu23366kWbNmDfr27YvRo0fj6aefRt++ffHWW28p6++880589NFHKC0txbBhw5CcnIzp06djxowZShqdToecnBy8//77eOihh/DXv/4VOTk5GDBggKrfl4jIk8nauZnEcIeaiHfffRepqamYMWMG2rVrhxkzZiA1NRXvvvuuy/dNpDZ3uObU4BHzQJG6ZJ4TR/T+SQyZz3m1Xa8PlNrt80Ufd55zYoj47n5+fjh58iTatWunLDt//jy6dOlytaO7Wvng7y7t/kXy1nudWzfhIyIi78LOzSSbtm3bIicnx6FVS05ODtq2bSswV0R0KxhAERGRqjgqF8lk3LhxSp+n+n2g4uPjxWaMiFrM4/tAeQNZ2ouSI/7uRETeb9myZYiPj8fixYvRpUsXLF68GPHx8RxAgsiDsQbKDXC0GDnxdyciksOyZcsYMBF5EdZAERERERERNREDKCIiIiIioiZiAEVEqmP/LyIiIvJU7ANFRKpj/y8iIiLyVKyBIiIiIiIiaiIGUERERETkVGyqLQaPuzrYhI+IiIiInIpNtcXgcVcHa6CIiIiIiIiaiAEUERERERFREzGAIiIiIiIiaiIGUEREErlRB2N2LiYiIro5BlBERBJJSUmB1WpV/uGNfcr/63c8dgWODkVERN6Ao/AREZEq1B4dqseGk7Besl13vd/6E9df10aDI2O7uCJbRETk4RhAERGRV7JesrU4QLtRcEVERHJjEz4iIiIiIqImYgBFRERERETURAygiIiIyGU4eAgReRsGUEREROQyNxr5UY3RH0Vi8EjknTiIBBEREZELqD3yJBGpgzVQRERERERETcQAioiIiIiIqInYhI+kcrOJNYHrz//CiTWJiIiIiAEUSYUTa/5XWloaTCbTddcbjUav7txNRERE1BIMoIgkxc7N6riVWk+ANZ9ERO6O93n5uHUAtWLFCmzduhWHDh1CmzZtMGDAALz66quIiIhQ0thsNqSnp+Odd96B1WpF//79sXz5coSHhytprFYr5s6di+3btwMAHnvsMSxduhR+fn5Kmn379iE5ORnff/89/P398eKLL2Lu3LnQaDRKmi1btuD111/Hv//9b/Ts2ROvvPIKnnzySRWOhHPxQidSz63UegLeV/NJRORteJ+Xj1sHULt27cLvf/97PPjgg7DZbHj99dfx29/+FoWFhfD39wcAZGZmYvXq1Vi9ejVCQ0OxdOlSxMbGYs+ePbj99tsBABMmTMDx48exadMmaDQavPTSS5g8eTI2btwIADhz5gxiY2MxaNAg7Ny5EwcPHsT06dPRrl07JCQkAACKioowfvx4pKSk4Mknn8TWrVvx4osvYseOHRgwYICYA9RCvNCJiIiIiFrGrQOozZs3O/z91ltvoVu3bti9ezdGjhwJm82G7OxszJo1CzExMQCA7OxshIaGIjc3F/Hx8SgpKYsabFwAAB7CSURBVMHf/vY3bN++HTqdDgDwxhtvYOTIkTh48CBCQ0OxadMmXLhwAdnZ2fD19UVERAR+/vlnvPnmm5gxYwY0Gg2ys7Pxm9/8BklJSQCAsLAwFBQUIDs7G+vWrVP3wBARERERALasIfW5dQB1rerqaly5ckVpenf06FFYLBZER0craXx9fTFo0CAUFhYiPj4eRUVF6NChgxI8AUBUVBTat2+PwsJChIaGoqioCHq9Hr6+vkqa4cOH4w9/+AOOHj2KHj16YM+ePZg0aZJDfoYPH463337bxd+aiIiIiK6HLWtIbR41D9S8efNw7733IjIyEgBgsVgAAIGBgQ7pAgMDUVZWBgAoKytDx44dHfoyaTQadOrUySFNY9uwr7Pv60b7ISIiIiIi7+cxNVAvv/wydu/eje3bt8PHx8dhXf3gCLg6sMS1AdO1bpbGZrM1WH6z/dDNsZqdiGQg+l53s/3zPktE1HIeEUClpKRg8+bN2Lp1K3r06KEs12q1AK7WEoWEhCjLy8vLldqioKAglJeXOwQ7NpsNFRUVDmmurUkqLy8H8N+aKK1W22iaa2ul6MZYzS4n0S+TRGoTfa/jnHfy4X2WZCPynHf7AMpoNGLz5s349NNP8etf/9phXffu3aHVapGXl4cHH3wQAFBTUwOz2YxFixYBACIjI1FdXY2ioiKlH1RRURHOnTun/B0ZGYnU1FTU1NTgtttuAwDk5eWhc+fO6N69OwBg4MCByMvLw0svvaTsPy8vz6FvFRE1TvTLJBGRt+N9lmQj8px36wAqKSkJGzduxF/+8hf4+fkpfZ7at2+PDh06QKPRYOrUqcjIyEBoaCh69+6N5cuXo3379jAYDACujpb3yCOPYPbs2cjMzITNZsPs2bMxYsQIhIaGAgAMBgNMJhOmTZuGpKQkHDp0CCtXrnSYB2rKlCkYNWoUVqxYgSeeeAKffvopCgoKlLmlmoOlRETqupVrjtcbERER1efWAdTatWsBQBmi3M5oNCIlJQUAMHPmTFy4cAHJycnKRLqbN29W5oACgDVr1sBoNGL06NEAgJEjR2Lp0qXK+jvvvBMfffQRkpKSMGzYMPj5+WH69OmYMWOGkkan0yEnJwdLlixBWloaevbsiZycnBbNAcVSIiJ1sTkTkTzY/4tIPbJWCrh1AGW1Wm+aRqPRICUlRQmoGuPv73/T4cbvuecefP755zdMExMT0yCYIyKixsn6YJWd6N+dBSYkE0++3gDPvebcOoAi8iaib3JsxkZqk/XBKjv+7kTq4fUmBgMoIpWIvsmxVJaIiFxFdCEhkZoYQBERERHRLRFdSEikpl+JzgAREREREZGnYABFRERERETURGzCR0TkxTh4CBERkXMxgCIi8mIcPISIiMi52ISPiIiIiIioiVgDRUREROQEHMqbSA4MoIiIiIicgEN5E8mBTfiIiIiIiIiaiAEUERERERFREzGAIiIiIiIiaiIGUERERERERE3EAIqIiIiIiKiJOAqfhObcuQXndn5yC59/CsA052WIiMjL8D5LROS9GEBJaMXpGCwc3fIH84r1J7DQifkh8mZ8kZYT77NERN6LARSpii+T4tzKsedxbzmZX6RFX+8yn/Myf3cimYi+z8qKARSpSuaXSdFu5djzuHsukS/Soq93mc95mb87kUxE32dlxQBKAJYWEJFa+CJNJAe+W4jDYy+GyOPOAMqFxowZAwDYuHGjw/IVy3dgX9c7lb+vXW//3PXWrxgUg30bv27x57FmGsZsv03Y50Wud0ZpvDt/P653/no8tkLezx+rAeK3yPn5YzUYs/22lp8/nvz5//fZG27/Bp+fE1qNczvr0FJz/uWDMWM6tHj/t5J/0c/HW/n8itOO7wbN/vz6E9gn6Pl69fNPYcyYluXfWS/SLX1+2AuqWvz5/1dY5W7PP3dfv+OxFQ4FhM3/fIxDIWH99Q2uzWswgCKp3MoD5lZL4+03+JzJjsvP7XzM4e/rrWcJFRF5ghU7W2HfwZYHQDuO1WBEV9fmkZzPGQHcCDh+vjn7dkYAQ81XP3htyfuNp77bMIAiUonMN3jRN9hb2X+Gdcst7ZuI1DUn+jISo08DaHoB1X8/64N9B12aPSKvUr+ZeEvebzz13YYBFBG5nOgb7K3s/9oqfiJyb7dSA7bilRjWfhHRTWmsVqtNdCZk47f+BKzxd0n5eebdMz/vyXm/1c97ct5v9fOenPdb/bwn5/1WP+/Jeb/VzzPvnvl5T877rX7ek/N+q58XuW/WQBERERE5AUdjI5IDAygiIhfzW3+i5Z9to3FiTtQn83cXraXHnse95UTPycPrTRwee7kwgGqmtWvX4o9//CMsFgv69OmDtLQ0DBo0SHS2mk3khS76JsOXCjFE/+6i3Kx5wK02QXBnMn93QOw5f6Pj6u3HXVayX28i8djLhwFUM2zevBnz5s1DRkYGoqKisHbtWsTFxWH37t3o2tVzep2KvNBF32RE7190ECEqeBR93O37aNHnPDh4s+N3b8HnvOCcl5nM57xIMhfOiiT6u8v4uzOAaobVq1fj+eefxwsvvAAAWLZsGf7+978jJycHr776arO2JfpkJ/WJfqESvX+R+N2vT9bv7s3f2x3IXPsm6/Nd5sJZkUR/d1l/dwZQTXTp0iXs3bsXCQkJDsujo6NRWFjYrG2JPtmJiIhcReZnnMzfnUgmvxKdAU9RUVGBuro6BAYGOiwPDAxEWVnZLW07LS0Nfn5+yj/Mvsfh77S0tFvavrvvXySR3130ced3b3z/3nzcb7Z/mb4773NyfHfReJ/ld+c15504D1QTlZaWIjw8HNu2bXMYNCI9PR0ffvgh9uzZIzB3nistLQ0mk+m6641GI1JSUlTMERGR88l8rxP53WU+7iQn0ee8LNc7A6gmunTpEjp37ox169bht7/9rbI8KSkJ+/fvx7Zt2wTmjoiIiIiI1MAmfE3Upk0b9OvXD3l5eQ7L8/LyoNPpBOWKiIiIiIjUxEEkmmH69OmYPHky+vfvD51Oh5ycHJw6dQrx8fGis0ZERERERCpgANUMo0ePRmVlJZYtWwaLxfL/t3fnQVGc+RvAHwR1FY9BNAPqDAoMp+ABCiHxihFQVlBgoylcXSNgUJboSpQohkDIggRJNB4oxFpU1BghgonCbimB4QhmN0YiuyLE4IEHioIiHhzz+yPr/ERUOjrdusnzqaJKXgafbqrm2++3++0e2NraYs+ePVAqlc9604iIiIiISAK8B4qIiIiIiEgg3gNFREREREQkEBsoIiIiIiIigdhAERERERERCcQGioiIiIiISCA2UERERERERAKxgSIiIiIiIhKIDRQREREREZFAbKCIiIiIiIgEYgNFREREREQkEBsoIiIiIiIigdhAERERERERCcQGioiIiIiISCA2UERERERERAKxgSIiIiIiIhKIDRQREREREZFAbKCIiIiIiIgEYgNFREREREQkEBsoIiIiIiIigdhAERERERERCcQG6jmTmpoKR0dHyOVyjB8/HsXFxZLkFhUVYdasWbC1tYVMJkN6erokuQCQlJSEiRMnQqFQwMLCAjNnzsS///1vSbJTUlLg5uYGhUIBhUKByZMnIzc3V5LsB61ZswYymQxvv/22JHlxcXGQyWTtvqysrCTJvufixYt48803YWFhAblcDhcXFxQWFoqe6+Dg0GHfZTIZXnvtNdGzAaC1tRWxsbHa97qjoyNiY2PR0tIiSf6NGzcQERGBYcOGwcTEBO7u7vjuu+9Eyeqstmg0GsTFxcHGxgYmJibw8vLCf/7zH0mys7Oz4evrCwsLC8hkMqjVap3kCslvbm5GVFQU3NzcMHDgQFhbWyMwMBBnz54VPRsAYmNjMXr0aAwcOBBmZmbw9vZGaWmpTrKF5N/vrbfegkwmwyeffCJJdkhISIf3/quvvqqTbCH5AFBVVYXZs2dDqVTC1NQU48aNQ0VFhejZD6t7MpkM4eHhT50tJL+xsRFvv/027OzsYGJiAmdnZ2zYsEGS7NraWoSEhMDGxgampqbw8/PDjz/+qJNsIfMYsWqdkGwxa11n+WLWOiH7LlatYwP1HMnMzERERASWLl2KgoICjBkzBn/4wx90dkB9nJs3b8LOzg7x8fHo0aOH6Hn3KywsxPz585Gbm4vs7GwYGBhg+vTpuHbtmujZAwcORHR0NPLz85GXl4dx48YhICAAx48fFz37ft9++y3S0tJgb28vaa5KpUJFRYX2S6qGHQDq6+vh4eEBjUaDPXv2oLS0FAkJCRgwYIDo2Xl5ee32Oz8/H3p6epg+fbro2QDw8ccfIzU1FatXr8aRI0cQHx+PlJQUJCUlSZIfFhaGw4cPY9OmTSguLsbEiRMxffp0nD9/XudZndWWtWvXYsOGDVi9ejUOHz6MAQMGYMaMGbhx44bo2U1NTRgzZgw++OCDp876pflNTU04duwYwsPDkZ+fj507d6Kmpgb+/v46aaQ723eVSoXExEQUFxcjJycHZmZm8Pf3R21t7VNnC8m/JysrC9999x1MTU11kis0e8KECe1qwOeffy5ZfnV1NTw8PGBmZobs7GyUlJQgMjIShoaGomffv88VFRXYvXs3AOis9nWWv3LlSvz9739HcnIySktLsXTpUkRHR2u3Q6xsjUaDgIAAnDp1Cunp6SgoKIBCoYCPjw9u3rz51NlC5jFi1Toh2WLWus7yxax1QvZdrFqnV19fr3mq/4F0ZtKkSbC3t8e6deu0Y6NGjYKPjw+ioqIk245BgwYhISEBAQEBkmXer7GxEUqlEunp6ZgyZYrk+UOGDEFUVBTmzZsnSV5DQwPGjx+PtWvXIiEhAXZ2dvjwww9Fz42Li9MevJ+FmJgYFBUVPbMrfvdLTEzEunXrcOLECfTs2VP0vJkzZ8LIyAjJycnasTfffBPXrl3DZ599Jmr2rVu3MHjwYGzbtg1eXl7a8fHjx2Py5MmIjIwULfvB2qLRaGBjY4OgoCDtGfBbt25BpVLh/fff1+l78HF1ra6uDhYWFti/fz/Gjh2rs0yh+fecOHECrq6uKCoq0unJFCHZ169fh1KpREZGBiZNmqSz7MflnzlzBh4eHti3bx/8/f0RHByMP//5z6Jnh4SE4OrVq6K/1x6VHxgYCD09PaSkpEie/aCwsDAUFxfjn//8pyT5L774IqZNm4YVK1Zox6ZOnQp7e3udHvcezK6qqoKzszPUajUcHBwAAG1tbbCyssK7776LOXPm6Cwb6DiPkbLWPW4OJUWtEzKHE6vWCcnWVa3jFajnxN27d/H999/jlVdeaTf+yiuv6HRZxf+CxsZGtLW1QSaTSZrb2tqKjIwM3Lx5E2PGjJEsd/HixfDx8cH48eMly7ynuroatra2cHR0xBtvvIHq6mrJsr/66is4OTlh3rx5sLS0xMsvv4wtW7ZAo5H2nI5Go8H27dsxc+ZMSZonAHB1dUVhYSFOnjwJ4OeDiVqtxuTJk0XPbmlpQWtrK373u9+1G+/Ro4fkzfTp06dx6dKldnWvR48ecHNz+83VPQDaM9FS1767d+8iLS0Nffr00U4uxdbS0oLAwECEh4fD2tpaksz7lZSUwNLSEk5OTggLC8Ply5clyW1ra0NOTg6sra3h5+cHCwsLTJw4EZmZmZLk3+/GjRvIzMzE3LlzJct0dXVFTk4Ozp07BwAoLS3F8ePHdd60P+jOnTsA0K7udenSBd27dxel7j04j5Gy1j2rOdQvyRer1nWWrctaZ/BUv006U1dXh9bW1g7LlwYMGKCzJRX/KyIiIuDg4CBZE1NeXg53d3fcvn0bhoaG2LFjh2RL6dLS0nDq1Cls3rxZkrz7OTs7Y+PGjVCpVLhy5Qo+/PBDuLu745tvvkG/fv1Ez6+ursann36KhQsXYvHixfjhhx+wfPlyAEBwcLDo+ffk5eXh9OnT+OMf/yhZ5uLFi9HY2AgXFxfo6+ujpaUF4eHhCAwMFD27d+/eGDNmDBITE2Frawu5XI69e/fiyJEjMDc3Fz3/fpcuXQKAh9a9CxcuSLotz9rdu3cRGRkJT09PDBo0SJLMnJwczJ8/H01NTTAxMcEXX3yBF154QZLsuLg4GBkZYf78+ZLk3e/VV1/FtGnTYGZmhjNnziA2Nhbe3t74+uuv0b17d1GzL1++jMbGRiQlJWHFihWIiopCQUEBgoKC0LNnT3h6eoqaf7+MjAzcuXMHr7/+umSZq1evxpIlSzBs2DAYGPw8BU1ISBB9v62srKBQKBATE4N169bB0NAQGzduRE1NjbYO6dKD8xgpa53Uc6hfmi9mrXtUthi1jg3Uc0ZPT6/d9xqNpsPYr9mKFSvwzTffICcnB/r6+pJkqlQqqNVqNDQ0IDs7GyEhIfjyyy9hZ2cnam5lZSViYmJw8OBBdOvWTdSsh3nwaoezszNGjBiBnTt3IjQ0VPT8trY2jBw5Urs8dfjw4Th16hRSU1MlbaDS0tIwatQoODo6SpaZmZmJ3bt3IzU1FTY2Nvjhhx8QEREBpVKp86UkD7N582YsWrQIdnZ20NfXx/Dhw+Hv749jx46Jnv0wv/W619LSguDgYDQ0NGDXrl2S5Y4dOxZqtRp1dXVIS0vDn/70J/zjH/+AiYmJqLmFhYXYuXOnzh/aIZSfn5/23/b29hgxYgQcHByQm5sLb29vUbPb2toA/Lxs7V6ddXR0xPfff4/U1FRJG6i0tDR4eXmhf//+kmVu3rwZpaWl2LVrFxQKBYqLi7Fq1SoolUqdPsjjQV27dsX27dsRGhqKoUOHQl9fHxMmTBDlqv/j5jFi17pnMYf6Jfli1rrHZYtR67iE7zlhbGwMfX39Dlebrly5IslN9c+Dd955BxkZGcjOzsaQIUMky+3WrRvMzc21k3kHBwds3LhR9NwjR46grq4OL774IoyNjWFsbIyioiKkpqbC2NhYu+RAKr169YKNjQ1OnTolSZ5cLu+wdMfKykq7tEMKly9fxoEDByRdwgIA7777LkJDQ+Hn5wd7e3vMmjULixYtwkcffSRJ/tChQ3HgwAHU1NSgvLwchw8fRnNzM8zMzCTJv0culwPAb7rutbS0YP78+SgvL0dWVpYkV3/vMTQ0hLm5OUaPHo3169eja9eu2LZtm+i5arUaFy9ehLW1tbb2nT17FlFRUaKfuHoYU1NTDBw4UJLaZ2xsDAMDg2de+8rKynD06FFJa9+tW7cQExOD6OhoTJkyBcOGDUNwcDB8fX119gTGxxkxYgQKCwtx+vRpVFRUICMjA1evXtVp3XvUPEaKWves5lBC88WsdZ1li1Hr2EA9J7p164YRI0YgLy+v3XheXh5cXFye0VZJZ/ny5di7dy+ys7Mlf5T2g9ra2nD37l3Rc7y8vFBcXAy1Wq39GjlyJPz8/KBWqyW/KnX79m1UVlZqC73YXF1dUVVV1W6sqqoKCoVCknwASE9PR/fu3eHr6ytZJvDzU4kePEOmr6+vPTstFUNDQ5iYmKC+vh6HDh3C1KlTJc03MzODXC5vV/du376NkpKS30Tda25uxrx581BeXo79+/dL9t57FKlqX2BgIIqKitrVPlNTUyxcuBBZWVmi5z+orq4OFy5ckOTv361bN4waNQqVlZXtxqWufWlpaVAqlZgwYYJkmc3NzWhubn7mta9v377o378/fvzxRxw9elRnde9x8xixa92znkN1li9mrXuSfddFreMSvufIokWLsGDBAjg5OcHFxQVbt27FxYsXJXkaXGNjo/bsW1tbG86dO4eysjIYGRmJXtTDw8Px2WefYceOHZDJZNq1woaGhujVq5eo2e+99x7c3d0xaNAgNDY2Yu/evSgsLMSePXtEzQX+//M47tezZ08YGRlJchb23hrkwYMHa++Bampqkmw9/MKFC+Hu7o7ExET4+vqirKwMW7ZswapVqyTJ12g02LZtG3x9fdG7d29JMu/x9PTExx9/DDMzM9jY2KCsrAwbNmzArFmzJMk/dOgQ2traoFKp8NNPP2HVqlVQqVSiPHmzs9oSEhKCNWvWQKVSwdLSEomJiTA0NIS/v7/o2deuXcPZs2fR0NAAAPjpp5/Qt29fyOVynRzgH5dvamqKuXPn4ujRo9i1axf09PS0ta9Pnz5P/XESj8vu27cv1q1bB09PT8jlctTV1SElJQXnz5/X2eOsO/vbP3jW3cDAAHK5HCqVStRsIyMjxMfHw9vbG3K5HGfOnEFMTAwGDBiA3//+90+d3Vm+QqFAWFgY5s2bBzc3N4wbNw5qtRqZmZk6+fxFIcfypqYmfP755wgLC9P5UtnO8l966SVER0fD0NAQCoUCRUVF2L17N6Kjo0XP3rdvH/r16welUony8nJERETAy8urw8O7nkRn8xg9PT3Rap2QOZSYta6z/JaWFtFqXWfZ169fF63W8THmz5nU1FSsXbsWly5dgq2tLf7617/ipZdeEj1XrVZj2rRpHcZff/11bNq0SdTsRz0tZfny5XjnnXdEzQ4JCYFarUZtbS369OkDe3t7hIWFif5EoEfx8vKS7DHmb7zxBoqLi1FXV4f+/fvD2dkZK1euhI2NjejZ9+Tm5iImJgZVVVUYPHgwgoKCsGDBAknufykoKIC3tzcOHToEJycn0fPud+PGDXzwwQf48ssvceXKFcjlcvj5+WHZsmUdno4nhi+++ALR0dE4f/48jIyM4O3tjcjISPTt21fnWZ3VFo1Gg/j4ePztb39DfX09nJyckJiYqJOTCJ1lp6enY9GiRR1+rqva87j8iIgIDB8+/KG/t2HDhqduZh+XvWbNGgQFBeFf//oXrl69in79+mHkyJFYunQpnJ2dnypXSP7DjikODg46e4z547KTkpIQEBCAsrIyNDQ0QC6XY+zYsVi5ciUGDx781Nmd5d/b9/T0dCQlJaGmpgbm5ub4y1/+opOTBkKyd+zYgbfeegvHjx/X6edvCcm/dOkSoqOjkZeXh2vXrkGhUGDOnDkIDQ196rrfWXZycjI++eQT1NbWQi6XY9asWVi2bJlOVnsImceIVeuEZItZ6zrLP336tGi1rrPspqYm0WodGygiIiIiIiKBeA8UERERERGRQGygiIiIiIiIBGIDRUREREREJBAbKCIiIiIiIoHYQBEREREREQnEBoqIiIiIiEggNlBEREREREQCsYEiIqJfpfT0dMhkMnz77bcP/fmSJUse+UGMREREj8IGioiIiIiISCA2UERERBLRaDS4ffv2s94MIiJ6CmygiIiI/mvbtm1wc3ODXC6HpaUlFixYgAsXLrR7jZeXF7y8vDr8blxcXIclgTKZDEuWLMG+ffvg5uaGF154ARkZGQCA/Px8TJkyBWZmZhg0aBCcnZ2xdOlS8XaOiIh0wuBZbwAREZGYrl+/jrq6ug7jD14J+uijjxAdHQ03NzfExMTg3LlzSElJQUlJCQoKCp74fqmSkhJkZWUhKCgIcrkcVlZWOHHiBF577TXY2dkhIiICPXv2RHV1NXJzc58og4iIpMMGioiIftX8/Pw6fU1dXR3i4+Px8ssvY9++fTAw+Pnw6OrqioCAAKxfvx6RkZFPlF9RUYH8/Hw4OjpqxzZt2oQ7d+5g7969MDY21o5HRUU9UQYREUmHDRQREf2qrV69GtbW1h3GN23apL3i8/XXX+POnTtYuHChtnkCfl6up1KpkJub+8QNlIuLS7vmCQB69+4NAPjqq68we/ZsdOnCFfVERP8r2EAREdGv2qhRozB69OgO41lZWdp/nzlzBgBgZWXV4XVWVlYoLCx84vwhQ4Z0GPPz88P27dsRFhaG9957D+PGjcPUqVMxY8YMdO3a9YmziIhIfDzlRURE9Bgajabd93p6eg99XWtr60PHe/To8dCxgwcPIjs7G7Nnz0ZlZSWCg4MxadIk3Lp16+k3moiIRMMGioiIfvOUSiUA4OTJkx1+VllZqf058POT9RoaGjq87t5VLKG6dOmCcePGISYmBkVFRVizZg3Kysqwf//+X7j1REQkJTZQRET0mzdhwgR0794dycnJ7a4kHTx4EJWVlfDw8NCOmZub4+TJk6itrdWOnT9/HgcOHBCcd/Xq1Q5jw4cPBwDU19c/yS4QEZFEeA8UERH95hkbGyMiIgLR0dHw8fHBtGnTUFNTgy1btkCpVCI0NFT72jlz5mD9+vWYMWMG5s6di4aGBmzduhUWFhY4duyYoLyEhAQUFhbCw8MDSqUS9fX12Lp1KwwNDeHp6SnWbhIRkQ6wgSIiIgKwZMkSGBsbIzk5GatWrUKvXr3g4+ODqKiodp8BZWlpiU8//RSxsbFYuXIlhg4divfffx+VlZWCG6ipU6fi3Llz2LVrF65cuYJ+/fph9OjRWLZsWbvlgkRE9PzRq6+v13T+MiIiIiIiIuI9UERERERERAKxgSIiIiIiIhKIDRQREREREZFAbKCIiIiIiIgEYgNFREREREQkEBsoIiIiIiIigdhAERERERERCcQGioiIiIiISCA2UERERERERAKxgSIiIiIiIhLo/wAe6nKprG8BLgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "demand_profile.boxplot(column=list(demand_profile.columns), grid=False, figsize=(12,8));\n",
    "plt.yticks(np.arange(0,250000,20000));\n",
    "plt.axhline(y=5000, ls=(0,(1,1)), alpha=0.7, color='k', label='Low \\u03C8')\n",
    "plt.axhline(y=70000, ls='--', alpha=0.7, color='m', label='Medium \\u03C8');\n",
    "plt.axhline(y=180000, ls='dotted', alpha=0.7, color='c', label='Medium \\u03C8');\n",
    "plt.xlabel('Hours');\n",
    "plt.ylabel('kWh/hr')\n",
    "plt.legend();\n",
    "plt.title('Hourly Demand Distribution');\n",
    "#plt.savefig('demand_dist.jpg', quality=95)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file C:\\Users\\garg1\\gurobi.lic\n",
      "Academic license - for non-commercial use only\n",
      "Changed value of parameter MIPGap to 0.01\n",
      "   Prev: 0.0001  Min: 0.0  Max: inf  Default: 0.0001\n",
      "Parameter OutputFlag unchanged\n",
      "   Value: 1  Min: 0  Max: 1  Default: 1\n"
     ]
    }
   ],
   "source": [
    "gurobi_env=Env(empty=False);\n",
    "gurobi_env.setParam('MIPGap',0.01);\n",
    "gurobi_env.setParam('OutputFlag',1);\n",
    "# gurobi_env.setParam('MIPFocus',3);\n",
    "# gurobi_env.setParam('ImproveStartTime',300);\n",
    "#gurobi_env.setParam('Cuts',3);\n",
    "#gurobi_env.setParam('CliqueCuts',2);\n",
    "#gurobi_env.setParam('FeasibilityTol',1e-3);\n",
    "#gurobi_env.setParam('Presolve',0);\n",
    "#gurobi_env.setParam('NodefileStart',0.5);\n",
    "#gurobi_env.setParam('Threads',8);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "times=range(nr_t)\n",
    "fuels=range(nr_fuel)\n",
    "customers=range(nr_cust)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* __With Emission Model Function__\n",
    "\n",
    "    * Create a user defined function to build & solve optimization model which considers the emission in the objection function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def with_emission(nr_t,nr_fuel,nr_cust,alpha_rate,fuel_supply,epsilon_rate,BIGM_price,caligraphy_n,beta_rate,\n",
    "                  demand_profile,demand,psi_rate,omega_rate,cost_xi,cost_gamma,BIGM_kkt1,BIGM_kkt2,BIGM_kkt3,\n",
    "                  cost_fuel,cost_ramp,cost_emission,data=data2):\n",
    "    \n",
    "    \"\"\"Function to Write Model with Carbon Emissions\"\"\"\n",
    "    \n",
    "    MCE=Model(env=gurobi_env, name='MCE')\n",
    "\n",
    "    #Defining Leaders' varfiable\n",
    "    PIl = MCE.addVars(times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"PIl\")\n",
    "    PIh = MCE.addVars(times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"PIh\")\n",
    "    X = MCE.addVars(times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"X\")\n",
    "    Y = MCE.addVars(fuels,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"Y\")\n",
    "    P = MCE.addVars(times,vtype=GRB.BINARY,name=\"P\")\n",
    "\n",
    "    #defining followers variable\n",
    "    Wl = MCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"Wl\")\n",
    "    Wh = MCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"Wh\")\n",
    "    Vplus = MCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"Vplus\")\n",
    "    Vminus = MCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"Vminus\")\n",
    "    U = MCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"U\")\n",
    "\n",
    "    #Defining KKT condition variables\n",
    "    MUa = MCE.addVars(customers,times,lb=-GRB.INFINITY,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"MUa\")\n",
    "    MUb = MCE.addVars(customers,lb=-GRB.INFINITY,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"MUb\")\n",
    "    MUc = MCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"MUc\")\n",
    "    MUd = MCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"MUd\")\n",
    "\n",
    "    #McCormic\n",
    "    rhoa = MCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhoa\")\n",
    "    rhob = MCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhob\")\n",
    "    rhoc = MCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhoc\")\n",
    "    rhod = MCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhod\")\n",
    "    rhoe = MCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhoe\")\n",
    "    rhof = MCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhof\")\n",
    "    rhog = MCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhog\")\n",
    "\n",
    "    MCE.update()\n",
    "    #Adding Constraints:\n",
    "\n",
    "    #leader's\n",
    "    MCE.addConstrs(\n",
    "        (Wl.sum('*',t) + Wh.sum('*',t) - X[t] == 0 for t in range(nr_t)), name=\"9a\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (quicksum(alpha_rate[s]*Y[s,t] for s in range(nr_fuel)) - X[t] == 0 for t in range(nr_t)), name=\"9b\")\n",
    "\n",
    "    #MCE.addConstrs((Y.sum(s,'*') <= fuel_supply[s] for s in range(nr_fuel)), name=\"9c\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (PIl[t] - PIl[t+1] <= BIGM_price*P[t] for t in range(nr_t-1)), name=\"9d\")\n",
    "    MCE.addConstrs(\n",
    "        (PIl[t] - PIl[t+1] >= -BIGM_price*P[t] for t in range(nr_t-1)), name=\"9d1\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (PIh[t] - PIh[t+1] <= BIGM_price*P[t] for t in range(nr_t-1)), name=\"9e\")\n",
    "    MCE.addConstrs(\n",
    "        (PIh[t] - PIh[t+1] >= -BIGM_price*P[t] for t in range(nr_t-1)), name=\"9e1\")\n",
    "\n",
    "    MCE.addConstr(\n",
    "        (quicksum(P[t] for t in range(nr_t)) <= caligraphy_n), name=\"9f\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (PIl[t] <= PIh[t] for t in range(nr_t)), name=\"9g\")\n",
    "\n",
    "    #Followers\n",
    "    MCE.addConstrs(\n",
    "        (Wl[k,t]+Wh[k,t]+Vminus[k,t]+U[k,t]-Vplus[k,t] == demand_profile.iloc[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"10a\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (quicksum(Wl[k,t] + Wh[k,t] + U[k,t] for t in range(nr_t)) == demand[k] for k in range(nr_cust)), name=\"10b\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (Wl[k,t] <= psi_rate for k in range(nr_cust) for t in range(nr_t)), name=\"10c\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (Vplus[k,t] <= omega_rate.iloc[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"10d\")\n",
    "    \n",
    "    MCE.update()\n",
    "\n",
    "    #KKT Conditions\n",
    "    \n",
    "    #Stationary Constraints\n",
    "    MCE.addConstrs(\n",
    "        (MUa[k,t] + MUb[k] - MUc[k,t] <= PIl[t] for k in range(nr_cust) for t in range(nr_t)), name=\"11a\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (MUa[k,t] + MUb[k] <= PIh[t] for k in range(nr_cust) for t in range(nr_t)), name=\"11b\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (-MUa[k,t] - MUd[k,t] <= cost_xi.iloc[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11c\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (MUa[k,t] <= 0 for k in range(nr_cust) for t in range(nr_t)), name=\"11d\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (MUa[k,t] + MUb[k] <= cost_gamma[t] for k in range(nr_cust) for t in range(nr_t)), name=\"11e\")\n",
    "    \n",
    "    #Complementary Slackness Condition Constraints\n",
    "    MCE.addConstrs(\n",
    "        (psi_rate - Wl[k,t] <= psi_rate*(1-rhoa[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11g\")\n",
    "    MCE.addConstrs(\n",
    "        (MUc[k,t] <= BIGM_kkt1*rhoa[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11g1\")\n",
    "\n",
    "    MCE.addConstrs(\n",
    "        (omega_rate.iloc[k,t]-Vplus[k,t] <= omega_rate.iloc[k,t]*(1-rhob[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11h\")\n",
    "    MCE.addConstrs(\n",
    "        (MUd[k,t] <= BIGM_kkt1*rhob[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11h1\")\n",
    "    \n",
    "    MCE.addConstrs(\n",
    "        (PIl[t] - MUa[k,t] - MUb[k] + MUc[k,t] <= BIGM_kkt2*(1-rhoc[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11i\")\n",
    "    MCE.addConstrs(\n",
    "        (Wl[k,t] <= psi_rate*rhoc[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11i1\")\n",
    "    \n",
    "    MCE.addConstrs(\n",
    "        (PIh[t] - MUa[k,t] - MUb[k]<= BIGM_kkt2*(1-rhod[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11j\")\n",
    "    MCE.addConstrs(\n",
    "        (Wh[k,t] <= (demand_profile.iloc[k,t]+omega_rate.iloc[k,t]-psi_rate)*rhod[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11j1\")\n",
    "    \n",
    "    MCE.addConstrs(\n",
    "        (cost_xi.iloc[k,t] + MUa[k,t] + MUd[k,t]  <= BIGM_kkt1*(1-rhoe[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11k\")\n",
    "    MCE.addConstrs(\n",
    "        (Vplus[k,t] <= omega_rate.iloc[k,t]*rhoe[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11k1\")\n",
    "    \n",
    "    MCE.addConstrs(\n",
    "        (0 - MUa[k,t] <= BIGM_kkt1*(1-rhof[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11l\")\n",
    "    MCE.addConstrs(\n",
    "        (Vminus[k,t] <= demand_profile.iloc[k,t]*rhof[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11l1\")\n",
    "    \n",
    "    MCE.addConstrs(\n",
    "        (cost_gamma[t] - MUa[k,t] - MUb[k] <= BIGM_kkt1*(1-rhog[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11m\")\n",
    "    MCE.addConstrs(\n",
    "        (U[k,t] <= (demand_profile.iloc[k,t]+omega_rate.iloc[k,t])*rhog[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11m1\")\n",
    "\n",
    "    #setting up the objective\n",
    "\n",
    "    MCE.setObjective(\n",
    "        quicksum(demand_profile.iloc[k,t]*MUa[k,t] - psi_rate*MUc[k,t] - omega_rate.iloc[k,t]*MUd[k,t] - cost_gamma[t]*U[k,t]\n",
    "                     -cost_xi.iloc[k,t]*Vplus[k,t] for k in range(nr_cust) for t in range(nr_t)) +\n",
    "        quicksum(demand[k]*MUb[k] for k in range(nr_cust)) -\n",
    "        quicksum(cost_fuel[s]*Y[s,t] for s in range(nr_fuel) for t in range(nr_t))-\n",
    "        quicksum(cost_emission*beta_rate[s]*Y[s,t] for s in range(nr_fuel) for t in range(nr_t)), GRB.MAXIMIZE)\n",
    "\n",
    "    MCE.update()\n",
    "    \n",
    "    MCE.reset(0)\n",
    "    print('--------------------------\\nStarted Solving MCE Model')\n",
    "    MCE.optimize()\n",
    "    \n",
    "    MCE_obj = MCE.objVal\n",
    "    print('Optimality Gap: {:.2f}'.format(MCE.MIPGap))\n",
    "    \n",
    "    #--------Capturing Details of Leaders----------------\n",
    "    print('\\n\\tGetting Results\\n--------------------------')\n",
    "    \n",
    "    #First create a list of all variables which are required to be output\n",
    "    pil = [PIl[t] for t in range(nr_t)];\n",
    "    pih = [PIh[t] for t in range(nr_t)];\n",
    "    x = [X[t] for t in range(nr_t)];\n",
    "    f1 = [Y[0,t] for t in range(nr_t)]; #Coal usage\n",
    "    f2 = [Y[1,t] for t in range(nr_t)]; #Gas usage\n",
    "    f3 = [Y[2,t] for t in range(nr_t)]; #Oil usage\n",
    "    w_l = [];\n",
    "    w_h = [];\n",
    "    u = []\n",
    "    vplus = []\n",
    "    vminus = []\n",
    "    for t in range(nr_t):\n",
    "        w_l.append(sum(MCE.getAttr(GRB.Attr.X, [Wl[k,t] for k in range(nr_cust)])))\n",
    "        w_h.append(sum(MCE.getAttr(GRB.Attr.X, [Wh[k,t] for k in range(nr_cust)])))\n",
    "        u.append(sum(MCE.getAttr(GRB.Attr.X, [U[k,t] for k in range(nr_cust)])))\n",
    "        vplus.append(sum(MCE.getAttr(GRB.Attr.X, [Vplus[k,t] for k in range(nr_cust)])))\n",
    "        vminus.append(sum(MCE.getAttr(GRB.Attr.X, [Vminus[k,t] for k in range(nr_cust)])))\n",
    "    \n",
    "    df_leader = pd.DataFrame(index=range(nr_t));\n",
    "    df_leader.index.name = 'Day-hours'\n",
    "    \n",
    "    df_leader['Low_Price'] = MCE.getAttr(GRB.Attr.X, pil);\n",
    "    df_leader['High_Price'] = MCE.getAttr(GRB.Attr.X, pih);\n",
    "    df_leader['Elec_Generated'] = MCE.getAttr(GRB.Attr.X, x);\n",
    "    df_leader['Low_Price Demand'] = w_l\n",
    "    df_leader['High_Price_Demand'] = w_h\n",
    "    df_leader['Upward_Demand_Shifted'] = vplus\n",
    "    df_leader['Downward_Demand_Shifted'] = vminus\n",
    "    df_leader['Comp_Demand'] = u\n",
    "    df_leader['Coal_Used'] = MCE.getAttr(GRB.Attr.X, f1);\n",
    "    df_leader['Coal_Emissions (lbs)'] = df_leader['Coal_Used']*data.iloc[0,4];\n",
    "    df_leader['Gas_Used'] = MCE.getAttr(GRB.Attr.X, f2);\n",
    "    df_leader['Gas_Emissions (lbs)'] = df_leader['Gas_Used']*data.iloc[1,4];\n",
    "    df_leader['Oil_Used'] = MCE.getAttr(GRB.Attr.X, f3);\n",
    "    df_leader['Oil_Emissions (lbs)'] = df_leader['Oil_Used']*data.iloc[2,4];\n",
    "    \n",
    "    #--------Capturing Details of Followers----------------\n",
    "    cust = ['Cust'+str(j) for j in range(nr_cust)]\n",
    "    df_wh = pd.DataFrame(index=range(nr_t));\n",
    "    df_wh.index.name = 'Day-hours'\n",
    "    \n",
    "    df_wl = pd.DataFrame(index=range(nr_t));\n",
    "    df_wl.index.name = 'Day-hours'\n",
    "    \n",
    "    df_vplus = pd.DataFrame(index=range(nr_t));\n",
    "    df_vplus.index.name = 'Day-hours'\n",
    "    \n",
    "    df_vminus = pd.DataFrame(index=range(nr_t));\n",
    "    df_vminus.index.name = 'Day-hours'\n",
    "    \n",
    "    df_ukt = pd.DataFrame(index=range(nr_t));\n",
    "    df_ukt.index.name = 'Day-hours'\n",
    "    \n",
    "    for k in range(nr_cust):\n",
    "        df_wh[cust[k]] = MCE.getAttr(GRB.Attr.X, [Wh[k,t] for t in range(nr_t)])\n",
    "        df_wl[cust[k]] = MCE.getAttr(GRB.Attr.X, [Wl[k,t] for t in range(nr_t)])\n",
    "        df_vplus[cust[k]] = MCE.getAttr(GRB.Attr.X, [Vplus[k,t] for t in range(nr_t)])\n",
    "        df_vminus[cust[k]] = MCE.getAttr(GRB.Attr.X, [Vminus[k,t] for t in range(nr_t)])\n",
    "        df_ukt[cust[k]] = MCE.getAttr(GRB.Attr.X, [U[k,t] for t in range(nr_t)])\n",
    "\n",
    "    \n",
    "    return MCE_obj, df_leader.round(0), df_wh.round(0), df_wl.round(0), df_vplus.round(0), df_vminus.round(0), df_ukt.round(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* __Without Emission Model Function__\n",
    "\n",
    "    * Create a user defined function to build & solve optimization model which does not consider the emission in the objection function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def without_emission(nr_t,nr_fuel,nr_cust,alpha_rate,fuel_supply,epsilon_rate,BIGM_price,caligraphy_n,beta_rate,\n",
    "                  demand_profile,demand,psi_rate,omega_rate,cost_xi,cost_gamma,BIGM_kkt1,BIGM_kkt2,BIGM_kkt3,\n",
    "                  cost_fuel,cost_ramp,cost_emission,data=data2):\n",
    "    \n",
    "    \"\"\"Without Emission Model\"\"\"\n",
    "    \n",
    "    MWCE=Model(env=gurobi_env, name='MWCE')\n",
    "\n",
    "    #Defining Leaders' varfiable\n",
    "    PIl = MWCE.addVars(times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"PIl\")\n",
    "    PIh = MWCE.addVars(times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"PIh\")\n",
    "    X = MWCE.addVars(times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"X\")\n",
    "    Y = MWCE.addVars(fuels,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"Y\")\n",
    "    P = MWCE.addVars(times,vtype=GRB.BINARY,name=\"P\")\n",
    "\n",
    "    #defining followers variable\n",
    "    Wl = MWCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"Wl\")\n",
    "    Wh = MWCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"Wh\")\n",
    "    Vplus = MWCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"Vplus\")\n",
    "    Vminus = MWCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"Vminus\")\n",
    "    U = MWCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"U\")\n",
    "\n",
    "    #Defining KKT condition variables\n",
    "    MUa = MWCE.addVars(customers,times,lb=-GRB.INFINITY,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"MUa\")\n",
    "    MUb = MWCE.addVars(customers,lb=-GRB.INFINITY,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"MUb\")\n",
    "    MUc = MWCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"MUc\")\n",
    "    MUd = MWCE.addVars(customers,times,lb=0.0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name=\"MUd\")\n",
    "\n",
    "    #McCormic\n",
    "    rhoa = MWCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhoa\")\n",
    "    rhob = MWCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhob\")\n",
    "    rhoc = MWCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhoc\")\n",
    "    rhod = MWCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhod\")\n",
    "    rhoe = MWCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhoe\")\n",
    "    rhof = MWCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhof\")\n",
    "    rhog = MWCE.addVars(customers,times,vtype=GRB.BINARY,name=\"rhog\")\n",
    "\n",
    "    MWCE.update()\n",
    "    #Adding Constraints:\n",
    "\n",
    "    #leader's\n",
    "    MWCE.addConstrs(\n",
    "        (Wl.sum('*',t) + Wh.sum('*',t) - X[t] == 0 for t in range(nr_t)), name=\"9a\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (quicksum(alpha_rate[s]*Y[s,t] for s in range(nr_fuel)) - X[t] == 0 for t in range(nr_t)), name=\"9b\")\n",
    "\n",
    "    #MWCE.addConstrs((Y.sum(s,'*') <= fuel_supply[s] for s in range(nr_fuel)), name=\"9c\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (PIl[t] - PIl[t+1] <= BIGM_price*P[t] for t in range(nr_t-1)), name=\"9d\")\n",
    "    MWCE.addConstrs(\n",
    "        (PIl[t] - PIl[t+1] >= -BIGM_price*P[t] for t in range(nr_t-1)), name=\"9d1\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (PIh[t] - PIh[t+1] <= BIGM_price*P[t] for t in range(nr_t-1)), name=\"9e\")\n",
    "    MWCE.addConstrs(\n",
    "        (PIh[t] - PIh[t+1] >= -BIGM_price*P[t] for t in range(nr_t-1)), name=\"9e1\")\n",
    "\n",
    "    MWCE.addConstr(\n",
    "        (quicksum(P[t] for t in range(nr_t)) <= caligraphy_n), name=\"9f\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (PIl[t] <= PIh[t] for t in range(nr_t)), name=\"9g\")\n",
    "\n",
    "    #Followers\n",
    "    MWCE.addConstrs(\n",
    "        (Wl[k,t]+Wh[k,t]+Vminus[k,t]+U[k,t]-Vplus[k,t] == demand_profile.iloc[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"10a\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (quicksum(Wl[k,t] + Wh[k,t] + U[k,t] for t in range(nr_t)) == demand[k] for k in range(nr_cust)), name=\"10b\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (Wl[k,t] <= psi_rate for k in range(nr_cust) for t in range(nr_t)), name=\"10c\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (Vplus[k,t] <= omega_rate.iloc[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"10d\")\n",
    "    \n",
    "    MWCE.update()\n",
    "\n",
    "    #KKT Conditions\n",
    "    \n",
    "    #Stationary Constraints\n",
    "    MWCE.addConstrs(\n",
    "        (MUa[k,t] + MUb[k] - MUc[k,t] <= PIl[t] for k in range(nr_cust) for t in range(nr_t)), name=\"11a\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (MUa[k,t] + MUb[k] <= PIh[t] for k in range(nr_cust) for t in range(nr_t)), name=\"11b\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (-MUa[k,t] - MUd[k,t] <= cost_xi.iloc[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11c\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (MUa[k,t] <= 0 for k in range(nr_cust) for t in range(nr_t)), name=\"11d\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (MUa[k,t] + MUb[k] <= cost_gamma[t] for k in range(nr_cust) for t in range(nr_t)), name=\"11e\")\n",
    "    \n",
    "    #Complementary Slackness Condition Constraints\n",
    "    MWCE.addConstrs(\n",
    "        (psi_rate - Wl[k,t] <= psi_rate*(1-rhoa[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11g\")\n",
    "    MWCE.addConstrs(\n",
    "        (MUc[k,t] <= BIGM_kkt1*rhoa[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11g1\")\n",
    "\n",
    "    MWCE.addConstrs(\n",
    "        (omega_rate.iloc[k,t]-Vplus[k,t] <= omega_rate.iloc[k,t]*(1-rhob[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11h\")\n",
    "    MWCE.addConstrs(\n",
    "        (MUd[k,t] <= BIGM_kkt1*rhob[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11h1\")\n",
    "    \n",
    "    MWCE.addConstrs(\n",
    "        (PIl[t] - MUa[k,t] - MUb[k] + MUc[k,t] <= BIGM_kkt2*(1-rhoc[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11i\")\n",
    "    MWCE.addConstrs(\n",
    "        (Wl[k,t] <= psi_rate*rhoc[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11i1\")\n",
    "    \n",
    "    MWCE.addConstrs(\n",
    "        (PIh[t] - MUa[k,t] - MUb[k]<= BIGM_kkt2*(1-rhod[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11j\")\n",
    "    MWCE.addConstrs(\n",
    "        (Wh[k,t] <= (demand_profile.iloc[k,t]+omega_rate.iloc[k,t]-psi_rate)*rhod[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11j1\")\n",
    "    \n",
    "    MWCE.addConstrs(\n",
    "        (cost_xi.iloc[k,t] + MUa[k,t] + MUd[k,t]  <= BIGM_kkt1*(1-rhoe[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11k\")\n",
    "    MWCE.addConstrs(\n",
    "        (Vplus[k,t] <= omega_rate.iloc[k,t]*rhoe[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11k1\")\n",
    "    \n",
    "    MWCE.addConstrs(\n",
    "        (0 - MUa[k,t] <= BIGM_kkt1*(1-rhof[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11l\")\n",
    "    MWCE.addConstrs(\n",
    "        (Vminus[k,t] <= demand_profile.iloc[k,t]*rhof[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11l1\")\n",
    "    \n",
    "    MWCE.addConstrs(\n",
    "        (cost_gamma[t] - MUa[k,t] - MUb[k] <= BIGM_kkt1*(1-rhog[k,t]) for k in range(nr_cust) for t in range(nr_t)), name=\"11m\")\n",
    "    MWCE.addConstrs(\n",
    "        (U[k,t] <= (demand_profile.iloc[k,t]+omega_rate.iloc[k,t])*rhog[k,t] for k in range(nr_cust) for t in range(nr_t)), name=\"11m1\")\n",
    "\n",
    "    #setting up the objective\n",
    "\n",
    "    MWCE.setObjective(\n",
    "        quicksum(demand_profile.iloc[k,t]*MUa[k,t] - psi_rate*MUc[k,t] - omega_rate.iloc[k,t]*MUd[k,t] - cost_gamma[t]*U[k,t]\n",
    "                     -cost_xi.iloc[k,t]*Vplus[k,t] for k in range(nr_cust) for t in range(nr_t)) +\n",
    "        quicksum(demand[k]*MUb[k] for k in range(nr_cust)) -\n",
    "        quicksum(cost_fuel[s]*Y[s,t] for s in range(nr_fuel) for t in range(nr_t)), GRB.MAXIMIZE)\n",
    "\n",
    "    MWCE.update()\n",
    "    \n",
    "    MWCE.reset(0)\n",
    "    print(f'--------------------------\\nStarted Solving MWCE Model {psi_rate}')\n",
    "    MWCE.optimize()\n",
    "    \n",
    "    MWCE_obj = MWCE.objVal\n",
    "    \n",
    "    \n",
    "    #--------Capturing Details of Leaders----------------\n",
    "    \n",
    "    print('\\n\\tGetting Results\\n--------------------------')\n",
    "    \n",
    "    #First create a list of all variables which are required to be output\n",
    "    pil = [PIl[t] for t in range(nr_t)];\n",
    "    pih = [PIh[t] for t in range(nr_t)];\n",
    "    x = [X[t] for t in range(nr_t)];\n",
    "    f1 = [Y[0,t] for t in range(nr_t)]; #Coal usage\n",
    "    f2 = [Y[1,t] for t in range(nr_t)]; #Gas usage\n",
    "    f3 = [Y[2,t] for t in range(nr_t)]; #Oil usage\n",
    "    w_l = [];\n",
    "    w_h = [];\n",
    "    u = []\n",
    "    vplus = []\n",
    "    vminus = []\n",
    "    for t in range(nr_t):\n",
    "        w_l.append(sum(MWCE.getAttr(GRB.Attr.X, [Wl[k,t] for k in range(nr_cust)])))\n",
    "        w_h.append(sum(MWCE.getAttr(GRB.Attr.X, [Wh[k,t] for k in range(nr_cust)])))\n",
    "        u.append(sum(MWCE.getAttr(GRB.Attr.X, [U[k,t] for k in range(nr_cust)])))\n",
    "        vplus.append(sum(MWCE.getAttr(GRB.Attr.X, [Vplus[k,t] for k in range(nr_cust)])))\n",
    "        vminus.append(sum(MWCE.getAttr(GRB.Attr.X, [Vminus[k,t] for k in range(nr_cust)])))\n",
    "    \n",
    "    df_leader = pd.DataFrame(index=range(nr_t));\n",
    "    df_leader.index.name = 'Day-hours'\n",
    "    \n",
    "    df_leader['Low_Price'] = MWCE.getAttr(GRB.Attr.X, pil);\n",
    "    df_leader['High_Price'] = MWCE.getAttr(GRB.Attr.X, pih);\n",
    "    df_leader['Elec_Generated'] = MWCE.getAttr(GRB.Attr.X, x);\n",
    "    df_leader['Low_Price Demand'] = w_l\n",
    "    df_leader['High_Price_Demand'] = w_h\n",
    "    df_leader['Upward_Demand_Shifted'] = vplus\n",
    "    df_leader['Downward_Demand_Shifted'] = vminus\n",
    "    df_leader['Comp_Demand'] = u\n",
    "    df_leader['Coal_Used'] = MWCE.getAttr(GRB.Attr.X, f1);\n",
    "    df_leader['Coal_Emissions (lbs)'] = df_leader['Coal_Used']*data.iloc[0,4];\n",
    "    df_leader['Gas_Used'] = MWCE.getAttr(GRB.Attr.X, f2);\n",
    "    df_leader['Gas_Emissions (lbs)'] = df_leader['Gas_Used']*data.iloc[1,4];\n",
    "    df_leader['Oil_Used'] = MWCE.getAttr(GRB.Attr.X, f3);\n",
    "    df_leader['Oil_Emissions (lbs)'] = df_leader['Oil_Used']*data.iloc[2,4];\n",
    "    \n",
    "    #--------Capturing Details of Followers----------------\n",
    "    cust = ['Cust'+str(j) for j in range(nr_cust)]\n",
    "    df_wh = pd.DataFrame(index=range(nr_t));\n",
    "    df_wh.index.name = 'Day-hours'\n",
    "    \n",
    "    df_wl = pd.DataFrame(index=range(nr_t));\n",
    "    df_wl.index.name = 'Day-hours'\n",
    "    \n",
    "    df_vplus = pd.DataFrame(index=range(nr_t));\n",
    "    df_vplus.index.name = 'Day-hours'\n",
    "    \n",
    "    df_vminus = pd.DataFrame(index=range(nr_t));\n",
    "    df_vminus.index.name = 'Day-hours'\n",
    "    \n",
    "    df_ukt = pd.DataFrame(index=range(nr_t));\n",
    "    df_ukt.index.name = 'Day-hours'\n",
    "    \n",
    "    for k in range(nr_cust):\n",
    "        df_wh[cust[k]] = MWCE.getAttr(GRB.Attr.X, [Wh[k,t] for t in range(nr_t)])\n",
    "        df_wl[cust[k]] = MWCE.getAttr(GRB.Attr.X, [Wl[k,t] for t in range(nr_t)])\n",
    "        df_vplus[cust[k]] = MWCE.getAttr(GRB.Attr.X, [Vplus[k,t] for t in range(nr_t)])\n",
    "        df_vminus[cust[k]] = MWCE.getAttr(GRB.Attr.X, [Vminus[k,t] for t in range(nr_t)])\n",
    "        df_ukt[cust[k]] = MWCE.getAttr(GRB.Attr.X, [U[k,t] for t in range(nr_t)])\n",
    "    \n",
    "    return MWCE_obj, df_leader.round(0), df_wh.round(0), df_wl.round(0), df_vplus.round(0), df_vminus.round(0), df_ukt.round(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computational Experiments\n",
    "\n",
    "\n",
    "* Define the values of BIG-M and $\\Psi$ for which both models will be solved"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "BIGM_price=100; #Big M for capturing the price changes\n",
    "BIGM_kkt1=20;  #Big M for McCormic Constraints-MUd\n",
    "BIGM_kkt2=50;  #Big M for McCormic Constraints for Mua, MUb, MUc\n",
    "BIGM_kkt3=200;  #Big M for McCormic Constraints-MUa\n",
    "psi_labels = ['Low \\u03C8', 'Medium \\u03C8', 'High \\u03C8']\n",
    "psi_values = [5000, 70000,180000] #low, Medium, High Psi Values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* __Store the Results__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "Profits_Actual_MCE = {}\n",
    "Profits_Wghited_MCE = {}\n",
    "Emission_Cost_MCE = {}\n",
    "Retailer_Profits_MCE = {}\n",
    "\n",
    "Profits_Actual_MWCE = {}\n",
    "Profits_Wghited_MWCE = {}\n",
    "Emission_Cost_MWCE = {}\n",
    "Retailer_Profits_MWCE = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve for `Low` $\\Psi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Psi Value: 5000\n",
      "Discarded solution information\n",
      "--------------------------\n",
      "Started Solving MCE Model\n",
      "Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (win64)\n",
      "Optimize a model with 26615 rows, 18218 columns and 71168 nonzeros\n",
      "Model fingerprint: 0x70f749fa\n",
      "Variable types: 9794 continuous, 8424 integer (8424 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 2e+05]\n",
      "  Objective range  [1e-02, 3e+06]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e-02, 3e+06]\n",
      "Presolve removed 6274 rows and 2576 columns\n",
      "Presolve time: 0.52s\n",
      "Presolved: 20341 rows, 15642 columns, 58457 nonzeros\n",
      "Variable types: 8256 continuous, 7386 integer (7386 binary)\n",
      "\n",
      "Root relaxation: objective 1.537215e+08, 7662 iterations, 0.28 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 1.5372e+08    0 1107          - 1.5372e+08      -     -    5s\n",
      "     0     0 1.5353e+08    0  631          - 1.5353e+08      -     -    7s\n",
      "     0     0 1.5352e+08    0  681          - 1.5352e+08      -     -    8s\n",
      "     0     0 1.5352e+08    0  676          - 1.5352e+08      -     -    8s\n",
      "     0     0 1.5352e+08    0  677          - 1.5352e+08      -     -    8s\n",
      "     0     0 1.5347e+08    0  401          - 1.5347e+08      -     -   10s\n",
      "     0     0 1.5347e+08    0  509          - 1.5347e+08      -     -   10s\n",
      "     0     0 1.5347e+08    0  497          - 1.5347e+08      -     -   11s\n",
      "     0     0 1.5347e+08    0  525          - 1.5347e+08      -     -   11s\n",
      "     0     0 1.5347e+08    0  510          - 1.5347e+08      -     -   11s\n",
      "     0     0 1.5345e+08    0  267          - 1.5345e+08      -     -   15s\n",
      "     0     0 1.5345e+08    0  299          - 1.5345e+08      -     -   39s\n",
      "     0     0 1.5345e+08    0  288          - 1.5345e+08      -     -   40s\n",
      "     0     0 1.5345e+08    0  207          - 1.5345e+08      -     -   42s\n",
      "     0     0 1.5345e+08    0  238          - 1.5345e+08      -     -   42s\n",
      "     0     0 1.5345e+08    0  245          - 1.5345e+08      -     -   42s\n",
      "     0     0 1.5344e+08    0  229          - 1.5344e+08      -     -   44s\n",
      "     0     0 1.5344e+08    0  204          - 1.5344e+08      -     -   46s\n",
      "     0     2 1.5344e+08    0  204          - 1.5344e+08      -     -   55s\n",
      "   250   293 1.5343e+08   49  126          - 1.5344e+08      -  47.8   60s\n",
      "   917  1036 1.5235e+08  224  342          - 1.5344e+08      -  30.8   65s\n",
      "  1469  1564 1.5343e+08   29  143          - 1.5344e+08      -  25.7   70s\n",
      "  1940  2023 1.5333e+08  105   90          - 1.5344e+08      -  40.1   75s\n",
      "  2195  2119 1.5223e+08  217  963          - 1.5344e+08      -  44.3   82s\n",
      "  2196  2120 1.5319e+08    6  341          - 1.5344e+08      -  44.2   85s\n",
      "  2200  2123 1.5341e+08   34  236          - 1.5344e+08      -  44.2   92s\n",
      "  2204  2125 1.5223e+08  193  234          - 1.5344e+08      -  44.1   96s\n",
      "  2207  2127 1.5217e+08  389  214          - 1.5344e+08      -  44.0  102s\n",
      "  2210  2129 1.5339e+08  139  246          - 1.5344e+08      -  44.0  107s\n",
      "  2213  2131 1.5341e+08   89  322          - 1.5344e+08      -  43.9  111s\n",
      "  2217  2134 1.5318e+08   19  301          - 1.5343e+08      -  43.8  116s\n",
      "  2221  2137 1.5217e+08  380  203          - 1.5343e+08      -  43.7  168s\n",
      "  2224  2139 1.5222e+08  321  247          - 1.5343e+08      -  43.7  170s\n",
      "  2227  2141 1.5202e+08  386  236          - 1.5343e+08      -  43.6  175s\n",
      "  2234  2145 1.5223e+08  187  184          - 1.5343e+08      -  43.5  181s\n",
      "  2238  2148 1.5223e+08  162  123          - 1.5342e+08      -  43.4  185s\n",
      "  2242  2151 1.5340e+08   16  154          - 1.5342e+08      -  43.3  190s\n",
      "  2245  2153 1.5222e+08  240  177          - 1.5342e+08      -  43.3  195s\n",
      "H 2245  2044                    1.534220e+08 1.5342e+08  0.00%  43.3  206s\n",
      "\n",
      "Cutting planes:\n",
      "  Learned: 23\n",
      "  Gomory: 248\n",
      "  Cover: 5\n",
      "  Implied bound: 66\n",
      "  Projected implied bound: 13\n",
      "  Clique: 32\n",
      "  MIR: 838\n",
      "  Flow cover: 1852\n",
      "  Network: 1\n",
      "  RLT: 38\n",
      "  Relax-and-lift: 1151\n",
      "\n",
      "Explored 2245 nodes (149425 simplex iterations) in 206.59 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 1.53422e+08 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-02)\n",
      "Best objective 1.534220296879e+08, best bound 1.534234854356e+08, gap 0.0009%\n",
      "Optimality Gap: 0.00\n",
      "\n",
      "\tGetting Results\n",
      "--------------------------\n",
      "Discarded solution information\n",
      "--------------------------\n",
      "Started Solving MWCE Model 5000\n",
      "Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (win64)\n",
      "Optimize a model with 26615 rows, 18218 columns and 71168 nonzeros\n",
      "Model fingerprint: 0xa00084c1\n",
      "Variable types: 9794 continuous, 8424 integer (8424 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 2e+05]\n",
      "  Objective range  [1e-02, 3e+06]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e-02, 3e+06]\n",
      "Presolve removed 6274 rows and 2576 columns\n",
      "Presolve time: 0.23s\n",
      "Presolved: 20341 rows, 15642 columns, 58457 nonzeros\n",
      "Variable types: 8256 continuous, 7386 integer (7386 binary)\n",
      "\n",
      "Root relaxation: objective 2.483412e+08, 7595 iterations, 0.11 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 2.4834e+08    0  881          - 2.4834e+08      -     -    2s\n",
      "     0     0 2.4815e+08    0  616          - 2.4815e+08      -     -   14s\n",
      "     0     0 2.4814e+08    0  591          - 2.4814e+08      -     -   15s\n",
      "     0     0 2.4814e+08    0  583          - 2.4814e+08      -     -   15s\n",
      "     0     0 2.4810e+08    0  428          - 2.4810e+08      -     -   16s\n",
      "     0     0 2.4810e+08    0  547          - 2.4810e+08      -     -   17s\n",
      "     0     0 2.4809e+08    0  544          - 2.4809e+08      -     -   17s\n",
      "     0     0 2.4809e+08    0  558          - 2.4809e+08      -     -   17s\n",
      "     0     0 2.4809e+08    0  568          - 2.4809e+08      -     -   17s\n",
      "     0     0 2.4807e+08    0  306          - 2.4807e+08      -     -   19s\n",
      "     0     0 2.4807e+08    0  305          - 2.4807e+08      -     -   20s\n",
      "     0     0 2.4807e+08    0  314          - 2.4807e+08      -     -   20s\n",
      "     0     0 2.4807e+08    0  270          - 2.4807e+08      -     -   22s\n",
      "     0     0 2.4807e+08    0  291          - 2.4807e+08      -     -   22s\n",
      "     0     0 2.4807e+08    0  299          - 2.4807e+08      -     -   23s\n",
      "     0     0 2.4806e+08    0  308          - 2.4806e+08      -     -   24s\n",
      "     0     0 2.4806e+08    0  307          - 2.4806e+08      -     -   25s\n",
      "     0     2 2.4806e+08    0  275          - 2.4806e+08      -     -   33s\n",
      "    55    88 2.4806e+08   17  203          - 2.4806e+08      -  93.2   35s\n",
      "   432   618 2.4689e+08   59  452          - 2.4806e+08      -  55.4   40s\n",
      "  1200  1301 2.4685e+08  285  172          - 2.4806e+08      -  35.4   45s\n",
      "  1669  1768 2.4805e+08   54   96          - 2.4806e+08      -  37.7   50s\n",
      "  2017  2157 2.4772e+08  134  190          - 2.4806e+08      -  47.1   55s\n",
      "  2207  2159 2.4684e+08  247  968          - 2.4806e+08      -  47.5   61s\n",
      "  2215  2165 2.4803e+08   70  290          - 2.4806e+08      -  47.3   66s\n",
      "  2219  2167 2.4805e+08   28  300          - 2.4806e+08      -  47.2   70s\n",
      "  2220  2168 2.4806e+08    6  264          - 2.4806e+08      -  47.2   80s\n",
      "  2224  2171 2.4772e+08  165  273          - 2.4805e+08      -  47.1   88s\n",
      "  2226  2172 2.4683e+08  115  290          - 2.4805e+08      -  47.1   91s\n",
      "  2227  2173 2.4665e+08  196  252          - 2.4805e+08      -  47.0   96s\n",
      "  2230  2175 2.4691e+08   43  265          - 2.4805e+08      -  47.0  100s\n",
      "  2235  2178 2.4803e+08   19  266          - 2.4805e+08      -  46.9  105s\n",
      "  2240  2181 2.4684e+08  260  378          - 2.4805e+08      -  46.8  110s\n",
      "  2244  2184 2.4772e+08  143  191          - 2.4805e+08      -  46.7  115s\n",
      "  2248  2187 2.4772e+08   79  229          - 2.4804e+08      -  46.6  121s\n",
      "  2249  2187 2.4768e+08  238  156          - 2.4804e+08      -  46.6  166s\n",
      "  2253  2190 2.4716e+08  264   97          - 2.4804e+08      -  46.5  170s\n",
      "  2258  2193 2.4684e+08  220   67          - 2.4804e+08      -  46.4  175s\n",
      "  2261  2195 2.4696e+08  272   61          - 2.4804e+08      -  46.3  180s\n",
      "  2268  2201 2.4734e+08  149  968          - 2.4804e+08      -  61.5  189s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  2269  2202 2.4685e+08  123  410          - 2.4804e+08      -  61.5  193s\n",
      "  2273  2205 2.4771e+08  101   67          - 2.4804e+08      -  61.3  196s\n",
      "  2275  2206 2.4699e+08  182   65          - 2.4804e+08      -  61.3  200s\n",
      "  2276  2207 2.4802e+08   14   66          - 2.4804e+08      -  61.3  207s\n",
      "H 2277  2096                    2.480417e+08 2.4804e+08  0.00%  61.2  214s\n",
      "\n",
      "Cutting planes:\n",
      "  Learned: 8\n",
      "  Gomory: 169\n",
      "  Cover: 7\n",
      "  Implied bound: 14\n",
      "  Projected implied bound: 9\n",
      "  Clique: 7\n",
      "  MIR: 652\n",
      "  Flow cover: 1553\n",
      "  Network: 1\n",
      "  RLT: 4\n",
      "  Relax-and-lift: 480\n",
      "\n",
      "Explored 2277 nodes (173684 simplex iterations) in 214.71 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 2.48042e+08 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-02)\n",
      "Best objective 2.480416660706e+08, best bound 2.480430623515e+08, gap 0.0006%\n",
      "\n",
      "\tGetting Results\n",
      "--------------------------\n"
     ]
    }
   ],
   "source": [
    "psi_rate = psi_values[0]\n",
    "\n",
    "print(f'Psi Value: {psi_rate}')\n",
    "MCE_obj,low_MCE_leader,low_MCE_wh,low_MCE_wl,low_MCE_vplus,low_MCE_vminus,low_MCE_ukt = with_emission(nr_t,nr_fuel,nr_cust,alpha_rate,fuel_supply,epsilon_rate,BIGM_price,caligraphy_n,beta_rate,\n",
    "                  demand_profile,demand,psi_rate,omega_rate,cost_xi,cost_gamma,BIGM_kkt1,BIGM_kkt2,BIGM_kkt3,\n",
    "                  cost_fuel,cost_ramp,cost_emission)\n",
    "\n",
    "Profits_Actual_MCE['Low'] = MCE_obj\n",
    "\n",
    "MWCE_obj,low_MWCE_leader,low_MWCE_wh,low_MWCE_wl,low_MWCE_vplus,low_MWCE_vminus,low_MWCE_ukt = without_emission(nr_t,nr_fuel,nr_cust,alpha_rate,fuel_supply,epsilon_rate,BIGM_price,caligraphy_n,beta_rate,\n",
    "                  demand_profile,demand,psi_rate,omega_rate,cost_xi,cost_gamma,BIGM_kkt1,BIGM_kkt2,BIGM_kkt3,\n",
    "                  cost_fuel,cost_ramp,cost_emission)\n",
    "\n",
    "Profits_Actual_MWCE['Low'] = MWCE_obj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_MCE_leader.to_csv('low_MCE_leader.csv')\n",
    "low_MCE_wh.to_csv('low_MCE_wh.csv')\n",
    "low_MCE_wl.to_csv('low_MCE_wl.csv')\n",
    "low_MCE_vplus.to_csv('low_MCE_vplus.csv')\n",
    "low_MCE_vminus.to_csv('low_MCE_vminus.csv')\n",
    "low_MCE_ukt.to_csv('low_MCE_ukt.csv')\n",
    "\n",
    "low_MWCE_leader.to_csv('low_MWCE_leader.csv')\n",
    "low_MWCE_wh.to_csv('low_MWCE_wh.csv')\n",
    "low_MWCE_wl.to_csv('low_MWCE_wl.csv')\n",
    "low_MWCE_vplus.to_csv('low_MWCE_vplus.csv')\n",
    "low_MWCE_vminus.to_csv('low_MWCE_vminus.csv')\n",
    "low_MWCE_ukt.to_csv('low_MWCE_ukt.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve for `Medium` $\\Psi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Psi Value: 70000\n",
      "Discarded solution information\n",
      "--------------------------\n",
      "Started Solving MCE Model\n",
      "Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (win64)\n",
      "Optimize a model with 26615 rows, 18218 columns and 71168 nonzeros\n",
      "Model fingerprint: 0xb452fb7e\n",
      "Variable types: 9794 continuous, 8424 integer (8424 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 2e+05]\n",
      "  Objective range  [1e-02, 3e+06]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e-02, 3e+06]\n",
      "Presolve removed 9400 rows and 5721 columns\n",
      "Presolve time: 0.21s\n",
      "Presolved: 17215 rows, 12497 columns, 49586 nonzeros\n",
      "Variable types: 6381 continuous, 6116 integer (6116 binary)\n",
      "\n",
      "Root relaxation: objective 1.537224e+08, 6898 iterations, 0.10 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 1.5372e+08    0  965          - 1.5372e+08      -     -    2s\n",
      "     0     0 1.5353e+08    0  530          - 1.5353e+08      -     -    3s\n",
      "     0     0 1.5352e+08    0  567          - 1.5352e+08      -     -    3s\n",
      "     0     0 1.5352e+08    0  562          - 1.5352e+08      -     -    3s\n",
      "     0     0 1.5348e+08    0  385          - 1.5348e+08      -     -    4s\n",
      "     0     0 1.5347e+08    0  516          - 1.5347e+08      -     -    4s\n",
      "     0     0 1.5347e+08    0  480          - 1.5347e+08      -     -    4s\n",
      "     0     0 1.5347e+08    0  503          - 1.5347e+08      -     -    5s\n",
      "     0     0 1.5347e+08    0  503          - 1.5347e+08      -     -    5s\n",
      "     0     0 1.5347e+08    0  504          - 1.5347e+08      -     -    5s\n",
      "     0     0 1.5345e+08    0  254          - 1.5345e+08      -     -    6s\n",
      "     0     0 1.5345e+08    0  268          - 1.5345e+08      -     -   20s\n",
      "     0     0 1.5345e+08    0  274          - 1.5345e+08      -     -   20s\n",
      "     0     0 1.5345e+08    0  272          - 1.5345e+08      -     -   20s\n",
      "     0     0 1.5345e+08    0  205          - 1.5345e+08      -     -   21s\n",
      "     0     0 1.5345e+08    0  217          - 1.5345e+08      -     -   21s\n",
      "     0     0 1.5345e+08    0  253          - 1.5345e+08      -     -   22s\n",
      "     0     0 1.5345e+08    0  221          - 1.5345e+08      -     -   22s\n",
      "     0     0 1.5345e+08    0  202          - 1.5345e+08      -     -   23s\n",
      "     0     2 1.5345e+08    0  201          - 1.5345e+08      -     -   27s\n",
      "   467   586 1.5343e+08   59  111          - 1.5344e+08      -  17.6   30s\n",
      "  1402  1490 1.5240e+08  319  239          - 1.5344e+08      -  57.6   35s\n",
      "  2012  2055 1.5215e+08  489  187          - 1.5344e+08      -  63.1   40s\n",
      "  2431  2290 1.5231e+08  154  887          - 1.5344e+08      -  68.8   46s\n",
      "  2437  2294 1.5200e+08  336  198          - 1.5344e+08      -  68.6   69s\n",
      "  2438  2295 1.5228e+08  237  200          - 1.5344e+08      -  68.6   70s\n",
      "  2443  2298 1.5231e+08  246  212          - 1.5344e+08      -  68.5   75s\n",
      "  2449  2302 1.5225e+08  332  160          - 1.5343e+08      -  68.3   80s\n",
      "  2455  2306 1.5341e+08   49  142          - 1.5342e+08      -  68.1   85s\n",
      "  2460  2310 1.5229e+08  207   72          - 1.5342e+08      -  68.0   90s\n",
      "  2464  2312 1.5239e+08  332   96          - 1.5342e+08      -  67.9   95s\n",
      "  2471  2318 1.5231e+08  247  887          - 1.5342e+08      -  74.9  103s\n",
      "  2472  2319 1.5342e+08   17  131          - 1.5342e+08      -  74.8  105s\n",
      "  2479  2324 1.5224e+08  475   64          - 1.5342e+08      -  74.6  110s\n",
      "  2483  2326 1.5225e+08  314   50          - 1.5342e+08      -  74.5  116s\n",
      "  2487  2329 1.5139e+08  345   51          - 1.5342e+08      -  74.4  120s\n",
      "H 2490  2214                    1.534220e+08 1.5342e+08  0.00%  74.3  126s\n",
      "\n",
      "Cutting planes:\n",
      "  Learned: 12\n",
      "  Gomory: 236\n",
      "  Implied bound: 3\n",
      "  MIR: 411\n",
      "  Flow cover: 1193\n",
      "  RLT: 5\n",
      "  Relax-and-lift: 343\n",
      "\n",
      "Explored 2490 nodes (212440 simplex iterations) in 126.50 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 1.53422e+08 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-02)\n",
      "Best objective 1.534220296879e+08, best bound 1.534231105524e+08, gap 0.0007%\n",
      "Optimality Gap: 0.00\n",
      "\n",
      "\tGetting Results\n",
      "--------------------------\n",
      "Discarded solution information\n",
      "--------------------------\n",
      "Started Solving MWCE Model 70000\n",
      "Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (win64)\n",
      "Optimize a model with 26615 rows, 18218 columns and 71168 nonzeros\n",
      "Model fingerprint: 0x4aae7e3c\n",
      "Variable types: 9794 continuous, 8424 integer (8424 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 2e+05]\n",
      "  Objective range  [1e-02, 3e+06]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e-02, 3e+06]\n",
      "Presolve removed 9400 rows and 5721 columns\n",
      "Presolve time: 0.30s\n",
      "Presolved: 17215 rows, 12497 columns, 49586 nonzeros\n",
      "Variable types: 6381 continuous, 6116 integer (6116 binary)\n",
      "\n",
      "Root relaxation: objective 2.483420e+08, 6772 iterations, 0.18 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 2.4834e+08    0  983          - 2.4834e+08      -     -    3s\n",
      "     0     0 2.4815e+08    0  556          - 2.4815e+08      -     -    5s\n",
      "     0     0 2.4814e+08    0  603          - 2.4814e+08      -     -    6s\n",
      "     0     0 2.4814e+08    0  594          - 2.4814e+08      -     -    6s\n",
      "     0     0 2.4814e+08    0  598          - 2.4814e+08      -     -    6s\n",
      "     0     0 2.4810e+08    0  357          - 2.4810e+08      -     -    8s\n",
      "     0     0 2.4809e+08    0  512          - 2.4809e+08      -     -    8s\n",
      "     0     0 2.4809e+08    0  485          - 2.4809e+08      -     -    8s\n",
      "     0     0 2.4809e+08    0  510          - 2.4809e+08      -     -    9s\n",
      "     0     0 2.4809e+08    0  511          - 2.4809e+08      -     -    9s\n",
      "     0     0 2.4807e+08    0  282          - 2.4807e+08      -     -   11s\n",
      "     0     0 2.4807e+08    0  286          - 2.4807e+08      -     -   40s\n",
      "     0     0 2.4807e+08    0  299          - 2.4807e+08      -     -   40s\n",
      "     0     0 2.4807e+08    0  290          - 2.4807e+08      -     -   40s\n",
      "     0     0 2.4807e+08    0  245          - 2.4807e+08      -     -   42s\n",
      "     0     0 2.4807e+08    0  234          - 2.4807e+08      -     -   42s\n",
      "     0     0 2.4807e+08    0  241          - 2.4807e+08      -     -   42s\n",
      "     0     0 2.4807e+08    0  243          - 2.4807e+08      -     -   42s\n",
      "     0     0 2.4806e+08    0  221          - 2.4806e+08      -     -   43s\n",
      "     0     0 2.4806e+08    0  205          - 2.4806e+08      -     -   44s\n",
      "     0     2 2.4806e+08    0  205          - 2.4806e+08      -     -   50s\n",
      "   651   798 2.4686e+08  181  294          - 2.4806e+08      -  33.0   55s\n",
      "  1565  1716 2.4806e+08    9  222          - 2.4806e+08      -  38.0   60s\n",
      "  2300  2240 2.4681e+08   83  890          - 2.4806e+08      -  41.9   68s\n",
      "  2303  2242 2.4670e+08  369  211          - 2.4806e+08      -  41.8   70s\n",
      "  2305  2244 2.4683e+08  345  241          - 2.4806e+08      -  41.8   77s\n",
      "  2308  2246 2.4686e+08  133  253          - 2.4805e+08      -  41.7   80s\n",
      "  2311  2248 2.4658e+08  296  195          - 2.4805e+08      -  41.7   85s\n",
      "  2316  2251 2.4682e+08  379   81          - 2.4804e+08      -  41.6   90s\n",
      "  2321  2254 2.4684e+08  216   82          - 2.4804e+08      -  41.5   95s\n",
      "  2326  2258 2.4670e+08  394   50          - 2.4804e+08      -  41.4  100s\n",
      "  2332  2263 2.4758e+08    5  890          - 2.4804e+08      -  48.4  108s\n",
      "  2333  2264 2.4657e+08  294  167          - 2.4804e+08      -  48.3  110s\n",
      "  2337  2267 2.4804e+08   28   73          - 2.4804e+08      -  48.3  115s\n",
      "  2339  2268 2.4684e+08   70   69          - 2.4804e+08      -  48.2  121s\n",
      "  2342  2270 2.4684e+08  257   45          - 2.4804e+08      -  48.2  126s\n",
      "  2343  2271 2.4687e+08  129   44          - 2.4804e+08      -  48.1  130s\n",
      "H 2344  2157                    2.480417e+08 2.4804e+08  0.00%  48.1  136s\n",
      "\n",
      "Cutting planes:\n",
      "  Learned: 19\n",
      "  Gomory: 188\n",
      "  Cover: 1\n",
      "  Implied bound: 4\n",
      "  Projected implied bound: 1\n",
      "  MIR: 349\n",
      "  Flow cover: 1364\n",
      "  RLT: 11\n",
      "  Relax-and-lift: 501\n",
      "\n",
      "Explored 2344 nodes (141100 simplex iterations) in 136.08 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 2.48042e+08 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-02)\n",
      "Best objective 2.480416660706e+08, best bound 2.480427545044e+08, gap 0.0004%\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\tGetting Results\n",
      "--------------------------\n"
     ]
    }
   ],
   "source": [
    "psi_rate = psi_values[1]\n",
    "\n",
    "print(f'Psi Value: {psi_rate}')\n",
    "MCE_obj,med_MCE_leader,med_MCE_wh,med_MCE_wl,med_MCE_vplus,med_MCE_vminus,med_MCE_ukt = with_emission(nr_t,nr_fuel,nr_cust,alpha_rate,fuel_supply,epsilon_rate,BIGM_price,caligraphy_n,beta_rate,\n",
    "                  demand_profile,demand,psi_rate,omega_rate,cost_xi,cost_gamma,BIGM_kkt1,BIGM_kkt2,BIGM_kkt3,\n",
    "                  cost_fuel,cost_ramp,cost_emission)\n",
    "\n",
    "Profits_Actual_MCE['Medium'] = MCE_obj\n",
    "\n",
    "MWCE_obj,med_MWCE_leader,med_MWCE_wh,med_MWCE_wl,med_MWCE_vplus,med_MWCE_vminus,med_MWCE_ukt = without_emission(nr_t,nr_fuel,nr_cust,alpha_rate,fuel_supply,epsilon_rate,BIGM_price,caligraphy_n,beta_rate,\n",
    "                  demand_profile,demand,psi_rate,omega_rate,cost_xi,cost_gamma,BIGM_kkt1,BIGM_kkt2,BIGM_kkt3,\n",
    "                  cost_fuel,cost_ramp,cost_emission)\n",
    "\n",
    "Profits_Actual_MWCE['Medium'] = MWCE_obj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "med_MCE_leader.to_csv('med_MCE_leader.csv')\n",
    "med_MCE_wh.to_csv('med_MCE_wh.csv')\n",
    "med_MCE_wl.to_csv('med_MCE_wl.csv')\n",
    "med_MCE_vplus.to_csv('med_MCE_vplus.csv')\n",
    "med_MCE_vminus.to_csv('med_MCE_vminus.csv')\n",
    "med_MCE_ukt.to_csv('med_MCE_ukt.csv')\n",
    "\n",
    "med_MWCE_leader.to_csv('med_MWCE_leader.csv')\n",
    "med_MWCE_wh.to_csv('med_MWCE_wh.csv')\n",
    "med_MWCE_wl.to_csv('med_MWCE_wl.csv')\n",
    "med_MWCE_vplus.to_csv('med_MWCE_vplus.csv')\n",
    "med_MWCE_vminus.to_csv('med_MWCE_vminus.csv')\n",
    "med_MWCE_ukt.to_csv('med_MWCE_ukt.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve for `High` $\\Psi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Psi Value: 180000\n",
      "Discarded solution information\n",
      "--------------------------\n",
      "Started Solving MCE Model\n",
      "Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (win64)\n",
      "Optimize a model with 26615 rows, 18218 columns and 71168 nonzeros\n",
      "Model fingerprint: 0x1f391335\n",
      "Variable types: 9794 continuous, 8424 integer (8424 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 2e+05]\n",
      "  Objective range  [1e-02, 3e+06]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e-02, 3e+06]\n",
      "Presolve removed 10674 rows and 6303 columns\n",
      "Presolve time: 0.38s\n",
      "Presolved: 15941 rows, 11915 columns, 45900 nonzeros\n",
      "Variable types: 6029 continuous, 5886 integer (5886 binary)\n",
      "Found heuristic solution: objective 1.128395e+08\n",
      "\n",
      "Root relaxation: objective 1.537216e+08, 6672 iterations, 0.16 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 1.5372e+08    0  917 1.1284e+08 1.5372e+08  36.2%     -    3s\n",
      "     0     0 1.5353e+08    0  548 1.1284e+08 1.5353e+08  36.1%     -   17s\n",
      "     0     0 1.5353e+08    0  577 1.1284e+08 1.5353e+08  36.1%     -   18s\n",
      "     0     0 1.5353e+08    0  573 1.1284e+08 1.5353e+08  36.1%     -   18s\n",
      "     0     0 1.5353e+08    0  576 1.1284e+08 1.5353e+08  36.1%     -   18s\n",
      "     0     0 1.5348e+08    0  381 1.1284e+08 1.5348e+08  36.0%     -   19s\n",
      "     0     0 1.5347e+08    0  485 1.1284e+08 1.5347e+08  36.0%     -   19s\n",
      "     0     0 1.5347e+08    0  456 1.1284e+08 1.5347e+08  36.0%     -   19s\n",
      "     0     0 1.5347e+08    0  472 1.1284e+08 1.5347e+08  36.0%     -   19s\n",
      "     0     0 1.5347e+08    0  472 1.1284e+08 1.5347e+08  36.0%     -   19s\n",
      "     0     0 1.5345e+08    0  277 1.1284e+08 1.5345e+08  36.0%     -   20s\n",
      "     0     0 1.5345e+08    0  294 1.1284e+08 1.5345e+08  36.0%     -   20s\n",
      "     0     0 1.5345e+08    0  294 1.1284e+08 1.5345e+08  36.0%     -   21s\n",
      "     0     0 1.5345e+08    0  242 1.1284e+08 1.5345e+08  36.0%     -   21s\n",
      "     0     0 1.5345e+08    0  235 1.1284e+08 1.5345e+08  36.0%     -   22s\n",
      "     0     0 1.5345e+08    0  235 1.1284e+08 1.5345e+08  36.0%     -   22s\n",
      "     0     0 1.5345e+08    0  221 1.1284e+08 1.5345e+08  36.0%     -   22s\n",
      "     0     0 1.5345e+08    0  221 1.1284e+08 1.5345e+08  36.0%     -   23s\n",
      "     0     2 1.5345e+08    0  221 1.1284e+08 1.5345e+08  36.0%     -   24s\n",
      "    31    36 1.5344e+08    8  179 1.1284e+08 1.5344e+08  36.0%   109   38s\n",
      "   211   299 1.5340e+08   31  131 1.1284e+08 1.5344e+08  36.0%  35.0   40s\n",
      "   613   764 1.5228e+08  102  329 1.1284e+08 1.5344e+08  36.0%  53.6   45s\n",
      "  1245  1385 1.5344e+08    6  204 1.1284e+08 1.5344e+08  36.0%  67.4   51s\n",
      "  1571  1777 1.5343e+08   59   70 1.1284e+08 1.5344e+08  36.0%  76.0   56s\n",
      "  2174  2101 1.5343e+08   31  221 1.1284e+08 1.5344e+08  36.0%  74.0   65s\n",
      "  2176  2102 1.5239e+08  121  880 1.1284e+08 1.5343e+08  36.0%  73.9   94s\n",
      "  2177  2103 1.5229e+08   61  755 1.1284e+08 1.5343e+08  36.0%  73.9   95s\n",
      "  2187  2110 1.5236e+08  182  169 1.1284e+08 1.5343e+08  36.0%  73.5  100s\n",
      "  2194  2114 1.5235e+08  223  110 1.1284e+08 1.5343e+08  36.0%  73.3  105s\n",
      "  2199  2118 1.5343e+08   36  126 1.1284e+08 1.5343e+08  36.0%  73.1  115s\n",
      "  2205  2122 1.5342e+08    7  129 1.1284e+08 1.5342e+08  36.0%  72.9  120s\n",
      "  2209  2124 1.5239e+08  112   78 1.1284e+08 1.5342e+08  36.0%  72.8  125s\n",
      "  2211  2126 1.5224e+08  305   77 1.1284e+08 1.5342e+08  36.0%  72.7  130s\n",
      "  2214  2128 1.5217e+08   92   82 1.1284e+08 1.5342e+08  36.0%  72.6  137s\n",
      "  2216  2129 1.5172e+08  327   63 1.1284e+08 1.5342e+08  36.0%  72.6  140s\n",
      "  2221  2132 1.1361e+08   27   48 1.1284e+08 1.5342e+08  36.0%  72.4  145s\n",
      "  2225  2136 1.5341e+08   26  880 1.1284e+08 1.5342e+08  36.0%  81.0  176s\n",
      "  2232  2141 1.5338e+08   14   54 1.1284e+08 1.5342e+08  36.0%  80.7  180s\n",
      "  2236  2147 1.5342e+08   29   38 1.1284e+08 1.5342e+08  36.0%  86.7  185s\n",
      "  2238  2150 1.5342e+08   30   87 1.1284e+08 1.5342e+08  36.0%  86.7  208s\n",
      "  2246  2155 1.5234e+08   31 1041 1.1284e+08 1.5342e+08  36.0%  88.8  213s\n",
      "  2266  2169 1.5342e+08   34   23 1.1284e+08 1.5342e+08  36.0%  88.3  232s\n",
      "  2298  2193 1.5337e+08   38   38 1.1284e+08 1.5342e+08  36.0%  88.6  275s\n",
      "  2579  2387 1.5244e+08   59  403 1.1284e+08 1.5342e+08  36.0%  94.3  280s\n",
      "  3182  2572 1.5242e+08  117  338 1.1284e+08 1.5342e+08  36.0%  86.1  285s\n",
      "  3551  2788 1.5240e+08  153  291 1.1284e+08 1.5342e+08  36.0%  88.9  290s\n",
      "  3684  2859 1.5240e+08  174  265 1.1284e+08 1.5342e+08  36.0%  93.0  295s\n",
      "  3777  2856 1.5239e+08  185  258 1.1284e+08 1.5342e+08  36.0%  96.5  306s\n",
      "  4000  3025 1.5239e+08  207  195 1.1284e+08 1.5342e+08  36.0%   102  310s\n",
      "  4380  3225 1.5239e+08  268  116 1.1284e+08 1.5342e+08  36.0%   110  316s\n",
      "  4474  3219 infeasible  288      1.1284e+08 1.5342e+08  36.0%   112  353s\n",
      "  4540  3221 1.5245e+08   49  383 1.1284e+08 1.5342e+08  36.0%   114  357s\n",
      "  4682  3424 1.5215e+08  164  273 1.1284e+08 1.5342e+08  36.0%   113  360s\n",
      "  5276  3624 1.5039e+08  407  523 1.1284e+08 1.5342e+08  36.0%   109  366s\n",
      "  5509  3675     cutoff  449      1.1284e+08 1.5342e+08  36.0%   111  456s\n",
      "  5673  3991 1.5297e+08   78  576 1.1284e+08 1.5342e+08  36.0%   112  505s\n",
      "  5962  4006 1.5245e+08  216  448 1.1284e+08 1.5342e+08  36.0%   108  510s\n",
      "  6176  4030 1.3551e+08   56  894 1.1284e+08 1.5342e+08  36.0%   113  516s\n",
      "  6251  4065 1.2363e+08   60  627 1.1284e+08 1.5342e+08  36.0%   116  531s\n",
      "  6316  4247 1.1389e+08   61  821 1.1284e+08 1.5342e+08  36.0%   116  536s\n",
      "H 6326  1089                    1.534220e+08 1.5342e+08  0.00%   116  536s\n",
      "\n",
      "Cutting planes:\n",
      "  Learned: 34\n",
      "  Gomory: 273\n",
      "  Cover: 28\n",
      "  Implied bound: 12\n",
      "  Clique: 9\n",
      "  MIR: 561\n",
      "  Flow cover: 2305\n",
      "  Network: 1\n",
      "  RLT: 9\n",
      "  Relax-and-lift: 389\n",
      "\n",
      "Explored 6551 nodes (798031 simplex iterations) in 536.36 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 2: 1.53422e+08 1.12839e+08 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-02)\n",
      "Best objective 1.534220296879e+08, best bound 1.534220296879e+08, gap 0.0000%\n",
      "Optimality Gap: 0.00\n",
      "\n",
      "\tGetting Results\n",
      "--------------------------\n",
      "Discarded solution information\n",
      "--------------------------\n",
      "Started Solving MWCE Model 180000\n",
      "Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (win64)\n",
      "Optimize a model with 26615 rows, 18218 columns and 71168 nonzeros\n",
      "Model fingerprint: 0x9b4dec49\n",
      "Variable types: 9794 continuous, 8424 integer (8424 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 2e+05]\n",
      "  Objective range  [1e-02, 3e+06]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e-02, 3e+06]\n",
      "Presolve removed 10674 rows and 6303 columns\n",
      "Presolve time: 0.34s\n",
      "Presolved: 15941 rows, 11915 columns, 45900 nonzeros\n",
      "Variable types: 6029 continuous, 5886 integer (5886 binary)\n",
      "Found heuristic solution: objective 2.074591e+08\n",
      "\n",
      "Root relaxation: objective 2.483412e+08, 6488 iterations, 0.12 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 2.4834e+08    0  936 2.0746e+08 2.4834e+08  19.7%     -    2s\n",
      "     0     0 2.4815e+08    0  544 2.0746e+08 2.4815e+08  19.6%     -   16s\n",
      "     0     0 2.4815e+08    0  583 2.0746e+08 2.4815e+08  19.6%     -   16s\n",
      "     0     0 2.4815e+08    0  584 2.0746e+08 2.4815e+08  19.6%     -   16s\n",
      "     0     0 2.4810e+08    0  415 2.0746e+08 2.4810e+08  19.6%     -   17s\n",
      "     0     0 2.4809e+08    0  475 2.0746e+08 2.4809e+08  19.6%     -   17s\n",
      "     0     0 2.4809e+08    0  467 2.0746e+08 2.4809e+08  19.6%     -   17s\n",
      "     0     0 2.4809e+08    0  484 2.0746e+08 2.4809e+08  19.6%     -   17s\n",
      "     0     0 2.4809e+08    0  484 2.0746e+08 2.4809e+08  19.6%     -   17s\n",
      "     0     0 2.4807e+08    0  252 2.0746e+08 2.4807e+08  19.6%     -   18s\n",
      "     0     0 2.4807e+08    0  257 2.0746e+08 2.4807e+08  19.6%     -   18s\n",
      "     0     0 2.4807e+08    0  267 2.0746e+08 2.4807e+08  19.6%     -   18s\n",
      "     0     0 2.4807e+08    0  189 2.0746e+08 2.4807e+08  19.6%     -   19s\n",
      "     0     0 2.4807e+08    0  201 2.0746e+08 2.4807e+08  19.6%     -   20s\n",
      "     0     0 2.4807e+08    0  207 2.0746e+08 2.4807e+08  19.6%     -   20s\n",
      "     0     0 2.4806e+08    0  225 2.0746e+08 2.4806e+08  19.6%     -   21s\n",
      "     0     0 2.4806e+08    0  224 2.0746e+08 2.4806e+08  19.6%     -   22s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     0     2 2.4806e+08    0  212 2.0746e+08 2.4806e+08  19.6%     -   23s\n",
      "    23    29 2.4806e+08    5  215 2.0746e+08 2.4806e+08  19.6%   235   25s\n",
      "    28    33 2.4687e+08    6  480 2.0746e+08 2.4806e+08  19.6%   445   61s\n",
      "   174   282 2.4805e+08   27  126 2.0746e+08 2.4806e+08  19.6%   181   65s\n",
      "   725   865 2.4805e+08   20  125 2.0746e+08 2.4806e+08  19.6%  68.3   70s\n",
      "  1253  1561 2.4686e+08   54  438 2.0746e+08 2.4806e+08  19.6%  63.5   76s\n",
      "  1620  1697 2.4685e+08  169  307 2.0746e+08 2.4806e+08  19.6%  55.3  172s\n",
      "  1756  1727 2.4680e+08  194  324 2.0746e+08 2.4806e+08  19.6%  53.9  208s\n",
      "H 1759  1727                    2.075687e+08 2.4806e+08  19.5%  53.8  208s\n",
      "  1822  2124 2.4684e+08  194  281 2.0757e+08 2.4806e+08  19.5%  55.3  212s\n",
      "  2284  2125 2.1252e+08  246  224 2.0757e+08 2.4806e+08  19.5%  54.5  231s\n",
      "  2286  2126 2.4700e+08  117  883 2.0757e+08 2.4806e+08  19.5%  54.5  249s\n",
      "  2287  2127 2.4701e+08  131  720 2.0757e+08 2.4806e+08  19.5%  54.5  250s\n",
      "  2298  2134 2.1286e+08  157  256 2.0757e+08 2.4805e+08  19.5%  54.2  255s\n",
      "  2304  2138 2.1241e+08  468  104 2.0757e+08 2.4804e+08  19.5%  54.1  260s\n",
      "  2312  2144 2.1290e+08   78   59 2.0757e+08 2.4804e+08  19.5%  53.9  284s\n",
      "  2313  2144 2.4681e+08  157   59 2.0757e+08 2.4804e+08  19.5%  53.9  285s\n",
      "  2319  2153 2.4804e+08   16   52 2.0757e+08 2.4804e+08  19.5%  61.9  293s\n",
      "  2327  2159 2.4800e+08   17   54 2.0757e+08 2.4804e+08  19.5%  62.1  306s\n",
      "  2347  2172 2.4804e+08   20   38 2.0757e+08 2.4804e+08  19.5%  63.9  330s\n",
      "  2359  2187 2.4696e+08   21 1029 2.0757e+08 2.4804e+08  19.5%  64.9  340s\n",
      "  2386  2231 2.4804e+08   24   46 2.0757e+08 2.4804e+08  19.5%  70.2  370s\n",
      "  2423  2263 2.4804e+08   29   41 2.0757e+08 2.4804e+08  19.5%  75.0  378s\n",
      "  2473  2416 2.4706e+08   36  419 2.0757e+08 2.4804e+08  19.5%  83.0  402s\n",
      "  2650  2361 2.4706e+08   60  403 2.0757e+08 2.4804e+08  19.5%  91.2  443s\n",
      "  2654  2521 2.3885e+08   60  417 2.0757e+08 2.4804e+08  19.5%  91.6  448s\n",
      "  2817  2470 2.4705e+08   89  370 2.0757e+08 2.4804e+08  19.5%  95.2  502s\n",
      "  2821  2637 2.4658e+08   89  373 2.0757e+08 2.4804e+08  19.5%  95.2  506s\n",
      "  2991  2736 2.4588e+08  116  328 2.0757e+08 2.4804e+08  19.5%  96.4  510s\n",
      "  3277  2867 2.4701e+08  163  278 2.0757e+08 2.4804e+08  19.5%   100  531s\n",
      "  3385  2915 2.4656e+08  183  390 2.0757e+08 2.4804e+08  19.5%   105  537s\n",
      "  3479  2964 2.4691e+08  196  242 2.0757e+08 2.4804e+08  19.5%   109  542s\n",
      "  3573  3005 2.4689e+08  206  285 2.0757e+08 2.4804e+08  19.5%   114  547s\n",
      "  3686  3044 2.4689e+08  221  271 2.0757e+08 2.4804e+08  19.5%   120  552s\n",
      "  3774  3291 2.4688e+08  240  239 2.0757e+08 2.4804e+08  19.5%   122  558s\n",
      "  4139  3415 2.4655e+08  291  294 2.0757e+08 2.4804e+08  19.5%   129  564s\n",
      "  4514  3514 2.2976e+08  343  536 2.0757e+08 2.4804e+08  19.5%   134  568s\n",
      "  4854  3668 2.4803e+08   33   92 2.0757e+08 2.4804e+08  19.5%   135  572s\n",
      "  5156  3571 2.4747e+08   97  524 2.0757e+08 2.4804e+08  19.5%   134  606s\n",
      "  5274  3949 2.4644e+08  144  428 2.0757e+08 2.4804e+08  19.5%   136  612s\n",
      "  5684  4210 2.4543e+08  283  575 2.0757e+08 2.4804e+08  19.5%   135  616s\n",
      "  6156  4464 2.4654e+08  240  265 2.0757e+08 2.4804e+08  19.5%   132  622s\n",
      "* 6504  1285              31    2.480417e+08 2.4804e+08  0.00%   127  622s\n",
      "\n",
      "Cutting planes:\n",
      "  Learned: 57\n",
      "  Gomory: 247\n",
      "  Cover: 7\n",
      "  Implied bound: 51\n",
      "  Clique: 49\n",
      "  MIR: 684\n",
      "  Flow cover: 2130\n",
      "  RLT: 19\n",
      "  Relax-and-lift: 644\n",
      "\n",
      "Explored 6735 nodes (872930 simplex iterations) in 622.23 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 3: 2.48042e+08 2.07569e+08 2.07459e+08 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-02)\n",
      "Best objective 2.480416660706e+08, best bound 2.480424216463e+08, gap 0.0003%\n",
      "\n",
      "\tGetting Results\n",
      "--------------------------\n"
     ]
    }
   ],
   "source": [
    "psi_rate = psi_values[2]\n",
    "\n",
    "print(f'Psi Value: {psi_rate}')\n",
    "MCE_obj,high_MCE_leader,high_MCE_wh,high_MCE_wl,high_MCE_vplus,high_MCE_vminus,high_MCE_ukt = with_emission(nr_t,nr_fuel,nr_cust,alpha_rate,fuel_supply,epsilon_rate,BIGM_price,caligraphy_n,beta_rate,\n",
    "                  demand_profile,demand,psi_rate,omega_rate,cost_xi,cost_gamma,BIGM_kkt1,BIGM_kkt2,BIGM_kkt3,\n",
    "                  cost_fuel,cost_ramp,cost_emission)\n",
    "\n",
    "Profits_Actual_MCE['High'] = MCE_obj\n",
    "\n",
    "MWCE_obj,high_MWCE_leader,high_MWCE_wh,high_MWCE_wl,high_MWCE_vplus,high_MWCE_vminus,high_MWCE_ukt = without_emission(nr_t,nr_fuel,nr_cust,alpha_rate,fuel_supply,epsilon_rate,BIGM_price,caligraphy_n,beta_rate,\n",
    "                  demand_profile,demand,psi_rate,omega_rate,cost_xi,cost_gamma,BIGM_kkt1,BIGM_kkt2,BIGM_kkt3,\n",
    "                  cost_fuel,cost_ramp,cost_emission)\n",
    "\n",
    "Profits_Actual_MWCE['High'] = MWCE_obj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "high_MCE_leader.to_csv('high_MCE_leader.csv')\n",
    "high_MCE_wh.to_csv('high_MCE_wh.csv')\n",
    "high_MCE_wl.to_csv('high_MCE_wl.csv')\n",
    "high_MCE_vplus.to_csv('high_MCE_vplus.csv')\n",
    "high_MCE_vminus.to_csv('high_MCE_vminus.csv')\n",
    "high_MCE_ukt.to_csv('high_MCE_ukt.csv')\n",
    "\n",
    "high_MWCE_leader.to_csv('high_MWCE_leader.csv')\n",
    "high_MWCE_wh.to_csv('high_MWCE_wh.csv')\n",
    "high_MWCE_wl.to_csv('high_MWCE_wl.csv')\n",
    "high_MWCE_vplus.to_csv('high_MWCE_vplus.csv')\n",
    "high_MWCE_vminus.to_csv('high_MWCE_vminus.csv')\n",
    "high_MWCE_ukt.to_csv('high_MWCE_ukt.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "({'Low': 248041666.07059878,\n",
       "  'High': 248041666.07059884,\n",
       "  'Medium': 248041666.0705971},\n",
       " {'Low': 153422029.68793482,\n",
       "  'High': 153422029.68793482,\n",
       "  'Medium': 153422029.68793476})"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Profits_Actual_MWCE, Profits_Actual_MCE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "({'Low': 164567350.0, 'Medium': 164567350.0, 'High': 164567360.0},\n",
       " {'Low': 89496170.0, 'Medium': 89496180.0, 'High': 89496160.0})"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Emission_Cost_MWCE, Emission_Cost_MCE"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}