NEUP_MaxProfit_V23.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30 10:11:34.979920 start\n",
      "Running initial pass\n",
      "Building Model\n",
      "Params & Vars established\n",
      "WARNING: DEPRECATED: Using __getitem__ to return a set value from its\n",
      "    (ordered) position is deprecated.  Please use at()  (deprecated in 6.1,\n",
      "    will be removed in (or after) 7.0) (called from <ipython-\n",
      "    input-4-74cc1116e232>:346)\n",
      "WARNING: DEPRECATED: Using __getitem__ to return a set value from its\n",
      "    (ordered) position is deprecated.  Please use at()  (deprecated in 6.1,\n",
      "    will be removed in (or after) 7.0) (called from <ipython-\n",
      "    input-4-74cc1116e232>:353)\n",
      "WARNING: DEPRECATED: Using __getitem__ to return a set value from its\n",
      "    (ordered) position is deprecated.  Please use at()  (deprecated in 6.1,\n",
      "    will be removed in (or after) 7.0) (called from <ipython-\n",
      "    input-4-74cc1116e232>:360)\n",
      "WARNING: DEPRECATED: Using __getitem__ to return a set value from its\n",
      "    (ordered) position is deprecated.  Please use at()  (deprecated in 6.1,\n",
      "    will be removed in (or after) 7.0) (called from <ipython-\n",
      "    input-4-74cc1116e232>:383)\n",
      "WARNING: DEPRECATED: Using __getitem__ to return a set value from its\n",
      "    (ordered) position is deprecated.  Please use at()  (deprecated in 6.1,\n",
      "    will be removed in (or after) 7.0) (called from <ipython-\n",
      "    input-4-74cc1116e232>:393)\n",
      "Model Built\n",
      "Solving...\n",
      "-6267671.557167619\n",
      "iMSR was chosen for operation\n",
      "2.0 modules used\n",
      "____________________________0___________________\n",
      "-732860.8772478672\n",
      "____________________________1___________________\n",
      "-662746.8294017536\n",
      "____________________________2___________________\n",
      "-331729.1227186456\n",
      "2022-08-30 10:12:06.883876 end\n",
      "0:00:31.903956\n",
      "Profit                        -1.050190e+06\n",
      "Annualized Capital Cost        7.641326e+05\n",
      "Fixed Cost                     1.970137e+05\n",
      "Variable Cost                  1.012380e+04\n",
      "Fuel Cost                      2.142601e+05\n",
      "Start-Up Cost                  6.000000e+04\n",
      "TES Annualized Capital         0.000000e+00\n",
      "Electric Wholesale Revenues    1.953402e+05\n",
      "Feasibility                    0.000000e+00\n",
      "dtype: float64\n",
      "2022-08-30 10:12:06.891882 start\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\"\\nA = main(hrcount = int(87),sites = 'Base4_MSR.csv',lmps = 'SPP_SOUTH_2021.csv',NGCost = 8, MaxMod=2, TES = 200, naics = 324110, EP=10, startHr = 0, dem = 327)\\ne = datetime.now()\\nprint(e,'end')\\nd = e-s\\nprint(d)\\nQuickGraphOutputs(A.Output,L,'DAM_test2',hr = 72, TES = True)\\nWaterfallPlot(A,'TestFull', 'SPP 2021, 327 MWht, 300 C, Full year')\\n\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAosAAACcCAYAAAAXpuEdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXxMV//A8c9kIxshUrEkrcjWUkUQsYTYpYLyNK0wiqpWn1JRKmksUxpiKbVT+pSmS2r79WltXdTylBJLaEsUSYVUSRBLQibb/f0xMrJNMols+L5fr7zM3HPPud+bXMl3zrn3HJWiKApCCCGEEEIUwaSqAxBCCCGEENWXJItCCCGEEMIgSRaFEEIIIYRBkiwKIYQQQgiDJFkUQgghhBAGSbIohBBCCCEMkmRRCCGEEEIYJMmiEEIIIYQwyMxQwaefflrmRkeOHFnmukIIIYQQovpQGVrBxdPTs2wNqlTExsY+UFBCCCGEEKJ6MNizCLBhwwZatGhhdGPHjx9nyJAhDxyUEEIIIYSoHgzes9i9e3dq165dqsbs7Ozo1q3bAwclhBBCCCGqB4PD0HmFhoby5ptv4uTkVKgsPj6e+fPns3LlygoJUAghhBBCVB2Dw9CnTp0iN4/85ptv8Pb25tatW4X22717NwcOHKi4CIUQQgghRJUxmCz+5z//YevWrahUKkDXu1hQbjI5cODACgpPCCGEEEJUJYPD0Ldv3yY2NhZFUXjllVeYPn06rq6u+fYxMTGhVq1auLm56ZNKIYQQQgjx6DCYLMbHx+Pi4gJAdHQ0zZo1w9raulKDE0IIIYQQVcvgMLS/vz9PPPEEPj4++Pj44OzsLMmiEEIIIcRjxmDP4vbt2zl27BjHjh3jzz//JCcnhyeffFKfPHp7e5d6ah0hhBBCCPFwMWrqnDt37nD8+HGOHTvG0aNHOXHiBOnp6Xh6euLj40P79u3p3LlzZcQrhBBCCCEqkVHJYkE5OTn88ssvbNy4kd27d5OdnS1L/AkhhBBCPIKKXe4vr+vXr3PgwAF+/fVXjhw5woULF7C2tqZz58506NChImMUQgghhBBVxGDPYkZGBocPH2b//v0cOHCAM2fOYGJiQsuWLfHx8aFDhw60aNECU1PTyo65SMePH6dGjRrl1p5Wqy3X9sSjSa4TURK5RkRJHuQa0Wq1tGzZUv8+MzOTxMRE0tPTyys88RioWbMmjRs3xtzcvMhyg8nic889R2ZmJh4eHnh7e+Pj40O7du2wtLSs0IDLKjY2lqeffrraticeTXKdiJLINSJK8iDXSMG6f/31F7a2ttjb28v8x8IoiqJw7do1bt++TZMmTYrcx+AwtFarpVatWjg5OeHk5ISzs/MDJ4rZ2dlMnTqVv/76C1NTU+bMmYOiKISEhKBSqXBzc2PGjBmYmJiwYcMGoqKiMDMzY+zYsfj5+T3QsYUQQohHXXp6Ok899ZQkisJoKpUKe3t7kpOTDe5jMFnMHX4+cOAAa9as4YMPPqB+/fp06NCBDh064OPjg729fakC2r17NwBRUVEcOnRInyxOmDABb29vpk+fzq5du2jZsiWRkZFs3rwZrVZLUFAQHTt2xMLColTHE0IIIR43kiiK0irpmjGYLNrb2xMQEEBAQAAAcXFx+gdcNBoNaWlpuLu76xNHX1/fEoPp0aMHXbt2BeDSpUvUq1ePPXv20K5dOwB8fX3Zv38/JiYmtGrVCgsLCywsLHB2dub06dO0aNHCYNtarbZcn8hOT08v1J5l/G/U/vU7KP0D5OIR5ZCdTUo1uW9XVE9yjYhc17sPJaNB4WG+ov7eCFGdGP00dNOmTWnatClqtZrs7GxiYmKIiooiMjKSdevWGX2hm5mZMWXKFH788UeWLFnC7t279RmttbU1t2/fJjU1FVtbW30da2trUlNTi223Ro0aFX7PovKfUPj9V6jjUG7HEQ+3zMxMgzcECwFyjYj77F5+C1URf6ce9J7F6iQiIoKTJ0+SnJxMeno6Tk5O1KlTh6FDhxIVFcWiRYsqLZZu3bqxY8cOgw8PZWVlsWrVKvbu3avfJyAggJdeeqnMx9yyZQtLlizhlVdeITU1lRUrVrBnzx7q168PwLVr1/D19WXWrFkcPnyYzp074+/vD0Dfvn3x8fFh+vTpAEyZMoWePXuyfv16NBoNTZs2zXeszz//nGHDhhEWFsb333/P119/XWifXD///DPLly/HzMyMwYMHExgYWKrzMjpZ/Oeff4iJieH48eMcP36c06dPk5OTw7PPPou3t3epDjp37lwmTZpEYGAgWq1Wvz0tLY1atWphY2NDWlpavu15k8eqoNxNhd/2w6A3UL0+q0pjEdVHnDy8IEog14h4nISEhAC6pCk+Pp5JkyYBcOjQoaoMq0iLFi0iJyeHqKgoTE1NSUtL4/XXX6dNmzYGky5j9OvXj5EjR7J06VKeeuopduzYwYgRIwDd6ngNGjQAoFOnThw9ehR/f38uXryIs7Mz0dHR+nZiYmKYNm0a69evL/I4K1euZNiwYYSHh3PhwgWD8WRmZjJnzhw2bdqEpaUlQ4YMwc/PDwcH4zu+DCaLuYlh7r/JycmYmprSrFkz2rdvz/jx4/Hy8irVQy/ffPMNV65c4fXXX8fS0hKVSkXz5s05dOgQ3t7e7Nu3j/bt29OiRQs++ugjtFotGRkZxMXF4e7ubvRxKsSxvZCVCe16VG0cQgghhJGUic+Xa3uqhdvKXDchIYHRo0dz/fp1/Pz8GDduHH/++ScffPABAHZ2dsyePZtTp06xYMECzM3NCQwM5JNPPqFNmzacOXOGJk2aYG9vz5EjR7CwsODjjz/m2rVraDQatFotN27c4N///jc9ehT/tzorK4sdO3bwww8/6KcAtLa2JjIyUj/a+eGHH3L48GEURWHEiBH07dsXtVqNp6cnZ8+eJTU1lcWLF9OoUSODx/H392fnzp36ZHH37t36B3Z9fHxYu3YtAHv27KFbt278/PPPnDt3jho1alC/fn1sbGwAWL58OVevXuXu3bssXLiQrVu3cvPmTTQaDRqNpthzjYuLw9nZWb9Es5eXF0eOHKFv377F1svLYLI4ZMgQTE1NefrppwkICMDb25s2bdpgZWVldOMF9erVi9DQUIYOHUpWVhbvvfceTZs2Zdq0aSxcuBAXFxd69+6NqakparWaoKAgFEUhODi46ucpO/QjWNlCc5+qjUMIIYR4CGm1WlasWEF2djZdu3Zl3LhxTJs2jdmzZ+Pq6srGjRtZu3YtHTp0QKvVsnHjRgCWLFlCv3798PLyok+fPoSGhhIcHMywYcM4d+4cKSkpjBw5Em9vb44dO8bSpUtLTBZTUlKoXbs2Zma6NOjLL79kx44dpKWl0b9/f5o0aUJiYiJRUVFotVoCAwPp2LEjAC1atCAsLIxFixaxbds2xowZY/A49erVw9LSkosXL5KTk4Ojo6M+n6lbty4qlYrbt2+zb98+Zs6cSVZWFvv27aN27dr5llHu0qULAwYMYOnSpezcuZOxY8fy+eefl5goAmW6ta8gg8niqlWraNOmjT6rLQ9WVlYsXry40PbPP/+80LbAwMBSj6lXFEVRIPpH8OqKylyeyBZCCPFweJCewPLm5uamn9UkN0mLi4vj/fffB3TDpbnz/BWc769Zs2YA1KpVSz9EXKtWLbRaLQ4ODqxcuZJNmzahUqnIysoqMRY7Oztu3LhBdnY2pqamBAUFERQUxFdffcXVq1c5c+YMJ0+eRK1WA7qeyEuXLgHwzDPPAODo6MjVq1dLPNbzzz/Ptm3byMrKIiAggP379+vLfHx8OHDgACkpKTRo0ABfX1/mzZuHtbW1vjcSoHnz5oAu+TTmmHmVx619BpNFBwcHEhISjG4o9wf5SPrrFCT/DcOnVHUkQgghxEOpqOlZmjRpwty5c2nYsCFHjx7Vz/VnYmJSYt1cixcv5sUXX6RLly5s3ryZ//u//ysxFnNzc3r16sVHH31EcHAwJiYmaLVaTpw4QePGjXFxccHb25tZs2aRk5PDihUraNy4cSnPWKd3796MGjUKa2tr3nzzzXzJYseOHVm2bJl+VhgnJydu3LjB1atX8fT0LLZdA2uqFNK0aVMSEhK4ceMGVlZWHDlyhFdffbVU52AwWRw8eHC+H05uUAV/YIqioFKpqt0TWeXq0I+6f9v1LFR0AYgB7ubZNhjIffZxE5D3M44H0Oreay26H4BMqvFoUoAbwEUg8d42/3v/3gR2FNi/A+B87/VBIBNoDDQCpD/70XUb3e+Cevfe/wUUfBSgH5A7xvMdkAaoAAfADXCq+DCFqBAajYYpU6aQnZ0NQHh4OElJSaVqo0+fPoSHh7N69WoaNGhASkqKUfUmT57M2rVrGTp0KGZmZqSmptKjRw9GjhyJpaUl0dHRBAUFcefOHXr06FHmkVZbW1scHR1xcnIqlAR7eXlx8uRJ3n77bf02T09Po4aJmzZtyqRJk1iwYEGx+5mbmxMSEsKrr76KoigMHjxY/3S2sQwu99eqVSvu3r2Lp6cnffv2pWXLlsVm9rlZcVWpyOX+lGB/SLvFjx//Qji6X9yT7+33ITCpQN1bQG4Hry2Q90c+8V4dgIXAu0BDdEmBE/ApkHtX6BHAEWiAJJTVUSawLT4eCxcXfUL4AtD6XvkQICrP/p5A7keqWOCZAu1FAbkTNrwEbMhT1h749d7rZHTXSe410/jel0zOUv3cBvbGxVGjaVM6ANb3tqvRfchMRPfBYTiQ+7zjZ8ArBdqJA1zuvXZBl1DmUt+rA7ATWEb+6+JfeY4rqqfyXO5PlpesWgWfAq8sarW6yOl1SqO4a8dgz+LBgwfZu3cvO3fuZPXq1dja2tK7d2/69u1Lq1atDFV75CipN+CPQ/DyBNai+wWf907Kl4Eu3P/UD/eTPYCjQE6e93XyvG4PhHK/5+kUkPtseQbQDl3vlCm6hDEMeONe+UHuJw+g63n6973XWmBFgfPwBbzuvd4LHJO6Rtc9ClxH9zPKAL68V3YFeMHFRV/PBGjC/WTxJaAt9/9w5/0c58L9xDFXwzyv5wGv3jvmRfJfU2eAgjdEbAVyn3mcA1xC1yNZ497x/3Wv7AKwuUDdwHv7Anx9r24uqWt83SeBQfde70T3u+EmwL1f3kfRjSqo0P2fdgP80F0bbfK0M4DC10bewa8f0X1QyUF3DdrlKUsD/kbXM5l7V1N/7ieLHkBN7n/Q6HfvCyAauD84pjOO+38klgDZecraAp2kbqnrisq1a9cu1q1bV2j78OHD6dmz8Ghhedi6dSv29vaMHDmyQtovKCwsjNjYWPbv31/kAy/lcq6KEdLT05Xvv/9emThxotKqVSvF19dXCQ8PV44ePWpM9Upx6tSpCmkvZ8//KTndais5v/+qtFQUxb9cj2JYhqIo2xRFWaUoylRFUV5RFOW7POX/VhSFPF+2ecpuFihDUZQP85QHS91S1zVRFKWxoii+iqLk3CvLUhRl0cWLyq+KoiQqipKpVJ6biqKcUhTlB0VRPlEU5VKeskBFUWor98+nW56yn5TC57svT3mnAmVS1/i6nfKUxSqK8paiKHMVRZmfmKjsUxTltnL/2qkMdxRFOZvnmDmKorypKEqAoigtFUWppyjKjDz7z1YKn296nnKLAmUhUrdMdYvyIH+/CtYt77+F4vFR3LVjcBjakIyMDPbt28f333/P3r17sbS0pHfv3rz33nsPlrU+oIoahlbmvQkHtqNsPkdtUzNGAYWf5658d9D1LuRV+96/Crqh8LxqoutlAt39lRlSt1R1rSm6G746D/mkobtf1pT7Pd9Z97bnlffcUsnfIyJ1ja9rSdH3llbna0RB19MJkI6upz2vWnnKbxYoq4Hu/4rULV3dosgwtKgOyjQMbYiFhQXdu3fHzs4OW1tbNmzYQGRkZJUnixVBycmB6J/AqxsmpmZconDSUVWKm+1Sxf1kqCiW3B/ulrplr1vdFXWfmhnFn29xt29LXePrPizy3oVek+ITmuK+F1LX+LpCPIyMThazs7P59ddf+eGHH9i1axcpKSk8++yzBAcH06tXr4qMseqc+w1SkqC97vyqdsFBIYQQQojKV2yymDvk/MMPP7Bnzx5SU1Np1aoVr7/+Or169cLR0bGy4qwah34AlQradCcbeSJZCCGEEI8fE0MFb7/9Nt7e3kyYMIHk5GSCg4PZt28fX3zxBcOHD3/0E0XQza/o0RpVHQemo3siOaekOkIIIcRj6tChQ/j4+KBWq/Vf48ePB3TTu8TFxZXLcS5evMj48eMJDAxk+PDhjBkzhrNnzz5Qm2q1mn/961+cO3eObt26MXr06Hzln376KR4eHmRkZNCxY0f93JAxMTF4enry+++/A7plDTt37szFixeLXInu0qVL/PzzzyQkJJS4Wl1OTg7Tp0/npZdeQq1Wl2qxlPJksGfx+++/x8TEhBYtWlCjRg327t3L3r17i9xXpVKxcuXKCguyKpim3oDTR0D9LgDn0N2jZDC7FkIIIQTt27dn0aJFFdb+3bt3GTt2LLNmzdJP5ffbb78xc+ZMIiMjH6jtuXPn6ucqvHLlCtevX6du3boA7N27l9q1a2NhYYGnpyexsbE0b96cvXv30qdPH/bu3cuzzz5LTEwMbdu2NTg39cGDB4mPj6dbt24sXLiQiRMnGoznp59+IiMjg6+//prjx48TERFRJfmWwWSxbdu2+td51xR8XFjHHgJFAe/egC5ZdK3akIQQQohS6VrO7e0phzZu375NWFiYfqWVqVOn4uHhwcaNG/nqq6/Iycmhe/fujBs3rsj6u3fvpn379vnmfG7RogWffaabnv6ff/5h2rRpaLVaatSowaxZs8jOzuadd97B0dGRixcv8uyzz+rXpDakd+/e7Ny5k6CgIOLi4nB2dtb3Xnbs2JEjR47QvHlzDh48yPz583nnnXd46623iI6OpnPnzgBcv36dN998k+TkZDw8PHj//ff5+OOPSU9Pp1WrVnh4eBQbw9GjR/VttWzZkj/++MOI73D5M5gsPmh2/rCzPnkA7BzAvSUKcBboWNVBCSGEENXcwYMHUavV+vddunTJN6S7atUq2rdvT1BQEOfPnyc0NJRly5axZs0avv32WywsLIiIiCAtLQ1r68LzOiQmJuLs7Kx/P3bsWFJTU0lKSmL9+vXMnTsXtVpNly5d+PXXX1mwYAHBwcGcP3+eTz75BEtLS3r06EFycjIODg4Gz6Nfv35MmzaNoKAgvv32WwICAti1axcAHTp0YPny5fTr1w9LS0ucnJxQFIXr169z+PBhhgwZglarJTU1lTlz5mBra0vPnj25ceMGY8aMIT4+nu7du5OYmGjw+ACpqan5lhk0NTUlKysLM7NST2bzQAwerWvXrqxatarEhazzio2NZezYsezZs6c8YqsySnY2NqcOQqfnUZmYkIluxYyqXdBQCCGEKJ09VXDMkoahz5w5w8GDB9mxYwcAt27d4uLFi7i5uVGzpm5SouKm43N0dMzXw5Y7LBsYGEhWVhZnzpxh9erVrF27FkVRMDfXLYbq7OysT7wcHBzQagvOlplfgwYNAF1P5bFjx5gwYYK+zMPDg/Pnz/O///1P3/PXqVMnDh06REZGBg4ODiQmJuLk5ETt2rrJlOzt7bl7926xxyzIxsYm3+huTk5OpSeKUEyyePnyZU6fPq2/gdMYZ86c4cqVK+USWJW6cgHF1BS8dVPmmKNbak8IIYQQD8bFxYX+/fsTEBDAtWvX2LhxI87OzsTHx5ORkYGFhQXjx48nLCyM+vXrF6rfvXt31qxZw/Hjx2nZsiUACQkJXL58GZVKhYuLC6NGjaJ169bExcVx+PBhAIP3EBbH39+fiIgIWrVqla++SqXSD53PnDkTAF9fX+bNm0e7du3y7VeQiYkJOTnGPS7bunVrdu/ejb+/P8ePH8fd3b3U51Aeik1PQ0NDS9WYoihl+mFUN6qGTTgXvhVPz+LvJRBCCCFEfgWHoQHWrFmjf/3GG28QFhbGhg0bSE1N5a233qJu3bq89tprDBs2DJVKhZ+fX5GJIoC1tTUrV67kww8/ZMGCBfph2VmzZtGoUSOmTJmCRqNBq9WSnp5OWFjZu3v69OlDeHg433zzTaGyjh07snTpUlxddU80tGjRgvj4eIKDg4tt093dnZUrV9KsWTOee+65Yvft2bMn+/fv5+WXX0ZRFGbPnl3mc3kQBpf7i46OLnOjebPqylJRy/0BbAbOAKVLncXjQJbWEiWRa0SURJb7qx7UajUajUb/NHRlSExMZOLEiWzYsKHSjmlImZb7q4qEr7raABxDkkUhhBCismg0miLnZVyzZo3+3sbyNmXKFCIiIvS9hRUpISGByZMnA1VzrqVR+XdJPoRk2hwhhBCicmk0mko9XmXPAvPkk09Wix5FY8gc0yVQkGRRCCGEEI8vSRZLcBOwANyqOhAhhBBCiCpg1DD0pUuXaNiwYUXHUi3ZAcnImtBCCCGEeDwZ1bM4dOhQvv3224qOpVqTLlghhBBCPI6M6lnMzMykVq1aFR1LtXQZsAIez7MXQgghjBcREcHJkydJTk4mPT0dJycn6tSpwy+//EKzZs3y7btu3ToyMzPRaDQkJSWhUqmwsbFBo9FQp04dg8c4e/Ys8+fP5+7du9y5c4cuXbowbtw4oqOjiYqKKnb1mIr0888/s3z5cszMzBg8eDCBgYFVEkdFMCpZHDduHBqNhhEjRuDi4oK9vX2hfQpeBI+KYCAaKPxAuxBCCCHyCgkJAWDLli3Ex8czadIkEhMTuXz5cpFPG0dFRVGvXj0iIiIAXQK5fPlypk6dWmT7t27dYuLEiSxdupSnnnqK7Oxs3n77baKionBxcam4EytBZmYmc+bMYdOmTVhaWjJkyBD8/PyKXXv6YWJUsjhjxgwA/Q8z7yotuau2xMbGlthOZmYm7733Hn///TcZGRmMHTsWV1dXQkJCUKlUuLm5MWPGDExMTNiwYQNRUVGYmZkxduxY/Pz8ynJ+D+ws8nCLEEKIh9Tx0+XbXkvPcm2uUaNGbNq0idatW9OuXTvUajUG1goBYNeuXXh7e/PUU08BYGpqyty5czE3NycmJoaEhARGjx7N9evX8fPz0/c4Llu2DID09HT9/mPHjsXOzg5fX1+8vb15//33sba2xt7enho1ahAREUFkZCRbt25FpVLh7+/P8OHDi4wrLi4OZ2dn/TrQXl5eHDlyhL59+5br96uqGJUsfvbZZ+VysG+//RY7Ozvmz59PSkoKL7zwAp6enkyYMAFvb2+mT5/Orl27aNmyJZGRkWzevBmtVktQUBAdO3bEwsKiXOIwVu60Oe0r9ahCCCHEo+XcuXP5lgBs1qwZISEhdO3alYyMDDZt2kRoaCju7u5MnToVD4+il9tNSkrCyckp3zZra2v9a61Wy4oVK8jOzqZr166MGzdOP2xdv359Vq1axc6dOwkICCA5OZnNmzdjYWHBCy+8wLx583Bzc2PRokVcuXKFc+fOsX37dr788ktUKhUjRoygU6dORfZgpqamYmtrmy+m1NTUB/22VRtGJYt5V3PJysoiJSWFOnXqYGZWujm9+/TpQ+/evfXvTU1NOXnypL59X19f9u/fj4mJCa1atcLCwgILCwucnZ05ffo0LVq0MNi2Vqs1qnfTWOnp6fx65gw33d2xvXyZ2JSUcmtbPDrS09PL9boTjx65RkRJKvQaKeeewLJydXUtchg6JiYGHx8fevXqRXZ2Nv/9738JDQ1ly5YtRbbTsGFDTp06lW/bxYsXuXz5MgBubm76jqXcHKV+/fqEh4djZWXFlStXaN26NQCNGzfW75uUlISbm24c0cvLi+3bt3PmzBkuXbrEiBEjALh58yYXLlwoMlm0sbEhLS1N/z4tLS1f8viwMzrbO378OB999BFHjx4lOzubjRs38umnn9KoUaMSF83OlZv9p6amMn78eCZMmMDcuXP1w9rW1tbcvn27TBl6jRo1yn1taDd3d74A2jg64u7oWG5ti0eHrMMqSiLXiCjJg64N/TDbtm0b1tbWBAcHY2pqioeHR7GjiH5+fqxevZohQ4bg7OxMZmYmERERdOjQAVdX13y3yeWaOnUqP/30EzY2NkyZMkU/zG1icn+eE0dHR86dO4erqysnTpwAwMXFBVdXV9auXYtKpWLdunW4u7sXGVfTpk1JSEjgxo0bWFlZceTIEV599dUH+dZUK0Yli7/++iuvvfYabdq0ITg4mPnz5wPg7u7ORx99hJ2dHSNHjjTqgP/88w///ve/CQoKIiAgQN8W6DLxWrVqVZsM3QYIqvSjCiGEEI+WgsPQALNnz2bChAnMmjWLAQMGYGlpiZWVFeHh4QbbsbGxISIigqlTp6IoCmlpafj5+REUFER0dHSRdQYMGEBgYCC1atWiXr16JCUlFdpnxowZvPfee1hZWWFubk79+vXx9PTEx8eHIUOGkJGRQYsWLahfv36RxzA3NyckJIRXX30VRVEYPHiwwX0fSooRBg0apEyePFlRFEXJzMxUPDw8lD/++ENRFEVZvHix0rt3b2OaUZKTk5U+ffooBw4c0G97/fXXlYMHDyqKoijTpk1Ttm3bpiQlJSn9+vVT0tPTlVu3bim9e/dW0tPTi2371KlTRsVgrPJuTzya5DoRJZFrRJTkQa6RgnXleiubzz//XLl27ZqiKIqycOFCZenSpVUcUeUr7toxqmfx7Nmz+qHmgl283t7erF271qjEdNWqVdy6dYsVK1awYsUKAMLCwvjggw9YuHAhLi4u9O7dG1NTU9RqNUFBQSiKQnBwMDVq1ChNDlwuZgJ3gIhKP7IQQgjx+NJoNMTFFZ60bs2aNdSsWbPcj2dvb8+oUaOwsrLC1tZWP/tLQbt27WLdunWFtg8fPpyePXuWe1zVhVHJor29PefOnaNTp06FyuLi4oqcd7EoU6dOLXLupM8//7zQtsDAwCqf0PI7wPC0oEIIIYSoCBqNplKP16dPH/r06VPift27d6d79+6VEFH1YlSyOHDgQBYvXoytrS2+vr4AZGdnc+DAAZYtW8bgwYMrNMiqoKCbY3FoVQcihBBCCDhq1HgAACAASURBVFGFjEoW33rrLS5fvkxYWJh+GPrll19GURR69uzJ+PHjKzTIqnDD1JSbgGtVByKEEEIIUYWMShZNTU2ZM2cOo0ePJjo6mps3b2Jra4uXlxeentVjDqfydtvEhHZA86oORAghhBCiChmVLAYHB9OlSxc6d+7MkCFDKjqmasE5M5NDVR2EEEIIIUQVMyl5F7hz5w4zZ86kc+fODB48mCVLlnDixIli128UQgghxOMlIiICtVpNnz596Nq1K2q1mvHjx3Po0CGjF/AoL926dUOr1Rosz8rKYtmyZbz44osMGzaMYcOG8fXXXz/QMbds2ULXrl359NNPWbp0KU8//TRXrlzRl1+7do1mzZqxZcsWQkND2b59u76sb9++zJw5U/9+ypQp/PTTT6jV6iKfDM99ODgsLIw2bdoUuU9ed+/e5eWXXy5xv6IYlSyuXr2a6OhoPv30Uzp27MjevXsZMmQIPj4+TJo0ia1bt5b6wNXdiZo1OVXybkIIIYS4JyQkhMjISMaMGUO/fv2IjIxkyZIlVR1WkRYtWkRaWhpRUVF8/vnnrF69mu+++65MyVRe/fr10y9U8tRTT7Fjxw592fbt22nQoAEAnTp14ujRo4BuyUJnZ+d8E4vHxMTQvn17g8dZuXIlAOHh4SWuAPT7778zdOhQLl68WKZzMnq5PzMzM7y9vfH29mbixIkcOHCAJUuWsHXrVrZt20a/fv3KFEB1NdvREUfgh6oORAghhCgjZeLz+Tf0CkLVRzfPh7LzC/jhy3zFqoXbjKpbFgkJCYwePZrr16/j5+fHuHHj+PPPP/nggw8AsLOzY/bs2Zw6dYoFCxZgbm5OYGAgn3zyCW3atOHMmTM0adIEe3t7jhw5goWFBR9//DHXrl1Do9Gg1Wq5ceMG//73v+nRo0exsWRlZbFjxw5++OEHTE1NAd3SwpGRkfoHeT/88EMOHz6MoiiMGDGCvn37olar8fT05OzZs6SmprJ48WIaNWpk8Dj+/v7s3LlTv7707t278fPzA8DHx0c/T/WePXvo1q0bP//8M+fOnaNGjRrUr18fGxsbAJYvX87Vq1e5e/cuCxcuZOvWrdy8eRONRmPUNEMZGRksX76cd999t8R9i2JUspiRkcFvv/3G0aNHOXr0KDExMaSmpuLk5MSgQYNo165dmQ5enSVYWFB4VkkhhBBClIVWq2XFihVkZ2fTtWtXxo0bx7Rp05g9ezaurq5s3LiRtWvX0qFDB7RaLRs3bgRgyZIl9OvXDy8vL/r06UNoaCjBwcEMGzaMc+fOkZKSwsiRI/H29ubYsWMsXbq0xGQxJSWF2rVrY2amS4O+/PJLduzYQVpaGv3796dJkyYkJiYSFRWFVqslMDCQjh07AtCiRQvCwsJYtGgR27ZtY8yYMQaPU69ePSwtLbl48SI5OTk4OjrqFxmpW7cuKpWK27dvs2/fPmbOnElWVhb79u2jdu3adO7cWd9Oly5dGDBgAEuXLmXnzp2MHTuWzz//3Oj5KL28vIzazxCjkkUvLy+ysrJwdXWlXbt2DBgwgLZt2/LEE0880MGrq+vALVNTmTZHCCHEQy1vT2Ghsj5DoZiewuLqloWbmxsWFhYA+iQtLi6O999/H4DMzEyaNGkCoP83V7NmzQCoVasWTZs21b/WarU4ODiwcuVKNm3ahEqlIisrq8RY7OzsuHHjBtnZ2ZiamhIUFERQUBBfffUVV69e5cyZM5w8eVK/nnVWVhaXLl0C4JlnngHA0dGRq1evlnis559/nm3btpGVlUVAQAD79+/Xl/n4+HDgwAFSUlJo0KABvr6+zJs3D2tra31vJEDz5rq5WerVq2fUMcubUfcsdu7cmVq1apGQkMCff/7JmTNnOHv2LHfu3Kno+KrE2Xv/ulVpFEIIIcSjo+BywaBLCufOnUtkZCSTJ0+mS5cuAJiYmJRYN9fixYsZMGAA8+fPx9vb26iHb83NzenVqxcfffQROTk5gK7n88SJE6hUKlxcXPD29iYyMpL169fTt29fGjduXJrT1evduze7du3iyJEjeHt75yvr2LEj69ev14/QOjk5cePGDRISEkqcmrAyHzI2qmcxdx3n06dPc+jQIaKjo/n6669JTU3F09OTtm3bMmXKlAoNtDI9C3z11190KvDJRgghhBDlR6PRMGXKFLKzswHdwxpJSUmlaqNPnz6Eh4ezevVqGjRoQEpKilH1Jk+ezNq1axk6dChmZmakpqbSo0cPRo4ciaWlJdHR0QQFBXHnzh169Oihv3+wtGxtbXF0dMTJyalQEuzl5cXJkyd5++239ds8PT1JTU0tsd2mTZsyadIkFixYUKa4SkOllCE1TU9P5/Dhw3z11Vf8/PPPqFQqYmNjKyI+o8XGxpb4NFBVticeTXKdiJLINSJK8iDXSMG6cr1VrS1bthAfH8+kSZMq9bhqtRqNRqMfoi+L4q4dox9wiYmJ4dChQxw8eJDff/8dlUpF69atmTx5sn69aCGEEEKI6mLXrl2sW7eu0Pbhw4fTs2fPCjnm1q1bsbe310+fU9HCwsKIjY1l//79RT7wUh7nalTP4nPPPUdGRgYNGzakc+fO+Pr60r59e6ysrB7o4OVJehZFVZDrRJRErhFREulZFNXBA/csTpw4kc6dO+Pi4lKugQkhhBBCiOrNqGTxlVdeIT09nS+++ILo6GhSU1Oxs7PDy8uLgQMHVqseRiGEEEIIUX6Mmjrn2rVrDBo0iPDwcC5cuIBKpSIuLo4PPviAF154gWvXrlV0nEIIIYQQogoY1bM4d+5c0tPT+e677/I9aRMXF8eYMWNYsGABc+bMqbAghRBCCCFE1TCqZ3Hv3r1MmDCh0CPZTZs25e2332bPnj0VEZsQQgghHiKHDh3Cx8cHtVqt/xo/fjygm94lLi6uXI5z8eJFxo8fT2BgIMOHD2fMmDGcPXu25IrFUKvV/Otf/+LcuXN069aN0aNH5yv/9NNP8fDwICMjg44dO+rnhoyJicHT05Pff/8d0E3u3blzZy5evEhgYGCh41y6dImff/6ZhIQEAgMDi9ynoBMnTuhXk6kKRvUsgm5ZHUPb7969W24BCSGEEOLh1b59exYtWlRh7d+9e5exY8cya9YsWrVqBcBvv/3GzJkziYyMfKC2586dq+8Yu3LlCtevX6du3bqAruOsdu3aWFhY4OnpSWxsLM2bN2fv3r306dOHvXv38uyzzxITE0Pbtm0Nrjpz8OBB4uPj6datGwsXLmTixInFxrRmzRq+/fZbLC0tH+jcHoRRPYvPPvssn332mT6LzpWVlcW6det49tlnKyQ4IYQQQpRd1wJf6/KUrSui3Ni6D+L27duMHz9e3/P4559/ArBx40YGDRrEwIEDWbp0qcH6u3fvpn379vpEEaBFixZ89tlnAPzzzz+MHj0atVrN6NGj+eeff0hMTOSll17i7bffZtCgQcyYMaPEOHv37s3OnTsB3W13zs7OmJubA7pl+o4cOQLokr933nmHffv2ARAdHU3nzp0BuH79Om+++SYvvvgiU6dOJTs7m48//pitW7eya9cuo75fzs7OxX4/KoNRPYuTJk3i5ZdfpmfPnnTv3l2/kPWuXbu4evUq69evr+g4hRBCCPEQOHjwYL4h0y5duuQb0l21ahXt27cnKCiI8+fPExoayrJly/Q9aBYWFkRERJCWloa1tXWh9hMTE3F2dta/Hzt2LKmpqSQlJbF+/Xrmzp2LWq2mS5cu/PrrryxYsIDg4GDOnz/PJ598gqWlJT169CA5ORkHBweD59GvXz+mTZtGUFAQ3377LQEBAfoEr0OHDixfvpx+/fphaWmJk5MTiqJw/fp1Dh8+zJAhQ9BqtaSmpjJnzhxsbW3p2bMnN27cYMyYMcTHx9O9e3cSExNL/H727t3bqP0qklHJoqenJ1999RUrVqxg69at3Lp1i9q1a+Pl5cUbb7xBs2bNKjpOIYQQQpTSnmLKRtz7Kkvd4pQ0DH3mzBkOHjzIjh07ALh16xYXL17Ezc2NmjVrAvDee+8ZrO/o6Mgff/yhf79y5UoAAgMDycrK4syZM6xevZq1a9eiKIq+N9DZ2Vm/vrODgwNarbbY82jQoAGg66k8duwYEyZM0Jd5eHhw/vx5/ve//+l7ETt16sShQ4fIyMjAwcGBxMREnJycqF27NgD29vYP7W17Rt+z+PTTT1d5N6gQQgghHm4uLi7079+fgIAArl27xsaNG3F2diY+Pp6MjAwsLCwYP348YWFh1K9fv1D97t27s2bNGo4fP07Lli0BSEhI4PLly6hUKlxcXBg1ahStW7cmLi6Ow4cPAxi8h7A4/v7+RERE0KpVq3z1VSoVHh4ebNy4kZkzZwLg6+vLvHnzaNeuXb79CjIxMSEnJ6fUsVSlYpPF06dPExUVxaVLl3BycuLFF1/E09OzsmITQgghxEOm4DA06B7SyPXGG28QFhbGhg0bSE1N5a233qJu3bq89tprDBs2DJVKhZ+fX5GJIoC1tTUrV67kww8/ZMGCBWRlZWFmZsasWbNo1KgRU6ZMQaPRoNVqSU9PJywsrMzn0qdPH8LDw/nmm28KlXXs2JGlS5fi6uoK6O6bjI+PJzg4uNg23d3dWblyJc2aNeO5554rc2yVyeDa0L/++iuvvfYaOTk51KlTh5SUFFQqFeHh4QwcOLDMBzxx4gQLFiwgMjKShIQEQkJCUKlUuLm5MWPGDExMTNiwYQNRUVGYmZkxduxY/Pz8SmxX1oYWVUGuE1ESuUZESWRt6OpBrVaj0WgKTRNYkRITE5k4cSIbNmyotGMaUqa1oVesWIGnpyfLli3D0dGRmzdvMnnyZD788MMyJ4sFH/+eM2cOEyZMwNvbm+nTp7Nr1y5atmxJZGQkmzdvRqvVEhQURMeOHbGwsCjTMYUQQgjx8NFoNEXOy7hmzRr9vY3lbcqUKUREROh7CytSQkICkydPBqrmXEvDYLJ4+vRp5s2bh6OjIwC1a9dm0qRJDBgwgEuXLtGwYcNSHyz38e93330XgJMnT+rH9n19fdm/fz8mJia0atUKCwsLLCwscHZ25vTp07Ro0aIs5yeEEEKIh5BGo6nU4z3oHI2l9eSTT1aLHkVjGEwW09LS9BNR5nryySdRFIVbt26VKVks+Pi3oij6mz+tra25ffs2qamp2Nra6vextrYmNTW1xLa1Wi2xsbGljsmQ9PT0cm1PPJrkOhElkWtElESuEVHdGUwWc3JyCj3FY2am273g5NxlZWJyf07wtLQ0atWqhY2NDWlpafm2500eDalRo4bcsygqnVwnoiRyjYiSPOg9i0JUNKNWcKkozzzzDIcOHQJg3759tGnThhYtWnD06FG0Wi23b98mLi4Od3f3qgxTCCGEEOKxVezUOVu3buXo0aP697m9jd999x3R0dH67SqVihEjRpT64FOmTGHatGksXLgQFxcXevfujampKWq1mqCgIBRFITg4mBo1apS6bSGEEEII8eCKTRZz11ksaN26dfnelyZZbNy4sf6GziZNmvD5558X2icwMJDAwECj2qsIyp3b1IreiWJniarBU1UWhxBCCPEwiYiI4OTJkyQnJ5Oeno6TkxN16tThl19+KbTa27p168jMzESj0ZCUlIRKpcLGxgaNRkOdOnUMHuPs2bPMnz+fu3fvcufOHbp06cK4ceOIjo4mKiqq2NVjKtrdu3cZOXIk4eHhlToFT0Ur9mnox1babRp+9j7UNIWXxld1NEIIIcRDISQkBIAtW7YQHx/PpEmTSExM5PLly0U+bRwVFUW9evWIiIgAdAnk8uXLmTp1apHt37p1i4kTJ7J06VKeeuopsrOzefvtt4mKisLFxaXiTswIv//+OzNmzODKlStVGkdFMHq5v8eJyqEhdxs2pebhnyRZFEII8fA6XqDjx7Ge7gvg8lXdV14t86zSVlzdctKoUSM2bdpE69atadeuHWq1GgNrhQCwa9cuvL29eeqppwAwNTVl7ty5mJubExMTQ0JCAqNHj+b69ev4+fnpexyXLVsG6J48z91/7Nix2NnZ4evri7e3N++//z7W1tbY29tTo0YNIiIiiIyMZOvWrahUKvz9/Rk+fLjB2DIyMli+fLl+esBHiSSLBqQ9056aezei3E1FZWlT1eEIIYQQD61z587lWwKwWbNmhISE0LVrVzIyMti0aROhoaG4u7szdepUPDw8imwnKSkJJyenfNusra31r7VaLStWrCA7O5uuXbsybtw4/bB1/fr1WbVqFTt37iQgIIDk5GQ2b96MhYUFL7zwAvPmzcPNzY1FixZx5coVzp07x/bt2/nyyy/1t9t16tTJYA+ml5dXOXynqidJFg1Ie9ob+91fw7nf4Vmfqg5HCCGEKL28PYUFldRTWFzdUnJ1dS1yGDomJgYfHx969epFdnY2//3vfwkNDWXLli1FttOwYUNOnTqVb9vFixe5fPkyAG5ubvoV33Kn+6tfvz7h4eFYWVlx5coVWrduDeieocjdNykpCTc3N0CX9G3fvp0zZ85w6dIl/TMZN2/e5MKFC1U+3F0VqnTqnOrsjmsr2BKHShJFIYQQokJs27aNtWvXArohZQ8Pj2KX9/Xz8+N///sfFy5cACAzM5OIiAjOnDkDUGh+aICpU6cye/ZsIiIieOKJJ/TD3HnnenZ0dOTcuXMAnDhxAgAXFxdcXV357LPPiIyMZNCgQY/tVH5G9Szm5OTk+6Y+FkzNUNnUruoohBBCiIdewWFogNmzZzNhwgRmzZrFgAEDsLS0xMrKivDwcIPt2NjYEBERwdSpU1EUhbS0NPz8/AgKCso3pV9eAwYMIDAwkFq1alGvXj2SkpIK7TNjxgzee+89rKysMDc3p379+nh6euLj48OQIUPIyMigRYsW1K9f/8G+EQ8plVLcnaT3dOzYkf79+zNgwAA8PcuvW7o8lfcqCbntKZfOg01tVLUMP8YvHl+yOocoiVwjoiQPuoJL3rpyvZXNF198Qd++falbty6LFi3C3Nyct956q6rDqlTFXTtG9Syq1Wq2bt3KunXrcHV1ZcCAAQQEBDzyGbbyz3lQt4Rx82Hga1UdjhBCCPFY0Gg0xMXFFdq+Zs0aatasWe7Hs7e3Z9SoUVhZWWFra6ufyqegXbt2FZprGmD48OH07Nmz3OOqLozqWcwVGxvLt99+y/bt20lOTqZdu3YMHDiQnj175nsaqSpUWM+iuhU4u6MK/7rc2haPDvkUL0oi14goifQsiuqguGunVDciPv3000yZMoU9e/YQGRlJVlYWoaGhdOrUiXfffZc//vijXAKuVtr2gOP/Q8nQVnUkQgghhBCVrtRPrRw+fBiNRsP48eOJiYmhY8eOBAcHc+PGDQIDA/n0008rIs6q07Y7pN+BkwerOhIhhBBCiEpn1D2LJ0+eZOvWrezYsYPLly/j5ubGqFGj6N+/Pw4ODoBuvD4kJIQVK1YwcuTICg26UrXsDPP/C83aV3UkQgghhBCVzqhkcfDgwdStW5d+/foxcOBAnnnmmSL38/T05MaNG+UaYFVTWVpD6y5VHYYQQgghRJUwKllcuXIlvr6+mJqaFrvfiBEj9DOdP2oURYHsLFRm5lUdihBCCFEtRUREcPLkSZKTk0lPT8fJyYk6derwyy+/0KxZs3z7rlu3jszMTDQaDUlJSahUKmxsbNBoNNSpU/R0dYcOHWLChAm4urqiKApZWVkMHz4cf39/YmNj2bVrl8Epb7p168aOHTuoUaNGuZ/3o85gspi3h7BVq1bcvn272Ibs7OzKL6pqRrn6D4ztCqOmQl91ifsLIYQQj6OQkBAAtmzZQnx8PJMmTSIxMZHLly8XudxfVFQU9erV009Vs27dOpYvX87UqVMNHqN9+/YsWrQIgLS0NNRqNU2aNOHpp5+WJ8EriMFksX379kUum2NIbGxsuQRULdk7gkoFR36WZFEIIcRDo2vXrvlG/datW8e6devYs2fPA+1TXho1asSmTZto3bo17dq1Q61WU4oZ/bC2tuall15i586d3Lp1i6ioKBYtWkRISAgXLlxAq9Xy6quv4u/vD+jmb0xMTARg2bJlmJqaEhYWxu3bt0lJSeHFF18kKCgItVpNnTp1uHXrFitWrCAkJISkpCQaNGjA4cOH+eWXX/jzzz/54IMPAF2H2ezZs7G1tS3X7091YTBZnD17dqmSxUeZSqVCadsd9m9Dyc5GVcJwvBBCCCHuK7jcX7NmzQgJCaFr165kZGSwadMmQkNDcXd3Z+rUqXh4eBjdtr29PSdPntS/T01N5dChQ2zevBmA/fv368sGDx5MmzZtCAkJYf/+/Tz55JM8//zz9OrViytXrqBWqwkKCgIgICCAnj17sn79eho3bsySJUuIi4ujX79+AEybNo3Zs2fj6urKxo0bWbt2LcHBwQ/0faquDCaLgwYN0r+eN28e/v7+NG/evFKCqpba9oCdX8Cfx+CZtlUdjRBCCFGivL2DUPSzBWXZp7RcXV2LHIaOiYnBx8eHXr16kZ2dzX//+19CQ0PZsmWL0W1funQJR0dH/XsbGxumTZvGtGnTSE1NpX///vqy3DymXr16pKenU69ePdavX88PP/yAjY0NWVlZ+n2bNGkCQFxcHL6+vgA0bdqUunXr6re///77AGRmZur3fxQZ9YDL119/TadOnSo6luqtdVd4JVQ3JF3BFEWBm9cg+W9waITKrp5u+9dL4ND3kHwJUpLByRXVyj26svOn4d/d8jcUsgpVZ91/EmX2GNi/9X6Zs7vULUvd60nw5QJwaAwODbG8k4Xi1BCVTW0qU+4wTW7vv/LDV/B3vO6aSb4Ez3VENWyyruzHKPhoYv4G/nMIVX0nXfmItrp6uboHopr4kdQtbd392+CzCKjXEJ5orPt3wKv6XZXMDFTmFlQmJTsb/o7TXRPJf+u+vPxQ3fvAqyx+B3746n6F2vaovvxdV5aSDMOey9/g2Nmo+o3QlS8Pge2fSd3S1P0gClUrX8R927Ztw9ramuDgYExNTfHw8MDCwvj/J6mpqWzcuJHFixeTnJwMQFJSEidPnmT58uVotVq6dOnCgAEDAAqNmP7nP/+hZcuWBAUFcfDgQfbu3asvy93X3d2dmJgYevTowYULF0hJSQF0yeTcuXNp2LAhR48e1R//UWRUstiyZUt+/PFHvL29S3wi+lGlqlUHhk8BQInZB3G/6xMGHBqhcmioK1MUSLuVv7JFTVQWuqevFO1dSLsNV+/98k6/g6r7i7qypER4J0D3iz3z3ooxk5dDn6G617dTIDsHPFpD3Sd0f4xy1aoD/V8ln4Z5PuW07Qb2edbyrvOE1C1L3auX4Puv4I7uga8nAep+BR36AqB8PEP3x9mhIdRrBE2eRuXdS1eWlamb4D0vS2tUprr/hsrdVMjO1m3PzoLMDFT1GujKzp6Abz6GpL91MST9Dav2gpOrbv+vPoLEs1DXUXds0zxP7Tu5Fz5fqzz31fR8CVJv3n/v3krqlqVuDUuwbwhJiXAyWvf/tV+eOWeHNEfJztL/zqDrC6h6vgzc+xBy5zbUttfdHw1gXev+h4G8xwSwqIHKQrc+rnL+NJw+cv/aSL2Javo63X53U2FkuwLnY3N/dKRtD6hpdb8s32vLwt+LJnkeHmjRAfLODiF1S67r0IjHVcFhaNDd7jZhwgRmzZrFgAEDsLS0xMrKivDw8GLbOnjwIGq1GhMTE7Kzsxk3bhwuLi76ZM3BwYHk5GQGDhyIlZUVo0aNwsys6HTHz88PjUbDd999h52dHaampmRkZOTb51//+hchISEMHTqUhg0b6p+m1mg0TJkyhex7v7dLivthZtTa0NOnT2fLli1YW1vz5JNP6rtg9Y2oVKxcubLCgjRGRa0NXRRlRShsznO+1rVQfXtBV5Z6AwY8lb/Cm7NRDX5TV748BLasKrru3TRY+Pa9nolGuj8qHl76RFRUH0raLUj+m4tHD+LUoz+q2va67cvehWP7dH+0025B6y6o5v9XV3bkZ5gyKH9DS75H1cxbVz6uF5yKvl/WtgeqiE26smN7Ye7Y+4mGQ0MY/Ob9Hq8bV8GmtkztVI0o6XeghiWnT5/G09MTvlqkSyRze/i6vYjq5bd1+365ED6Zmb+B75P1P0+lVz3dB4hcQe+genWaruyzubB+jm67nYPu2lj2Eyozc92H1583Qb0GuuumXgN9kimqD1kbuno7duwYd+7coVOnTpw/f57Ro0fz008/VXVY5a64a8eonsW//vqLVq3uf4JOS0srn8geVmNnw7B3ITlR1wt44+r9MouaMLbAp4tnO9x/3fF5aNT0/i/vPImgytIawtZWcPCiPKisa4F1LdLuok8UAVRvzdO/VtJuwd08/1ecXAtfG/Wd779+4XXoohsqwcQUnN3vt9u6C3x9ynA8925VENWHKk/Pk0qlgqCJhnfuHKD7XXDzWp4G8qzG+vosUHLuv/f0uv+630jo8dK9RDD//HEqlQrujVwI8TDRaDTExcUV2r5mzRpq1qzcDzxOTk5MnDiRZcuWkZWVxfTp0yv1+NWBUT2LD4PK7FkUIpdcJ6Ikco2IkkjPoqgOHrhn8dKlSyXu07ChDJUKIYQQVU1RFJn6TpRKSf2GRiWL3bp1K/HCe6Qn5RZCCCEeAjVr1uTatWvY29tLwiiMoigK165dK3Z436hkcdmyZYW23blzhyNHjrBr1y7mzp1b9igNyMnJQaPR8Oeff2JhYcEHH3zAk08+We7HEUIIIR4VjRs3JjEx8ZGexkWUv5o1a9K4cWOD5UYliz169Chye//+/alXrx5ffPFFuc/D+NNPP5GRkcHXX3/N8ePHiYiIqPInroUQQojqzNzc/JGeHFpUDVONRqN5kAZUKhWrVq3ijTfeKKeQdDZs2EC7du1wd3fH0dGRuXPnMmrUKIP7X716FQcHh3I7/qBex8MaAAAAB3NJREFUfcm4cAnPGrXIuPAPS9d+wvjJk2jXrjMZF/7B9K+/6T4gQL/Pqes3+fjjNbzzzjsEPeNFxoV/yLjwD90HBJByKQmrOk9w+fJVvlixiikTJxLYvI1+nw6DBnD58lVq13bA5twF/rNqNcGTJhHYvA2XL1/lr+QU+vXz53biFVqY25Bx4R/9Pq3b+3L58lUuX77Ki337SMyVHHOtO1ks++Thivlh/D4/zDG/P3cu3t4PV8wP4/f5YYh51Btjyv3vV3n/7ROiKEb1LBZn+/btFbJwdmpqKjY2Nvr3pqamZGVlGZxYU6vVVsh9k7n3fCqKgqJAdnYOOTkKee8Fvb89B0WBgveJKopCdnaO/nXednNf5+TkkJ2dU8R2Rb89N4a8ctuVmCVmibl6xgzKQxfzw/h9fhhiNvQ3Kj09Xe77F9WaUVPnBAQEFNqmKApXr17l5s2bBAcHM2ZM0Z+YymrOnDk899xz+Pv7A+Dr68u+ffsM7i9T54iqINeJKIlcI6Ik5Tl1jhAVwaiexWbNmhX5VJWNjQ2+vr507ty53ANr3bo1u3fvxt/fn+PHj+Pu7l5yJSGEEEIIUa6MShYjIiIqOo5Cevbsyf79+3n55ZdRFIXZs2dXegxCCCGEEI87o+9ZTE9PZ/PmzURHR5OamoqdnR1eXl76hbrLm4mJCTNnzix5RyGEEEL8f3t3F9Jk+8Bx/LtmFmQvTBGCooVUmgn2or2AREkUdVQHMQIP6iCMkVpB0UZtpqsJ0eGWClGmkYJKSB0E1TqIDjIZVGzFyoKik4kxEvFl7jmQZ//Hajz7w2P3sN/n7Lp2M34H496P++W6RGZNWs8sDg0NUVVVxcePH1m3bh25ublEo1HevXvHypUruXPnDrm5uf/2NbMqGAyyYMGCfz9QRERkjhgbG6O0tNToGDLHpVUWz549S39/P62trRQUFCTn379/z/HjxykvL+fKlSuzGlREREREfr956Rz09OlT6urqZhRFgIKCAmprawkEArORTUREREQMllZZBFiyZEnK+dHR0f8skIiIiIhkjrTKYklJCW1tbcTj8Rnzk5OT3Lx5k5KSklkJJyIiIiLGSuuZxXA4jM1mw2KxUFlZSV5eHtFolEePHhGNRrl16xYbN278HXlFRERE5DdKqyzC9CrxPp+P/v5+YrEYS5cuZfPmzVRXV1NcXDzbOUVERETEAGmXRRERERH586RclPvNmzf/1xfNlauLU1NTuN1u3r59S3Z2No2NjaxatcroWGKwiYkJHA4HX758YXx8nBMnTlBZWQlAX18f7e3tdHZ2GpxSMsHQ0BCHDh3ixo0bjI+P43K5MJvNWK1WPB4P8+al/V6hzEE9PT309vYC02skhkIhurq6aGhowGw2k52dTVNTE3l5eQYnFfmflFcWCwsLf7kf9I8SiQQmk4lQKPSfhzPCw4cPefz4MV6vl2AwSHNzM36/3+hYYrDu7m7C4TBOp5Ph4WEOHjxIIBAgFArh9XoZHR2lq6vL6JhisImJCerq6ohEIvh8Pq5du8bhw4fZuXMnZ86c4cCBA+zevdvomJIh6uvrKSwspK+vD6fTSVFREXfv3mVwcJDz588bHU8kKeWVxba2thnjeDzO0aNHqa+vZ/Xq1bMezCgvX76koqICgNLSUl6/fm1wIskE+/btY+/evcmx2WxmeHiYq1ev4nA4uHDhgoHpJFM0NTVhs9loaWkBoKioiG/fvpFIJBgZGSErK+0dVmWOe/XqFZFIBJfLxa5du8jPzwem/2u1G5lkmpRnrvLy8hnjv5fN2bBhw5y55fwr379/JycnJzk2m81MTk7qJP+HW7RoETD9+6ipqaG2than04nD4dCJXYDp24sWi4WKiopkWbRarVy6dAm/38/ixYvZunWrwSklUzQ3N2O32wGSRXFgYID29nY6OjqMjCbyEzWgH+Tk5DAyMpIcT01NqSgKAF+/fsVut3PkyBGsViufPn3C7XYzNjZGJBLB4/HgdDqNjikG6e7uxmQy8fz5c0KhEOfOnSMcDtPb28uaNWvo6OjA6/XicrmMjioGi8VifPjwgW3btiXnHjx4gN/vp6WlBYvFYmA6kZ+pBf1g06ZNPHnyhP379xMMBlm7dq3RkSQDRKNRjh07xsWLF9m+fTsA9+/fB+Dz58+cPn1aRfEP98+rQVVVVbjdbux2e/JORX5+PgMDA0bFkwzy4sULduzYkRzfu3ePzs5Obt++zbJlywxMJvJrKos/2LNnD8+ePcNms5FIJLh8+bLRkSQDXL9+nVgshs/nw+fzAdDa2srChQsNTiaZrLGxkVOnTpGVlcX8+fNpaGgwOpJkgMHBQVasWAFMP+Ll8XhYvnw5J0+eBKCsrIyamhojI4rMkPY6i/F4nOLiYnp6eli/fv1s5xIRERGRDJCyLFZXV/80FwgE2LJly4wXQABMJpOWlxERERGZg1Lehv7nSx5/KysrS/mZiIiIiMw92u5PRERERFLSvlMiIiIikpLKooiIiIikpLIoIiIiIimpLIqIiIhISiqLIiIiIpLSX/x1Y04xCMzoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAACfCAYAAADtTjDeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZdr48e/MJDNJJr1BSCAhtITeITQrUuzoirKiu6yu+rJr3RVEBFdQ5HWXfbeJ5bduAURYwbWggICKFANEAgETpAbSe5lJpp/fH0MCgZBJY2ZC7s91cREmM+fcQ07mPk+7H5WiKApCCCGE6LDUng5ACCGEEG0jyVwIIYTo4CSZCyGEEB2cJHMhhBCig5NkLoQQQnRwksyFEEKIDs7H3Se0Wq3Mnz+f3Nxc1Go1S5YswcfHh/nz56NSqejTpw+LFy9GrZb7DCGEEKI53J7Mv/nmG2w2Gx988AG7d+/m//7v/7BarTz99NOMGTOGRYsWsX37diZPntzkcdLT09HpdO0Wl9lsbtfjiWuDXBfiUnJNiEu565owm80MHTq00e+5PZn37NkTu92Ow+HAYDDg4+NDeno6o0ePBmDSpEns3r3bZTLX6XQkJye3W1yZmZntejxxbZDrQlxKrglxKXddE5mZmVf8ntuTeUBAALm5uUybNo3y8nLeeust9u/fj0qlAkCv11NdXe3yOGazuck31lImk6ldjyeuDXJdiEvJNSEu5Q3XhNuT+T//+U8mTJjAc889R35+Pg8//DBWq7X++0ajkeDgYJfHkZa5cAe5LsSl5JoQl/KGlrnbZ5kFBwcTFBQEQEhICDabjf79+5OamgrAzp07GTlypLvDEq3w6aE8/r33jKfDEEKITs/tLfOf/exnLFiwgFmzZmG1WnnmmWcYOHAgL730EitWrCAxMZEpU6a4OyzRQg6HwrLPnXeJD6UkeDYYIYTo5NyezPV6PX/6058ue3z16tXuDkW0wcFzFeRVmlCrwGp34KuRpYRCCOEp8gksWmXT4XwAHAoUVJo8HI0QQnRuksxFizkcCp9n5BOu1wKQW1Hr4YiEEKJzc9nNvnTp0lYffOHCha1+rfBe358tp6DKxDM39+WP234kp1ySuRBCeJLLZL569WoiIyPRarXNPqjFYqG0tFSS+TXqs8P5aH3UzE6J54/bfiRXkrkQQnhUsybAvfnmmwwePLjZB01PT+f+++9vdVDCe9nPd7Hf0C+KcL2W6CAduRU1ng5LCCE6NZdj5g899BBRUVEtOmh0dDQPPfRQq4MS3uvAmTKKqs3cNrgbALFh/jJmLoQQHuYymS9YsICYmJgWHbRbt24sWLCg1UEJ77UpIx8/XzU3JkUDEBvq36xu9k8P5fHsuvSrHZ4QQnRKLVpnXlJSwiuvvEJqairV1dUoinLZczxdn1ZcPc4u9gJuTIpGr3NeOnFhAWw9WojDoaBWq6742s8z8tn6QyG//8mQJp8nhBCi5VqUzF966SVSU1OZMWMGcXFx9ZujiM4h9XQpJQYztw7qVv9YbJg/FruDYoOZLsF+V3ztqWIjdodCZa2VMH3zJ1MKIYS3WPn1ScICfLl/dA9Ph3KZFiXz7777joULF3LPPfdcrXiEF9t0OB9/X019FztAXKg/ADnltVdM5g6HwulSIwClRrMkcyFEh7R231lUKrwymbeoaIxer6dLly5XKxbhxWx2B5uPFHBTcjT+Wk3947FhzmTe1CS43IpaLDYHAKUGy9UNVAghrpLiajPZpTWcK/O+FTwtSuY/+clPeO+997BY5AO5s0k9XUap0cJtgxtOhoytb5lf+eI+VWKs/7rUKNeOEKLjMZpt1FrtAOw5WeLhaC7nspv98ccfr/9aURS+++47Jk2aRFJSEn5+DbtVVSoVK1eubP8ohcd9djgfvVbD9f2iGzyu1/kQFuDb5Iz208WG+q8lmQshOqISg7n+610nSpk5yru62l0mc6PR2ODfI0aMAMBut1/2PXFtstodbD6Sz839u+Dnq7ns+67Wmp8qMaLXajBa7JRe9AshhBAdRXG187MrKkjHnhMlLlfwuJvLZL5q1Sp3xCG8WFp2OeU1VqYNbLzeQGyoP6eKr3xjd6rYSK/oQM6W1VAmLXMhRAdU1zK/fXA33tt9mmOF1STHBHs4qgtcjplPmDCB3/72t2zcuJG8vDx3xCS8TPq5CgDG9Axv9PuxoQHkVtQ2WncA4FSxgcRIPRF6rUyAE0J0SMXnP7vuGuZcmrv7hHeNm7tsmd988818//33bNq0CUVR6N69O2PHjiUlJYUxY8YQHt74B7y4dhw6V0F8RMAVl5TFhflTY7FTXmOt3xa1Tq3FTl6licSoQPIqTA3GnYQQoqMorjajUkH/mGASo/TsOlHCIxMTPR1WPZfJ/OWXXwbAYDCQlpbGwYMHSUtL45NPPsFisdC3b19SUlIYO3Yso0aNIiAg4GrHLNzs0LkKRiRc+aatfnlaee1lyfz0+ZnsPSP1ZOZXcaLIcNnrhRDC25UYzIQHaPHRqJnQO5IP03Kw2BxofVq0KOyqaXbRmMDAQK677jquu+46AGw2Gz/88ANpaWmkpqayevVqADIyMlwe6+2332bHjh1YrVYeeOABRo8ezfz581GpVPTp04fFixejVnvHf1BnV1RtIq/SxJy4kCs+p255Wm5FDYMued6pEmfyTozSExGoJfW0dLMLITqekmozUUE6AMb1iuTfe7NJP1fB6CsMP7pbqzJmUVERmzZt4j//+Q8ffvghu3btwt/fn+uvv97la1NTUzl48CBr165l1apVFBQUsGzZMp5++mnef/99FEVh+/btrQlLXAWHz1UCMLR76BWfExd2oQrcpeomxvWM1BOu11FeY8HuaHxsXQghvFWJwUxkoDOZpyRGoFbBLi8aN29Wy7y2tpZ9+/axe/du9uzZw8mTJ/H19WXYsGHccccdpKSkMGjQoGbVat+1axd9+/Zl7ty5GAwGnn/+edavX8/o0aMBmDRpErt372by5MlNHsdsNrfrpi4mk0k2iWnE9vQy1CrQVOeTmVnY6HMURcHfR0XGqVwyIxuOiaefLCJKryH75HFshkoUBfalHyXU//Ilbt5IrgtxKbkmOqe8MgPJ0X71P/s+ETq2HT7LtDi7V1wTLpP57NmzSU9Px263k5SUxHXXXccLL7zAyJEj0el0LT5heXk5eXl5vPXWW+Tk5PDEE0+gKEr9jYBer6e6utrlcXQ6HcnJyS0+/5VkZma26/GuFTl7UunXNZhhgwY0+bweESXUqPwv+z8s3V5K35hQkpOTOWnNg9RSImLj6dsl6GqG3W7kuhCXkmuic6o0Z9OrW1T9z/6mbBVv7zxFXM/e5Jw+4ZZroqkbBpfJfP/+/YSGhvLggw8yZcoU+vTp06ZgQkNDSUxMRKvVkpiYiE6no6CgoP77RqOR4GDvWbvXmSmKwuGcSqYP6uryubFhl+9rrigKp4qN3DUsFoAIvfPmr8Rg7jDJXAgh6kq51o2ZA0zoHcmbX59k3+kyujXxWndxOWa+YsUKJk+ezEcffcQdd9zBxIkTef755/nvf/9LUVFRi084YsQIvv32WxRFobCwkNraWlJSUkhNTQVg586djBw5suXvRLS77NIaKmutDIm78nh5ndjQy6vAFRvMVJttJEbpAYgIdM50l8IxQoiOpK76W92YOcDw+DB0PmqvGTd32TKfPn0606dPByA7O5s9e/awd+9eli1bRlVVFYmJiYwbN45x48YxevRo9Hp9k8e74YYb2L9/P/feey+KorBo0SLi4uJ46aWXWLFiBYmJiUyZMqV93p1ok0M5zmIxQ5qY/FYnNsyfylor1SYrQX6+AJw+P/ktMSoQgIjzy9akcIwQoiOpq48ReVHL3M9Xw6iEcPacKOW+3lGeCq1ei/Yzj4+PJz4+ngceeABFUcjIyCA1NZW0tDSefPJJoHlL055//vnLHqtb2ia8R/q5Cvx81fSJDnT53AvL02pJ6upM5nW7pSVGOm/wQgO0qFSy2YoQomOpS+ZRgQ3niY3vHcnyzVmU1YZ5IqwGWpTM65SWlnLw4EHS09M5fPgwmZmZKIrC4MGD2zs+4UGHcyoZFBuCj8b1Csa4iwrHJHV1znk4VWxA66Om2/lEr1GrCA/QymYrQogOpb6bPahhUazxvSMASM+vZbzbo2rIZTK32+1kZWWRnp5en8Bzc3PRaDQMHDiQ0aNH89hjjzFixIjLtkQVHZfV7uBIbiWzx8Y36/n1VeAuGjc/VWykZ4QezUU7C4XrtTJmLoToUIoNFlQqCA9omMwHdAshxN+X9Pwr7xrpLi6T+YgRIzCbzajVavr378/UqVMZPXo0I0eOlNKt17BjBdWYbY5mjZcDROp1aH3UDWa0ny4x0q9rw1nrEYGy2YoQomMpMZiJ0Gsv66XUqFWM6xXBgdPFDZZYe4LLZP7AAw8wduxYRowYQWCg67FTcW2on/zWjJnsAGq1ithQ//oqcFa7g7NlNUwd2HBZW4ReR1ZBVfsGK4QQV1FxtbnBTPaLjesdyRdHCjhdYqyf7OsJLpP5vHnz3BGH8DKHz1USFuBL93D/Zr8mNtSfnPPd7GfLarA5lMsu7ohArUyAE0J0KBeXcr3UsPO9l1kF1d6dzIcPH97sg6lUKtLS0toUkPAOh3IqGNI9tEXdRnFh/mzLdNYeOFW/LK3hUsVwvZaKGitWuwPfZkysE0IITysxmEmIb3zZdcL51Tp1O0R6istkXlNTA0BSUhKTJk1Cq218T2tx7TCabfxYWM2UAa4rv10sNtSfEoMZk9XO6brd0iIb/gJEnL+7La+xEB0kEyaFEN5NUZTz3eyN575AnQ9hfhqyS708ma9atYotW7awZcsW1qxZww033MDUqVMlsV/DjuRW4lBgSPcrb3vamLoZ7XkVtZwqNhKu1xJ6yezPyIsKx0gyF0J4O6PFjsnquGI3O0C3YF/OlNS4MarLuUzmo0aNYtSoUbz44oukpaWxefNmXnnlFYxGI9dffz3Tp09n4sSJktivIXWT3wY3c/JbnbrCMTnlzmR+aascnN3sICVdhRAdQ0kjpVwv1S3Ih0NFnm2ZN3vQUqVSMXLkSBYuXMjOnTt5++23iYiIYMmSJaSkpPDcc8+xbdu2qxmrcJNDOZXEhfk3efE25uK15qdKDJeNl8OFbvYSKRwjhOgA6qu/BTXdMi+qNmM029wV1mVaPQNp5MiRLFiwgC+++IIHH3yQLVu28Otf/7o9YxMecuhcRbPXl1+sa7AfGrWKrPwqSgyWRmd2RkjLXAjRgTS2ycqlugU7S1hnl3quq71V5VwNBgNff/01W7du5dtvv8VmszF+/HhuueWW9o5PuFmJwUxOeS0PpTSv8tvFfDRqugb78e1x5y5CPRvpZg/x90WjVknhGCFEh3Bhk5UrDyV3C3Im8zOlRvp388wW3s1O5uXl5Wzfvp0vv/ySvXv3olarmTBhAr/73e+48cYbpaDMNeJwC4vFXCo2zJ99p8sA6NVIN7tarSIsQNaaCyE6hmKDBbXKWfDqSupa5mc8OKPdZTJfs2YNW7duJS0tDa1Wy6RJk1i+fDnXX389/v7NLygiOoZD5ypRq2BgbMtmsteJC/VnH84yhz3CG1+XGRkom60IITqG4moz4Xptgz0mLhXgqyYyUMcZD641d5nMlyxZglqtZtiwYUycOBGdTkdBQQEffPDBZc9VqVT87Gc/uxpxCjc5lFNB3y5B6HWtGoGp3z2te5g/Wp/Gp2TIZitCiI6iqepvF0uICOCMt4+ZOxwO0tLSXFZ3k2TesSmKQkZOJTcmRbf6GHUz2hsbL68TEajjSG5lq88hhBDuUmIwNzmTvU5CpJ5vjxe7IaLGuUzmWVlZ7ohDeIG8ShOlRguD4lrXxQ4QG+rcSa+pGsUReulmF0J0DMXVZhIirtw4qZMQEcCHaWZqLDYCtK3r2WwLl0vTfvWrX7Fu3TpycnLa9cSlpaVcd911nDx5kuzsbB544AFmzZrF4sWLcTgc7Xou0TwZOc7WcmvHywESIp3J/NKtTy8WoddSZbJhscnPWQjhvRRFOd/N7rooWl2Ndk8tT3OZzEeNGsW2bdu47bbbmDJlCkuWLOGrr76qr9neGlarlUWLFuHn5yznuWzZMp5++mnef/99FEVh+/btrT62aL0juZVo1Cr6x7R+aUVcWAAb/2ccdw+LveJzwgNlrbkQwvvVlXJtVjf7+da7pybBuUzmDz/8MO+++y779u1j8eLFaLVaVqxYwZgxY5g9ezbvvPMOR48ebdFJly9fzv333090tHNs9ujRo4wePRqASZMmsWfPnla8FdFWGbmV9IkOxM9X06bjDO8R1uSOaHVLPEqN0tUuhPBezSkYUyc+wtkr6alJcM3u2NdqtYwbN45x48Yxb948ioqK2LVrF7t27eK9995DrVaTkpLCH/7whyaPs3HjRsLDw5k4cSLvvPMO4OzKqNtqU6/XU11d7TIes9lMZmZmc8N3yWQytevxOhpFUUg/W8ro2ICr/v9QXWIC4OAPJ1BXBlzVc7VVZ78uxOXkmug8jhQ6P6tqygrJzLxyXjKZTOScPkGon4b0k7lkdnF/r2OrR+mjo6OZMWMGM2bMcM6Czshg165dLl+3YcMGVCoVe/fuJTMzk3nz5lFWVlb/faPRSHCw625enU5HcnJya8O/TGZmZrser6PJq6il0nSaCQPiSU5OuKrn0kUaYHMegRFdSU6+cne8N+js14W4nFwTnccZWz6Qx/D+fUhuorJb3TXRu0s5FTbVVbs+mrqJbFYyLysrQ6vVNqjyVl5ezueff05hYSGJiYlMnz6dwYMHuzzWmjVr6r+ePXs2L7/8Mm+88QapqamMGTOGnTt3Mnbs2OaEJdpRRm7bJ781V103u2y2IoTwZsXN2GTlYvERenad8MzytCbHzB0OB0uWLGHChAl8++239Y9nZWUxdepUli1bxoYNG3jhhRe49dZbOXPmTKuCmDdvHn/5y1+YOXMmVquVKVOmtOo4ovWO5Dorv7Vl8ltzBfv74KNWyQQ4IYRXK6k2o1Zd2LrZlZ6RARRWOZenuVuTLfP169ezbt06nn32WcaMGVP/+Pz58wkNDeXTTz8lOjqas2fPMnfuXF555RXee++9Zp981apV9V+vXr26FeGL9uKc/BaEv7Ztk9+aQ6VSERGolc1WhBBerdhgcVnK9WLxEReWpyW7oWF0sSu2zHNzc/n4448ZP348Q4YM4eTJk+zfv5/NmzeTlZXFjTfeSHZ2Nvv376ewsJBbbrmFAwcOcODAAfLy8tz5HkQbKYrCkdxKt3Sx1wnX62SzFSGEV2tuKdc6PevXmrt/edoVW+aZmZnk5OQQERHBxo0b6x/PyspCo9FQUVHR4PHCwkKsVisbNmzg5ptvplu3blc3ctFuCqpMlBgsDIp1351kZKBWlqYJIbxacXXzSrnWqVuedrrE/cvTrpjMb775Zr777jtSU1NZuHAher2empoa7rvvPsaMGcOyZcsaPP/ZZ58lISHhsseF96ur/NaWMq4tFa7XcrbMc5sSCCGEKyUGc5P7TFwqyM+XyECtd7XMAX7+85/z8ccfM3nyZJKSkvjxxx+prKzktddeq3/OP/7xDzZu3MiJEydYunTpVQ9YtL8Lk9/cl8wj9DoZMxdCeK26Uq4taZmDc9zcE/uaNzmbPTY2lo0bNzJt2jRUKhUTJ05k7dq1DZagHTlyhICAAFasWME999xz1QMW7S8jt5Le0YFumfxWJyJQi8Fsw2S1u+2cQgjRXM7PJ0ez6rJfLCFCzxlv6mavk5OTwzPPPNNgjfnFGqv4ZjAYyMjIICUlpe0RiqtKURQycquY1DfSreeN0F+oz94t1N+t5xZCCFdKzvcctmQCHDh3T9vwvYlai92tDSSXtdnnzJnDqVOnWnTQkydPMmfOnFYHJdynsMpMicHMIDfOZIcL6zalq10I4Y1KWlgwpk583Yz2Mvd2tbtsmSuKwnvvvUdkZPNbbiUlJW0KSrhPXeU3dyfziEDZbEUI4b1assnKxXrW755WQ1JX960QcpnMu3XrxuHDh1t84JiYmFYFJNwro27yWxN1h6+GCGmZCyG8WF3LvKXJPD6ybvc0L2uZ79ixwx1xCA85kltJr6hAArSt3nOnVSJkT3MhhBdraSnXOsF+vkTo3b88zeWYubi2ZeRWur2LHSBQ54NWo6ZEutmFEF6o2GAmXK9rdinXi8VHBHC6xMta5uLaVVhlorja7NYyrnXq6rOXSTe7EJ2S2WZn548l5FXUUlhlorDKTFG1iaIqMwmRAbw4vT89zldU84TiakuLl6XVSYjUs/dkaTtH1DRJ5p2YJyq/XSxcr5X67EJ0MnaHwkcHc/njlz+SW1ELgI9aRXSQjuhgP7qHB7DreAmT//gNT97Uh0cnJqL1cX8ncmsKxtRJiNCz8ftcty5Pk2TeiWXkVqJy07anjYkIlM1WrjWVtVZ2nyjB5lDQqFRo1KBRq9GonVX/kmKC0Pm4b+2t8B6KovDlD4W8seUYx4sMDIoNYeldAxkUF0J4gBb1Rd3ZBZUmfvfpUd7YcoyPDuby6l0DGZMY4dZ4i6vNJLaglOvFEs6/7mxZDf26BrVnWFckybwTq5v8ptd55jKI1Gs5VWzwyLlF+1EUhcM5laxJzeaTQ3mYrI4rPtdXo6Jf1yAGxYYyOC6EwXEh9I8JRqVq+bik6DgOnCnj1c8zOXi2gsQoPSt/OpypA7te8efeNcSPlQ+OYEdWIYs+PsrMd77j3hFxPDoxkb5dAq/69VJXyjWy1S3zug1XjN6VzNesWcOaNWvIy8sjLi6OmTNn8tOf/hS1WubPdWQZuZWM7+3eym8XC9drZTZ7B2Yw2/j0UB5rUrM5kluFv6+Gu4bGcu+IOEIDfLE7nF2qdoeCzeGgoNLE4dxKMnIq2XQ4j7X7zgLQK0rPQykJzBgeS5Cfr4fflWhPDofCym9O8oetx4gO8mP5PYO4Z3gcPprm5Y4bk7qQkhjJX3Yc552dp/gwLYdeUXpuHdyNWwfFXLXEbjDbMNtaXsq1zoV9zd03Cc5lMl+zZg1LliyhZ8+e3HDDDWRnZ/Paa6+Rl5fHvHnz3BGjuAqKqkwUeWjyW52IQB01Frvbyx4KqDJZOVlkoNZip8Zip8Zqx2SxU2u1E+TnQ9cQP2JC/Oka7Ff/s8mvrOXAmXLSsp1/fsivwu5QSOoaxJI7B3DnsFiCXSTjaYOc9ScUReFsWQ2pp8pYs+8siz85yv9uzuKeEXE8lBJP7+jmtWZOlxgprjYzKiFMWvdeprLGyrPr09meVcQdQ7qxbMagVvUC+ms1PD81iZ+P78nmowVsOpzHX3cc58/bj9MrSs+M4XH8fHxCuy6vrSvl2tox8xB/X8L1Ws6Uuq9Gu8t3v379eu644w6WL19e/8vyhz/8gdWrV/Ob3/wGjUY+hDuaylorb359EnB/5beL1ReOMZqJ0zY+a7W42syxgmqOFVZzvLCayf27cFNyF3eGeU0oMZjZf7qM1NNl7DtdRmZBFYrSvNeG+Pui81FTdL4ilr+vhqHdQ3niul7ckBTF8B4tT6QqlYr4CD3xEXruG9Wd9HMV/HvvGT7Yd45/781mWI9QhsSFMig2hIGxIfSK0uOjUWO22dl3uowdWUV8fay4fvnPkLgQnp+a5NGeps6m2mTFz1eDbyOt7CO5lTyxJs059n3HAB5KiW/zzVZUkI7ZY+OZPTaeomoTW44U8NnhfN7Ycox/7z3Dc7f0457hca1aSnap1lZ/u1hCRABn3Lg8zWUyz87OZv78+Q1+ELNmzeLdd9/l3LlzJCQktOiEVquVBQsWkJubi8Vi4YknnqB379715+jTpw+LFy+WLvyroLLWyj93n+Hvu05RZbJx6+AYhvUI9Vg8dYVjSg0W4sKcydzuUPjvwVw2HswhK7+6wQQ5rUbNf9Nz+e/c8a0qk2i1O9h0OB+D2YaPWoVarcJHrUKjVuHvqyFMryUsQEu4XkuIv2+7fCi0haIonCw2UFlrpWuIP9FBukY/OC9lttnJyq8m/VwFh85VkH6uglPnP1T8fNUM7xHGUzf1YWC3EAL9fAjQagjQavDz1eDvq6HKZKOg0kRBVS35lSYKKk0YzDYGxYYwIj6M5JjgZsXREkO7hzK0+1BenJ7MB/vPsT2zkHX7z/HPPWfq4+4dHcipYiM1Fjs6HzUpvSL4+fgEtBo1f9lxgp/+v1Qm9I7k+an9GBzXtuu61mKnuNpM93D/Dt3iVxSF40UGwvXaNiUmcN4Q7jtdRuqpUlJPl5FVUI1Wo6ZPl0D6xwTTv1sw/WOCOVls5OVPjxIeoOWDX6YwIj6snd7NBdFBfsxOSWB2SgIHzpSxdFMmz394mH/sPsPCW5Mb3NSZrHZOFBn4Ib+K3PJaVCrQqJy//yoVqFUqTFY7BpMNg9lGtcnG2TJni7ot/2cPj0ug2mRr83ttLpfJ3GQyodc3nNEXFRUFQE1Ny7sQPvnkE0JDQ3njjTcoLy/n7rvvJikpiaeffpoxY8awaNEitm/fzuTJk1t8bNG4S5P45P5dnB/mHmyVw4XKSmVGCw6HwuajBaz48kdOFBnoFaXnpuRo+nUNpl+XIPp1DUJRFG79yy6eWP09n/xqfIvGV612B0+uPcgXRwqa9XyVCkL9fekT7sut5f5M6htFQkTAVf9gLzda2HWihG+PF/Pt8RLyK00NYooM1BET4kd0kB8+ahUKCooCdY3somozmXlVWOzOSWiRgTqGdg/hvlHdGd0znIHdQlwu84kI1NGzlbN42yoiUMfcG3oz94be2B0Kp4oNHMmr5EhuFVkFVdw9LJYbk6IZ1yuywdDMXcNiWZN6lr99dYI7/rqb6YO6cu+IOPpEBxEb6t9gpvSlqkxWjhdWk5FTyeHcSo7kVnKiyIBDgZ+O6cErdw70+I1daxRVm1iw8QjbMgsBiAzUktQ1mKSuzt+nYT3C6B3d+G6YdfIra1mbepbPjxRwosg5WdXfV8PIhDCmDYyhxmrjh7wqdmQV8Z+0nPrXTegdyZ/uH1q/B8PVNDIhnI/+ZxyfHs5n+RdZ/PT/pXJd36l99skAABOASURBVCjCAnz5Ib+Kk8VG7A7X3VABWg2BOh+C/HwI9PPl1sExJEa1/vfgzqGxrX5ta6gUpenOtqSkJNavX99gD3O73c6AAQPYsGEDAwYMaNEJjUYjiqIQGBhIeXk59957LxaLhZ07d6JSqdi2bRu7d+9m8eLFTR4nPT0dna79LhSTyYSfn1+7Hc9bVNTaeWpTLkVGGyndA5g1JIzeEVf/F6w5Cqqt/HzjOab0CeJEqZmTZRa6h/jy0LAwxvfQN5o4Mwpqmb81n5TuAbx4fZdmJVebQ+F/dxbxbbaRR0aGc0NiIHYHOBSl/m+TTaHKbKfK5KDSbKfKbKesxs6h/BryDc4917sG+jAiNoCB0X6EB2gI9XP+CdSpUatUOBSF0ho7+dVW8qutFFTbKKmxYbYpmO0KZpvj/N8KCqA+30LQqJ1fm6wKp8stKECgVs3QGH+Gd/MnMsCH0hrnsUpq7JQYbZTV2nEoCioAlYq6/4UgnZq+ETr6RupIivIjMkDToVuWLWW0OPjohwo2Hq2k1ub8aNP5qOge7EuPUC1Reh8qTHaKjTZKjDaKa2zUWi98BIb5aegdoaV3hI5qs4PPjlUxvkcAz0+KRnu+N8LbPysUReHr00ZW7ivBbFO4f3Aofj5qTpebOVNuIbvCisXufM/dQ3wZ30PP+Hg9vcK1qM5fx4fya/nsWBXfnatBUWBojD9DY/wZ1NWPPhE6fC65uVEUhfJaOyfLLJhsDsb10HvkBshid/BJZhXrj1Sg06hIDNeRGKYlMVxLzzAtMUG+qFXgUOr+OG+GfTWqNsXrzmsiOTm50cfbNGOgNR8Sda18g8HAk08+ydNPP91gPF6v11NdXe3yODqd7opvqjUyMzPb9XjewGZ3MPvv+6g0O/jP4ymMSgj3dEgN9DDbYOM5thyvpkd4ACvuS+bOobFN/lIlJ0Ol5hSvfp7JnlI/HpmY2OQ5bHYHT61L59tsIwtvTXb5/EtlZmbiH9WDnceL2fljMV+dLGXTsaoGz/FRqwgN0FJVa61vEQNo1Cq6BOkI0Png7+uDn5+aYF9nd7ZaxfnZ3g5sDgWHoqBWqbhrZDgT+0QyOC60Q7YGvcHIIfCC2caxgiqOFxo4XuT8k1VYzVenDUTodXQL9SM5LpgbQ/yJCfEjMSqQQbEhdAnWNfhcG7brNEs++4HXdlfxzkMjCfH39erPiuJqMwv/m8GWo0UM6xHK738yhF5RDVvfdofC6RIje0+W8MWRAtYfKeWDjAq6h/szsU8Ue0+WcrrESLhey2PX9WLW6B50D29eJbbxV+NNtdCQgbDwXsWtN7HuuiYyMzOv+L1mJfPly5cTFHT57NLXXnuNwMALF4pKpWLlypUuj5efn8/cuXOZNWsWt99+O2+88Ub994xGI8HBnilicq15Y+sx9p4q5fc/GeJ1iRyc3VpP39yHqCAd943s3uxx2Ecm9iQtu5xlX2QxOC6U0T0bf282u4Nn1h9i0+F8Xpze8kReJyFST0Kkc/mU2WbndImRMoOFEqOFkmozpUYzZUYLwf6+9AgPoEd4APHhemJC/dp9bFk0T6DOhxHx4YyIb3htOBxKk13ul/rFhJ5EBmr5zX8OMfPtvfx7zuj2DrVd2OwOPjqYy2ufZ2K02FkwPYlfTEhs9IZQo1bROzqQ3tGBzE5JoMxo4csfCvjiSAH/OXCOwXGhPHVTH6YN6tphC/x0pt6oOi6T+ahRowBnkm3O466UlJQwZ84cFi1aREpKCgD9+/cnNTWVMWPGsHPnTsaOHduiY4rLbT5SwNvfnOKnY3pw74g4T4fTKJVKxdM3923V6974yWDu+Otu5r7/PZuenEB0UMMuLpvdwbPrD/HpoTxemJbEo5Nal8gvpfPRuHWPYtG+WpLI69w5NJawAC2Pr05jxso9vDQpAm9pl1tsDjZ+n8ObX5/kbFkNQ7uH8vufDG720j5wzl2ZOaoHM0f1QFHc26IV7cflmHl7W7p0KV988QWJiRc+XF988UWWLl2K1WolMTGRpUuXulzy1t7dGt7cddZSp4oN3PHX3fSKDmT9Y2M77N21K8cKqrnrb7vp2yWQMYkRVNVaqTJZqay1kl9p4lSxkfnTknj8ul6tPse1dF2ItjmcU8HP/7GfUqOFIXEh3JzsXCaZHBPk9gRostpZt/8cb31zkvxKE0PiQvj1jX24KTlakrEHuLOb/UrncXsyby+SzBtXY7Fx1992U1xt5rMnJxIb6u/pkK6qj9Nzef7Dw4BzPXSwvy/Bfj6E+PsyuX9XZo3p0abjXyvXhWgf+ZW1vLX5IIdKFA7lVKAoEBvqz03J0fzP9b3pGnL1J0Fl5FQy51/764vl/PrGPkzsEylJ3IO8IZm77Ga//fbbW3SyTz/9tEXPF00zmm1UmawE+Prgr9U0uaxIURTmb8jgeJGBf88Zfc0ncnB2gd42uJtMFhNuERPiz/2Dw/hdcjJF1Sa+yipiW2YR6/afY8vRAv7+8KiruuSzymRl7vvf46tW8cEvxzLWzZuPCO/lMpkfP34clUpFv379pIVyle3IKiQzv5ozJUayS2s4XWqsr0RUx1fjLHASoPW5aImFcn6ZlUJ5jZXf3NKXiX2iPPQu3E8SufCE6CC/+rHmrIIq5vxjPz95ay9/fmAYk/u3f5VCRVFYsDGD3Ipa1v1yLCO9cFKr8ByXyXzZsmVs2bKF3bt3U1tby9SpU5k2bRpJSUnuiK/TOHi2nDn/PABAl2Ad8RF6bugXRXyEnrAALbVWO7UWG8bztcyNZmdlIbVKhVp9/m+ViviIAOaM7+nJtyJEp5PUNZj/zh3PI/8+wC9XHWDhrf2ZMz6hXbu+1x84x2eH8/ntlH6SyMVlXCbzu+++m7vvvhuDwcD27dvZvHkz//jHP4iJiZHE3o42HynAV6Ni7ws3tbnsohDC/aKD/fjgl2N5Zl06Sz77gTMlRhbf3r/ZO4Q15XhhNYs/Ocr43hFtmtAprl3NLhoTGBjInXfeyZ133lmf2Ldu3crMmTPp2rUrU6dOZerUqdIV3wqK4ixlOq5XpCRyITqwAK0PK386gtc3Z/HOzlPklNew8sER+Pk2vaLkTImRzUcLuHVQzGUFWkxWO796/yB6rQ9/vG+oDCuJRrXqlrEusf/tb39jz549TJ06lb///e/MmDGjvePrFLIKqskurWHKgK6eDkUI0UZqtYoF05NZctdAvv6xmMdWpWGy2q/4/OOF1dz71l5e/yKLSW98xZx/7mdHVmF9PfEln/3AscJq/nDfEKKDvbeMrPCsVpdzPXjwIFu3buXLL78kJyeH3r17c8stt7RnbJ3GlqMFqFRclUkzQgjPmD02Hq1GxbwNGTy+Oo23Z4+4rObDD3lVPPj3VDRqFWseGUPqqVLW7j/HnH8eIC7Mn4l9Ilm77xyPTUrk+n7RHnonoiNodjJ3OBykpqby5Zdfsm3bNoqKikhOTuaee+7hlltuoVcvGcdprc1HChgZH0ZUkHSxC3EtmTmqBw4FXtiYwROrv2flg8PrE3pGTiUP/j0Vf18N7z86hsSoQMb3juTXN/Vh69FCVn+Xzdp95xjSPZTnbunn4XcivJ3LZP7NN9+wZcsWduzYQWVlJQMHDuShhx5iypQpdO/e3R0xXtOyS41kFVSz8FaZayDEteiB0T1QFFjwUQb/s/p73nxwOEfzqnj4vX0E+/my9tGx9Ii4ME7uq1Fz6+AYbh0cw7myGkIDfF1uWyuEy2T+2GOPoVarGTZsGJMnT6Zbt26AsxJNYzu4SFd7y2w56txfW8bLhbh2zRrTA4eisPC/R3j4vX1k5FQSGaRjzSNjiAu78o5kzd2tTIhmdbM7HA7S0tJIS0tr8nkqlarJLdrE5bYcLWRAt2D5pRXiGvfg2HgUReGlj4+SGKXn/UfGuqX8q+gcXCbz7du3uyOOTqmoykRadjnPTW75zmFCiI5ndkoC/buF0CtKT2iA1tPhiGuIy2QeGxvbrAOVlpaSm5vb7OcL2PJDIQBTBkoXuxCdxYj4ME+HIK5B7Tar4rPPPmPmzJntdbhOYevRAhIj9fSJDvR0KEIIITowmSLpIZU1VvaeLGXKwK6ydaEQQog2kWTuIduzCrE5FJnFLoQQos0kmXvI5iMFxIT4Mfgq7n0shBCic/CaZO5wOFi0aBEzZ85k9uzZZGdnezqkq6bGYmPn8WJu6d8FtWyaIIQQoo1czmZfunRpsw7U1vXl27Ztw2KxsG7dOtLT03n99ddZuXJlm47prXb+WIzJ6pBZ7EIIIdqFy2S+Y8eOZh8sJiam1YGkpaUxceJEAIYOHcqRI0dafayW2n2ihGc+PItak+eW8xnMNsICfBmdEO6W8wkhhLi2tWsybwuDwUBg4IUlWhqNBpvNho9P4yGazeZ2qzZXW2VlaFctKnXTew63H1+Gxfhz/MdjbjqfaC2TySRVDUUDck2IS3nDNdHqLVDbW2BgIEajsf7fDofjiokcQKfTkZzcPpuTJAOxwZntdjxx7cjMlOtCNCTXhLiUu66Jpm4YvGYC3PDhw9m5cycA6enp9O0rJU6FEEKI5vCalvnkyZPZvXs3999/P4qi8Nprr3k6JCGEEKJD8JpkrlareeWVVzwdhhBCCNHhqBRFUTwdRGukp6ej0+k8HYYQQgjhFmazmaFDhzb6vQ6bzIUQQgjh5DUT4IQQQgjROpLMhRBCiA5OkrkQQgjRwUkyF0IIITo4SeZCCCFEByfJXAghhOjgvKZojKc4HA5efvlljh07hlarZenSpcTHx3s6LOFmVquVBQsWkJubi8Vi4YknnuCmm24C4NNPP2X16tWsW7fOw1EKTyktLWXGjBm89957WCwWFi9ejEajISEhgVdffRW1WtpFncXGjRv56KOPgAsbfq1fv54lS5ag0WjQarUsX76cyMhIt8bV6a/Ai/dRf+6553j99dc9HZLwgE8++YTQ0FDef/993n33XZYsWQI4Nzb48MMPkXIMnZfVamXRokX4+fkB8Ne//pW5c+eydu1aLBYLX3/9tWcDFG41Y8YMVq1axapVqxgwYAALFy7k1Vdf5aWXXmLVqlVMnjyZd9991+1xdfpk7sl91IX3mDp1Kk899VT9vzUaDeXl5fz+979nwYIFHoxMeNry5cu5//77iY6OBiA5OZmKigoURcFoNDa5u6O4dmVkZHDixAlmzpzJihUr6ndNs9vtHqlO2umT+ZX2URedi16vJzAwEIPBwJNPPslTTz3Fiy++yIIFC9Dr9Z4OT3jIxo0bCQ8Pr7/hB+q71qdNm0ZpaSljxozxYITCU95++23mzp0LUH+j9/3337N69Wp+9rOfuT2eTn9L2dJ91MW1Kz8/n7lz5zJr1iwSEhLIzs7m5Zdfxmw2c+LECV599VVefPFFT4cp3GjDhg2oVCr27t1LZmYm8+bNIysri48++og+ffqwZs0aXn/9dRYvXuzpUIUbVVVVcerUKcaOHVv/2Oeff87KlSt55513CA8Pd3tMnT5rDR8+nK+++orp06fLPuqdWElJCXPmzGHRokWkpKQAsGnTJgBycnJ49tlnJZF3QmvWrKn/evbs2bz88svMnTu3vjcvOjqa77//3lPhCQ/Zv38/48aNq//3xx9/zLp161i1ahWhoaEeianTJ3PZR10AvPXWW1RVVfHmm2/y5ptvAvDuu+/WT3oSos7SpUt55pln8PHxwdfXt36ypOg8Tp8+TVxcHOAcI3/11VeJiYnh17/+NQCjRo3iySefdGtMsmuaEEII0cF1+glwQgghREcnyVwIIYTo4CSZCyGEEB2cJHMhhBCig5NkLoQQQnRwksyF6GRmz57NY4891uj3UlNT6devHxkZGW6OSgjRFpLMhRBCiA5OkrkQQgjRwUkyF0JcUVZWFo888gijR49m9OjR/Pa3v6WkpKT++/Pnz+e2225r8Jpt27bRr18/cnJyAGe3/ksvvcQvfvELhg8fzvLly936HoToDCSZC9EJKYqCzWa77I/D4ah/TmZmJjNnzsRqtfL666+zYMECDhw4wIMPPkhNTU2Lzrdx40bi4uL485//zLRp09r77QjR6XX62uxCdEbffPMNAwYMaPI5b775JuHh4bz77rtotVoABg4cyO23386GDRuYPXt2s8+n1+tZuHAhvr6+bYpbCNE4SeZCdEIjRozghRdeuOzxo0eP1m/nuX//fm677bb6RA7Qu3dv+vXrx/79+1uUzHv06CGJXIirSJK5EJ1QUFAQgwYNuuzxi7vPq6qqiIiIuOw5ERERGAyGFp2vseMIIdqPjJkLIRoVEhJCaWnpZY+XlJTU79msUqkajLMDGI1Gt8QnhLhAkrkQolEjRoxg+/btWCyW+sdOnjzJjz/+yPDhwwHnWHhpaWmDhJ6Wlub2WIXo7CSZCyEa9fjjj1NcXMyjjz7KV199xSeffMKjjz5KbGwsd911FwCTJk2ioqKC3/3ud6SmpvKnP/2Jbdu2eThyITofSeZCiEYNHDiQf/3rX9hsNp566ileffVVRo4cydq1awkMDAScyfyZZ55h+/bt/PKXvyQzM5PXX3/dw5EL0fmoFEVRPB2EEEIIIVpPWuZCCCFEByfJXAghhOjgJJkLIYQQHZwkcyGEEKKDk2QuhBBCdHCSzIUQQogOTpK5EEII0cFJMhdCCCE6uP8Pxq8uu21u3LMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## V23 Streamlining for Greatlakes\n",
    "\n",
    "## V22 Improving the day-by-day DAM suboptimizations into a whole year rather than full year.\n",
    "# Adds in initial states into the formulation (lb, ub # mods, reactor gen, elec gen, one/off state)\n",
    "# Takes a first pass with reduced hours to find the number of modules needed \n",
    "# Then fixes the number of modules for each daily optimization (need to add a lower bound of number lb = ub)\n",
    "# \n",
    "\n",
    "## V21 Adding in subsets of the year start date and number of hours\n",
    "\n",
    "## V20 Adding in the hourly profiles. will do this using the naics code\n",
    "\n",
    "## V19 Making the script more able to handle the per facility format. Demands in particular are now not a csv which is brough in, but a single hourly value. 06/13/2022\n",
    "\n",
    "## V18 Dropping the Avoided cost of ff in the objective function. Add in more post process for comparing and analizing revenue breakdown 06/09/2022\n",
    "\n",
    "## V17 More cleanup\n",
    "\n",
    "## V16 Cleaning up some constraintes based upon the writing. Making the total charge and discharge for the TES for TES scale not for the modules.\n",
    "\n",
    "## V15 Making TES paramters optimization -- Lied, did not do that. Just updated some things\n",
    "\n",
    "## V14 Removing Fast solve as it was far slower. And inputing TES without optimization\n",
    "\n",
    "## V13 This is a short-cut for resolving results more quickly, by brute force\n",
    "\n",
    "## V12 Adding the fast solve functionality with solely a float value number of modules option. \n",
    "# Adding this possibility on top of current layout to provide layered effect.\n",
    "# Can run fast, choose best 3, then run full analysis.\n",
    "\n",
    "## V11 Adding in thermal penalties for coolant to working fluid transfer. 01/21/22\n",
    "# This is primarily done in the Generator spreadsheet 'sites.csv'\n",
    "# The Temperature coefficient is 0.75 for all systems currently and are adjusted to occur prior to importing.\n",
    "# Energy loss to the heat transfer are adjusted to the maximum thermal capacity BUT the electric capacity already accounts for this in the efficiceny therefore it is only applied to the pTCap value.\n",
    "\n",
    "## V10 Fixing the MSL for thermal baseline issue.\n",
    "# This means that the vTGen is now the total thermal generation from the 'boiler' aspect of the cycle.\n",
    "# vGen remains unchanged\n",
    "# The thermal demand for the facility is satisfied by (vTGen-vGen/pEff)\n",
    "# This allows for the system to have a more specific link between the thermal requirements of the nuclear reactor itself which has a distinct ramping and range\n",
    "# and for the electricicty to behave a certain way without tryign to directly couple them in vTGenabovemin which would provide complications if the \n",
    "# vTGen were simply the heat to industry.\n",
    "# new value is put into the results vHeatGen which is the MWht for the supplied industrial heating.\n",
    "# this opens the question is the pTVOC is solely the thermal cycle and if the pVOC is just to operate the turbine and generate electricty outside of the thermal behavior.\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import cplex\n",
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import math\n",
    "\n",
    "import plotly.graph_objects as go\n",
    "import plotly.io as pio\n",
    "import plotly.express as px\n",
    "\n",
    "import pyomo.environ as pe\n",
    "import pyomo.opt\n",
    "from pyomo.opt import SolverStatus\n",
    "\n",
    "from datetime import datetime\n",
    "s = datetime.now()\n",
    "print(s,'start')\n",
    "\n",
    "\n",
    "import Facility_Processing_2015 as FP\n",
    "\n",
    "class EconomicAssessment_OneGen:\n",
    "    def __init__(self):\n",
    "        \n",
    "        self.Generator_IDs = None             # Names of the egenrators being evaluated\n",
    "        \n",
    "        #### sets ####\n",
    "        self.G = None                         # Generator indexes [0,1,2....]\n",
    "        self.T = None                         # Timestep indexes [0,1,2....8759]\n",
    "        self.M = None #list(range(6))              # Maximum number of modules for any generator\n",
    "        #self.YR = None                        # Number of years\n",
    "        \n",
    "        #### paramteres ####\n",
    "        \n",
    "        # Environmental paramters/System\n",
    "        self.pLMP = None                      # LMP electric $MWhe [t] or [t,yr]\n",
    "        #self.pTLMP = None                     # LMP heat $/MWht [g,t] or [t,yr]\n",
    "        self.pTempdemanded = None             # Minimum temperature which the reactor must supply for the given TLMP single value\n",
    "        self.pTEnergyDemand = None            # Thermal energy required per hour in MWt [t]\n",
    "        \n",
    "        # Engineering parameters for generators\n",
    "        self.pOnoroffinitial = None           # initial state, default, 1 = On [g]\n",
    "        self.pStartupfixedcost = None         # $/MW [g]\n",
    "        self.pVOC = None                      # Variable O&M ($/MWhe) [g]\n",
    "        self.pTVOC = None                     # Variable O&M ($/MWht) [g]\n",
    "        self.pCap = None                      # Maximum Electric capacity MWe == pTCap*ThermalEff [g]\n",
    "        self.pTCap = None                     # Maximum Thermal capacity MWt [g]\n",
    "        self.pThermaleff = None               # Thermal efficiency [g]\n",
    "        self.pThermaTransfereff = None        # V11 Thermal transsfer efficicnecy between coolant and workign fluid.\n",
    "        self.pMSL = None                      # Minimum stable thermal load (MWt) [g]\n",
    "        self.pMSLTurb = None                  # Minimum stable trubine load minumum core + turbine power [g]\n",
    "        self.pMDT = None                      # Minimum down time (hours) [g]\n",
    "        self.pGenabovemininitial = None       # Initial generation above the min (MWe) [g]\n",
    "        self.pRGenInitial = None              # Initial generation of reactor (MWth) [g]\n",
    "        self.pRamprate = None                 # Maximum ramp (up and down) (MW/hr) [g]\n",
    "        \n",
    "        self.pCAPEX = None                    # Capital costs of the generator $/MWe [g]\n",
    "        self.pFOPEX = None                    # Fixed O&M costs of the generator $/MWe [g]\n",
    "        self.pModulesize = None               # Module size in MWe [g]\n",
    "        self.pOutlettemp = None               # SMR outlet temperature [g]\n",
    "        self.pWorkingTemp = None              # SMR working fluid temp \n",
    "        self.pMaxMod = None                   # maximum number of modules [g]\n",
    "        self.pMinMod = None                   # minimum number of modules [g] V222 71522    \n",
    "        \n",
    "        self.pIR = 0.07                       # Discount/ Interest Rate 7% (OMB Circular)\n",
    "        self.pLifetime = 30                   # Lifetime, years\n",
    "        \n",
    "        # V14 Additions MV 3/1/2022\n",
    "        ## TES Parameters ####\n",
    "        self.pTESCapacity = None              # Storage thermal capacity (MWh-thermal)\n",
    "        self.pTESChargeEff = None             # Charging Efficiency (%) \n",
    "        self.pTESDischargeEff = None          # Discharging Efficicency (%) \n",
    "        self.pTESLossRate = None              # Rate at which the stored thermal energy is lost (%/hr)\n",
    "        self.pTESPowerCap = None              # Power in or out maximum (MWt)\n",
    "        \n",
    "        ## TES Variables #### \n",
    "        self.vTESSOC = None                   # TES State of Charge (MWht)\n",
    "        self.vTESSOC_Charge = None            # TES Charge amount (MWht)\n",
    "        self.vTESSOC_Discharge = None         # TES Discharge amount (MWht)\n",
    "        # V15\n",
    "        self.vTESCapacity = None              # DV of the total storage (MWht)\n",
    "        self.vTESPowerCap = None              # DV of CHarge/Discharge Power (MWt)\n",
    "        # V22 \n",
    "        self.pTESInitSOC = None               # Start SOC of TES (MWht)\n",
    "        self.pTES_Ch0 = None\n",
    "        self.pTES_Dis0 = None\n",
    "        #### variables ####\n",
    "        \n",
    "        # Float variables\n",
    "        self.vGen = None                      # Electric generation (MWhe) [g,t] --> [t]\n",
    "        self.vTGen = None                     # Heat generation (MWht) [g,t] --> [t]\n",
    "        self.vGenabovemin = None              # Electric generation above the MSLTurb (MWhe) [m,t]\n",
    "        self.vTGenabovemin = None             # Thermal Generation above MSL (MWht) [m,t]\n",
    "        \n",
    "        # Binary variables\n",
    "        self.vTurnon = None                   # Does the generator turn on (1) or not (0) [m,t] --> [t]\n",
    "        self.vTurnoff = None                  # Gnerator turns off(1) or not (0) at time t [m,t] --> [t]\n",
    "        self.vOnoroff = None                  # Generator state at time t [m,t] --> [t]\n",
    "        self.vEOnoroff = None                 # Electricity generation on (1) or off (0) [m,t] --> [t]\n",
    "        self.vTOnoroff = None                 # Heat generation on (1) or off (0) [m,t] --> [t]\n",
    "        \n",
    "        # results/objective variables  --These are not explicit model components and ar calculated from other variable values\n",
    "        self.rEProfit = None                  # Profits from vGen*LMP [m,t]\n",
    "        self.rTProfit = None                  # Profits from vTgen*TLMP [m,t]\n",
    "        self.rProfit = None                   # Hourly profits ($) [m,t]\n",
    "        \n",
    "        self.vTotalProfits = None             # Total Profits ($) [-] single value\n",
    "        self.rFeasible = None                 # is this configurations feasible     \n",
    "        \n",
    "        self.vGenerator = None                # Index of the generator chosen for the\n",
    "        \n",
    "        #V18 \n",
    "        self.rACC = None                      # Annualized capital cost\n",
    "        self.rFOMC = None                     # Fixed cost\n",
    "        self.rVC = None                       # Variable cost\n",
    "        self.rFC = None                       # Fuel cost\n",
    "        self.rEP = None                       # Electricity profit\n",
    "        self.rAVC = None                      # Avoided cost of natural gas\n",
    "        self.rSUC = None                      # Startup costs\n",
    "        self.rTES = None                      # TES CAPEX\n",
    "        self.RampCost = 1.0\n",
    "        \n",
    "        \n",
    "    def BuildModel(self, DAM = False):\n",
    "        pTS = []\n",
    "        for i in range(len(self.G)):\n",
    "            if self.pOutlettemp[i] >= self.pTempdemanded:\n",
    "                pTS.append(1)\n",
    "            else:\n",
    "                pTS.append(0)\n",
    "            \n",
    "        self.pThermalSat = pTS\n",
    "        \n",
    "        if not DAM:\n",
    "            print('Building Model')\n",
    "        model = pe.ConcreteModel()\n",
    "        \n",
    "        #### SETS ####     \n",
    "        \n",
    "        model.G = pe.Set(initialize = self.G, doc = \"Generator IDs\")\n",
    "        model.T = pe.Set(initialize = self.T, doc = \"Timestep\")\n",
    "        model.M = pe.Set(initialize = self.M, doc = \"Modules\")\n",
    "        #model.YR = pe.Set(initialize = self.YR, doc = 'Years of model runs-- CURRENTLY NOT USED')\n",
    "        \n",
    "        \n",
    "        #### PARAMETERS ####\n",
    "        \n",
    "        model.pLMP = pe.Param(model.T,initialize = {self.T[i]: self.pLMP[i] for i in range(len(self.T))} , doc = \"LMP Electric $/MWhe\")\n",
    "        #model.pTLMP = pe.Param(model.G, model.T,initialize = self.pTLMP, doc = \"LMP Thermal $/MWht\")\n",
    "        model.pTEnergyDemand = pe.Param(model.T,initialize = {self.T[i]: self.pTEnergyDemand[i] for i in range(len(self.T))} , doc = \"Energy demand MWht\")\n",
    "        \n",
    "        model.pOnoroffinitial = pe.Param(model.M,initialize = {self.M[i]: self.pOnoroffinitial[i] for i in range(len(self.M))}, doc = \"Inital state, 1 on, 0 off\")\n",
    "        model.pStartupfixedcost = pe.Param(model.G,initialize = {self.G[i]: self.pStartupfixedcost[i] for i in range(len(self.G))}, doc = \"Startup costs, $500 default\")\n",
    "        #model.pVOC = pe.Param(model.G,initialize = {self.G[i]: self.pVOC[i] for i in range(len(self.G))}, doc = \"Variable O&M Electric $/MWhe\")\n",
    "        model.pTVOC = pe.Param(model.G,initialize = {self.G[i]: self.pTVOC[i] for i in range(len(self.G))}, doc = \"Variable O&M Thermal $/MWhe\")\n",
    "        model.pCap = pe.Param(model.G,initialize = {self.G[i]: self.pCap[i] for i in range(len(self.G))} , doc = \"Maximum Electric Capacity MWe\")\n",
    "        model.pTCap = pe.Param(model.G,initialize = {self.G[i]:self.pTCap[i] for i in range(len(self.G))}, doc = \"Maximum Thermal Capacity MWt\")\n",
    "        model.pThermaleff = pe.Param(model.G,initialize = {self.G[i]: self.pThermaleff[i] for i in range(len(self.G))}, doc = \"Thermal efficiency (MWe/MWt)\")\n",
    "        model.pThermalTransfereff = pe.Param(model.G,initialize = {self.G[i]: self.pThermaTransfereff[i] for i in range(len(self.G))}, doc = \"Thermal transfer efficiency (MW -coolant/MWt- working fluid)\")\n",
    "        model.pMSL = pe.Param(model.G,initialize = {self.G[i]: self.pMSL[i] for i in range(len(self.G))}, doc = \"Minimum Stable Load Thermal MWt\")\n",
    "        model.pMSLTurb = pe.Param(model.G,initialize = {self.G[i]: self.pMSLTurb[i] for i in range(len(self.G))}, doc = \"Minimum Stable Load of Turbine MWe\")\n",
    "        model.pMDT = pe.Param(model.G,initialize = {self.G[i]: self.pMDT[i] for i in range(len(self.G))} , doc = \"Mimimum Down time (hours)\")\n",
    "        #model.pGenabovemininitial = pe.Param(model.M,initialize = {self.M[i]: self.pGenabovemininitial[i] for i in range(len(self.M))}, doc = \"Initial Generation above the Minimun MWe\")\n",
    "        model.pRGenInitial = pe.Param(model.M,initialize = {self.M[i]: self.pRGenInitial[i] for i in range(len(self.M))}, doc = \"Initial Generation MWt\")\n",
    "        model.pRamprate = pe.Param(model.G,initialize = {self.G[i]: self.pRamprate[i] for i in range(len(self.G))}, doc = \"Ramp rate up and down MW/hr\")        \n",
    "        model.pFOPEX = pe.Param(model.G,initialize = {self.G[i]: self.pFOPEX[i] for i in range(len(self.G))}, doc = \"FOPEX pparamter $/MWe\")\n",
    "        model.pCAPEX = pe.Param(model.G,initialize = {self.G[i]: self.pCAPEX[i] for i in range(len(self.G))}, doc = \"CAPEX values$/MWe\")\n",
    "        model.pWorkingTemp = pe.Param(model.G,initialize = {self.G[i]: self.pWorkingTemp[i] for i in range(len(self.G))}, doc = \"working temperature which the reactor can supply in Degrees C\")\n",
    "        model.pThermalSat  = pe.Param(model.G,initialize = {self.G[i]: self.pThermalSat[i] for i in range(len(self.G))}, doc = \"SMR can supply the temperature\")\n",
    "        model.pMaxMod = pe.Param(model.G, initialize = {self.G[i]: self.pMaxMod[i] for i in range(len(self.G))}, doc = \"maximum number of modules\")\n",
    "        #V22 MV 7/15/2022\n",
    "        model.pMinMod = pe.Param(model.G, initialize = {self.G[i]: self.pMinMod[i] for i in range(len(self.G))}, doc = \"maximum number of modules\")\n",
    "        model.pTESInitSOC = pe.Param(model.G, initialize = {self.G[i]: self.pTESInitSOC[i] for i in range(len(self.G))}, doc = \"initial SOC\")\n",
    "        if DAM:\n",
    "            model.pTES_Ch0 = pe.Param(model.G, initialize = {self.G[i]: self.pTES_Ch0[i] for i in range(len(self.G))}, doc = \"initial charge DAM\")\n",
    "            model.pTES_Dis0 = pe.Param(model.G, initialize = {self.G[i]: self.pTES_Dis0[i] for i in range(len(self.G))}, doc = \"initial discharge DAM\")\n",
    "\n",
    "        \n",
    "        # V14 Additions MV 3/1/2022\n",
    "        model.pTESCapacity = pe.Param(model.G,initialize = {self.G[i]: self.pTESCapacity[i] for i in range(len(self.G))}, doc = \"TES max storage\")\n",
    "        model.pTESChargeEff = pe.Param(model.G,initialize = {self.G[i]: self.pTESChargeEff[i] for i in range(len(self.G))}, doc = \"TES Charge Efficiency\")\n",
    "        model.pTESDischargeEff = pe.Param(model.G,initialize = {self.G[i]: self.pTESDischargeEff[i] for i in range(len(self.G))}, doc = \"TES Discharge Efficiency\")\n",
    "        model.pTESLossRate = pe.Param(model.G,initialize = {self.G[i]: self.pTESLossRate[i] for i in range(len(self.G))}, doc = \"TES Loss Rate if just sitting\")\n",
    "        model.pTESPowerCap = pe.Param(model.G,initialize = {self.G[i]: self.pTESPowerCap[i] for i in range(len(self.G))}, doc = \"Maximum power in or out of stroage\")\n",
    "        \n",
    "        #### VARIABLES ####\n",
    "        model.vGenerator = pe.Var(model.G,within=pe.Binary,doc = \"Choice of generator\")\n",
    "        model.vModuleS = pe.Var(model.G,model.M,within = pe.Binary, doc = \"choice of modules for full optimization\")\n",
    "        \n",
    "        model.vGen = pe.Var(model.G,model.M,model.T,within=pe.NonNegativeReals, doc = \"Electric generation (MWhe)\")\n",
    "        model.vTGen = pe.Var(model.G,model.M,model.T,within=pe.NonNegativeReals, doc = \" Heat generation (MWht) \")\n",
    "        model.vRGen = pe.Var(model.G,model.M,model.T,within=pe.NonNegativeReals, doc = \" Total Reactor generation (MWht) \")\n",
    "        \n",
    "        model.vGenabovemin = pe.Var(model.G,model.M,model.T,within=pe.NonNegativeReals, doc = \"Electric generation above the MSLTurb (MWhe)\")\n",
    "        model.vRGenabovemin = pe.Var(model.G,model.M,model.T,within=pe.NonNegativeReals, doc = \"Thermal Generation above MSL (MWht)\")\n",
    "        \n",
    "        model.vTurnon = pe.Var(model.G,model.M,model.T,within=pe.Binary, doc = \"Does the generator turn on (1) or not (0)\")\n",
    "        model.vTurnoff = pe.Var(model.G, model.M,model.T,within=pe.Binary, doc = \"Generator turns off(1) or not (0) at time t\")\n",
    "        model.vOnoroff = pe.Var(model.G,model.M,model.T,within=pe.Binary, doc = \"Generator state at time t\")\n",
    "        model.vEOnoroff = pe.Var(model.G,model.M,model.T,within=pe.Binary, doc = \"Electricity generation on (1) or off (0)\")\n",
    "        model.vTOnoroff = pe.Var(model.G,model.M,model.T,within=pe.Binary, doc = \"Heat generation on (1) or off (0)\")\n",
    "        \n",
    "        # V14 Additions MV 3/1/2022\n",
    "        model.vTESSOC = pe.Var(model.G,model.T,within=pe.NonNegativeReals, doc = \"TES State of Charge\")\n",
    "        model.vTESSOC_Charge = pe.Var(model.G,model.M,model.T,within=pe.NonNegativeReals, doc = \"TES SOC charge\")\n",
    "        model.vTESSOC_Discharge = pe.Var(model.G,model.T,within=pe.NonNegativeReals, doc = \"TES SOC charge\")\n",
    "        \n",
    "        # V15 \n",
    "        #model.vTESCapacity = pe.Var(model.G,within=pe.NonNegativeReals, doc = \"DV of TES Capacity\")\n",
    "        #model.vTESPowerCap = pe.Var(model.G,within=pe.NonNegativeReals, doc = \"DV of TES Power Charge/Discharge\")\n",
    "        \n",
    "        if not DAM:\n",
    "            print('Params & Vars established')\n",
    "        \n",
    "        \n",
    "        #### OBJECTIVES ####\n",
    "        \n",
    "        def ERevenues(model):\n",
    "            return sum(sum(sum(model.vGen[g,m,t]*model.pLMP[t]for g in model.G) for m in model.M) for t in model.T)\n",
    "        \n",
    "        #def TRevenues(model):\n",
    "        #    return sum(sum((sum(model.vTGen[g,m,t] for m in model.M)+ model.vTESSOC_Discharge[g,t]*model.pTESDischargeEff[g])*model.pTLMP[g,t] for g in model.G) for t in model.T)\n",
    "        \n",
    "        def RampCost(model):\n",
    "            return sum(sum(sum((((model.vRGenabovemin[g,m,model.T.prevw(t)])-(model.vRGenabovemin[g,m,t]))**2)*self.RampCost for t in model.T) for m in model.M) for g in model.G)\n",
    "        \n",
    "        def Costs(model):\n",
    "            return (sum(model.pStartupfixedcost[g]*sum(sum(model.vTurnon[g,m,t] for t in model.T) for m in model.M) for g in model.G) \n",
    "                    + sum(sum(sum((model.vRGen[g,m,t]*model.pTVOC[g]) for g in model.G) for m in model.M) for t in model.T) \n",
    "                    + sum(sum(model.pCap[g]*(model.pCAPEX[g]*model.vModuleS[g,m]*(self.pIR/(1-((1+self.pIR)**(-1*self.pLifetime))))+model.pFOPEX[g]*model.vModuleS[g,m]) for m in model.M) for g in model.G)*(len(self.T)/8760)\n",
    "                    + sum(model.vGenerator[g]*model.pTESCapacity[g]*self.pTESCapex*(self.pIR/(1-((1+self.pIR)**(-1*self.pLifetime))))*(len(self.T)/8760) for g in model.G))\n",
    "        \n",
    "        \n",
    "        def Obj_Profit(model):\n",
    "            return ERevenues(model)-Costs(model)-RampCost(model) # +TRevenues(model)\n",
    "        model.Obj_Profit = pe.Objective(rule = Obj_Profit, sense=pe.maximize, doc = \"Maximize the profits by balancing thermal and electric generation\")\n",
    "                \n",
    "\n",
    "        #### CONSTRAINTS ####\n",
    "        \n",
    "        model.GeneratorChooser = pe.Constraint(expr = sum(model.vGenerator[g] for g in model.G)<=1)\n",
    "        model.ModChooserS = pe.Constraint(expr = sum(sum(model.vModuleS[g,m] for g in model.G)for m in model.M)>=1)\n",
    "\n",
    "        def ThermalEquality(model,t):\n",
    "            return sum(sum(model.vTGen[g,m,t] for m in model.M) + model.vTESSOC_Discharge[g,t]*model.pTESDischargeEff[g] for g in model.G) == model.pTEnergyDemand[t]\n",
    "        model.ThermalEquality = pe.Constraint(model.T, rule = ThermalEquality, doc = \"Limits the thermal output (and therefore the 'revenue from heat')\")\n",
    "        \n",
    "        def ModuleLimit(model,g):\n",
    "            return sum(model.vModuleS[g,m] for m in model.M) <= model.pMaxMod[g]\n",
    "        model.ModuleLimit = pe.Constraint(model.G,rule = ModuleLimit, doc = \"limits the number of available modules\")\n",
    "        \n",
    "        def ModuleLimitLower(model,g):\n",
    "            return sum(model.vModuleS[g,m] for m in model.M) >= model.pMinMod[g]\n",
    "        model.ModuleLimitLower = pe.Constraint(model.G,rule = ModuleLimitLower, doc = \"limits the lower number of available modules\")\n",
    "        \n",
    "        def ModGenS(model,g,m):\n",
    "            return model.vModuleS[g,m] <= model.vGenerator[g]\n",
    "        model.ModGenS = pe.Constraint(model.G,model.M,rule = ModGenS, doc = \"cannot build modules for the non-optimal gens Slow\")\n",
    "        \n",
    "        def ModStatus(model,g,m,t):\n",
    "            return model.vOnoroff[g,m,t] <= model.vModuleS[g,m]\n",
    "        model.ModStatus = pe.Constraint(model.G,model.M, model.T, rule = ModStatus, doc = \"The status of the plant must be the previous state plus the turnon/off\")\n",
    "        \n",
    "        def ECapacityLimit(model,g,m,t):\n",
    "            return model.vGen[g,m,t] <= model.pCap[g]\n",
    "        model.ECapacityLimit = pe.Constraint(model.G,model.M,model.T,rule = ECapacityLimit, doc = \"The sum of the electric generation must not surpass maximum generation\")\n",
    "        \n",
    "        def RGenEquality(model,g,m,t):\n",
    "            return (model.vTGen[g,m,t] + model.vGen[g,m,t]/model.pThermaleff[g] + \n",
    "                    model.vTESSOC_Charge[g,m,t]/model.pTESChargeEff[g]) == model.vRGen[g,m,t]*model.pThermalTransfereff[g] # Updated V14 MV 03012022\n",
    "        model.RGenEquality = pe.Constraint(model.G,model.M,model.T,rule = RGenEquality, doc = \"electricity and thermal must equal ttoal\")\n",
    "        \n",
    "        def EStatus(model,g,m,t): \n",
    "            return model.vEOnoroff[g,m,t] <= model.vOnoroff[g,m,t]\n",
    "        model.Estatus = pe.Constraint(model.G,model.M,model.T, rule = EStatus, doc = \"The electric generation cannot be on when the overall status is off\")\n",
    "        \n",
    "        def Estatus2(model,g,m,t):\n",
    "            return model.vGen[g,m,t] <= model.pCap[g]*model.vEOnoroff[g,m,t]\n",
    "        model.Estatus2 = pe.Constraint(model.G,model.M,model.T,rule = Estatus2, doc = \"The sum of the electric generation must not surpass maximum generation\")\n",
    "        \n",
    "        def RStatus(model,g,m,t):\n",
    "            return model.vRGen[g,m,t]*model.pThermaleff[g]<= model.pCap[g]*model.vOnoroff[g,m,t]\n",
    "        model.RStatus = pe.Constraint(model.G,model.M,model.T,rule = RStatus, doc = \"The sum of the electric and thermal generation must not surpass maximum generation\")\n",
    "        \n",
    "        def TAllow(model,g,m,t):\n",
    "            return model.vOnoroff[g,m,t] <= model.pThermalSat[g]\n",
    "        model.TAllow = pe.Constraint(model.G,model.M,model.T, rule = TAllow, doc = \"Must satisfty heat constraint\")\n",
    "        \n",
    "        def RGenAbove(model,g,m,t):\n",
    "            return model.vRGenabovemin[g,m,t]  == ((model.vRGen[g,m,t])-(((model.pMSL[g])/model.pThermaleff[g])*model.vOnoroff[g,m,t]))\n",
    "        model.RGenAbove = pe.Constraint(model.G,model.M,model.T, rule = RGenAbove, doc = \"This defines the generation above the minimum stable load\")\n",
    "        \n",
    "        def EGenAbove(model,g,m,t):\n",
    "            return model.vGenabovemin[g,m,t]  == (model.vGen[g,m,t]-(model.pMSLTurb[g]*model.vEOnoroff[g,m,t]))\n",
    "        model.EGenAbove = pe.Constraint(model.G,model.M,model.T, rule = EGenAbove, doc = \"This defines the generation above the minimum stable load\")\n",
    "        \n",
    "        def RRamprateUp(model,g,m,t):\n",
    "            if t == model.T[1]: \n",
    "                return (model.pRGenInitial[m])-(model.vRGenabovemin[g,m,t]) >= -1*model.pRamprate[g]/(model.pThermaleff[g])\n",
    "                #return pe.Constraint.Feasible # V11C THis means that we can start iur reactor at any output.\n",
    "            return (model.vRGenabovemin[g,m,model.T.prev(t)])-(model.vRGenabovemin[g,m,t]) >= -1*model.pRamprate[g]/(model.pThermaleff[g])\n",
    "        model.RRamprateUp = pe.Constraint(model.G,model.M,model.T, rule = RRamprateUp, doc = \"The hourly change must not exceed the ramprate between time steps\")\n",
    "        \n",
    "        def RRamprateDn(model,g,m,t):\n",
    "            if t == model.T[1]:\n",
    "                return ((model.pRGenInitial[m])-(model.vRGenabovemin[g,m,t]))*(model.pThermaleff[g]) <= model.pRamprate[g]\n",
    "                #return pe.Constraint.Feasible # V11C THis means that we can start iur reactor at any output.\n",
    "            return (model.vRGenabovemin[g,m,model.T.prev(t)]-model.vRGenabovemin[g,m,t])*(model.pThermaleff[g]) <= model.pRamprate[g]\n",
    "        model.RRamprateDn = pe.Constraint(model.G,model.M,model.T, rule = RRamprateDn, doc = \"The hourly change must not exceed the ramprate between time steps\")\n",
    "        ''''''\n",
    "        def Status(model,g,m,t):\n",
    "            if t == model.T[1]:\n",
    "                #return pe.Constraint.Feasible\n",
    "                return model.vOnoroff[g,m,t] == (model.pOnoroffinitial[m]+ model.vTurnon[g,m,t] - model.vTurnoff[g,m,t])\n",
    "            return model.vOnoroff[g,m,t] == (model.vOnoroff[g,m,model.T.prev(t)] + model.vTurnon[g,m,t] - model.vTurnoff[g,m,t])\n",
    "        model.Status = pe.Constraint(model.G,model.M, model.T, rule = Status, doc = \"The status of the plant must be the previous state plus the turnon/off\")\n",
    "        \n",
    "        def SingleGenStatus(model,g,m,t):\n",
    "            return model.vOnoroff[g,m,t] <= model.vGenerator[g]\n",
    "        model.SingleGenStatus = pe.Constraint(model.G,model.M, model.T, rule = SingleGenStatus, doc = \"The status of the plant must be the previous state plus the turnon/off\")\n",
    "        '''\n",
    "        def TurnOn(model,g,m,t):\n",
    "            if t == model.T[1]:\n",
    "                return model.vTurnon[g,m,t] == (model.vOnoroff[g,m,t]-model.pOnoroffinitial[g] + model.vTurnoff[g,m,t])\n",
    "            return model.vTurnon[g,m,t] == (model.vOnoroff[g,m,t]-model.vOnoroff[g,m,model.T.prev(t)] + model.vTurnoff[g,m,t])\n",
    "        model.TurnOn = pe.Constraint(model.G,model.M,model.T, rule = TurnOn, doc = \"Ensure turn on variable goes from off to on states\")\n",
    "        \n",
    "        def TurnOff(model,g,m,t):\n",
    "            if t == model.T[1]:\n",
    "                return model.vTurnoff[g,m,t] == (model.pOnoroffinitial[g] - model.vOnoroff[g,m,t] + model.vTurnon[g,m,t])\n",
    "            return model.vTurnoff[g,m,t] == model.vOnoroff[g,m,model.T.prev(t)] - model.vOnoroff[g,m,t] + model.vTurnon[g,m,t]\n",
    "        model.TurnOff = pe.Constraint(model.G,model.M,model.T, rule = TurnOff, doc = \"Ensure turn off variable goes from on to off states\")\n",
    "        '''\n",
    "        def MinDownTime(model,g,m,t):\n",
    "            if t == model.T[1]:\n",
    "                return pe.Constraint.Feasible\n",
    "            if t <= model.pMDT[g]:\n",
    "                return sum((1-model.vOnoroff[g,m,model.T.prev(t,td+1)]) for td in range(t)) >= t*model.vTurnon[g,m,t]\n",
    "            return sum((1-model.vOnoroff[g,m,model.T.prev(t,td+1)]) for td in range(model.pMDT[g])) >= model.pMDT[g]*model.vTurnon[g,m,t]\n",
    "        model.MinDownTime = pe.Constraint(model.G,model.M,model.T, rule = MinDownTime, doc = \"We must enforce minimum down time\")\n",
    "        \n",
    "        # V14 Updates\n",
    "        \n",
    "        def TES_SOC(model,g,t):\n",
    "            if t == model.T[1]:\n",
    "                if DAM == True:\n",
    "                    return model.vTESSOC[g,t] == (model.pTESInitSOC[g]*(1-model.pTESLossRate[g]) \n",
    "                                      - model.pTES_Dis0[g] + model.pTES_Ch0[g])\n",
    "                else:\n",
    "                    return model.vTESSOC[g,t] == model.pTESInitSOC[g]\n",
    "            return model.vTESSOC[g,t] == (model.vTESSOC[g,model.T.prev(t)]*(1-model.pTESLossRate[g]) \n",
    "                                          - model.vTESSOC_Discharge[g,model.T.prev(t)]\n",
    "                                          + sum(model.vTESSOC_Charge[g,m,model.T.prev(t)] for m in model.M))\n",
    "        model.TES_SOC = pe.Constraint(model.G,model.T, rule = TES_SOC, doc = \"SOC is equal to the in/out plus the previous step with losses.\")\n",
    "        \n",
    "        def TES_SOC_Max(model,g,t):\n",
    "            return model.vTESSOC[g,t] <= model.pTESCapacity[g]\n",
    "        model.TES_SOC_Max = pe.Constraint(model.G,model.T, rule = TES_SOC_Max, doc = \"SOC has to be below the capacity limit\")\n",
    "        \n",
    "        def TES_SOC_Power(model,g,t):\n",
    "            return model.vTESSOC_Discharge[g,t] <= model.pTESPowerCap[g]\n",
    "        model.TES_SOC_Power = pe.Constraint(model.G, model.T, rule = TES_SOC_Power, doc = \"TES cannot surpass the power out max\")\n",
    "        \n",
    "        def TES_SOC_Power2(model,g,t):\n",
    "            return sum(model.vTESSOC_Charge[g,m,t]for m in model.M) <= model.pTESPowerCap[g]\n",
    "        model.TES_SOC_Power2 = pe.Constraint(model.G, model.T, rule = TES_SOC_Power2, doc = \"TES cannot surpass the power out max\")\n",
    "        \n",
    "        self.model = model\n",
    "        if not DAM:\n",
    "            print('Model Built')\n",
    "        \n",
    "        \n",
    "    def SolveModel(self, solver='cplex_direct', DAM = False):\n",
    "        self.BuildModel(DAM = DAM)\n",
    "        if not DAM:\n",
    "            print('Solving...')\n",
    "        opt = pyomo.opt.SolverFactory(solver,tee = True) #iterlim = 1,\n",
    "        opt.options.mipgap = 0.1\n",
    "        opt.options.optimalitytarget = 3\n",
    "        \n",
    "\n",
    "        results = opt.solve(self.model, tee = False, logfile=\"CPLEX_test_log.log\")\n",
    "        # print('>>Solver status is {} and solver termination condition is {}'.format(results.solver.status,results.solver.termination_condition))\n",
    "        \n",
    "        if results.solver.termination_condition =='optimal':\n",
    "            print(pe.value(self.model.Obj_Profit))\n",
    "            self.vTotalProfits = pe.value(self.model.Obj_Profit)\n",
    "            self.rFeasible = True\n",
    "            self.ResultData(DAM = DAM)\n",
    "            self.ResultOutput(DAM = DAM)\n",
    "            \n",
    "        else:\n",
    "            self.rFeasible = False\n",
    "            self.vTotalProfits= None\n",
    "            self.rACC = 0\n",
    "            self.rFOMC = 0\n",
    "            self.rVC = 0\n",
    "            self.rFC = 0\n",
    "            self.rEP = 0\n",
    "            self.rAVC = 0\n",
    "            self.rSUC = 0\n",
    "            self.rTES = 0\n",
    "            self.ModCount = 0\n",
    "            \n",
    "    def ResultData(self, DAM = False):\n",
    "        self.vGen = np.array([[[pe.value(self.model.vGen[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        self.vTGen = np.array([[[pe.value(self.model.vTGen[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        self.vRGen = np.array([[[pe.value(self.model.vRGen[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        \n",
    "        self.vTESSOC_Discharge = np.array([[pe.value(self.model.vTESSOC_Discharge[g,t]) for g in self.G] for t in self.T])\n",
    "        self.vTESSOC_Charge = np.array([[[pe.value(self.model.vTESSOC_Charge[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        self.vTESSOC = np.array([[pe.value(self.model.vTESSOC[g,t]) for g in self.G] for t in self.T])\n",
    "        \n",
    "        # Binary variables\n",
    "        self.vTurnon = np.array([[[pe.value(self.model.vTurnon[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        self.vTurnoff = np.array([[[pe.value(self.model.vTurnoff[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        self.vOnoroff = np.array([[[pe.value(self.model.vOnoroff[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        self.vEOnoroff = np.array([[[pe.value(self.model.vEOnoroff[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        #self.vTOnoroff = np.array([[[pe.value(self.model.vTOnoroff[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        self.vGenabovemin = np.array([[[pe.value(self.model.vGenabovemin[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        self.vRGenabovemin = np.array([[[pe.value(self.model.vRGenabovemin[g,m,t]) for g in self.G] for m in self.M] for t in self.T])\n",
    "        \n",
    "        self.vTotalProfits = pe.value(self.model.Obj_Profit)\n",
    "        self.vGenerator = [pe.value(self.model.vGenerator[g]) for g in self.G]\n",
    "        self.vModuleS = np.array([[pe.value(self.model.vModuleS[g,m]) for g in self.G] for m in self.M])\n",
    "        self.vModDAM = [i[0] for i in self.vModuleS]\n",
    "        self.ModCount = int(self.vModuleS.sum())\n",
    "        \n",
    "    def ResultOutput(self, DAM = False):\n",
    "        self.vRGen_PD = pd.DataFrame(self.vRGen[:,:,self.vGenerator.index(1)])\n",
    "        self.vGen_PD = pd.DataFrame(self.vGen[:,:,self.vGenerator.index(1)])\n",
    "        self.vTGen_PD = pd.DataFrame(self.vTGen[:,:,self.vGenerator.index(1)])\n",
    "        self.vTurnon_PD = pd.DataFrame(self.vTurnon[:,:,self.vGenerator.index(1)])\n",
    "        self.vOnoroff_PD = pd.DataFrame(self.vOnoroff[:,:,self.vGenerator.index(1)])\n",
    "        \n",
    "        self.vTESSOC_Discharge = pd.DataFrame(self.vTESSOC_Discharge[:,self.vGenerator.index(1)])\n",
    "        self.vTESSOC_Charge = pd.DataFrame(self.vTESSOC_Charge[:,:,self.vGenerator.index(1)])\n",
    "        \n",
    "        self.rHeatDelivered = self.vTGen_PD  - (self.vGen_PD/self.pThermaleff[self.vGenerator.index(1)])\n",
    "        \n",
    "        self.rERevenue = self.vGen_PD.multiply(self.pLMP, axis = 0)\n",
    "        if not DAM:\n",
    "            print(self.Generator_IDs[self.vGenerator.index(1)],\"was chosen for operation\")\n",
    "            print(sum(sum(self.vModuleS)),\"modules used\")\n",
    "        \n",
    "        self.Output = pd.DataFrame()\n",
    "        m = [i for i, x in enumerate(self.vModuleS[:,0].tolist()) if x == 1]\n",
    "        for x in m:\n",
    "            self.Output['Reactor_Gen [MWht]_'+str(x)] = self.vRGen_PD[x]\n",
    "            self.Output['Thermal_Gen [MWht]_'+str(x)] = self.vTGen_PD[x]\n",
    "            self.Output['Elec_Gen [MWhe]_'+str(x)] = self.vGen_PD[x]\n",
    "            self.Output['Elec_Gen [MWht]_'+str(x)] = self.Output['Elec_Gen [MWhe]_'+str(x)]/float(self.pThermaleff[0])\n",
    "            self.Output['TES_Charge_'+str(x)] = self.vTESSOC_Charge[x]\n",
    "        self.Output['TES_Disharge'] = self.vTESSOC_Discharge[0]\n",
    "        self.Output['SOC'] = self.vTESSOC[:,self.vGenerator.index(1)]\n",
    "        self.Output\n",
    "        \n",
    "        self.rACC = float(sum(self.pCap)*sum(self.pCAPEX)*sum(self.vModuleS)*(self.pIR/(1-((1+self.pIR)**(-1*self.pLifetime))))*(len(self.T)/8760))\n",
    "        self.rFOMC = float(sum(self.pCap)*sum(self.pFOPEX)*sum(self.vModuleS)*(len(self.T)/8760))\n",
    "        self.rVC = float(self.vRGen_PD.values.sum()*self.pVOC.values.sum()*self.pThermaleff[self.vGenerator.index(1)])\n",
    "        self.rFC = float(self.vRGen_PD.values.sum()*self.FuelCost.values.sum()*sum(self.pThermaleff))\n",
    "        self.rEP = float(self.rERevenue.values.sum())\n",
    "        self.rAVC = float((self.vTGen_PD*self.NGprice).values.sum())\n",
    "        self.rSUC = float(sum(self.pStartupfixedcost)* sum(sum(sum(self.vTurnon))))\n",
    "        self.rTES = float(sum(self.pTESCapacity)*self.pTESCapex*(self.pIR/(1-((1+self.pIR)**(-1*self.pLifetime))))*(len(self.T)/8760))\n",
    "        \n",
    "        \n",
    "def ParamsVarsPD(base,sites,lmps, dem, TES, naics, EP, heatLMP = 4, minTemp = 500, MaxMod = 2, hrProfiles = None,startHr = 0, hrcount=8760, DAM = False, DAM_data = None, init_run = False):\n",
    "    \n",
    "    # Environmental paramters/System\n",
    "    base.G = list(range(len(sites['Sites'])))\n",
    "    base.T = list(range(hrcount))\n",
    "    base.M = list(range(MaxMod)) #will need to be fixed eventually!!!!!!!!!!!!!!!!!!!!!!\n",
    "    base.Generator_IDs = sites['Sites'].tolist()\n",
    "    \n",
    "    \n",
    "\n",
    "    base.NGprice = heatLMP\n",
    "    HLMP = pd.DataFrame()\n",
    "    for x in base.G:\n",
    "        HLMP[x] = [heatLMP]*len(base.T)\n",
    "    base.pTLMP_PD = HLMP\n",
    "    HLMP = np.array(HLMP)\n",
    "    HLMP_Dict ={(base.G[i],base.T[j]):HLMP[j,i] for i in range(len(base.G)) for j in range(len(base.T))}\n",
    "    base.pTLMP = HLMP_Dict\n",
    "    \n",
    "    base.pLMP = (lmps['Settlement Point Price']).tolist()[startHr:startHr+hrcount]\n",
    "    \n",
    "    \n",
    "    base.pTempdemanded = minTemp\n",
    "    \n",
    "    if naics== None:\n",
    "        base.pTEnergyDemand = [dem]*len(base.pLMP)\n",
    "    \n",
    "    if naics != None:\n",
    "        ###########################################################3333\n",
    "        TotalHeat = dem*8760\n",
    "        TotalHours = sum(hrProfiles)\n",
    "        CorrectiveFactor = TotalHeat/TotalHours\n",
    "        NormalizedHour = [hrProfiles[i] * CorrectiveFactor for i in range(len(hrProfiles))]\n",
    "        base.pTEnergyDemand = NormalizedHour[startHr:startHr+hrcount]\n",
    "    if init_run:\n",
    "        base.pTEnergyDemand[::6]\n",
    "    \n",
    "    \n",
    "    # Engineering parameters for generators\n",
    "    base.pCAPEX = (sites['CAPEX $/kWe']*1000).tolist()   \n",
    "    base.pFOPEX = (sites['FOPEX $/kWe']*1000).tolist()   \n",
    "    base.pOutlettemp = sites['Outlet Temp (C)'].tolist()\n",
    "    base.pWorkingTemp = sites['Working Fluid Temp (C)'].tolist()\n",
    "        \n",
    "    base.pStartupfixedcost = sites['Startupfixedcost in $'].tolist()\n",
    "    \n",
    "    \n",
    "    \n",
    "    base.pCap = sites['Power in MWe'].tolist() \n",
    "    base.pTCap = (sites['Power in MWt']).tolist() \n",
    "    base.pThermaleff = (sites['Thermal Efficiency']).tolist() \n",
    "    base.pThermaTransfereff = sites['Thermal Transfer Efficiency'].tolist() #V11\n",
    "    base.pMSL = sites['MSL in MWe'].tolist()\n",
    "    base.pMSLTurb = sites['MSL_turb in MWe'].tolist()\n",
    "    \n",
    "        \n",
    "    \n",
    "    #base.pVOC = sites['Effective VOC in $/MWh-e'].tolist() # V18  \n",
    "    #base.pTVOC = sites['Effective HVOC in $/MWh-t'].tolist() # V18\n",
    "    \n",
    "    base.pVOC = sites['VOM in $/MWh-e']\n",
    "    base.FuelCost = sites['FC in $/MWh-e']\n",
    "    \n",
    "    base.pTVOC = ((base.pVOC+base.FuelCost)*base.pThermaleff).tolist()\n",
    "    \n",
    "    base.pMDT = sites['MDT in hours'].tolist()\n",
    "    diff = [base.pCap[i] - base.pMSL[i] for i in range(len(base.pMSL))]\n",
    "    base.pGenabovemininitial = [diff[i]/1.2 for i in range(len(diff))] # This value is currently not being used as it was too often limiting in the conditions. \n",
    "    base.pRamprate = sites['Ramp Rate (MW/hr)'].tolist()\n",
    "    \n",
    "    CapValue = TES\n",
    "    base.pTESCapacity = [CapValue]*len(base.G)\n",
    "    base.pTESChargeEff = [0.99]*len(base.G) # MIT report - >98% RTE\n",
    "    base.pTESDischargeEff = [0.99]*len(base.G) # MIT report - >98% RTE\n",
    "    base.pTESLossRate = [0.01]*len(base.G)\n",
    "    base.pTESPowerCap = [CapValue/EP]*len(base.G)\n",
    "    base.pTESCapex = 20000 # $/MWh-th generally from table 3 of https://aip.scitation.org/doi/pdf/10.1063/1.4984433 + MIT report For Nitrate Salt Storage\n",
    "\n",
    "    base.pTESInitSOC = [0]*len(base.G)\n",
    "    \n",
    "    base.YR = None\n",
    "    \n",
    "    # V22\n",
    "    if DAM:\n",
    "        base.pMaxMod = [DAM_data[0]]*len(base.G)\n",
    "        base.pMinMod = [DAM_data[0]]*len(base.G) \n",
    "        base.pTESInitSOC = [DAM_data[3]]*len(base.G)\n",
    "        base.pRGenInitial = []\n",
    "        for n in range(DAM_data[0]):\n",
    "            d = DAM_data[1][n]-base.pMSL[0]\n",
    "            if d > 0:\n",
    "                base.pRGenInitial.append(d)\n",
    "            else:\n",
    "                base.pRGenInitial.append(0)\n",
    "        base.pOnoroffinitial = DAM_data[2]\n",
    "        \n",
    "        base.pTES_Ch0 = [DAM_data[4]]*len(base.G)\n",
    "        base.pTES_Dis0 = [DAM_data[5]]*len(base.G)\n",
    "        \n",
    "    else:\n",
    "        base.pMaxMod = [MaxMod]*len(base.G) # THis can and likely will need to be changed, but for now it will work nicely.\n",
    "        base.pMinMod = [0]*len(base.G) \n",
    "        base.pRGenInitial = [0]*MaxMod\n",
    "        base.pOnoroffinitial = [0]*len(base.M) #sites['Onoroffinitial'].tolist()\n",
    "        \n",
    "def ImportData(sites , lmps , hrcount = 8760, startHr = 0):\n",
    "    Sites_pd = pd.read_csv(sites)\n",
    "    LMP_pd = pd.read_csv(lmps)\n",
    "    \n",
    "    LMP_pd = LMP_pd.iloc[startHr:startHr+hrcount]    \n",
    "    return Sites_pd, LMP_pd\n",
    "\n",
    "def NGTempCostCurve(Temp,NG_Cost = 3.0, AHF_Coeffs = [0,-0.00038,0.90556]):\n",
    "    HHV = 54 # MJ/kg\n",
    "    Density = 0.68 # kg/m3\n",
    "    cfTom3 = 35.31 # Unit conversion\n",
    "    AHF = AHF_Coeffs[0]*(Temp^2) + AHF_Coeffs[1]*(Temp) + AHF_Coeffs[2] # avaialble Heat fraction - Deep Patel Equation\n",
    "    \n",
    "    HHV = HHV*Density*(1/cfTom3)*(1/1000000)*(1000) # returns  TJ/thousand cf \n",
    "    Cost = NG_Cost*(1/HHV)*(1/AHF)*(1/277.778) # returns the Cost in $/MWh\n",
    "    return Cost\n",
    "\n",
    "def QuickGraphOutputs(Output,LMP,name,hr = 72,demand = False, TES = True):\n",
    "    \n",
    "    labels = Output.columns.tolist()\n",
    "    labels_2 = []\n",
    "    \n",
    "    if TES:\n",
    "        for i in labels:\n",
    "            if 'Reactor' in i:\n",
    "                pass\n",
    "            elif 'SOC' in i:\n",
    "                pass\n",
    "            elif 'MWhe' in i:\n",
    "                pass\n",
    "            else:\n",
    "                labels_2.append(i)\n",
    "        OG = Output[labels_2]\n",
    "\n",
    "        t = 0\n",
    "        td = hr\n",
    "        dashes = []\n",
    "        colors = ['xkcd:orangered','xkcd:cyan','pink']\n",
    "        pal = []\n",
    "        l = len(OG.columns.tolist())\n",
    "        d1 = 3\n",
    "        d2 = 0\n",
    "        while len(dashes)<l-1:\n",
    "            dashes+=[(d1, d2),(d1, d2),(d1, d2)]\n",
    "            pal+=colors\n",
    "            d2+=2\n",
    "        dashes+=[(1,2)]\n",
    "        pal+=['xkcd:black']\n",
    "        sns.set_style(\"whitegrid\")\n",
    "        ticks_24 = list(range(0,td+1,24))\n",
    "\n",
    "        fig = plt.figure()\n",
    "        plt.figure(figsize = (8,2))\n",
    "        g = sns.lineplot(data=OG.loc[t:td,:], palette=pal, dashes=dashes)\n",
    "        g.set_ylabel('Hourly Power [MWt]', fontsize=16)\n",
    "        #g.lines[0].set_linestyle(\"--\")\n",
    "        plt.xticks(ticks_24)\n",
    "        plt.legend(bbox_to_anchor=(1.01, 1), loc='upper left', borderaxespad=0)\n",
    "        g.figure.savefig(str(name)+'_Reactors.png', dpi = 1000, bbox_inches = \"tight\")\n",
    "\n",
    "        fig2 = plt.figure()\n",
    "        plt.figure(figsize = (8,2))\n",
    "        g2 = sns.lineplot(data = LMP.loc[t:td+1,'Settlement Point Price'])\n",
    "        g2.set_ylabel('LMP [$/MWh]', fontsize=16)\n",
    "        g2.set_xlabel('Hour', fontsize=16)\n",
    "        plt.xticks(ticks_24)\n",
    "        g2.figure.savefig(str(name)+'_LMPS.png', dpi = 1000, bbox_inches = \"tight\")\n",
    "    else:\n",
    "        for i in labels:\n",
    "            if 'Reactor' in i:\n",
    "                pass\n",
    "            elif 'SOC' in i:\n",
    "                pass\n",
    "            elif 'MWhe' in i:\n",
    "                pass\n",
    "            elif 'TES' in i:\n",
    "                pass\n",
    "            else:\n",
    "                labels_2.append(i)\n",
    "        OG = Output[labels_2]\n",
    "\n",
    "        t = 0\n",
    "        td = hr\n",
    "        dashes = []\n",
    "        colors = ['xkcd:orangered','xkcd:cyan']\n",
    "        pal = []\n",
    "        l = len(OG.columns.tolist())\n",
    "        d1 = 3\n",
    "        d2 = 0\n",
    "        while len(dashes)<l-1:\n",
    "            dashes+=[(d1, d2),(d1, d2)]\n",
    "            pal+=colors\n",
    "            d2+=2\n",
    "        sns.set_style(\"whitegrid\")\n",
    "        ticks_24 = list(range(0,td+1,24))\n",
    "\n",
    "\n",
    "        fig = plt.figure()\n",
    "        plt.figure(figsize = (16,4))\n",
    "        g = sns.lineplot(data=OG.loc[t:td,:], palette=pal, dashes=dashes)\n",
    "        g.set_ylabel('Hourly Power [MWt]', fontsize=16)\n",
    "        #g.lines[0].set_linestyle(\"--\")\n",
    "        plt.xticks(ticks_24)\n",
    "        plt.legend(bbox_to_anchor=(1.01, 1), loc='upper left', borderaxespad=0)\n",
    "        g.figure.savefig(str(name)+'_Reactors.png')\n",
    "\n",
    "        fig2 = plt.figure()\n",
    "        plt.figure(figsize = (16,4))\n",
    "        g2 = sns.lineplot(data=LMP.loc[t:td+1,'Settlement Point Price'])\n",
    "        g2.set_ylabel('LMP [$/MWh]', fontsize=16)\n",
    "        g2.set_xlabel('Hour', fontsize=16)\n",
    "        plt.xticks(ticks_24)\n",
    "        g2.figure.savefig(str(name)+'_LMPS.png')\n",
    "    \n",
    "def WaterfallPlot(Data,name, gtitle,hours = 8760):\n",
    "    C = 8760/hours\n",
    "    if isinstance(Data,list):\n",
    "        y = [-1*Data[0]*C,-1*Data[1]*C,-1*Data[2]*C,-1*Data[3]*C,-1*Data[4]*C,Data[5]*C, ((Data[5])-((Data[0]+Data[1])+Data[2]+Data[3]+Data[4]))*C]\n",
    "    else:\n",
    "        y = [-1*A.rACC*C,-1*A.rFOMC*C,-1*A.rVC*C,-1*A.rFC*C,-1*A.rSUC*C,A.rEP*C, ((A.rEP)-((A.rACC+A.rFOMC)+A.rVC+A.rFC+A.rSUC))*C] # ,A.rAVC*C\n",
    "    \n",
    "    fig = go.Figure(go.Waterfall(\n",
    "        orientation = \"v\",\n",
    "        measure = [\"relative\", \"relative\", \"relative\", \"relative\", \"relative\", \"relative\", \"total\"],\n",
    "        x = [\"Annualized Capital Cost\", \"Fixed Cost\", \"Variable Cost\", \"Fuel Cost\", \"Start-Up Cost\", \"Electric Wholesale Revenues\", \"Net Revenue\"],\n",
    "        textposition = \"outside\",\n",
    "        y = y,\n",
    "        connector = {\"line\":{\"color\":\"rgb(63, 63, 63)\"}}\n",
    "    )) #\"TES Annualized Capital\",  , \"relative\",    , \"Natural Gas Cost\"\n",
    "\n",
    "    fig.update_layout(\n",
    "            title = gtitle,\n",
    "            showlegend = True,\n",
    "            width=400, height=400,\n",
    "            paper_bgcolor='rgba(0,0,0,0)',\n",
    "            plot_bgcolor='rgba(0,0,0,0)'\n",
    "    )\n",
    "\n",
    "    fig.write_html(name+\".html\")\n",
    "    fig.write_image(name+\".png\")\n",
    "\n",
    "\n",
    "def main(Temp = 300, NGCost = 4, hrcount = 8760, dem = 100, lmps = 'lmps.csv',sites = 'sites_3.csv', genInd = None, TES = 0, MaxMod = 2, naics = None, year = 2019, EP = 10, startHr = 0, DAM = False, DAM_data = None, init_run = False):\n",
    "    if isinstance(lmps,(str)):\n",
    "        S,L= ImportData(lmps = lmps, sites = sites)\n",
    "    else:\n",
    "        S,L = sites,lmps\n",
    "    Assessment = EconomicAssessment_OneGen()\n",
    "    heatPrice = NGTempCostCurve(Temp,NG_Cost = NGCost)\n",
    "    if genInd == None:\n",
    "        pass\n",
    "    else:\n",
    "        S = S.loc[S.index == genInd]\n",
    "    \n",
    "    if naics == None:\n",
    "        hrProfiles = [1]*8760\n",
    "    else:\n",
    "        hrProfiles = FP.QuickProfile(year, naics)\n",
    "    s = datetime.now()\n",
    "    # print(s,'Profiles complete')\n",
    "    ParamsVarsPD(Assessment,S,L,dem, TES, naics, EP, heatLMP = heatPrice, minTemp = Temp, MaxMod =MaxMod, hrProfiles = hrProfiles, startHr = startHr, hrcount=hrcount, DAM = DAM, DAM_data = DAM_data, init_run = init_run)\n",
    "    Assessment.SolveModel(DAM = DAM)\n",
    "    return Assessment\n",
    "\n",
    "def DAMmain(Temp = 300, NGCost = 4, dem = 100, lmps = 'SPP_SOUTH_2021.csv',sites = 'Base4_MSR.csv', genInd = None, TES = 0, InitMaxMod = 2, naics = 324110, year = 2019, EP = 10, startHr = 0, days = 1):\n",
    "    # do the rough version\n",
    "    S,L= ImportData(lmps = lmps, sites = sites)\n",
    "    L2 = L.iloc[::12,:].reset_index()\n",
    "    print('Running initial pass')\n",
    "    hr = len(L2.index)\n",
    "    A = main(Temp = Temp, NGCost = NGCost, hrcount = hr, dem = dem, lmps = L2,sites = S, genInd = genInd, TES = TES, MaxMod = InitMaxMod, naics = naics, year = year, EP = EP, startHr = startHr, DAM = False, DAM_data = None, init_run = True)\n",
    "    \n",
    "    MaxMod = A.ModCount\n",
    "    InitGen = [0]*MaxMod\n",
    "    InitStat = [0]*MaxMod\n",
    "    InitSOC = 0\n",
    "    Ch0 = 0\n",
    "    Dis0 = 0\n",
    "    Output = pd.DataFrame()\n",
    "    hP = []\n",
    "    hACC = []\n",
    "    hFOMC = []\n",
    "    hVC = []\n",
    "    hFC = []\n",
    "    hSUC = []\n",
    "    hrTES = []\n",
    "    hEP = []\n",
    "    F = []\n",
    "    if A.ModCount > 0:\n",
    "        try: \n",
    "            for x in range(days):\n",
    "                if x == days-1:\n",
    "                    hrc = 24\n",
    "                else:\n",
    "                    hrc = 48\n",
    "                print(\"____________________________{}___________________\".format(x))\n",
    "                A = main(Temp = Temp, NGCost = NGCost, hrcount = hrc, dem = dem, lmps = lmps,sites = sites, genInd = genInd, TES = TES, MaxMod = MaxMod, naics = naics, year = year, EP = EP, startHr = 24*x, DAM = True, DAM_data = [MaxMod,InitGen,InitStat, InitSOC, Ch0, Dis0], init_run = False)\n",
    "\n",
    "                InitGen = A.vRGen_PD.iloc[23,:].tolist()\n",
    "                InitStat = A.vOnoroff_PD.iloc[-1,:].tolist()\n",
    "                Output = pd.concat((Output,A.Output.iloc[:24,:]))\n",
    "                InitSOC = float(A.vTESSOC.sum(axis=1)[23])\n",
    "                Ch0 = float(A.vTESSOC_Charge.sum(axis = 1)[23])\n",
    "                Dis0 = float(A.vTESSOC_Discharge.sum(axis = 1)[23])\n",
    "\n",
    "                vRGen = A.vRGen.sum(axis = 1)\n",
    "                vGen = A.vGen.sum(axis = 1)\n",
    "                TO = A.vTurnon.sum(axis = 1)[:24]\n",
    "                EP_List = [A.pLMP[k]*A.vGen[k] for k in range(24)]\n",
    "\n",
    "                rVC = float(sum(vRGen[:24])*A.pVOC*A.pThermaleff)\n",
    "                rFC = float(sum(vRGen[:24])*A.FuelCost*A.pThermaleff)\n",
    "                rSUC = float(sum(TO)*A.pStartupfixedcost)\n",
    "                rEP = float(sum(sum(EP_List)))        \n",
    "\n",
    "                hACC.append(A.rACC*(24/hrc))\n",
    "                hFOMC.append(A.rFOMC*(24/hrc))\n",
    "                hVC.append(rVC)\n",
    "                hFC.append(rFC)\n",
    "                hSUC.append(rSUC)\n",
    "                hrTES.append(A.rTES*(24/hrc))\n",
    "                hEP.append(rEP)\n",
    "                hP.append(rEP-(A.rACC*(24/hrc)+A.rFOMC*(24/hrc)+A.rTES*(24/hrc)+rVC+rFC+rSUC))\n",
    "                F.append(0)\n",
    "        except:\n",
    "            print('Failure - DAM Market')\n",
    "            F.append(2)\n",
    "            hACC.append(0)\n",
    "            hFOMC.append(0)\n",
    "            hVC.append(0)\n",
    "            hFC.append(0)\n",
    "            hSUC.append(0)\n",
    "            hrTES.append(0)\n",
    "            hEP.append(0)\n",
    "            hP.append(0)\n",
    "    else:\n",
    "        F.append(1)\n",
    "        hACC.append(0)\n",
    "        hFOMC.append(0)\n",
    "        hVC.append(0)\n",
    "        hFC.append(0)\n",
    "        hSUC.append(0)\n",
    "        hrTES.append(0)\n",
    "        hEP.append(0)\n",
    "        hP.append(0)\n",
    "                \n",
    "    df = {'Profit':hP,\"Annualized Capital Cost\":hACC, \"Fixed Cost\":hFOMC, \"Variable Cost\":hVC, \"Fuel Cost\":hFC, \n",
    "          \"Start-Up Cost\":hSUC, \"TES Annualized Capital\":hrTES, \"Electric Wholesale Revenues\":hEP, \"Feasibility\":F}\n",
    "    Finances = pd.DataFrame(df)\n",
    "    \n",
    "    Output = Output.reset_index(drop= True)\n",
    "    \n",
    "    Fl = Finances.sum().tolist()\n",
    "    WF = [Fl[1],Fl[2],Fl[3],Fl[4],Fl[5],Fl[7]] # This is missing the TES portion right now with Fl[6]\n",
    "    \n",
    "    return Output, L, Finances, WF\n",
    "\n",
    "\n",
    "A,L,F,WF = DAMmain(dem = 327, days = 3, InitMaxMod = 2, TES= 0)\n",
    "QuickGraphOutputs(A,L,'DAM_test1',hr = 72, TES = True)\n",
    "#WaterfallPlot(WF,'TestDAM', 'SPP 2021, 327 MWht, 300 C, DAM')\n",
    "e = datetime.now()\n",
    "print(e,'end')\n",
    "d = e-s\n",
    "print(d)\n",
    "print(F.sum())\n",
    "s = datetime.now()\n",
    "print(s,'start')\n",
    "#print(A)\n",
    "'''\n",
    "A = main(hrcount = int(87),sites = 'Base4_MSR.csv',lmps = 'SPP_SOUTH_2021.csv',NGCost = 8, MaxMod=2, TES = 200, naics = 324110, EP=10, startHr = 0, dem = 327)\n",
    "e = datetime.now()\n",
    "print(e,'end')\n",
    "d = e-s\n",
    "print(d)\n",
    "QuickGraphOutputs(A.Output,L,'DAM_test2',hr = 72, TES = True)\n",
    "WaterfallPlot(A,'TestFull', 'SPP 2021, 327 MWht, 300 C, Full year')\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'SPP_NG_0818.csv'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-5-57e44c0280b6>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m    141\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mpdfBestPer\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    142\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 143\u001b[1;33m \u001b[0mSPP_Facilities\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mread_csv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'SPP_NG_0818.csv'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    144\u001b[0m \u001b[0mSPP_Facilities\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mSPP_Facilities\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mSPP_Facilities\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'Temp_degC'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m<=\u001b[0m\u001b[1;36m800\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    145\u001b[0m \u001b[0mSPP_Facilities_Reduced\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mSPP_Facilities\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msample\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m5\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreplace\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrandom_state\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m21\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\envs\\NEUP\\lib\\site-packages\\pandas\\io\\parsers.py\u001b[0m in \u001b[0;36mread_csv\u001b[1;34m(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, squeeze, prefix, mangle_dupe_cols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, dialect, error_bad_lines, warn_bad_lines, delim_whitespace, low_memory, memory_map, float_precision)\u001b[0m\n\u001b[0;32m    684\u001b[0m     )\n\u001b[0;32m    685\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 686\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0m_read\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfilepath_or_buffer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkwds\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    687\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    688\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\envs\\NEUP\\lib\\site-packages\\pandas\\io\\parsers.py\u001b[0m in \u001b[0;36m_read\u001b[1;34m(filepath_or_buffer, kwds)\u001b[0m\n\u001b[0;32m    450\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    451\u001b[0m     \u001b[1;31m# Create the parser.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 452\u001b[1;33m     \u001b[0mparser\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mTextFileReader\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfp_or_buf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    453\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    454\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mchunksize\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0miterator\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\envs\\NEUP\\lib\\site-packages\\pandas\\io\\parsers.py\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, f, engine, **kwds)\u001b[0m\n\u001b[0;32m    944\u001b[0m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0moptions\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"has_index_names\"\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mkwds\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"has_index_names\"\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    945\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 946\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_make_engine\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mengine\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    947\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    948\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\envs\\NEUP\\lib\\site-packages\\pandas\\io\\parsers.py\u001b[0m in \u001b[0;36m_make_engine\u001b[1;34m(self, engine)\u001b[0m\n\u001b[0;32m   1176\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m_make_engine\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mengine\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m\"c\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1177\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mengine\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m\"c\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1178\u001b[1;33m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_engine\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mCParserWrapper\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0moptions\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1179\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1180\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mengine\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m\"python\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\envs\\NEUP\\lib\\site-packages\\pandas\\io\\parsers.py\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, src, **kwds)\u001b[0m\n\u001b[0;32m   2006\u001b[0m         \u001b[0mkwds\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"usecols\"\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0musecols\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2007\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 2008\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_reader\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mparsers\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mTextReader\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   2009\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0munnamed_cols\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_reader\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0munnamed_cols\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2010\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mpandas\\_libs\\parsers.pyx\u001b[0m in \u001b[0;36mpandas._libs.parsers.TextReader.__cinit__\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mpandas\\_libs\\parsers.pyx\u001b[0m in \u001b[0;36mpandas._libs.parsers.TextReader._setup_parser_source\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'SPP_NG_0818.csv'"
     ]
    }
   ],
   "source": [
    "def PerFacilityRuns(facilities,region,year,outDir, TESval = 0, MaxMod = 2, startHr = 0,days = 365):\n",
    "\n",
    "    try: \n",
    "        os.mkdir(outDir) \n",
    "    except OSError as error: \n",
    "        print('Directory Already Exists, Overwriting')\n",
    "\n",
    "    if region =='ERCOT':\n",
    "        if year == 2020:\n",
    "            LMPs = ['ERCOT_2020_HubAVG.csv']\n",
    "        if year == 2021:\n",
    "            print('Will fail until new data is entered')\n",
    "            LMPs = ['ERCOT_2021_HubAVG.csv']\n",
    "        if year == 'all':\n",
    "            print('Will fail until new data is entered')\n",
    "            LMPs = ['ERCOT_2020_HubAVG.csv','ERCOT_2021_HubAVG.csv']\n",
    "\n",
    "    if region =='SPP':\n",
    "        if year == 2020:\n",
    "            print('Will fail until new data is entered')\n",
    "            LMPs = ['SPP_SOUTH_2020.csv']\n",
    "        if year == 2021:\n",
    "            LMPs = ['SPP_SOUTH_2021.csv']\n",
    "        if year == 'all':\n",
    "            print('Will fail until new data is entered')\n",
    "            LMPs = ['ERCOT_2020_HubAVG.csv','ERCOT_2021_HubAVG.csv']\n",
    "\n",
    "    pdf2 = pd.DataFrame()\n",
    "    pdfBestPer = pd.DataFrame()\n",
    "    for LMP_Val in LMPs:\n",
    "        T = facilities['Temp_degC'].tolist()\n",
    "        \n",
    "        facID = facilities['FACILITY_ID'].tolist()\n",
    "\n",
    "        Ds = facilities['Thermal MWh/hr'].tolist()\n",
    "        \n",
    "        NAICS = facilities['FINAL_NAICS_CODE'].tolist()\n",
    "\n",
    "        S,L = ImportData(lmps = LMP_Val, sites = 'Base4.csv')\n",
    "        j = 0\n",
    "        \n",
    "        while j < len(T):\n",
    "            print(\"####### \",str(j),\"of \",str(len(T)),' ##########')\n",
    "            G = []\n",
    "            P = []\n",
    "            LMP_List = []\n",
    "            D_list = []\n",
    "            T_list = []\n",
    "            F_list = []\n",
    "            ACC = []\n",
    "            FOMC = []\n",
    "            VC = []\n",
    "            FC = []\n",
    "            SUC = []\n",
    "            rTES = []\n",
    "            EP = []\n",
    "            AVC = []\n",
    "            t = T[j]\n",
    "            D = Ds[j]\n",
    "            NAICS_j = NAICS[j]\n",
    "            N = []\n",
    "            DEMTYPE = []\n",
    "            GI = 0\n",
    "            while GI < len(S.index):\n",
    "                print('Checkpoint 2')\n",
    "                D_list.append(D)\n",
    "                T_list.append(t)\n",
    "                F_list.append(facID[j])\n",
    "                A,L,F,WF = DAMmain(days = days, sites = 'Base5.csv', lmps = LMP_Val,dem = D, genInd = GI, Temp = t, TES = TESval, InitMaxMod = MaxMod, startHr = 0, year = year, naics = NAICS[j])\n",
    "                LMP_List.append(LMP_Val)\n",
    "                G.append(S.loc[S.index==GI]['Sites'].tolist()[0])\n",
    "                print(WF)\n",
    "                TotalP = WF[5]-(WF[0]+WF[1]+WF[2]+WF[3]+WF[4])\n",
    "                P.append(TotalP)\n",
    "                if TotalP != None:\n",
    "                    A.to_csv(outDir+'/'+str(facID[j])+'_Op_'+str(S.loc[S.index==GI]['Sites'].tolist()[0])+'.csv') \n",
    "                    F.to_csv(outDir+'/'+str(facID[j])+'_Fin_'+str(S.loc[S.index==GI]['Sites'].tolist()[0])+'.csv') \n",
    "                    N.append((len(A.columns.tolist())-2)/5)\n",
    "                else:\n",
    "                    N.append(0)\n",
    "                \n",
    "                ACC.append(-1*WF[0])\n",
    "                FOMC.append(-1*WF[1])\n",
    "                VC.append(-1*WF[2])\n",
    "                FC.append(-1*WF[3])\n",
    "                SUC.append(-1*WF[4])\n",
    "                #rTES.append(-1*WF[5])\n",
    "                EP.append(WF[5])\n",
    "\n",
    "                DEMTYPE.append('hourly')\n",
    "                \n",
    "                ## THIS IS SOLELY FOR THE TEST OF THE INDUSTRIAL DEMAND HOURLY CURVES\n",
    "                '''\n",
    "                A,L,F,WF = DAMmain(days = days, sites = 'Base5.csv', lmps = LMP_Val,dem = D, genInd = GI, Temp = t, TES = TESval, InitMaxMod = MaxMod, startHr = 0, year = year, naics = None)\n",
    "                LMP_List.append(LMP_Val)\n",
    "                D_list.append(D)\n",
    "                T_list.append(t)\n",
    "                F_list.append(facID[j])\n",
    "                G.append(S.loc[S.index==GI]['Sites'].tolist()[0])\n",
    "                print(WF)\n",
    "                TotalP = WF[5]-(WF[0]+WF[1]+WF[2]+WF[3]+WF[4])\n",
    "                P.append(TotalP)\n",
    "                if TotalP != None:\n",
    "                    A.to_csv(outDir+'/'+str(facID[j])+'_Op_'+str(S.loc[S.index==GI]['Sites'].tolist()[0])+'.csv') \n",
    "                    F.to_csv(outDir+'/'+str(facID[j])+'_Fin_'+str(S.loc[S.index==GI]['Sites'].tolist()[0])+'.csv') \n",
    "                    N.append((len(A.columns.tolist())-2)/5)\n",
    "                else:\n",
    "                    N.append(0)\n",
    "                ACC.append(-1*WF[0])\n",
    "                FOMC.append(-1*WF[1])\n",
    "                VC.append(-1*WF[2])\n",
    "                FC.append(-1*WF[3])\n",
    "                SUC.append(-1*WF[4])\n",
    "                #rTES.append(-1*WF[5])\n",
    "                EP.append(WF[5])\n",
    "                \n",
    "                DEMTYPE.append('flat')\n",
    "                '''\n",
    "                \n",
    "                GI+=1 ##############################3\n",
    "                \n",
    "            # For the flat versus hourly\n",
    "            \n",
    "                \n",
    "            df = {'LMPType':LMP_List,'Facility ID':F_list,'Temperature Req':T_list,'Generator':G,'Modules':N,\n",
    "                  'Facility Demand':D_list,'Profit':P,\"Annualized Capital Cost\":ACC, \n",
    "                  \"Fixed Cost\":FOMC, \"Variable Cost\":VC, \"Fuel Cost\":FC, \"Start-Up Cost\":SUC, \n",
    "                  \"Electric Wholesale Revenues\":EP,'Demand type':DEMTYPE} #  \"Natural Gas Cost\":AVC}\"TES Annualized Capital\":rTES, \n",
    "\n",
    "            print(df)\n",
    "            j+=1\n",
    "            pdf = pd.DataFrame(df)\n",
    "            \n",
    "            prof_best = pdf['Profit'].max()\n",
    "            \n",
    "            pdfBestPer = pd.concat((pdfBestPer,pdf.loc[pdf['Profit']==prof_best]))\n",
    "            \n",
    "            pdf2 = pd.concat((pdf2,pdf))\n",
    "    pdfBestPer.to_csv(outDir+'/BestPer.csv')   \n",
    "    pdf2.to_csv(outDir+'/PerFacilityTest.csv')\n",
    "    return pdfBestPer\n",
    "\n",
    "SPP_Facilities = pd.read_csv('SPP_NG_0818.csv')\n",
    "SPP_Facilities = SPP_Facilities.loc[SPP_Facilities['Temp_degC']<=800]\n",
    "SPP_Facilities_Reduced = SPP_Facilities.sample(n=5, replace=False, random_state=21)\n",
    "SPP_Facilities_Reduced = SPP_Facilities_Reduced.sample(n=1, replace=False, random_state=1)\n",
    "SPP_Facilities_Reduced = pd.concat((SPP_Facilities_Reduced,SPP_Facilities.sample(n=3, replace=False, random_state=8))).reset_index()\n",
    "\n",
    "\n",
    "SPP_628 = PerFacilityRuns(SPP_Facilities_Reduced,'SPP',2021,'SPP_NG_Hourly830', TESval = 0, MaxMod = 12, days = 365)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>index</th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>COUNTY</th>\n",
       "      <th>COUNTY_FIPS</th>\n",
       "      <th>FACILITY_ID</th>\n",
       "      <th>FINAL_NAICS_CODE</th>\n",
       "      <th>FUEL_TYPE</th>\n",
       "      <th>MECS_NAICS</th>\n",
       "      <th>REPORTING_YEAR</th>\n",
       "      <th>STATE</th>\n",
       "      <th>Temp_degC</th>\n",
       "      <th>Industry</th>\n",
       "      <th>Label</th>\n",
       "      <th>Natural_gas</th>\n",
       "      <th>Total</th>\n",
       "      <th>MMTCO2E</th>\n",
       "      <th>Thermal MWh</th>\n",
       "      <th>Thermal MWh/hr</th>\n",
       "      <th>Fuel_group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>43</td>\n",
       "      <td>0</td>\n",
       "      <td>WASHINGTON</td>\n",
       "      <td>31177</td>\n",
       "      <td>1006425</td>\n",
       "      <td>311221</td>\n",
       "      <td>Natural Gas (Weighted U.S. Average)</td>\n",
       "      <td>311221</td>\n",
       "      <td>2015</td>\n",
       "      <td>NE</td>\n",
       "      <td>149</td>\n",
       "      <td>Wet Corn Milling</td>\n",
       "      <td>Wet Corn Milling_1006425</td>\n",
       "      <td>7147.487371</td>\n",
       "      <td>7147.487371</td>\n",
       "      <td>0.359827</td>\n",
       "      <td>1.985415e+06</td>\n",
       "      <td>226.645519</td>\n",
       "      <td>NG_Type</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>144</td>\n",
       "      <td>0</td>\n",
       "      <td>MERRICK</td>\n",
       "      <td>31121</td>\n",
       "      <td>1000344</td>\n",
       "      <td>325193</td>\n",
       "      <td>Other Biomass Gases</td>\n",
       "      <td>325193</td>\n",
       "      <td>2015</td>\n",
       "      <td>NE</td>\n",
       "      <td>450</td>\n",
       "      <td>Ethyl Alcohol</td>\n",
       "      <td>Ethyl Alcohol_1000344</td>\n",
       "      <td>2404.614616</td>\n",
       "      <td>2413.811169</td>\n",
       "      <td>0.121512</td>\n",
       "      <td>6.705036e+05</td>\n",
       "      <td>76.541511</td>\n",
       "      <td>NG_Type</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>KAY</td>\n",
       "      <td>40071</td>\n",
       "      <td>1001920</td>\n",
       "      <td>324110</td>\n",
       "      <td>Natural Gas (Weighted U.S. Average)</td>\n",
       "      <td>324110</td>\n",
       "      <td>2015</td>\n",
       "      <td>OK</td>\n",
       "      <td>600</td>\n",
       "      <td>Petroleum Refineries</td>\n",
       "      <td>Petroleum Refineries_1001920</td>\n",
       "      <td>497.482221</td>\n",
       "      <td>13373.313726</td>\n",
       "      <td>0.748174</td>\n",
       "      <td>3.714812e+06</td>\n",
       "      <td>424.065336</td>\n",
       "      <td>NG_Type</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>25</td>\n",
       "      <td>0</td>\n",
       "      <td>RENO COUNTY</td>\n",
       "      <td>20155</td>\n",
       "      <td>1004961</td>\n",
       "      <td>322130</td>\n",
       "      <td>Natural Gas (Weighted U.S. Average)</td>\n",
       "      <td>322130</td>\n",
       "      <td>2015</td>\n",
       "      <td>KS</td>\n",
       "      <td>275</td>\n",
       "      <td>Paperboard Mills</td>\n",
       "      <td>Paperboard Mills_1004961</td>\n",
       "      <td>16.272375</td>\n",
       "      <td>16.272375</td>\n",
       "      <td>0.000819</td>\n",
       "      <td>4.520108e+03</td>\n",
       "      <td>0.515994</td>\n",
       "      <td>NG_Type</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   index  Unnamed: 0       COUNTY  COUNTY_FIPS  FACILITY_ID  FINAL_NAICS_CODE  \\\n",
       "0     43           0   WASHINGTON        31177      1006425            311221   \n",
       "1    144           0      MERRICK        31121      1000344            325193   \n",
       "2      7           0          KAY        40071      1001920            324110   \n",
       "3     25           0  RENO COUNTY        20155      1004961            322130   \n",
       "\n",
       "                             FUEL_TYPE  MECS_NAICS  REPORTING_YEAR STATE  \\\n",
       "0  Natural Gas (Weighted U.S. Average)      311221            2015    NE   \n",
       "1                  Other Biomass Gases      325193            2015    NE   \n",
       "2  Natural Gas (Weighted U.S. Average)      324110            2015    OK   \n",
       "3  Natural Gas (Weighted U.S. Average)      322130            2015    KS   \n",
       "\n",
       "   Temp_degC              Industry                         Label  Natural_gas  \\\n",
       "0        149      Wet Corn Milling      Wet Corn Milling_1006425  7147.487371   \n",
       "1        450         Ethyl Alcohol         Ethyl Alcohol_1000344  2404.614616   \n",
       "2        600  Petroleum Refineries  Petroleum Refineries_1001920   497.482221   \n",
       "3        275      Paperboard Mills      Paperboard Mills_1004961    16.272375   \n",
       "\n",
       "          Total   MMTCO2E   Thermal MWh  Thermal MWh/hr Fuel_group  \n",
       "0   7147.487371  0.359827  1.985415e+06      226.645519    NG_Type  \n",
       "1   2413.811169  0.121512  6.705036e+05       76.541511    NG_Type  \n",
       "2  13373.313726  0.748174  3.714812e+06      424.065336    NG_Type  \n",
       "3     16.272375  0.000819  4.520108e+03        0.515994    NG_Type  "
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SPP_Facilities = pd.read_csv('SPP_NG_0818.csv')\n",
    "SPP_Facilities = SPP_Facilities.loc[SPP_Facilities['Temp_degC']<=800]\n",
    "SPP_Facilities_Reduced = SPP_Facilities.sample(n=5, replace=False, random_state=21)\n",
    "SPP_Facilities_Reduced = SPP_Facilities_Reduced.sample(n=1, replace=False, random_state=1)\n",
    "SPP_Facilities_Reduced = pd.concat((SPP_Facilities_Reduced,SPP_Facilities.sample(n=3, replace=False, random_state=8))).reset_index()\n",
    "SPP_Facilities_Reduced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>LMPType</th>\n",
       "      <th>Facility ID</th>\n",
       "      <th>Temperature Req</th>\n",
       "      <th>Generator</th>\n",
       "      <th>Modules</th>\n",
       "      <th>Facility Demand</th>\n",
       "      <th>Profit</th>\n",
       "      <th>Annualized Capital Cost</th>\n",
       "      <th>Fixed Cost</th>\n",
       "      <th>Variable Cost</th>\n",
       "      <th>Fuel Cost</th>\n",
       "      <th>Start-Up Cost</th>\n",
       "      <th>TES Annualized Capital</th>\n",
       "      <th>Electric Wholesale Revenues</th>\n",
       "      <th>NatGCost_multiplier</th>\n",
       "      <th>NG_100Cost_1USD</th>\n",
       "      <th>NGBE_Cost</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2</td>\n",
       "      <td>SPP_SOUTH_2021.csv</td>\n",
       "      <td>1001920</td>\n",
       "      <td>220</td>\n",
       "      <td>iMSR</td>\n",
       "      <td>6.0</td>\n",
       "      <td>1621.824519</td>\n",
       "      <td>-4.935116e+06</td>\n",
       "      <td>-3.183886e+06</td>\n",
       "      <td>-820890.410959</td>\n",
       "      <td>-41930.382589</td>\n",
       "      <td>-887414.617108</td>\n",
       "      <td>-180000.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>179043.712294</td>\n",
       "      <td>4.211594</td>\n",
       "      <td>683046.667613</td>\n",
       "      <td>7.225152</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>SPP_SOUTH_2021.csv</td>\n",
       "      <td>1004291</td>\n",
       "      <td>600</td>\n",
       "      <td>iMSR</td>\n",
       "      <td>1.0</td>\n",
       "      <td>287.375556</td>\n",
       "      <td>-8.468796e+05</td>\n",
       "      <td>-5.306477e+05</td>\n",
       "      <td>-136815.068493</td>\n",
       "      <td>-6741.415445</td>\n",
       "      <td>-142675.316473</td>\n",
       "      <td>-30000.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>5.109159</td>\n",
       "      <td>146824.747899</td>\n",
       "      <td>5.767962</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>SPP_SOUTH_2021.csv</td>\n",
       "      <td>1004381</td>\n",
       "      <td>220</td>\n",
       "      <td>iMSR</td>\n",
       "      <td>1.0</td>\n",
       "      <td>295.909093</td>\n",
       "      <td>-8.513165e+05</td>\n",
       "      <td>-5.306477e+05</td>\n",
       "      <td>-136815.068493</td>\n",
       "      <td>-6941.599884</td>\n",
       "      <td>-146912.019938</td>\n",
       "      <td>-30000.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4.211594</td>\n",
       "      <td>124624.900671</td>\n",
       "      <td>6.831030</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2</td>\n",
       "      <td>SPP_SOUTH_2021.csv</td>\n",
       "      <td>1004381</td>\n",
       "      <td>600</td>\n",
       "      <td>iMSR</td>\n",
       "      <td>6.0</td>\n",
       "      <td>1763.173827</td>\n",
       "      <td>-5.081633e+06</td>\n",
       "      <td>-3.183886e+06</td>\n",
       "      <td>-820890.410959</td>\n",
       "      <td>-42193.099600</td>\n",
       "      <td>-892974.759940</td>\n",
       "      <td>-180000.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>38314.517518</td>\n",
       "      <td>5.109159</td>\n",
       "      <td>900833.585376</td>\n",
       "      <td>5.641034</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>SPP_SOUTH_2021.csv</td>\n",
       "      <td>1007942</td>\n",
       "      <td>220</td>\n",
       "      <td>iMSR</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1439.181765</td>\n",
       "      <td>-4.235596e+06</td>\n",
       "      <td>-2.653238e+06</td>\n",
       "      <td>-684075.342466</td>\n",
       "      <td>-33761.125363</td>\n",
       "      <td>-714520.457180</td>\n",
       "      <td>-150000.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4.211594</td>\n",
       "      <td>606124.951796</td>\n",
       "      <td>6.987991</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>102</th>\n",
       "      <td>2</td>\n",
       "      <td>SPP_SOUTH_2021.csv</td>\n",
       "      <td>1002865</td>\n",
       "      <td>450</td>\n",
       "      <td>iMSR</td>\n",
       "      <td>1.0</td>\n",
       "      <td>275.890610</td>\n",
       "      <td>-8.465580e+05</td>\n",
       "      <td>-5.306477e+05</td>\n",
       "      <td>-136815.068493</td>\n",
       "      <td>-6726.899532</td>\n",
       "      <td>-142368.101704</td>\n",
       "      <td>-30000.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4.712701</td>\n",
       "      <td>130019.006150</td>\n",
       "      <td>6.511033</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>103</th>\n",
       "      <td>0</td>\n",
       "      <td>SPP_SOUTH_2021.csv</td>\n",
       "      <td>1002865</td>\n",
       "      <td>233</td>\n",
       "      <td>iPWR</td>\n",
       "      <td>1.0</td>\n",
       "      <td>125.250952</td>\n",
       "      <td>-5.516958e+05</td>\n",
       "      <td>-3.920736e+05</td>\n",
       "      <td>-101084.474886</td>\n",
       "      <td>-6009.555981</td>\n",
       "      <td>-67391.911969</td>\n",
       "      <td>-30000.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>44881.058271</td>\n",
       "      <td>4.237059</td>\n",
       "      <td>53069.567087</td>\n",
       "      <td>10.395709</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>104</th>\n",
       "      <td>2</td>\n",
       "      <td>SPP_SOUTH_2021.csv</td>\n",
       "      <td>1005812</td>\n",
       "      <td>450</td>\n",
       "      <td>iMSR</td>\n",
       "      <td>1.0</td>\n",
       "      <td>255.112358</td>\n",
       "      <td>-8.353291e+05</td>\n",
       "      <td>-5.306477e+05</td>\n",
       "      <td>-136815.068493</td>\n",
       "      <td>-6220.274051</td>\n",
       "      <td>-131645.880008</td>\n",
       "      <td>-30000.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4.712701</td>\n",
       "      <td>120226.836471</td>\n",
       "      <td>6.947942</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>105</th>\n",
       "      <td>2</td>\n",
       "      <td>SPP_SOUTH_2021.csv</td>\n",
       "      <td>1005683</td>\n",
       "      <td>450</td>\n",
       "      <td>iMSR</td>\n",
       "      <td>1.0</td>\n",
       "      <td>229.260758</td>\n",
       "      <td>-7.900902e+05</td>\n",
       "      <td>-5.306477e+05</td>\n",
       "      <td>-136815.068493</td>\n",
       "      <td>-6889.084475</td>\n",
       "      <td>-145800.583828</td>\n",
       "      <td>-30000.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>60096.231025</td>\n",
       "      <td>4.712701</td>\n",
       "      <td>108043.749206</td>\n",
       "      <td>7.312688</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>106</th>\n",
       "      <td>2</td>\n",
       "      <td>SPP_SOUTH_2021.csv</td>\n",
       "      <td>1003954</td>\n",
       "      <td>450</td>\n",
       "      <td>iMSR</td>\n",
       "      <td>1.0</td>\n",
       "      <td>199.467109</td>\n",
       "      <td>-8.052576e+05</td>\n",
       "      <td>-5.306477e+05</td>\n",
       "      <td>-136815.068493</td>\n",
       "      <td>-4863.504435</td>\n",
       "      <td>-102931.207867</td>\n",
       "      <td>-30000.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4.712701</td>\n",
       "      <td>94002.892421</td>\n",
       "      <td>8.566306</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>107 rows × 18 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Unnamed: 0             LMPType  Facility ID  Temperature Req Generator  \\\n",
       "0             2  SPP_SOUTH_2021.csv      1001920              220      iMSR   \n",
       "1             2  SPP_SOUTH_2021.csv      1004291              600      iMSR   \n",
       "2             2  SPP_SOUTH_2021.csv      1004381              220      iMSR   \n",
       "3             2  SPP_SOUTH_2021.csv      1004381              600      iMSR   \n",
       "4             2  SPP_SOUTH_2021.csv      1007942              220      iMSR   \n",
       "..          ...                 ...          ...              ...       ...   \n",
       "102           2  SPP_SOUTH_2021.csv      1002865              450      iMSR   \n",
       "103           0  SPP_SOUTH_2021.csv      1002865              233      iPWR   \n",
       "104           2  SPP_SOUTH_2021.csv      1005812              450      iMSR   \n",
       "105           2  SPP_SOUTH_2021.csv      1005683              450      iMSR   \n",
       "106           2  SPP_SOUTH_2021.csv      1003954              450      iMSR   \n",
       "\n",
       "     Modules  Facility Demand        Profit  Annualized Capital Cost  \\\n",
       "0        6.0      1621.824519 -4.935116e+06            -3.183886e+06   \n",
       "1        1.0       287.375556 -8.468796e+05            -5.306477e+05   \n",
       "2        1.0       295.909093 -8.513165e+05            -5.306477e+05   \n",
       "3        6.0      1763.173827 -5.081633e+06            -3.183886e+06   \n",
       "4        5.0      1439.181765 -4.235596e+06            -2.653238e+06   \n",
       "..       ...              ...           ...                      ...   \n",
       "102      1.0       275.890610 -8.465580e+05            -5.306477e+05   \n",
       "103      1.0       125.250952 -5.516958e+05            -3.920736e+05   \n",
       "104      1.0       255.112358 -8.353291e+05            -5.306477e+05   \n",
       "105      1.0       229.260758 -7.900902e+05            -5.306477e+05   \n",
       "106      1.0       199.467109 -8.052576e+05            -5.306477e+05   \n",
       "\n",
       "        Fixed Cost  Variable Cost      Fuel Cost  Start-Up Cost  \\\n",
       "0   -820890.410959  -41930.382589 -887414.617108      -180000.0   \n",
       "1   -136815.068493   -6741.415445 -142675.316473       -30000.0   \n",
       "2   -136815.068493   -6941.599884 -146912.019938       -30000.0   \n",
       "3   -820890.410959  -42193.099600 -892974.759940      -180000.0   \n",
       "4   -684075.342466  -33761.125363 -714520.457180      -150000.0   \n",
       "..             ...            ...            ...            ...   \n",
       "102 -136815.068493   -6726.899532 -142368.101704       -30000.0   \n",
       "103 -101084.474886   -6009.555981  -67391.911969       -30000.0   \n",
       "104 -136815.068493   -6220.274051 -131645.880008       -30000.0   \n",
       "105 -136815.068493   -6889.084475 -145800.583828       -30000.0   \n",
       "106 -136815.068493   -4863.504435 -102931.207867       -30000.0   \n",
       "\n",
       "     TES Annualized Capital  Electric Wholesale Revenues  NatGCost_multiplier  \\\n",
       "0                      -0.0                179043.712294             4.211594   \n",
       "1                      -0.0                     0.000000             5.109159   \n",
       "2                      -0.0                     0.000000             4.211594   \n",
       "3                      -0.0                 38314.517518             5.109159   \n",
       "4                      -0.0                     0.000000             4.211594   \n",
       "..                      ...                          ...                  ...   \n",
       "102                    -0.0                     0.000000             4.712701   \n",
       "103                    -0.0                 44881.058271             4.237059   \n",
       "104                    -0.0                     0.000000             4.712701   \n",
       "105                    -0.0                 60096.231025             4.712701   \n",
       "106                    -0.0                     0.000000             4.712701   \n",
       "\n",
       "     NG_100Cost_1USD  NGBE_Cost  \n",
       "0      683046.667613   7.225152  \n",
       "1      146824.747899   5.767962  \n",
       "2      124624.900671   6.831030  \n",
       "3      900833.585376   5.641034  \n",
       "4      606124.951796   6.987991  \n",
       "..               ...        ...  \n",
       "102    130019.006150   6.511033  \n",
       "103     53069.567087  10.395709  \n",
       "104    120226.836471   6.947942  \n",
       "105    108043.749206   7.312688  \n",
       "106     94002.892421   8.566306  \n",
       "\n",
       "[107 rows x 18 columns]"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SPP100 = pd.read_csv('SPP_627_6/BestPer.csv')\n",
    "SPP100\n",
    "\n",
    "SPP100['NatGCost_multiplier'] = NGTempCostCurve(SPP100['Temperature Req'],NG_Cost = 1.0, AHF_Coeffs = [0,-0.00038,0.90556]) # find the $/MWh-th of Natural gas from the $/mmbtu so 1 $/mmbtu = This output $/MWhth\n",
    "SPP100['NG_100Cost_1USD'] = SPP100['NatGCost_multiplier']*SPP100['Facility Demand']*100\n",
    "SPP100['NGBE_Cost'] = -1*SPP100['Profit']/SPP100['NG_100Cost_1USD']\n",
    "SPP100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Type = ['iPWR','iMSR','HTGR','Microreactor']\n",
    "Size = [1000,500,100,50,10]\n",
    "SupCurve = pd.DataFrame()\n",
    "for s in Size:\n",
    "    sub_df_s = SPP100.loc[SPP100['Facility Demand']<=s]\n",
    "    SupCurve_s = pd.DataFrame()\n",
    "    for t in Type:\n",
    "        subdf_t = sub_df_s.loc[sub_df_s['Generator']==t]\n",
    "        NGBEs = subdf_t['NGBE_Cost'].unique()\n",
    "        SupCurve_t = pd.DataFrame()\n",
    "        DemC = []\n",
    "        NumC = []\n",
    "        for i in NGBEs:\n",
    "            subdf = subdf_t.loc[subdf_t['NGBE_Cost']<= i]\n",
    "            DemC.append(subdf['Facility Demand'].sum())\n",
    "            NumC.append(len(subdf.index))\n",
    "        SupCurve_t['Natural Gas Cost'] = NGBEs\n",
    "        SupCurve_t['Facility Demand'] = DemC\n",
    "        SupCurve_t['Facility Count'] = NumC\n",
    "        SupCurve_t['Reactor Type'] = t\n",
    "        SupCurve_t['Demand Size'] = s\n",
    "        SupCurve = pd.concat((SupCurve,SupCurve_t))\n",
    "    SupCurve = pd.concat((SupCurve_s,SupCurve))\n",
    "SupCurve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc = {'figure.figsize':(15,3)})\n",
    "sns.set_style(\"whitegrid\")\n",
    "sc = sns.lineplot(SupCurve['Natural Gas Cost'],SupCurve['Facility Demand'], drawstyle='steps-pre',style = SupCurve['Reactor Type'], hue = SupCurve['Demand Size'])#, hue = All_models50['Region'], style = All_models50['Region'])\n",
    "sc.set(ylabel='Facility Demand [MWh-th/hr]')\n",
    "sc.set(xlabel='Natural Gas Cost [$/mmbtu]')\n",
    "sc.set(xlim=(5,40))\n",
    "plt.show()\n",
    "#sc.figure.savefig('exampleSupplyCurve.png', dpi = 700)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc = {'figure.figsize':(15,3)})\n",
    "sns.set_style(\"whitegrid\")\n",
    "sc = sns.lineplot(x = SupCurve['Natural Gas Cost'],y = SupCurve['Facility Count'], drawstyle='steps-pre',style = SupCurve['Reactor Type'])#, hue = All_models50['Region'], style = All_models50['Region'])\n",
    "sc.set(ylabel='Facility Count')\n",
    "sc.set(xlim=(0,40))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SPP10 = pd.read_csv('SPP_627_6/PerFacilityTest.csv')\n",
    "#SPP10 = pd.read_csv('BestPer.csv')\n",
    "#print(SPP10)\n",
    "\n",
    "SPP10['NatGCost_multiplier'] = NGTempCostCurve(SPP10['Temperature Req'],NG_Cost = 1.0, AHF_Coeffs = [0,-0.00038,0.90556]) # find the $/MWh-th of Natural gas from the $/mmbtu so 1 $/mmbtu = This output $/MWhth\n",
    "SPP10['NG_100Cost_1USD'] = SPP10['NatGCost_multiplier']*SPP10['Facility Demand']*100\n",
    "SPP10['NGBE_Cost'] = -1*SPP10['Profit']/SPP10['NG_100Cost_1USD']\n",
    "SPP10\n",
    "\n",
    "Type = ['iPWR','iMSR','HTGR','Microreactor']\n",
    "SupCurve = pd.DataFrame()\n",
    "for t in Type:\n",
    "    subdf_t = SPP10.loc[SPP10['Generator']==t]\n",
    "    NGBEs = subdf_t['NGBE_Cost'].unique()\n",
    "    SupCurve_t = pd.DataFrame()\n",
    "    DemC = []\n",
    "    NumC = []\n",
    "    for i in NGBEs:\n",
    "        subdf = subdf_t.loc[subdf_t['NGBE_Cost']<= i]\n",
    "        DemC.append(subdf['Facility Demand'].sum())\n",
    "        NumC.append(len(subdf.index))\n",
    "    SupCurve_t['Natural Gas Cost'] = NGBEs\n",
    "    SupCurve_t['Facility Demand'] = DemC\n",
    "    SupCurve_t['Facility Count'] = NumC\n",
    "    SupCurve_t['Reactor Type'] = t\n",
    "    SupCurve = pd.concat((SupCurve,SupCurve_t))\n",
    "SupCurve_t\n",
    "\n",
    "sns.set(rc = {'figure.figsize':(15,3)})\n",
    "sns.set_style(\"whitegrid\")\n",
    "sc = sns.lineplot(x = SupCurve['Natural Gas Cost'],y = SupCurve['Facility Demand'], drawstyle='steps-post',style = SupCurve['Reactor Type'])#, hue = All_models50['Region'], style = All_models50['Region'])\n",
    "sc.set(ylabel='Facility Demand [MWh-th/hr]')\n",
    "sc.set(xlabel='Natural Gas Cost [$/mmbtu]')\n",
    "sc.set(xlim=(0,50))\n",
    "plt.show()\n",
    "plt.tight_layout()\n",
    "#sc.figure.savefig('exampleSupplyCurve_best.png', dpi = 700, bbox_inches = \"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc = {'figure.figsize':(15,3)})\n",
    "sns.set_style(\"whitegrid\")\n",
    "sc = sns.lineplot(x = SupCurve['Natural Gas Cost'],y = SupCurve['Facility Count'], drawstyle='steps-post',style = SupCurve['Reactor Type'])#, hue = All_models50['Region'], style = All_models50['Region'])\n",
    "sc.set(ylabel='Facility Count')\n",
    "sc.set(xlim=(0,50))\n",
    "plt.show()\n",
    "#sc.figure.savefig('exampleSupplyCurve_Num_best.png', dpi = 700, bbox_inches = \"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "facIDs = SPP10['Facility ID'].unique().tolist()\n",
    "inds = []\n",
    "for i in SPP10['Facility ID'].tolist():\n",
    "    inds.append(facIDs.index(i))\n",
    "SPP10['Indexer'] = inds\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.scatterplot(data=SPP10, x=\"Indexer\", y=\"NGBE_Cost\", style= \"Generator\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SPP10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}