{ "cells": [ { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import time\n", "from datetime import datetime\n", "import matplotlib.pyplot as plt\n", "from scipy import interpolate\n", "\n", "import sys\n", "import os\n", "current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n", "parent = os.path.dirname(current)\n", "sys.path.append(parent)\n", "from functions.pressure_conversion import pressure_conversion\n", "from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n", "from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class\n", "from Turbinen.Turbinen_class_file import Francis_Turbine\n", "from Regler.Regler_class_file import PI_controller_class\n", "from Kraftwerk.Kraftwerk_class_file import Kraftwerk_class" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "def interpolate_validation_data(x,y,t_vec):\n", " interpol_fun = interpolate.interp1d(x,y,'linear',fill_value='extrapolate')\n", " return(interpol_fun(t_vec))\n" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "ename": "KeyError", "evalue": "\"None of ['Timestamp'] are in the columns\"", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mKeyError\u001b[0m Traceback (most recent call last)", "\u001b[1;32mv:\\georg\\Documents\\Persönliche Dokumente\\Arbeit\\Kelag\\Coding\\Python\\DT_Slot_3\\Kelag_DT_Slot_3\\Validation Data\\Validation_Untertweng.ipynb Cell 3\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 3\u001b[0m validation_data_UT \u001b[39m=\u001b[39m pd\u001b[39m.\u001b[39mread_csv(\u001b[39m'\u001b[39m\u001b[39mValidation_data_UT.csv\u001b[39m\u001b[39m'\u001b[39m)\u001b[39m.\u001b[39mset_index(\u001b[39m'\u001b[39m\u001b[39mTimestamp\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[0;32m 7\u001b[0m validation_data \u001b[39m=\u001b[39m validation_data_UT\u001b[39m.\u001b[39mjoin(validation_data_TB,how\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mouter\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[1;32m----> 8\u001b[0m validation_data\u001b[39m.\u001b[39;49mset_index(\u001b[39m'\u001b[39;49m\u001b[39mTimestamp\u001b[39;49m\u001b[39m'\u001b[39;49m)\n\u001b[0;32m 9\u001b[0m validation_data\u001b[39m.\u001b[39msort_index()\n\u001b[0;32m 12\u001b[0m comp_val_OL_M1_p \u001b[39m=\u001b[39m pressure_conversion(validation_data[\u001b[39m'\u001b[39m\u001b[39mOL_T1_p\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m.\u001b[39mto_numpy(copy\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m),\u001b[39m'\u001b[39m\u001b[39mbar\u001b[39m\u001b[39m'\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mPa\u001b[39m\u001b[39m'\u001b[39m)\n", "File \u001b[1;32mc:\\Users\\georg\\anaconda3\\envs\\DT_Slot_3\\lib\\site-packages\\pandas\\util\\_decorators.py:311\u001b[0m, in \u001b[0;36mdeprecate_nonkeyword_arguments..decorate..wrapper\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m 305\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mlen\u001b[39m(args) \u001b[39m>\u001b[39m num_allow_args:\n\u001b[0;32m 306\u001b[0m warnings\u001b[39m.\u001b[39mwarn(\n\u001b[0;32m 307\u001b[0m msg\u001b[39m.\u001b[39mformat(arguments\u001b[39m=\u001b[39marguments),\n\u001b[0;32m 308\u001b[0m \u001b[39mFutureWarning\u001b[39;00m,\n\u001b[0;32m 309\u001b[0m stacklevel\u001b[39m=\u001b[39mstacklevel,\n\u001b[0;32m 310\u001b[0m )\n\u001b[1;32m--> 311\u001b[0m \u001b[39mreturn\u001b[39;00m func(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n", "File \u001b[1;32mc:\\Users\\georg\\anaconda3\\envs\\DT_Slot_3\\lib\\site-packages\\pandas\\core\\frame.py:5494\u001b[0m, in \u001b[0;36mDataFrame.set_index\u001b[1;34m(self, keys, drop, append, inplace, verify_integrity)\u001b[0m\n\u001b[0;32m 5491\u001b[0m missing\u001b[39m.\u001b[39mappend(col)\n\u001b[0;32m 5493\u001b[0m \u001b[39mif\u001b[39;00m missing:\n\u001b[1;32m-> 5494\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mKeyError\u001b[39;00m(\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mNone of \u001b[39m\u001b[39m{\u001b[39;00mmissing\u001b[39m}\u001b[39;00m\u001b[39m are in the columns\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[0;32m 5496\u001b[0m \u001b[39mif\u001b[39;00m inplace:\n\u001b[0;32m 5497\u001b[0m frame \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\n", "\u001b[1;31mKeyError\u001b[0m: \"None of ['Timestamp'] are in the columns\"" ] } ], "source": [ "# import validation data\n", "validation_data_TB = pd.read_csv('Validation_data_TB.csv').set_index('Timestamp')\n", "validation_data_UT = pd.read_csv('Validation_data_UT.csv').set_index('Timestamp')\n", "\n", "\n", "\n", "validation_data = validation_data_UT.join(validation_data_TB,how='outer')\n", "validation_data.sort_index()\n", "\n", "\n", "comp_val_OL_M1_p = pressure_conversion(validation_data['OL_T1_p'].to_numpy(copy=True),'bar','Pa')\n", "comp_val_OL_M1_LA = validation_data['OL_T1_LA'].to_numpy(copy=True)/100.\n", "comp_val_OL_M2_p = pressure_conversion(validation_data['OL_T2_p'].to_numpy(copy=True),'bar','Pa')\n", "comp_val_OL_M2_LA = validation_data['OL_T2_LA'].to_numpy(copy=True)/100.\n", "comp_val_t_vec_TB = validation_data.index.to_numpy(copy=True) \n", "comp_val_t_vec_TB -= np.min(comp_val_t_vec_TB)\n", "\n", "comp_val_OL_level = validation_data['TB-Pegel'].to_numpy(copy=True)\n", "comp_val_UL_M1_p = pressure_conversion(validation_data['UL_T1_p'].to_numpy(copy=True),'bar','Pa')\n", "comp_val_UL_M1_LA = validation_data['UL_T1_LA'].to_numpy(copy=True)/100.\n", "comp_val_UL_M2_p = pressure_conversion(validation_data['UL_T2_p'].to_numpy(copy=True),'bar','Pa')\n", "comp_val_UL_M2_LA = validation_data['UL_T2_LA'].to_numpy(copy=True)/100.\n", "comp_val_t_vec_UT = validation_data.index.to_numpy(copy=True)\n", "comp_val_t_vec_UT -= np.min(comp_val_t_vec_UT) \n", "\n" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "f280f339c1e94ba2a20ff9341df5ea96", "version_major": 2, "version_minor": 0 }, "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABhtElEQVR4nO3deXhTZcI+/jtN2zQtaVlK24QlLcra4gLMyyKyKhYRGWAYwGWozPDFH6IURwSUyqKlgA7jOzDaQR0GZMC+ruMMLixDGRkqQqFYAWVrS4WW2tImXdMlz++PNAdD9+Ukac79ua5cmuTJeZ4e2uTOsx2VEEKAiIiIiBTDy9UNICIiIiLnYgAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhgGQCIiIiKFYQAkIiIiUhhvVzegI7Narbh27Rp0Oh1UKpWrm0NERETNIIRAcXExDAYDvLyU2RfGANgG165dQ69evVzdDCIiImqF7Oxs9OzZ09XNcAkGwDbQ6XQAbL9AgYGBLm4NERERNYfZbEavXr2kz3ElYgBsA/uwb2BgIAMgERFRB6Pk6VvKHPgmIiIiUjBZA2B1dTVWrVqFiIgIaLVa9OnTB+vWrYPVapXKxMTEQKVSOdxGjBghPX/jxg08/fTT6N+/P/z9/dG7d28888wzMJlMTdb/xhtvICIiAn5+fhg6dCi++uorh+ebqpuIiIjIE8k6BLxx40YkJiZix44diIyMxIkTJ/DEE08gKCgIS5YskcpFR0dj+/bt0n1fX1/p/69du4Zr167htddew6BBg5CVlYUnn3wS165dwwcffNBg3UlJSYiNjcUbb7yBe+65B3/5y18wefJknD17Fr17925W3URERESeSCWEEHId/KGHHkJoaCjeeecd6bGZM2fC398f7777LgBbL1xRURE++eSTZh/3/fffx2OPPYbS0lJ4e9efYYcPH44hQ4bgzTfflB4bOHAgfvnLXyIhIaHVdf+c2WxGUFAQTCYT5wASERF1EPz8lnkIePTo0Th48CDOnz8PADh9+jSOHDmCBx980KFccnIyQkJC0K9fPyxYsAB5eXmNHtf+D9ZQ+KusrERqaiomTZrk8PikSZNw9OjRNtVNRERE1NHJOgS8fPlymEwmDBgwAGq1GjU1NYiPj8fcuXOlMpMnT8asWbNgNBqRkZGBuLg4TJgwAampqdBoNHWOWVBQgJdffhkLFy5ssN78/HzU1NQgNDTU4fHQ0FDk5ua2um6LxQKLxSLdN5vNLTofRERERO5A1gCYlJSEXbt2Yffu3YiMjERaWhpiY2NhMBgwb948AMDs2bOl8lFRURg2bBiMRiP27t2LGTNmOBzPbDZjypQpGDRoEFavXt1k/bcu7xZCODzWkroBICEhAWvXrm3eD09ERETkpmQdAl62bBlWrFiBOXPmYPDgwXj88cexdOlSaQ5effR6PYxGIy5cuODweHFxMaKjo9GpUyd8/PHH8PHxafAYwcHBUKvVDr19AJCXl1enV7A5ddutXLkSJpNJumVnZzd4LCIiIiJ3JWsALCsrq3ONPbVa7bANzK0KCgqQnZ0NvV4vPWY2mzFp0iT4+vri008/hZ+fX6P1+vr6YujQodi/f7/D4/v378eoUaNaVPfPaTQaadNnbv5MREREHZWsAXDq1KmIj4/H3r17kZmZiY8//hibN2/G9OnTAQAlJSV47rnnkJKSgszMTCQnJ2Pq1KkIDg6WyhQXF2PSpEkoLS3FO++8A7PZjNzcXOTm5qKmpkaqa+LEidi6dat0/9lnn8Xbb7+Nv/71rzh37hyWLl2KK1eu4Mknn2x23URERESeSNY5gFu2bEFcXBwWLVqEvLw8GAwGLFy4EC+99BIAW29geno6du7ciaKiIuj1eowfPx5JSUnS9flSU1Nx7NgxAMDtt9/ucPyMjAyEh4cDAC5duoT8/HzpudmzZ6OgoADr1q1DTk4OoqKi8Nlnn8FoNDa7biIiIiJPJOs+gJ6O+wgRERF1PPz85rWAiYiIiBSHAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYRgAiYiIiBSGAZCIiIhIYWQNgNXV1Vi1ahUiIiKg1WrRp08frFu3DlarVSoTExMDlUrlcBsxYoT0/I0bN/D000+jf//+8Pf3R+/evfHMM8/AZDI1Wf8bb7yBiIgI+Pn5YejQofjqq68cnhdCYM2aNTAYDNBqtRg3bhzOnDnTfieAiIiIyA3JGgA3btyIxMREbN26FefOncOmTZvw6quvYsuWLQ7loqOjkZOTI90+++wz6blr167h2rVreO2115Ceno6//e1v+OKLL/Db3/620bqTkpIQGxuLF198EadOncK9996LyZMn48qVK1KZTZs2YfPmzdi6dSuOHz+OsLAw3H///SguLm7fE0FERETkRlRCCCHXwR966CGEhobinXfekR6bOXMm/P398e677wKw9QAWFRXhk08+afZx33//fTz22GMoLS2Ft7d3vWWGDx+OIUOG4M0335QeGzhwIH75y18iISEBQggYDAbExsZi+fLlAACLxYLQ0FBs3LgRCxcubLIdZrMZQUFBMJlMCAwMbHb7iYiIyHX4+S1zD+Do0aNx8OBBnD9/HgBw+vRpHDlyBA8++KBDueTkZISEhKBfv35YsGAB8vLyGj2u/R+sofBXWVmJ1NRUTJo0yeHxSZMm4ejRowCAjIwM5ObmOpTRaDQYO3asVIaIiIjIE9WfoNrJ8uXLYTKZMGDAAKjVatTU1CA+Ph5z586VykyePBmzZs2C0WhERkYG4uLiMGHCBKSmpkKj0dQ5ZkFBAV5++eVGe+jy8/NRU1OD0NBQh8dDQ0ORm5sLANJ/6yuTlZVV73EtFgssFot032w2N3EGiDzP/x3PhrmiCg9EhqFXV38AQGW1Fa8fOI9eXf3xq6E94aPm+jIiIncmawBMSkrCrl27sHv3bkRGRiItLQ2xsbEwGAyYN28eAGD27NlS+aioKAwbNgxGoxF79+7FjBkzHI5nNpsxZcoUDBo0CKtXr26yfpVK5XBfCFHnseaUsUtISMDatWubrJfIk20/molzOWb0DdVJAfBaUTneSL4EPx8vzPlFLxe3kIiImiLr1/Rly5ZhxYoVmDNnDgYPHozHH38cS5cuRUJCQoOv0ev1MBqNuHDhgsPjxcXFiI6ORqdOnfDxxx/Dx8enwWMEBwdDrVZLvXx2eXl5Uo9fWFgYADRa5lYrV66EyWSSbtnZ2Q3/8EQeylxeBQAI0t78G8wuLAMA9Ozi3+AXKCIich+yBsCysjJ4eTlWoVarHbaBuVVBQQGys7Oh1+ulx8xmMyZNmgRfX198+umn8PPza7ReX19fDB06FPv373d4fP/+/Rg1ahQAICIiAmFhYQ5lKisrcfjwYanMrTQaDQIDAx1uREpjqicA/lhYDgDo1UXrkjYREVHLyDoEPHXqVMTHx6N3796IjIzEqVOnsHnzZsyfPx8AUFJSgjVr1mDmzJnQ6/XIzMzECy+8gODgYEyfPh2Aredv0qRJKCsrw65du2A2m6W5d927d4darQYATJw4EdOnT8fixYsBAM8++ywef/xxDBs2DCNHjsS2bdtw5coVPPnkkwBsQ7+xsbFYv349+vbti759+2L9+vXw9/fHI488IudpIeqwqmusKLFUA7ilB/CGrQfQPiRMRETuTdYAuGXLFsTFxWHRokXIy8uDwWDAwoUL8dJLLwGw9Qamp6dj586dKCoqgl6vx/jx45GUlASdTgcASE1NxbFjxwAAt99+u8PxMzIyEB4eDgC4dOkS8vPzpedmz56NgoICrFu3Djk5OYiKisJnn30Go9EolXn++edRXl6ORYsWobCwEMOHD8e+ffukuonIkbmiWvr/QL+bbx/ZUg8gAyARUUcg6z6Ano77CJHS/FhYhl/++SisQuBk3P3S40/vOYX9Z3Px+uy7EB2lb+QIRESux89vBsA24S8Q0U1CCFgFoPbiIhAicm/8/JZ5CJiIlEOlUkHN7EdE1CFwt1YiIiIihWEAJCIiIlIYBkAiIiIihWEAJCIiIlIYBkAiapbLP5Vg7T/P4N2vs1zdFCIiaiMGQCJqlgt5Jdj+30x8kPqjq5tCRERtxABIRM1yo7QSABAc4OvilhARUVsxABJRs9gDYFcGQCKiDo8BkIiapaDEFgC7ddK4uCVERNRWDIBE1Cw3Si0AgG7sASQi6vAYAImoWUzlVQCAIK2Pi1tCRERtxQBIRM1SXFENAND58RLiREQdHQMgETXLzQDIHkAioo6OAZCImqW4wjYE3Ik9gEREHR4DIBE1S7HF1gPYScMASETU0TEAElGThBAotXAOIBGRp2AAJKImVVRZYRW2/w9gDyARUYfHAEhETSqp7f0DAH8ftQtbQkRE7YEBkIiaZB/+DfBVw8tL5eLWEBFRWzEAElGT7D2AXAFMROQZGACJqEn2AMj5f0REnoEBkIiaVMotYIiIPAoDIBE1SeoB9GUAJCLyBAyARNSkUksNAM4BJCLyFAyARNQkDgETEXkWBkAialKxtAiEewASEXkCBkAialIpVwETEXkUBkAiapJ0HWAGQCIij6ASQghXN6KjMpvNCAoKgslkQmBgoKubQySbiqoaFFdUQ+PjhUA/H1c3h4ioTfj5DfDrPBE1yc9HDT9eA5iIyGPIOgRcXV2NVatWISIiAlqtFn369MG6detgtVqlMjExMVCpVA63ESNGOBxn27ZtGDduHAIDA6FSqVBUVNRk3cXFxYiNjYXRaIRWq8WoUaNw/PhxhzLNqZuIiIjI08jaA7hx40YkJiZix44diIyMxIkTJ/DEE08gKCgIS5YskcpFR0dj+/bt0n1fX1+H45SVlSE6OhrR0dFYuXJls+r+3e9+h++++w7vvvsuDAYDdu3ahfvuuw9nz55Fjx49ml03EQFPvpuKzIJSvPzLKPwivKurm0NERG0kawBMSUnBtGnTMGXKFABAeHg49uzZgxMnTjiU02g0CAsLa/A4sbGxAIDk5ORm1VteXo4PP/wQ//jHPzBmzBgAwJo1a/DJJ5/gzTffxCuvvNLsuokIuPRTCS7klaCqxtp0YSIicnuyDgGPHj0aBw8exPnz5wEAp0+fxpEjR/Dggw86lEtOTkZISAj69euHBQsWIC8vr031VldXo6amBn5+fg6Pa7VaHDlyRNa6iTxReZXtSiBazgMkIvIIsvYALl++HCaTCQMGDIBarUZNTQ3i4+Mxd+5cqczkyZMxa9YsGI1GZGRkIC4uDhMmTEBqaio0Gk2r6tXpdBg5ciRefvllDBw4EKGhodizZw+OHTuGvn37trpui8UCi8Ui3Tebza1qH1FHU2EPgL4MgEREnkDWAJiUlIRdu3Zh9+7diIyMRFpaGmJjY2EwGDBv3jwAwOzZs6XyUVFRGDZsGIxGI/bu3YsZM2a0uu53330X8+fPR48ePaBWqzFkyBA88sgjOHnypFSmpXUnJCRg7dq1rW4TUUdVXskeQCIiTyLrEPCyZcuwYsUKzJkzB4MHD8bjjz+OpUuXIiEhocHX6PV6GI1GXLhwoU1133bbbTh8+DBKSkqQnZ2Nb775BlVVVYiIiGh13StXroTJZJJu2dnZbWojUUcghEAZh4CJiDyKrD2AZWVl8PJyzJhqtdphG5hbFRQUIDs7G3q9vl3aEBAQgICAABQWFuLLL7/Epk2bWl23RqNp9bA0UUdlqbbCvl08h4CJiDyDrD2AU6dORXx8PPbu3YvMzEx8/PHH2Lx5M6ZPnw4AKCkpwXPPPYeUlBRkZmYiOTkZU6dORXBwsFQGAHJzc5GWloaLFy8CANLT05GWloYbN25IZSZOnIitW7dK97/88kt88cUXyMjIwP79+zF+/Hj0798fTzzxRIvqJlI6+/w/ANwMmojIQ8jaA7hlyxbExcVh0aJFyMvLg8FgwMKFC/HSSy8BsPUGpqenY+fOnSgqKoJer8f48eORlJQEnU4nHScxMdFh7p19a5ft27cjJiYGAHDp0iXk5+dLZUwmE1auXIkff/wRXbt2xcyZMxEfHw8fH58W1U2kdPYVwD5qFXzUvHw4EZEn4LWA24DXEiQluPRTCSb+4TB0ft5IX/OAq5tDRNRm/PyWeQiYiDo++wpgf87/IyLyGAyARNSoCq4AJiLyOAyARNQo+xxALgAhIvIcDIBE1ChpE2iZhoCraqz4RfwB/CL+AL67apKlDiIiciTrKmAi6vjsPYByzgH8qdh2icVqK9ekERE5AwMgETVK7svAeXup8PmSewEA4d0CZKmDiIgcMQASUaPkngOoUqkwUK/MbRiIiFyFAZCIGlXuhFXAe765gm9/LMLDd/bAyNu6yVYPERHZcBEIETXKGfsAHrmQjz3fZOOHXLNsdRAR0U0MgETUKHsA9JMxAGp8bG9FFdVW2eogIqKbGACJqFHOGAK2zy+0bzpNRETyYgAkokY5JQB6245tYQ8gEZFTMAASUaOcMQfQzz4EzB5AIiKnYAAkokbZewA1ThkCZg8gEZEzMAASUaPKnNAD2C9Uhyl36HFHzyDZ6iAiopu4DyARNcoZQ8DRUWGIjgqT7fhEROSIAZCIGnVzEYh8bxeV1VbcKK1EjRDo0VkrWz1ERGTDIWAiapQzegD//X0eRiQcxDN7TslWBxER3cQASESNKqusBgBoZQyAgVpb76K5vEq2OoiI6CYGQCJqlH0RiJz7AAb6+QAAzBUMgEREzsAASEQNqrEKaXNmOYeApQBYXi1bHUREdBMDIBE16OcbM/v7yrcIxD4EXF5Vg6oa7gVIRCQ3rgImogbZh38BQOMt3/dFnZ8P4h4aBJ2fN4SQrRoiIqrFAEhEDSr/2fw/Ly+VbPWovVT47egI2Y5PRESOOARMRA0qq7LNyZNz/p/dv7+/jv87no38EovsdRERKR0DIBE1SOoBdEIAfOVf5/D8h9/i8k+lstdFRKR0DIBE1CBnbAJtF6CxzUgptXAlMBGR3BgAiahBztgD0C5AY6ujmAGQiEh2DIBE1KCyKucNAXdiDyARkdMwABJRg8or7YtA5N8wgEPARETOwwBIRA1y5iIQew9gCQMgEZHsGACJqEHSELAT5gBKAbCCAZCISG6yBsDq6mqsWrUKERER0Gq16NOnD9atWwer9ealnmJiYqBSqRxuI0aMcDjOtm3bMG7cOAQGBkKlUqGoqKjJuouLixEbGwuj0QitVotRo0bh+PHjDmWEEFizZg0MBgO0Wi3GjRuHM2fOtMvPTuQJXLIKuJIBkIhIbrIGwI0bNyIxMRFbt27FuXPnsGnTJrz66qvYsmWLQ7no6Gjk5ORIt88++8zh+bKyMkRHR+OFF15odt2/+93vsH//frz77rtIT0/HpEmTcN999+Hq1atSmU2bNmHz5s3YunUrjh8/jrCwMNx///0oLi5u2w9O5CHKnDgEHCANAdc0UZKIiNpK1pndKSkpmDZtGqZMmQIACA8Px549e3DixAmHchqNBmFhYQ0eJzY2FgCQnJzcrHrLy8vx4Ycf4h//+AfGjBkDAFizZg0++eQTvPnmm3jllVcghMDrr7+OF198ETNmzAAA7NixA6Ghodi9ezcWLlzYwp+WyPOU1w4B+/vIvwikU+02MCUVVbLXRUSkdLL2AI4ePRoHDx7E+fPnAQCnT5/GkSNH8OCDDzqUS05ORkhICPr164cFCxYgLy+vTfVWV1ejpqYGfn5+Do9rtVocOXIEAJCRkYHc3FxMmjRJel6j0WDs2LE4evRom+on8hTOHAIe2ScYf5p7NxZP6Ct7XURESifr1/rly5fDZDJhwIABUKvVqKmpQXx8PObOnSuVmTx5MmbNmgWj0YiMjAzExcVhwoQJSE1NhUajaVW9Op0OI0eOxMsvv4yBAwciNDQUe/bswbFjx9C3r+3DJTc3FwAQGhrq8NrQ0FBkZWXVe1yLxQKL5eZ1Ss1mc6vaR9RRlNXOx/NzQgDs3c0fvbv5y14PERHJHACTkpKwa9cu7N69G5GRkUhLS0NsbCwMBgPmzZsHAJg9e7ZUPioqCsOGDYPRaMTevXulodnWePfddzF//nz06NEDarUaQ4YMwSOPPIKTJ086lFOpVA73hRB1HrNLSEjA2rVrW90moo7GPgfQ3wmrgPNLLPjX6WtQe6nw+Mhw2esjIlIyWYeAly1bhhUrVmDOnDkYPHgwHn/8cSxduhQJCQkNvkav18NoNOLChQttqvu2227D4cOHUVJSguzsbHzzzTeoqqpCREQEAEhzDu09gXZ5eXl1egXtVq5cCZPJJN2ys7Pb1EYid+fMIeCfii1Y88+z+N+DbfvbJyKipskaAMvKyuDl5ViFWq122AbmVgUFBcjOzoZer2+XNgQEBECv16OwsBBffvklpk2bBgCIiIhAWFgY9u/fL5WtrKzE4cOHMWrUqHqPpdFoEBgY6HAj8mTlLrgUHDeCJiKSn6xDwFOnTkV8fDx69+6NyMhInDp1Cps3b8b8+fMBACUlJVizZg1mzpwJvV6PzMxMvPDCCwgODsb06dOl4+Tm5iI3NxcXL14EAKSnp0On06F3797o2rUrAGDixImYPn06Fi9eDAD48ssvIYRA//79cfHiRSxbtgz9+/fHE088AcA29BsbG4v169ejb9++6Nu3L9avXw9/f3888sgjcp4Wog5DuhKIE4aA7dvAVFRZUV1jhbea+9QTEclF1gC4ZcsWxMXFYdGiRcjLy4PBYMDChQvx0ksvAbD1Bqanp2Pnzp0oKiqCXq/H+PHjkZSUBJ1OJx0nMTHRYe6dfWuX7du3IyYmBgBw6dIl5OfnS2VMJhNWrlyJH3/8EV27dsXMmTMRHx8PHx8fqczzzz+P8vJyLFq0CIWFhRg+fDj27dvnUDeRkklzAJ1yLeCbIbPUUoMgfwZAIiK5qIQQwtWN6KjMZjOCgoJgMpk4HEwe6Y41X8JcUY0Dz47F7SGdZK+v34ufo7LGiv+umIAenbWy10dEysTPb14LmIgaUVFlm6/rjDmAwM1ewFLOAyQikhUDIBHVq7rGisoaWwB0xjYwANDJzzbUXFzBAEhEJCcGQCKqV0X1zdX6TusBrJ1ryB5AIiJ5MQASUb3sK4BVKkDj7Zy3CvtWMAyARETyYgAkonr9fAuYhq6O097sW8EUMwASEclK/r0diKhDkjaBdtL8PwB4blJ/PDn2NqesOCYiUjIGQCKqlz0A+jkxAA7uGeS0uoiIlIxDwERUL2kI2EkLQIiIyHkYAImoXhUuGAL+4rscLN59EruPXXFanURESsQASET1csUcwEs/leJf3+bg1JVCp9VJRKREDIBEVC/7ELCfE4eAA7W2a3WbyqucVicRkRIxABJRvW72ADrvbSKIAZCIyCkYAImoXq6YA8gASETkHAyARFQvV6wCDuS1gImInIIBkIjq5Yp9ANkDSETkHAyARFQvV6wCtgfAEks1qmusTquXiEhpGACJqF6umANoXwUMcBiYiEhODIBEVC9XzAH0UXvBv7Y+DgMTEcmHAZCI6lXmokvBcR4gEZH8GACJqF6umAMIMAASETkDAyAR1csVcwABINDPFgDNFQyARERyYQAkonpJ28A4eQiYl4MjIpIfAyAR1UtaBMIhYCIij+Pt6gYQkXuqqLLtw+fsADj3f3phbP/uiDQEOrVeIiIlYQAkonpJi0CcPAQ8LLyrU+sjIlIiBkAiqperhoCzCkpx+PxP6OLvi6l3GpxaNxGRUnAOIBHVIYRwybWAAeDsNTNe+scZ7Dia6dR6iYiUhAGQiOqwVN+8Di83giYi8jwMgERUh334FwD8vJ37NsFtYIiI5McASER12Id/fdVe8FY7922CPYBERPJjACSiOm7O/3P+W0SQvy0AWqqt0tVIiIiofTEAElEd9iFgf1/nbxTQydcbXirb/7MXkIhIHrIGwOrqaqxatQoRERHQarXo06cP1q1bB6v15gTzmJgYqFQqh9uIESMcjrNt2zaMGzcOgYGBUKlUKCoqclrdRErkqj0AAcDLS8V5gEREMpP16/3GjRuRmJiIHTt2IDIyEidOnMATTzyBoKAgLFmyRCoXHR2N7du3S/d9fX0djlNWVobo6GhER0dj5cqVTq2bSInsPYDO3gLGLkjrg6KyKgZAIiKZyBoAU1JSMG3aNEyZMgUAEB4ejj179uDEiRMO5TQaDcLCwho8TmxsLAAgOTnZ6XUTKZHUA+iCOYAA0FnrgywApjIGQCIiOcj67j569GgcPHgQ58+fBwCcPn0aR44cwYMPPuhQLjk5GSEhIejXrx8WLFiAvLy8Dl03UUdX4cIhYIBbwRARyU3WHsDly5fDZDJhwIABUKvVqKmpQXx8PObOnSuVmTx5MmbNmgWj0YiMjAzExcVhwoQJSE1NhUajcau6LRYLLBaLdN9sNre6fUTuzFWXgbPjVjBERPKSNQAmJSVh165d2L17NyIjI5GWlobY2FgYDAbMmzcPADB79mypfFRUFIYNGwaj0Yi9e/dixowZblV3QkIC1q5d2+o2EXUUrroMnB0DIBGRvGQNgMuWLcOKFSswZ84cAMDgwYORlZWFhIQEKYTdSq/Xw2g04sKFC25X98qVK/Hss89K981mM3r16tWmdhK5o5tzABkAiYg8kawBsKysDF5ejtMM1Wq1w1YstyooKEB2djb0er3b1a3RaNo0LE3UUVRUunYOYOfazaCLyipdUj8RkaeTdRHI1KlTER8fj7179yIzMxMff/wxNm/ejOnTpwMASkpK8NxzzyElJQWZmZlITk7G1KlTERwcLJUBgNzcXKSlpeHixYsAgPT0dKSlpeHGjRtSmYkTJ2Lr1q3tXjeRErm6B7BrgO2L1g2uAiYikoWsPYBbtmxBXFwcFi1ahLy8PBgMBixcuBAvvfQSAFuPXHp6Onbu3ImioiLo9XqMHz8eSUlJ0Ol00nESExMd5t6NGTMGALB9+3bExMQAAC5duoT8/Px2r5tIiVw9B7BbJ9t+nPnFliZKEhFRa6iEEMLVjeiozGYzgoKCYDKZEBgY6OrmELWb3//faXx48kesmDwAT469zen1p/9owtStRxAaqMGxF+5zev1E5Nn4+c1rARNRPSpcPARs7wEsKKkEv6MSEbU/BkAiqsPVcwC7+NsCYLVVoKx2QQoREbUfWecAElHH1C9UhxJLNcKC/FxSv5+PF3zUKlTVCJjKqxCg4VsVEVF74rsqEdWxYvIAl9avUqkQ6OeDgtJKmCuqYIDWpe0hIvI0HAImojrKKqtRXdPwnpnOYL8esLm82qXtICLyRAyARFTHxD8cxu0vfo5vfyxyWRtuBkDuBUhE1N4YAImojgoX7wMIAIF+thkqvBwcEVH7YwAkojoqqmzDv37eLgyA9h7ACgZAIqL2xgBIRA6EEKiotvcAuu4tIohzAImIZMMASEQOKmussO+9rHHpEDB7AImI5MIASEQO7MO/gGt7AAO1nANIRCQXBkAicmBfAKL2UsFX7cIA6MdVwEREcmEAJCIHP78OsEqlclk7grgIhIhINgyAROSgvMr1C0AAbgRNRCQnBkAiclBe6fo9AAHuA0hEJCcGQCJyUP6zIWBX4j6ARETyYQAkIgfSHEBf1wZA+xzAEks1rFbh0rYQEXkaBkAiclBeWXsVEBf3AOpqh4CFAIotnAdIRNSeGACJyIG7DAFrvNXSQhRuBUNE1L4YAInIgbsEQODmXoBcCEJE1L4YAInIgcVNtoEBuBcgEZFcXP8OT0Ruxb4NjKsXgQDcC5CISC4MgETk4OZG0G4QAGsXgnAOIBFR+2IAJCIHbjUHkEPARESyYAAkIgcVbhQApTmA7AEkImpXDIBE5MCt5gD62XsAOQeQiKg9MQASkQO3mgOo5fWAiYjkwABIRA4qqtzjSiAAh4CJiOTi7eoGEJF7cadFIH1DdZjzi16INAS6uilERB6FAZCIHEiLQHxdP0AwpHcXDOndxdXNICLyOK5/hyeielVU1eDC9WKcyzE7tV77IhB3GAKusQpczCvBqSuFEEK4ujlERB6DAZDITaVlF+H+P/4HT+0+6dR63WkIuKKqBvdtPozpbxyV5iYSEVHbMQASuanO/rYFEKYy5y6AuDkE7PoA6O+rhtpLBYCbQRMRtSdZA2B1dTVWrVqFiIgIaLVa9OnTB+vWrYPVevObfExMDFQqlcNtxIgRDsfZtm0bxo0bh8DAQKhUKhQVFbVL3UIIrFmzBgaDAVqtFuPGjcOZM2fa7ecnagv7ClhTeZVThz+lfQDdoAdQpVKhk8Y2VbmYAZCIqN3Iughk48aNSExMxI4dOxAZGYkTJ07giSeeQFBQEJYsWSKVi46Oxvbt26X7vr6+DscpKytDdHQ0oqOjsXLlynare9OmTdi8eTP+9re/oV+/fnjllVdw//3344cffoBOp2uHM0DUep21tr+DaqtAaWWNFITkJIRwqyFgwLYXoKm8CqZybgZNRNReZP1ESUlJwbRp0zBlyhQAQHh4OPbs2YMTJ044lNNoNAgLC2vwOLGxsQCA5OTkdqtbCIHXX38dL774ImbMmAEA2LFjB0JDQ7F7924sXLiw2XURycHPxwu+3l6orLaiqKzSKQGwqkbAWtvZqHGTAKjT+AAoZw8gEVE7knUIePTo0Th48CDOnz8PADh9+jSOHDmCBx980KFccnIyQkJC0K9fPyxYsAB5eXmy152RkYHc3FxMmjRJeo1Go8HYsWNx9OjReo9psVhgNpsdbkRyUalUDsPAzmDv/QPcpwdQ52cLvrwcHBFR+5G1S2H58uUwmUwYMGAA1Go1ampqEB8fj7lz50plJk+ejFmzZsFoNCIjIwNxcXGYMGECUlNTodFoZKs7NzcXABAaGurwutDQUGRlZdV7zISEBKxdu7bVbSLX+O/FfJRaqjGmX3e32NqkJTprffBTscVpC0HsC0DUXir4qFVOqbMpax6ORHWNQO9u/q5uChGRx5A1ACYlJWHXrl3YvXs3IiMjkZaWhtjYWBgMBsybNw8AMHv2bKl8VFQUhg0bBqPRiL1790pDs3LVDdh6WX5OCFHnMbuVK1fi2Wefle6bzWb06tWr1W0k5/jtjuOoqLLiq+fHo1fXxkPEP9KuoriiGmP7dW+yrDPYVwIXOasH8GcLQBr6O3C2gXpeBYSIqL3JGgCXLVuGFStWYM6cOQCAwYMHIysrCwkJCQ4h7Of0ej2MRiMuXLgga932OYe5ubnQ6/XS6/Ly8ur0CtppNJo29UqSawT4eqOiqhJllTVNlt1xNBMnrxThT3PvdnoAPJdjxp5vruC27p0wb1Q4ALhsCLij9ZQSEVHLyDoHsKysDF5ejlWo1WqHrVhuVVBQgOzsbIdQJkfdERERCAsLw/79+6XnKysrcfjwYYwaNapNdZN78dfYwkxpZeNzyIQQuHC9BADQL7ST7O261eWfSrEzJQufpedIjwXVrgQuctIQcLkbXQbObt+ZXPz2b8eRePiSq5tCRE70j7SrmLMtBXu/zWm6MLWYrD2AU6dORXx8PHr37o3IyEicOnUKmzdvxvz58wEAJSUlWLNmDWbOnAm9Xo/MzEy88MILCA4OxvTp06Xj5ObmIjc3FxcvXgQApKenQ6fToXfv3ujatSsAYOLEiZg+fToWL17crLpVKhViY2Oxfv169O3bF3379sX69evh7++PRx55RM7TQk4W4Gv7NS+zNN4DmGuuQLGlGmovFSKCA5zRNAelFltADfjZal9pM2gn9QBW2C8D5+0+PYC55goc/D4Pvt7uE0qpZUot1aiqsaKzv2/ThYlqfZp2DV9fvoERfboBaFunENUlawDcsmUL4uLisGjRIuTl5cFgMGDhwoV46aWXANh65NLT07Fz504UFRVBr9dj/PjxSEpKctiHLzEx0WHxxZgxYwAA27dvR0xMDADg0qVLyM/Pb3bdAPD888+jvLwcixYtQmFhIYYPH459+/ZxD0AnKLFU47q5AvogP/j7yru9iT1QNdUDeL629y+8mz80LghA9vb9PABOuUOP/mE6DHLSPLjS2gDo74QtZ5qrS21ouFFa6eKWUGt9dPJHrPvXWfx2dB+smDyg3Y6b/EMe7rk9GD5qfjnwNKbyKvznwk8AgCmDGf7kIOu7vE6nw+uvv47XX3+93ue1Wi2+/PLLJo+zZs0arFmzptEymZmZLaobsPUCNufY1P6OZ9zAE387joH6QHy+5F5Z6/KvvaSZvYetIReuFwMA+oXK8wWgstqK766ZUFhaiYkD684zlXoAf3YJtiG9u2BI7y6ytKc+ZZV12+BqXQNsAbCwjAGwo/r41FVU1QgEd2qfHkBLdQ1e/tdZ7Pr6Cv7fmD544cGB7XJcch/7z15HVY1A/1Ad+sr0nqx07vM1nxQlx1QBANAH+clel30IuLSJRSAllmpovL1ke7Mpq6zGjDdse0z+8Ep0nV5GqfetiR7REks1Pj75I/x9vTFzaM92bqO9De4TAG/2AHIj6I4oq6AUJ68UwUsFPHynoc3Hy75Rhqd2n8S3P5oAABpvr0Z3b+iI3v7qMnp28ceYfsGyj5C4q73fXgNgGwUheSjzN4tcLtdsC4BhTgiA9kUgZU30AMbe1w9PT+iLyuqGFym1RZDWBz5qFapqBApKKmHorHV4vqu/L/qH6mDo3Pg5yS+2IO4fZ9BJI2cAdJ+3hi4BtVvhlFV63Ae9EnxyyvZBfs/twQgJbNvf+/6z1/H7/0uDuaIaXfx98MfZd2Fc/5D2aKbbKLVUY/1n52AVwDcvTHSrv0VnMZVV4asLtildD3L4VzbK+81SsIt5xfjTwYsI0HgjYcZgl7Yl11QOANAH+uHdlEwcuZiPmUN6YlJkw5cEbK3m9gACtg2QtTL1fqlUKgR30iDHVIGfii11AuCCMX2wYEyfJo9jqQ2oGhkWRdhDsrv1AN7bNxhd/H1RWWN1yfxMah0hBD5JuwoAmH53j1Yfp6rGite+/AF/+c9lAMDdvTvjz48MqfM35AnOXDPDKoCwQL82B+aO6suzuai2CgwI0+H2EOfvyKAUDIAKUl5pxaenryE0UAOg/gD4+oHz+Dw9F7+7NwKzhsm3ybV9CDgsyA8nMgvx5ZnrGNwjSJYA2NweQGforrsZAFvLfrUOWQJglfv1APr5qPHub4e7uhnUCqd/NCEjvxRaHzUeaOXfdq6pAk/vOYnjmYUAgPn3RGDF5AEeuyr82x+LAACDewa5tiEuNKR3Fywc2we93WAzfk/mPu/yJLtgnW0uVUFJJaxWAS+vukNp+SUW/HC9GFdulMnallxpDqAWOj/b4otima71OmlQKHp28UekwfVXlAjuZNtIPL+k9QHQ3gMox2bNZdJWNO7Vy1ZcUYUSSzWCO2m44rMDGajXIfGxobhWVO6wur0lrhaV4dSVIug03tj0qzsw2cOHBO1zG+90swBYVWOVtqPqrPWBt4x/h7eHdMLKyVzYIzcGQAXpFmALH9VWgaLyKml15c+F6mxDDtdr5+jJ5U9z78aPheWINAQiNcv2zd4sUwAcauyKocaushy7pbrXBsD26AG8UVaJE5k30F2ngbFb++xbaB8ml2sYvLVGbfg3iiuq8e/fj0Wf7hwS6ig03mpER7WtV3+osStem3Un7uzV2SX7czpb+lVbALyjZ2fXNuQWP+QW46EtRwCAf4cegl+lFcTX2wtdajcWbiiAhNbOOclrQ0BpjoH6QNw/KBRdAnwRqLV9Dymu8PxVnkvu64uvnh/frLl+DbH3ABaVVeFXiSnY9XVWezVPuhZwgBsNAQNAp9reoxI3GMYn5/vl3T0UEf4qqmqkL3iDe7hXDyB5Hvd6lyfZBXfSoLCsCvklFvRH3e1Ougfaeqium+UNgD+n87OFUrmGgN1Je0xat39A2DXnGsfNZd+M2t16AKUAqIDfEVIuPx81UlZORH6JBV3qGaFxpageQcjcMMXVzaB2xACoMCP6dIOxm3+DqzztQ8B5Mg8B/1zPLlqMvj0YUT1cP0evJUxlVfjjgfPIvlGGt+cNc9r2JJZbtqkpb8cAWOauPYB+7AEk5bDPFSaSk3u9y5PsXv5lVKPP9+isxa+G9kRYoJ/T9lwb0adb7bUeOxY/Xy+8+3UWaqwC180Wp+xpCNTtAWzqEnctYb8SiDttAwPUHQI+l2PG0qQ0DAjT4fU5d9cpn3T8Cv5+7AoeukOP/zfmtgaP+9jbx6BSAS89NEi2DcCFEHjg9f8gLEiL12ffVe/cW6Ln3j+NKwVlWHp/P4y8rX3fD//9/XVs+uIH3D8oFL+f1L9dj00dFwOgwlVU1aCssgaBft7wVnshyN8Hr82609XNQqmlGpXVVnT292lxCK2oqoGl2gp/X7WsK0Y13mpEBAfgYl4JzuWanRYAb+0BbM8hYHe8EggA6G7pAcwxleP73GKo61nJDgCXfirFtz+a8Ivwhhf/VFZbcSyjAFU1QtYh7+tmC85fL8Gln0qln4PoViezCnE5vxRWIdr92GeumvF9brHTrilOHQMXgShcyuUCDHl5P6b9+b+uboqDf56+hrtf3o8FO0+0+LXvn8jGnWv34Zk9p2RomaMBYbZeox9yi2Wvy27SoFC8M28YFo+/HUD7zp0ss7jfPoAAsODePnhn3jDpGsr5xbbrAjc0VGafwhCia3go7XJ+CapqBHQab/SQcUPhy/klAIDeXf25hQ3Vq7rGKm29FS7DYpfva9+f+ofxmrp0E9+NFM5cu69TkNbHxS1xZF+E0pq5MPZeLK0M++TdyhUBsFdXf0wcGIp7+wYDuPlv2B7sw8n+brYP4N29u2DiwFApqP1U0vjvh30Ve0hgw78/3+fY/s0G6HVN9jJ/evoakn/Ia9UcxIz8UgBQxCpWap1rRRWotgr4entBL8PVP77PNQMABrAHkH6GAVDhispqN/b0d68AmFfcdA9OQ8qrnLeX3YAw2xvquRyz7HXdKqj236yonQKgEMJtt4G5lX0jbfvm5reSAqCu4Q9T+7+Z/d+wIUIIrP30DGK2H8fln0pa3NaMnxgA3dGXZ3Kxed8Prm4GACCjwPY7YuzqX+8G/YAtxD37f2nSZszNVVFVI30JGcAeQPoZ936XJ9mZ3LQH0P4B3r0V34alAOiEHkD7kMqln0pQVWN1yhDfB6k/4njGDdzZqzMA279heyzYqayxotpqm3/kbtvA3Cq/xDYE3L2BHkD7RuahjfQAnqvttR3YRK/IT8UWFJRWwksF9GvFQpHL7AF0Oxn5pVj095OosQro/HzatC9ne8is/R1paPhXCIGlSadxLseM3l39EXtfv2Yf+2JeCazC9iW/NV+oyXOxB7ADS/j8HKb86SucvFLY6mPYewCDtG1bmVhRVYOVH32LMZsOSStJ2+JmD07dN6xnk9Iw882jDQ67ljtxIUPPLlp00nijqkZI37LlduxyAZJOZCO3NuTUWEW7bI/y8+1k3G0RyK3yixseArZaBX49rBem3mmQNjavz30DQ/DQHXrc2avxDXfP1vYURgQHtOrye+P7d8e0uwzc2NeNRAQH4Lna1bDxn53DB6k/urQ9t4d0wpxf9MK4/t3rfV6lUklzfv96JAPmFmyaH6DxRsyocMwc0tNpW1VRx8AewA7s4vUSnLlmRvqPJgzp3aVVx2ioB3DOthRkFZRh6yN3N+syahpvLxy5mI/sG+X478UC3D8oVHpuy8EL2HUsC4+PMGLxhL7Nalee1INT9wP826smXMwrQUGpBahnM2t7kPFzQohRqVToH6ZDalYhzuWYW9VD1FL2N/8QnQa+3l6orLZdo9O+oXZr2S8D56v2cvvFCvmNzAH08lIh7qFBTR7jNyPD8ZuR4U2WO5fTvJ7Chjw+MhyPN6Mecq4nx/bBjVIL3voqA8s//BZBWh+H9y1nuuf2YNxze3CjZSZHhaFvSCdcyCvBzqOZzX4vjQgOwJqHI9ujmeRh3Ptdnho1yGD7QDp7rfXzz564Jxxb5t6N+waGODz+U7EFOaaKOluONESlUmFCf9sx/v19nsNzOeYKXDdbUFnT/O0N/jDrTvxh1p31DpvZLypfaql/+5N5o8Kx9ZG7cf9A57yZ93fyQhBzua23L1Br27Jn2+ND0cW/7XvLlbvpApD6NDUHsD3Z5wq2NgCSe1KpVHjhwYGYOaQnaqwCT+0+iWOXC1zdrAZ5eamweIKtF/DtIxncFJ3ajD2AHVhkbQA8k2Nq9TGiegQhqp6hqcTHhsJSbW3RlgTjB4RgR0oWDn2f5zAnzT5c171T8z+sRzXybbhTbUApbeANsKGfSS72idVXi8qdUp+9BzDQzxvj+oc0Ubr57IHa3wlzJ9uiqsaKwtqpC864YoI9AHIPNc+jUqmwceZgmMorceBcHn634wSSFo6Uvly7m4fuMOB/D1zA5fxSvJuShf9vXMObnBM1hT2AHdggvS3knM+1LUBoT31DdYjqESRdgaE5RvTpBq2PGrnmCmneFND4cF1r2Peoa+gKGO8cycCur7Nwo7SyXepryrS7euDEqvvwv/VckUIOUgBs54U70ibQLfg3d4WSimr06KyFv6+6XXo+G2O1Cmh8vOCjVrEH0EN5q72w9ZEh+J/wrujRRYvgFnxRdTZ1bS+g1kcty4bRpCzu/U5PjerZRQudxhvFlmpczCtx+QeUn48a99zeDQfO5eHQ93mINNgCqrRis51WoHWShoDrD4Cb9/2A0soajLytm1Muu+XsFdQDwwIRpPVBt3b+2eyLdwLcfAFIlwBf/HfFBKdcqtDLS4V/PX0vKqut8FFzAr2n8vNR4+2YYRDC/XZEuNXDdxowtl93dOP1gqmNGAA7MC8vFabdbYAQtkUY7mDKHXpovNUOQ7Dt3QM4LNy24OW27p3qPFdiqZYWM4TJsKGqO9j2m2GyHNd+3tx9Cxg7Z65o9HWTvy+ST2AbF1E5i7fai+GP2gUDYAf3yi8Hu7oJDqbf3RPT7+4p3bdaBWLv64v8kspGr8rQEo8ON+LR4cZ6n8s12VYP6zTe0mIRap5yqQeQ542IyNPxnZ5k5eWlwv8b47yJyvYNgNsrbCqJfRFIR+kBJCKi1uO4Bskqz1yBf39/XboWpdzsATAsyDOHf+Vkv4IKewCJiDwfAyDJ6ljGDcz/2wm89I8zTqnPfnWM0EauAUv1sy+qYQ8gEZHnYwAkWdm3YmnvFasNCfTzwUB9IG4LqbtAhBpn3wYmoANsBE1ERG3DsR6SlT0AdnFSAHxshBGPjah/gQg1zr4NjL8bDwH/6eAFpF81IWZUeJOXziIiooaxB5Bk5eweQGo9+zYw/m48BJyaVYj9Z68jp3a1NxERtQ4DIMnqRlltD6DMV2ygtiuvdP9FIParH6j5zkVE1CZ8GyVZ3ai9Ckg3N768Etl0hEUgNVZbAPRy4ibQRESeyH2/6pNHeGveMNwoqUTngI6xy76SSdvAuPEiEHsAVHsxABIRtYWsPYDV1dVYtWoVIiIioNVq0adPH6xbtw5Wq1UqExMTA5VK5XAbMWKEw3G2bduGcePGITAwECqVCkVFRU3WHR4eXue4KpUKTz31VIvqprbppPFG727+HeYyS0om9QD6uO/3QmkImD2ARERtIus7/caNG5GYmIgdO3YgMjISJ06cwBNPPIGgoCAsWbJEKhcdHY3t27dL9319HYcLy8rKEB0djejoaKxcubJZdR8/fhw1NTXS/e+++w73338/Zs2a5VCuqbqJlKK8A2wDM8TYBZ003rzSCxFRG8kaAFNSUjBt2jRMmTIFgK1Xbs+ePThx4oRDOY1Gg7CwsAaPExsbCwBITk5udt3du3d3uL9hwwbcdtttGDt2bIvqJlKKjrAKeOXkga5uAhGRR5B1CHj06NE4ePAgzp8/DwA4ffo0jhw5ggcffNChXHJyMkJCQtCvXz8sWLAAeXl57dqOyspK7Nq1C/Pnz4fqlqGjltRtsVhgNpsdbkSeolwKgO47BExERO1D1nf65cuXw2QyYcCAAVCr1aipqUF8fDzmzp0rlZk8eTJmzZoFo9GIjIwMxMXFYcKECUhNTYVG0z7DPJ988gmKiooQExPj8HhL605ISMDatWvbpU1E7qSy2orKGtvcXHfeBoaIiNqHSojaWdUyeO+997Bs2TK8+uqriIyMRFpaGmJjY7F582bMmzev3tfk5OTAaDTivffew4wZMxyeS05Oxvjx41FYWIjOnTs3ux0PPPAAfH198c9//rPRco3VDdh6AC0Wi3TfbDajV69eMJlMCAwMbHZ7iNyNqawKd67bBwA4/8pk+Hpzhygi8lxmsxlBQUGK/vyW9av+smXLsGLFCsyZMwcAMHjwYGRlZSEhIaHBAKjX62E0GnHhwoV2aUNWVhYOHDiAjz76qMmyTdWt0WjarVeSyJ2UVdlWAPuoVQx/REQKIOs7fVlZGby8HKtQq9UO28DcqqCgANnZ2dDr9e3Shu3btyMkJERaiNKY9q6bqKMotdjm/2l93HcBCBERtR9ZA+DUqVMRHx+PvXv3IjMzEx9//DE2b96M6dOnAwBKSkrw3HPPISUlBZmZmUhOTsbUqVMRHBwslQGA3NxcpKWl4eLFiwCA9PR0pKWl4caNG1KZiRMnYuvWrQ71W61WbN++HfPmzYO3t2NnZ3PrJlKCm1vAcP4fEZESyPpuv2XLFsTFxWHRokXIy8uDwWDAwoUL8dJLLwGw9Qamp6dj586dKCoqgl6vx/jx45GUlASdTicdJzEx0WHxxZgxYwDYevfsCzsuXbqE/Px8h/oPHDiAK1euYP78+XXa1ty6iZSgtNL9LwNHRETtR9ZFIJ6Ok0jJXZRYqnH5pxJ08fdFr67+LX79oe/z8MTfjmNwjyD88+nRMrSQiMh98POb1wIm8ghnr5nx67+koE9wAP793LgWv/5/Irpi39Ix4AXWiIiUgQGQyAMUV1QBAHR+rfuTDtB4o18opz4QESkFAyCRByix2ObwdWplAPxH2lUcy7iBCf1DcN+g0PZsGhERuSFu+EXkAcwVtQGwlat4v8m4gd3HruC7a6b2bBYREbkpBkAiD1BSGwB1fj6ter09QAa28vVERNSxMAASeYASi20OYGt7AM3lttcHahkAiYiUgHMAiTzAvFHhuH9QGLr4t7YHsG2LSIiIqGPhuz2RBwjR+SFE59fq1xdzCJiISFE4BExEPxsC5ndCIiIlYAAk8hB7v83B/L8dx1+PZLT4tfYhYPYAEhEpA7/uE3mIKzfK8O/v89AtwLdFr7NU16CiygqAi0CIiJSCPYBEHsKr9jpuNS28vLd9/p9KBehauYqYiIg6FgZAIg+hrk2ALcx/0vy/ThpveHnxasBERErAAEjkIVQqW3iztjABchNoIiLlYQAk8hD2zrtqawsDYDn3ACQiUhoGQCIP4eejBgBYqmpa9DppBTAXgBARKQYDIJGH8Pe1BcCyypYFQG4CTUSkPAyARB5CW9sDWN7SHkBuAk1EpDgMgEQewt/XFuDKW9gDyE2giYiUhwGQyENoWzkEbC6vHQLmHEAiIsVgACTyEK2dA3izB5BDwERESsEASOQh7AGwvLK6Ra+7OQeQPYBERErBAEjkIaQh4KoaiBZsBn1zI2j2ABIRKQUDIJGHsC8CEQKwVFub/TqpB5CLQIiIFIMBkMhD2LeBAVo2D5AbQRMRKQ8DIJGHUHup4Ott+5Mua8E8QG4ETUSkPAyARB7k5kKQ5vUAVtVYpd5CbgRNRKQcDIBEHsTfp2Vbwdh7/wCgk4YBkIhIKRgAiTxISzeDti8A6aTxhreabwdERErBd3wiDyJdDq6qeXMATeXcBJqISIkYAIk8SEt7AG+UVgIAugT4ytYmIiJyPwyARB6kpZeDswfArgyARESKImsArK6uxqpVqxAREQGtVos+ffpg3bp1sFpvblIbExMDlUrlcBsxYoTDcbZt24Zx48YhMDAQKpUKRUVFTdYdHh5e57gqlQpPPfWUVEYIgTVr1sBgMECr1WLcuHE4c+ZMu/38RM7W0lXAhWUMgERESiRrANy4cSMSExOxdetWnDt3Dps2bcKrr76KLVu2OJSLjo5GTk6OdPvss88cni8rK0N0dDReeOGFZtd9/Phxh2Pu378fADBr1iypzKZNm7B582Zs3boVx48fR1hYGO6//34UFxe34acmch2tj20uX4uHgP0ZAImIlETWmd8pKSmYNm0apkyZAsDWK7dnzx6cOHHCoZxGo0FYWFiDx4mNjQUAJCcnN7vu7t27O9zfsGEDbrvtNowdOxaArffv9ddfx4svvogZM2YAAHbs2IHQ0FDs3r0bCxcubHZdRO5C6gGs4hAwERE1TNYewNGjR+PgwYM4f/48AOD06dM4cuQIHnzwQYdyycnJCAkJQb9+/bBgwQLk5eW1azsqKyuxa9cuzJ8/HyqVCgCQkZGB3NxcTJo0SSqn0WgwduxYHD16tN7jWCwWmM1mhxuRO7k5BNy8VcAMgEREyiRrD+Dy5cthMpkwYMAAqNVq1NTUID4+HnPnzpXKTJ48GbNmzYLRaERGRgbi4uIwYcIEpKamQqPRtEs7PvnkExQVFSEmJkZ6LDc3FwAQGhrqUDY0NBRZWVn1HichIQFr165tlzYRyaGlq4A5B5CISJlkDYBJSUnYtWsXdu/ejcjISKSlpSE2NhYGgwHz5s0DAMyePVsqHxUVhWHDhsFoNGLv3r3S0GxbvfPOO5g8eTIMBkOd5+w9gnZCiDqP2a1cuRLPPvusdN9sNqNXr17t0kai9qD1adkikALOASQiUiRZA+CyZcuwYsUKzJkzBwAwePBgZGVlISEhQQqAt9Lr9TAajbhw4UK7tCErKwsHDhzARx995PC4fc5hbm4u9Hq99HheXl6dXkE7jUbTbr2SRHJo6TYwhRwCJiJSJFnnAJaVlcHLy7EKtVrtsA3MrQoKCpCdne0Qytpi+/btCAkJkRai2EVERCAsLExaHQzY5goePnwYo0aNape6iZxNW3slkLJmLAIRQqBviA59ugcwABIRKYysPYBTp05FfHw8evfujcjISJw6dQqbN2/G/PnzAQAlJSVYs2YNZs6cCb1ej8zMTLzwwgsIDg7G9OnTpePk5uYiNzcXFy9eBACkp6dDp9Ohd+/e6Nq1KwBg4sSJmD59OhYvXiy9zmq1Yvv27Zg3bx68vR1/VJVKhdjYWKxfvx59+/ZF3759sX79evj7++ORRx6R87QQyaYli0BUKhX+78mRcjeJiIjckKwBcMuWLYiLi8OiRYuQl5cHg8GAhQsX4qWXXgJg6w1MT0/Hzp07UVRUBL1ej/HjxyMpKQk6nU46TmJiosPiizFjxgCw9e7ZF3ZcunQJ+fn5DvUfOHAAV65ckQLnrZ5//nmUl5dj0aJFKCwsxPDhw7Fv3z6Huok6kpYsAqmuscIqAF9vXhCIiEhpVEII4epGdFRmsxlBQUEwmUwIDAx0dXOIkHKpAHPf+hq3dQ/Awd+Pa7Ts8cwbmP2XFAwzdmVPIBEpCj+/eS1gIo/i52P7k66oanierV1WQRmsAvDxrn/VOxEReS4GQCIPYh8CtlQ3PQR85UYZAKB31wBZ20RERO6HAZDIg/h52wJgc3oArxSUAgCM3fxlbRMREbkfBkAiD+LnYw+ATfcAZkk9gAyARERKI+sqYCJyLp2fNx4fYYSfjxesVgEvr4bn910tLAcA9OrCAEhEpDQMgEQeJEDjjZd/GdVkOatV4EbtVUC663h1GyIipeEQMJGHOZdjxonMG40OA5srqlBtte0AxauAEBEpDwMgkYd55K2v8avEFGTXzvGrT36Jrfcv0M+bG0ETESkQ3/mJPIy//XrAjVwNJL/EAgDo1onDv0RESsQASORh7JtBlzcyBGzvHezRWeuUNhERkXthACTyMPYewPJGegClTaC5ByARkSIxABJ5GPvVQBrrAcwqsAVAI/cAJCJSJAZAIg+jrd0MurE5gPZNoHkVECIiZWIAJPIw/vYewMrqBsvkF9sWgfA6wEREysSNoIk8jL0HsLEh4CPLxyO/pBKd/X2c1SwiInIjDIBEHsY+B7CxIWCVSsUrgBARKRiHgIk8zM0h4IYDIBERKRsDIJGHac4QMBERKRsDIJGH0TbjSiBERKRsDIBEHkZrvxIIAyARETWAi0CIPIx0JZDaIeAvz+TiurkCo27rhttDdK5sGhERuQn2ABJ5mJurgG37ACYdz8ZL/ziDE5mFrmwWERG5EfYAEnmYXl39MWWwHv1Cbb19+SW2TZ+7deK2L0REZMMASORh7urVGX9+dIh0v6CkEgAQ3MnXVU0iIiI3wwBI5IHyiitgqbKiu06Dn2p7AIPZA0hERLUYAIk80NxtX+PST6V4+zfDUFltBQB0Yw8gERHVYgAk8kDeXrb1XTdKbdf7raq2SquDiYiI+IlA5IG81SoAQEigBmkvTUJ1jdXFLSIiInfCbWCIPJC3ly0A1liF7b6af+pERHQTPxWIPJC6NgBW1QgXt4SIiNwRAyCRB/Kp7fGr5NAvERHVgwGQyAPZrwZSwesBExFRPWQNgNXV1Vi1ahUiIiKg1WrRp08frFu3DlbrzV6JmJgYqFQqh9uIESMcjrNt2zaMGzcOgYGBUKlUKCoqalb9V69exWOPPYZu3brB398fd911F1JTU1tUN1FH5F8bAO3XAyYiIvo5WVcBb9y4EYmJidixYwciIyNx4sQJPPHEEwgKCsKSJUukctHR0di+fbt039fXcb+ysrIyREdHIzo6GitXrmxW3YWFhbjnnnswfvx4fP755wgJCcGlS5fQuXNnh3JN1U3UEWl9bH/aZewBJCKiesgaAFNSUjBt2jRMmTIFABAeHo49e/bgxIkTDuU0Gg3CwsIaPE5sbCwAIDk5udl1b9y4Eb169XIId+Hh4XXKNVU3UUck9QBWVru4JURE5I5kHQIePXo0Dh48iPPnzwMATp8+jSNHjuDBBx90KJecnIyQkBD069cPCxYsQF5eXpvr/vTTTzFs2DDMmjULISEhuPvuu/HWW2/VKdeSui0WC8xms8ONyB1xCJiIiBojaw/g8uXLYTKZMGDAAKjVatTU1CA+Ph5z586VykyePBmzZs2C0WhERkYG4uLiMGHCBKSmpkKjaf21Sy9fvow333wTzz77LF544QV88803eOaZZ6DRaPCb3/ymVXUnJCRg7dq1rW4TkbNMGBCC7joN7ujZ2dVNISIiN6QSQsi2Udh7772HZcuW4dVXX0VkZCTS0tIQGxuLzZs3Y968efW+JicnB0ajEe+99x5mzJjh8FxycjLGjx+PwsLCOnP5buXr64thw4bh6NGj0mPPPPMMjh8/jpSUlBbXDdh6AC0Wi3TfbDajV69eMJlMCAwMbLQ9RERE5B7MZjOCgoIU/fktaw/gsmXLsGLFCsyZMwcAMHjwYGRlZSEhIaHBAKjX62E0GnHhwoU21a3X6zFo0CCHxwYOHIgPP/yw0dc0VrdGo2lTryQRERGRO5B1DmBZWRm8vByrUKvVDtvA3KqgoADZ2dnQ6/Vtqvuee+7BDz/84PDY+fPnYTQaZa+biIiIyJ3JGgCnTp2K+Ph47N27F5mZmfj444+xefNmTJ8+HQBQUlKC5557DikpKcjMzERycjKmTp2K4OBgqQwA5ObmIi0tDRcvXgQApKenIy0tDTdu3JDKTJw4EVu3bpXuL126FF9//TXWr1+PixcvYvfu3di2bRueeuqpFtVNRERE5HGEjMxms1iyZIno3bu38PPzE3369BEvvviisFgsQgghysrKxKRJk0T37t2Fj4+P6N27t5g3b564cuWKw3FWr14tANS5bd++XSpjNBrF6tWrHV73z3/+U0RFRQmNRiMGDBggtm3bJj3X3LobYzKZBABhMplafnKIiIjIJfj5LYSsi0A8HSeREhERdTz8/Oa1gImIiIgUhwGQiIiISGEYAImIiIgUhgGQiIiISGEYAImIiIgUhgGQiIiISGEYAImIiIgUhgGQiIiISGG8Xd2Ajsy+h7bZbHZxS4iIiKi57J/bSr4WBgNgGxQXFwMAevXq5eKWEBERUUsVFxcjKCjI1c1wCV4Krg2sViuuXbsGnU4HlUrl6uZ0CGazGb169UJ2drZiL78jF55befC8yofnVh48r00TQqC4uBgGgwFeXsqcDccewDbw8vJCz549Xd2MDikwMJBvTDLhuZUHz6t8eG7lwfPaOKX2/NkpM/YSERERKRgDIBEREZHCMACSU2k0GqxevRoajcbVTfE4PLfy4HmVD8+tPHheqTm4CISIiIhIYdgDSERERKQwDIBERERECsMASERERKQwDIBERERECsMASM2WkJAAlUqF2NjYep9fuHAhVCoVXn/9dYfHLRYLnn76aQQHByMgIAAPP/wwfvzxR4cyhYWFePzxxxEUFISgoCA8/vjjKCoqcihz5coVTJ06FQEBAQgODsYzzzyDysrKdvwJXaehc3vu3Dk8/PDDCAoKgk6nw4gRI3DlyhXpeZ7bptV3bktKSrB48WL07NkTWq0WAwcOxJtvvunwOp7butasWQOVSuVwCwsLk54XQmDNmjUwGAzQarUYN24czpw543AMnte6GjuvVVVVWL58OQYPHoyAgAAYDAb85je/wbVr1xyOwfNKLSaImuGbb74R4eHh4o477hBLliyp8/zHH38s7rzzTmEwGMQf//hHh+eefPJJ0aNHD7F//35x8uRJMX78eHHnnXeK6upqqUx0dLSIiooSR48eFUePHhVRUVHioYcekp6vrq4WUVFRYvz48eLkyZNi//79wmAwiMWLF8v1IztNQ+f24sWLomvXrmLZsmXi5MmT4tKlS+Jf//qXuH79ulSG57ZxDZ3b3/3ud+K2224Thw4dEhkZGeIvf/mLUKvV4pNPPpHK8NzWtXr1ahEZGSlycnKkW15envT8hg0bhE6nEx9++KFIT08Xs2fPFnq9XpjNZqkMz2tdjZ3XoqIicd9994mkpCTx/fffi5SUFDF8+HAxdOhQh2PwvFJLMQBSk4qLi0Xfvn3F/v37xdixY+sEwB9//FH06NFDfPfdd8JoNDoEwKKiIuHj4yPee+896bGrV68KLy8v8cUXXwghhDh79qwAIL7++mupTEpKigAgvv/+eyGEEJ999pnw8vISV69elcrs2bNHaDQaYTKZZPipnaOxczt79mzx2GOPNfhantvGNXZuIyMjxbp16xzKDxkyRKxatUoIwXPbkNWrV4s777yz3uesVqsICwsTGzZskB6rqKgQQUFBIjExUQjB89qQxs5rfb755hsBQGRlZQkheF6pdTgETE166qmnMGXKFNx33311nrNarXj88cexbNkyREZG1nk+NTUVVVVVmDRpkvSYwWBAVFQUjh49CgBISUlBUFAQhg8fLpUZMWIEgoKCHMpERUXBYDBIZR544AFYLBakpqa228/qbA2dW6vVir1796Jfv3544IEHEBISguHDh+OTTz6RyvDcNq6x39vRo0fj008/xdWrVyGEwKFDh3D+/Hk88MADAHhuG3PhwgUYDAZERERgzpw5uHz5MgAgIyMDubm5DudMo9Fg7Nix0vngeW1YQ+e1PiaTCSqVCp07dwbA80qt4+3qBpB7e++995CamooTJ07U+/zGjRvh7e2NZ555pt7nc3Nz4evriy5dujg8HhoaitzcXKlMSEhIndeGhIQ4lAkNDXV4vkuXLvD19ZXKdDSNndu8vDyUlJRgw4YNeOWVV7Bx40Z88cUXmDFjBg4dOoSxY8fy3Daiqd/bP/3pT1iwYAF69uwJb29veHl54e2338bo0aMB8Pe2IcOHD8fOnTvRr18/XL9+Ha+88gpGjRqFM2fOSD/PrT9vaGgosrKyAPC8NqSx89qtWzeHshUVFVixYgUeeeQRBAYGAuB5pdZhAKQGZWdnY8mSJdi3bx/8/PzqPJ+amor//d//xcmTJ6FSqVp0bCGEw2vqe31rynQUTZ1bq9UKAJg2bRqWLl0KALjrrrtw9OhRJCYmYuzYsQ0em+e28XML2ALg119/jU8//RRGoxH/+c9/sGjRIuj1+np7DO2Ufm4nT54s/f/gwYMxcuRI3HbbbdixYwdGjBgBoO7P25yflee14fP67LPPSs9VVVVhzpw5sFqteOONN5o8rtLPKzWOQ8DUoNTUVOTl5WHo0KHw9vaGt7c3Dh8+jD/96U/w9vZGcnIy8vLy0Lt3b+n5rKws/P73v0d4eDgAICwsDJWVlSgsLHQ4dl5envRNMywsDNevX69T/08//eRQ5tZvoIWFhaiqqqrzjbUjaOrcduvWDd7e3hg0aJDD6wYOHCitAua5rV9T57a0tBQvvPACNm/ejKlTp+KOO+7A4sWLMXv2bLz22msAeG6bKyAgAIMHD8aFCxekVau3/ry3njOe16b9/LzaVVVV4de//jUyMjKwf/9+qfcP4Hml1mEApAZNnDgR6enpSEtLk27Dhg3Do48+irS0NMTExODbb791eN5gMGDZsmX48ssvAQBDhw6Fj48P9u/fLx03JycH3333HUaNGgUAGDlyJEwmE7755hupzLFjx2AymRzKfPfdd8jJyZHK7Nu3DxqNBkOHDnXG6WhXTZ1bjUaDX/ziF/jhhx8cXnf+/HkYjUYAPLcNaerc1tTUoKqqCl5ejm9/arVa6nnluW0ei8WCc+fOQa/XIyIiAmFhYQ7nrLKyEocPH5bOB89r8/z8vAI3w9+FCxdw4MCBOsPCPK/UKk5fdkIdWn2rgH/u1lXAQti2J+jZs6c4cOCAOHnypJgwYUK92xPccccdIiUlRaSkpIjBgwfXuz3BxIkTxcmTJ8WBAwdEz549PWp7glvP7UcffSR8fHzEtm3bxIULF8SWLVuEWq0WX331lVSG57Z5bj23Y8eOFZGRkeLQoUPi8uXLYvv27cLPz0+88cYbUhme27p+//vfi+TkZHH58mXx9ddfi4ceekjodDqRmZkphLBtAxMUFCQ++ugjkZ6eLubOnVvvNjA8r44aO69VVVXi4YcfFj179hRpaWkOW8VYLBbpGDyv1FIMgNQirQmA5eXlYvHixaJr165Cq9WKhx56SFy5csWhTEFBgXj00UeFTqcTOp1OPProo6KwsNChTFZWlpgyZYrQarWia9euYvHixaKioqKdfjLXq+/cvvPOO+L2228Xfn5+4s4773TYp04IntvmuvXc5uTkiJiYGGEwGISfn5/o37+/+MMf/iCsVqtUhue2Lvu+fj4+PsJgMIgZM2aIM2fOSM9brVaxevVqERYWJjQajRgzZoxIT093OAbPa12NndeMjAwBoN7boUOHpGPwvFJLqYQQwpU9kERERETkXJwDSERERKQwDIBERERECsMASERERKQwDIBERERECsMASERERKQwDIBERERECsMASERERKQwDIBEREQN2LZtG8aNG4fAwECoVCoUFRU163VXr17FY489hm7dusHf3x933XUXUlNTpedjYmKgUqkcbiNGjJCez8zMrPO8/fb+++9L5cLDw+s8v2LFCoe2XLlyBVOnTkVAQACCg4PxzDPPoLKy0qFMeno6xo4dC61Wix49emDdunW4dZvgw4cPY+jQofDz80OfPn2QmJhY5+f+8MMPMWjQIGg0GgwaNAgff/xxnTJvvPEGIiIi4Ofnh6FDh+Krr75yeF4IgTVr1sBgMECr1WLcuHE4c+ZMM866oz//+c8YOHAgtFot+vfvj507d7b4GB7NtftQExERudbYsWPF9u3b633uj3/8o0hISBAJCQkCQJ0rZ9Tnxo0bwmg0ipiYGHHs2DGRkZEhDhw4IC5evCiVmTdvnoiOjna4tFtBQYH0fHV1tcNzOTk5Yu3atSIgIEAUFxdL5YxGo1i3bp1DuZ8/b7+82/jx48XJkyfF/v37hcFgcLi8m8lkEqGhoWLOnDkiPT1dfPjhh0Kn04nXXntNKnP58mXh7+8vlixZIs6ePSveeust4ePjIz744AOpzNGjR4VarRbr168X586dE+vXrxfe3t7i66+/lsq89957wsfHR7z11lvi7NmzYsmSJSIgIEBkZWVJZTZs2CB0Op348MMPRXp6unSllJ9fUrApb7zxhtDpdOK9994Tly5dEnv27BGdOnUSn376abOP4ekYAImISNEaC4B2hw4danYAXL58uRg9enSjZebNmyemTZvW/EYKIe666y4xf/58h8fqu/zmz3322WfCy8tLXL16VXpsz549QqPRCJPJJISwhaWgoCCHS74lJCQIg8EgXR7x+eefFwMGDHA49sKFC8WIESOk+7/+9a9FdHS0Q5kHHnhAzJkzR7r/P//zP+LJJ590KDNgwACxYsUKIYTtcoJhYWFiw4YN0vMVFRUiKChIJCYmSo8VFRWJBQsWiO7duwudTifGjx8v0tLSpOdHjhwpnnvuOYd6lixZIu65554Gz5XScAiYiIioHX366acYNmwYZs2ahZCQENx9991466236pRLTk5GSEgI+vXrhwULFiAvL6/BY6ampiItLQ2//e1v6zy3ceNGdOvWDXfddRfi4+MdhndTUlIQFRUFg8EgPfbAAw/AYrFIQ9IpKSkYO3YsNBqNQ5lr164hMzNTKjNp0iSHeh944AGcOHECVVVVjZY5evQoAKCyshKpqal1ykyaNEkqk5GRgdzcXIcyGo0GY8eOlcoIITBlyhTk5ubis88+Q2pqKoYMGYKJEyfixo0bAACLxQI/Pz+HerRaLb755hupvUrHAEhERNSOLl++jDfffBN9+/bFl19+iSeffBLPPPOMwxy0yZMn4+9//zv+/e9/4w9/+AOOHz+OCRMmwGKx1HvMd955BwMHDsSoUaMcHl+yZAnee+89HDp0CIsXL8brr7+ORYsWSc/n5uYiNDTU4TVdunSBr68vcnNzGyxjv99UmerqauTn5zdaxn6M/Px81NTUNFrG/t/Gyhw6dAjp6el4//33MWzYMPTt2xevvfYaOnfujA8++ACALXi+/fbbSE1NhRACJ06cwF//+ldUVVVJ7VU6b1c3gIiIyJnWr1+P9evXS/fLy8vx9ddfY/HixdJjn3/+Oe69995WHd9qtWLYsGFSHXfffTfOnDmDN998E7/5zW8AALNnz5bKR0VFYdiwYTAajdi7dy9mzJjhcLzy8nLs3r0bcXFxdepaunSp9P933HEHunTpgl/96ldSryAAqFSqOq8TQjg8fmsZUbsApD3K3PpYW8ukpqaipKRE+vnsysvLcenSJQBAXFwccnNzMWLECAghEBoaipiYGGzatAlqtRrEAEhERArz5JNP4te//rV0/9FHH8XMmTMdglePHj1afXy9Xo9BgwY5PDZw4EB8+OGHjb7GaDTiwoULdZ774IMPUFZWJoXHxthXEl+8eBHdunVDWFgYjh075lCmsLAQVVVVUi9bWFiY1LtmZx+ObqqMt7e3FMQaKmM/RnBwMNRqdaNlwsLCANh6AvV6fb1lrFYr9Ho9kpOT6/z8nTt3BmAb7v3rX/+Kv/zlL7h+/Tr0ej22bdsGnU6H4ODges+d0nAImIiIFKVr1664/fbbpZtWq0VISEidx1rrnnvuwQ8//ODw2Pnz52E0Ght8TUFBAbKzsx1Cj90777yDhx9+GN27d2+y7lOnTgGAdJyRI0fiu+++Q05OjlRm37590Gg0GDp0qFTmP//5j8PcwX379sFgMCA8PFwqs3//foe69u3bh2HDhsHHx6fRMvZha19fXwwdOrROmf3790tlIiIiEBYW5lCmsrIShw8flsoMGTIEubm58Pb2dvg3u/322+uEOx8fH/Ts2RNqtRrvvfceHnroIXh5MfoA4DYwRESkbI2tAs7JyRGnTp0Sb731lgAg/vOf/4hTp045bNkyYcIEsWXLFun+N998I7y9vUV8fLy4cOGC+Pvf/y78/f3Frl27hBBCFBcXi9///vfi6NGjIiMjQxw6dEiMHDlS9OjRo85WJxcuXBAqlUp8/vnnddp29OhRsXnzZnHq1Clx+fJlkZSUJAwGg3j44YelMvZtYCZOnChOnjwpDhw4IHr27OmwDUxRUZEIDQ0Vc+fOFenp6eKjjz4SgYGB9W4Ds3TpUnH27Fnxzjvv1NkG5r///a9Qq9Viw4YN4ty5c2LDhg0NbgPzzjvviLNnz4rY2FgREBAgMjMzpTIbNmwQQUFB4qOPPhLp6eli7ty5DtvAWK1WMXr0aHHnnXeKL774QmRkZIj//ve/4sUXXxTHjx8XQgjxww8/iHfffVecP39eHDt2TMyePVt07dpVZGRk1PvvrEQMgEREpGiNBcDVq1cLAHVuPy9vNBrF6tWrHV73z3/+U0RFRQmNRiMGDBggtm3bJj1XVlYmJk2aJLp37y58fHxE7969xbx588SVK1fq1L9y5UrRs2dPUVNTU+e51NRUMXz4cBEUFCT8/PxE//79xerVq0VpaalDuaysLDFlyhSh1WpF165dxeLFix22fBFCiG+//Vbce++9QqPRiLCwMLFmzRppCxi75ORkcffddwtfX18RHh4u3nzzzTptev/990X//v2Fj4+PGDBggPjwww/rlPnzn/8sjEaj8PX1FUOGDBGHDx92eN5qtYrVq1eLsLAwodFoxJgxY0R6erpDGbPZLJ5++mlhMBiEj4+P6NWrl3j00Uelc3j27Flx1113Ca1WKwIDA8W0adPE999/X6ctSqYS4patvomIiIjIo3EgnIiIiEhhGACJiIiIFIYBkIiIiEhhGACJiIiIFIYBkIiIiEhhGACJiIiIFIYBkIiIiEhhGACJiIiIFIYBkIiIiEhhGACJiIiIFIYBkIiIiEhhGACJiIiIFIYBkIiIiEhh/n/AU60mEkCmWgAAAABJRU5ErkJggg==", "text/html": [ "\n", "
\n", "
\n", " Figure\n", "
\n", " \n", "
\n", " " ], "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib widget\n", "plt.plot(validation_data_UT['TB-Pegel'])" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# define constants\n", "\n", " # for physics\n", "g = 9.81 # [m/s²] gravitational acceleration \n", "rho = 1000. # [kg/m³] density of water \n", "pUnit_calc = 'Pa' # [string] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", "pUnit_conv = 'mWS' # [string] for pressure conversion in print statements and plot labels\n", "\n", " # for KW OL \n", "OL_T1_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", "OL_T1_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", "OL_T1_closingTime = 10. # [s] closing time of turbine\n", "\n", "OL_T2_Q_nenn = 0.85/2 # [m³/s] nominal flux of turbine \n", "OL_T2_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", "OL_T2_closingTime = 10. # [s] closing time of turbine\n", "\n", " # for KW UL\n", "UL_T1_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", "UL_T1_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", "UL_T1_closingTime = 10. # [s] closing time of turbine\n", "\n", "UL_T2_Q_nenn = 0.85/2 # [m³/s] nominal flux of turbine \n", "UL_T2_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", "UL_T2_closingTime = 10. # [s] closing time of turbine\n", "\n", " # for PI controller\n", "Con_targetLevel = 6. # [m]\n", "Con_K_p = 0.1 # [-] proportional constant of PI controller\n", "Con_T_i = 1000. # [s] timespan in which a steady state error is corrected by the intergal term\n", "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n", "\n", " # for pipeline\n", "Pip_length = (535.+478.) # [m] length of pipeline\n", "Pip_dia = 0.9 # [m] diameter of pipeline\n", "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", "Pip_f_D = 0.014 # [-] Darcy friction factor\n", "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", " # derivatives of the pipeline constants\n", "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", "\n", " # for reservoir\n", "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", "Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n", "Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n", "Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n", "Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n", "Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n", "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", "\n", " # for general simulation\n", "flux_init = (OL_T1_Q_nenn+OL_T2_Q_nenn)/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", "simTime_target = np.max(comp_val_t_vec_UT) # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlCElEQVR4nO3dd3xUZdr/8c9Fr6EHAiGEEkpoCqFZcS0gFsTyLLursuou6urqurtSRFfsqLvr6toerLC6KhIULKiIBQuCQSUNAqEHQm8B0nP//pjD74kYENLOTOb7fr3yypn7nDNz3cyQb845M1fMOYeIiEgtvwsQEZHgoEAQERFAgSAiIh4FgoiIAAoEERHx1PG7gPJq3bq1i42N9bsMEZGQsmzZsp3OuTZlrQvZQIiNjSUpKcnvMkREQoqZbTjaOp0yEhERQIEgIiIeBYKIiAAKBBER8SgQREQEUCCIiIhHgSAiIoACQUQkZOQWFPPQ/BVk7TlUJfcfsh9MExEJJ1+v2cmkxBQ27j5EdItGXDW0U6U/hgJBRCSI7c8r5KH3V/Da0k3EtmrE6+OHMrRLqyp5LAWCiEiQWpC+jTvfTmFHTj7Xn9mF287pToO6tavs8RQIIiJBZueBfKbOS+Pd5Gx6tmvKc1cn0C+6eZU/rgJBRCRIOOeY+8MW7nknjYP5xfzl3O5cf2ZX6tWpnvf//OyjmNmLZrbdzFJLjbU0swVmttr73qLUuslmlmlmGWY2otT4QDNL8dY9YWbmjdc3sze88SVmFlvJcxQRCXpb9uZy3Ywk/vTGD8S2bsx7t5zGH8+Oq7YwgON72+nLwMgjxiYBC51zccBC7zZmFg+MBXp7+zxtZodPeD0DjAfivK/D93kdsMc51w14DHi4vJMREQk1JSWOV77ZwHmPLWLxml387cJ4Zt9wCnFtm1Z7LT97ysg5t6iM39pHA8O95RnAZ8BEb/x151w+sM7MMoHBZrYeiHDOLQYws5nAJcB8b5+p3n3NBp40M3POufJOSkQkFKzbeZCJicksXbeb07q15qFL+9KxZSPf6invNYS2zrlsAOdctplFeuMdgG9KbZfljRV6y0eOH95nk3dfRWa2D2gF7DzyQc1sPIGjDGJiYspZuoiIv4qKS3j+y3U8tmAV9erU4pHL+nFFQjTemXTfVPZF5bJm444xfqx9fjro3HRgOkBCQoKOIEQk5KRv2c/ExGRSNu/jvPi23HdJH9pGNPC7LKD8gbDNzKK8o4MoYLs3ngV0LLVdNLDFG48uY7z0PllmVgdoBuwuZ10iIkEpv6iYJz/J5JnP1tC8UV2e+vUARvVt5/tRQWnlvXw9DxjnLY8D5pYaH+u9c6gzgYvHS73TSzlmNtR7d9HVR+xz+L4uBz7R9QMRqUmWbdjDBU98yb8/yeTik9qz4LYzuaBfVFCFARzHEYKZvUbgAnJrM8sC7gamAbPM7DpgI3AFgHMuzcxmAelAEXCTc67Yu6sbCbxjqSGBi8nzvfEXgP94F6B3E3iXkohIyDtUUMSjH2bw8tfrad+sIS9fM4jhPSJ/fkefWKj+Mp6QkOCSkpL8LkNEpExfrt7JpDnJZO3J5ephnZgwsidN6vv/WWAzW+acSyhrnf/ViYjUIPsOFfLA++nMSsqiS+vGzLp+GIM7t/S7rOOiQBARqSQfpG7lrrmp7D5YwI3Du3Lr2XFV2oyusikQREQqaEdOoBndeynZxEdF8NJvB9GnQzO/yzphCgQRkXJyzjHnu83c+246uQXF3D6iB+PP6ELd2qH5xygVCCIi5bB5by53zEnh81U7GNipBQ9f1o9ukU38LqtCFAgiIiegpMTxypINPDx/JQ645+LeXDW0E7VqBddnCspDgSAicpzW7DjApMRkvl2/h9PjWvPgGH+b0VU2BYKIyM8oLC7huS/W8q+PV9Owbm3+fkV/LhvQIeg+aVxRCgQRkWNI3byPiYnJpG3Zz/l92nHP6N5ENg2OZnSVTYEgIlKGvMJi/v3Jap79fC0tGtXjmd8M4Py+UX6XVaUUCCIiR0hav5sJicms3XGQKwZGM+WCXjRvVM/vsqqcAkFExHMgv4hHP1jJzG820L5ZQ2ZeO5gzurfxu6xqo0AQEQE+X7WDO+aksGVfLuOGxXL7iB40DoJmdNUpvGYrInKEvYcKuO/dFSR+l0XXNo158/phJMSGRjO6yqZAEJGwNT8lm7vmprHnUAE3n9WNm3/RLaSa0VU2BYKIhJ3t+/P429w0PkjbSp8OEcy4dhC924deM7rKpkAQkbDhnGP2sizuezedvKISJo7sye9P70ydEG1GV9kUCCISFjbtPsQdb6XwxeqdDI5tybTL+tKlTWg3o6tsCgQRqdGKSxwzF6/n0Q8zMOC+0b35zZCa0YyusikQRKTGytyew8TEFJZt2MOZ3dvw4KV96dC8od9lBS0FgojUOIXFJfzv52t4YmEmjerX5p//058xJ9e8ZnSVTYEgIjVKStY+JiQmsyJ7Pxf0i2LqRb1p07S+32WFBAWCiNQIeYXF/Ovj1Tz3xVpaNa7H/141kBG92/ldVkhRIIhIyFuydheT5qSwbudBfpnQkTsu6EWzhnX9LivkKBBEJGTl5BXyyAcZ/OebDXRs2ZBXfzeEU7u19ruskKVAEJGQ9GnGdqbMSSF7fx7XntqZv47oTqN6+pFWEfrXE5GQsudgAfe9m86c7zcTF9mExBtPYUBMC7/LqhEUCCISEpxzvJeSzd1z09iXW8gtZ8dx01ldqV8nfJvRVTYFgogEvW3787jz7VQWpG+jX3QzXvndEHpFRfhdVo1ToY5OZnabmaWZWaqZvWZmDcyspZktMLPV3vcWpbafbGaZZpZhZiNKjQ80sxRv3ROmT4+ICIGjgje+3cg5//ycRat2cMeonsy58RSFQRUpdyCYWQfgFiDBOdcHqA2MBSYBC51zccBC7zZmFu+t7w2MBJ42s8PHes8A44E472tkeesSkZph465D/Ob5JUxMTCE+KoIP/3QG48/oqs6kVaiip4zqAA3NrBBoBGwBJgPDvfUzgM+AicBo4HXnXD6wzswygcFmth6IcM4tBjCzmcAlwPwK1iYiIai4xPHy1+v5+4cZ1K5lPDCmD78aFKNmdNWg3IHgnNtsZn8HNgK5wEfOuY/MrK1zLtvbJtvMIr1dOgDflLqLLG+s0Fs+cvwnzGw8gSMJYmJiylu6iASpVdtymDA7mR827eUXPSN5YEwfopqpGV11KXcgeNcGRgOdgb3Am2Z25bF2KWPMHWP8p4POTQemAyQkJJS5jYiEnoKiEp75bA1PfrqaJvXr8PjYk7i4f3s1o6tmFTlldA6wzjm3A8DM5gCnANvMLMo7OogCtnvbZwEdS+0fTeAUU5a3fOS4iISB5Zv2MjExmZVbc7i4f3vuviieVk3UjM4PFbk6sxEYamaNvHcFnQ2sAOYB47xtxgFzveV5wFgzq29mnQlcPF7qnV7KMbOh3v1cXWofEamhcguKefD9FYx5+iv2Hirk+asTeOJXJysMfFSRawhLzGw28B1QBHxP4HROE2CWmV1HIDSu8LZPM7NZQLq3/U3OuWLv7m4EXgYaEriYrAvKIjXY4jW7mDwnmfW7DvGrwTFMHtWTiAZqRuc3cy40T8UnJCS4pKQkv8sQkROwP6+QafNX8t8lG+nUqhEPXdqXU7qqGV11MrNlzrmEstbpk8oiUi0WrtjGlLdS2Z6Tx+9P78yfz+1Bw3pqOxFMFAgiUqV2HcjnnnfSmbd8Cz3aNuXZqwZyUsfmfpclZVAgiEiVcM4xb/kW7nknnZy8Qm47pzs3Du9KvTr6pHGwUiCISKXL3pfLnW+lsnDldvp3bM4jl/WjR7umfpclP0OBICKVpqTE8fq3m3jo/RUUlpRw5wW9uObUztRW24mQoEAQkUqxfudBJs1J5pu1uxnWpRXTLutLp1aN/S5LToACQUQqpKi4hJe+Ws8/FmRQt1Ytpl3al18O6qi2EyFIgSAi5bZy634mzk5medY+zunVlvsv6UO7Zg38LkvKSYEgIicsv6iYpz5dw9OfZtKsYV3+/auTubBflI4KQpwCQUROyPcb9zAxMZlV2w4w5uQO3HVhPC0b1/O7LKkECgQROS6HCor4x0erePGrdbSLaMCLv03gFz3b+l2WVCIFgoj8rK8zdzJpTgobdx/iyqExTBzZk6ZqRlfjKBBE5Kj25Rby0PsreP3bTXRu3ZjXxw9laJdWfpclVUSBICJl+ihtK3e+ncrOA/lcf2YXbjunOw3qqhldTaZAEJEf2Xkgn6nz0ng3OZue7Zry/LgE+kU397ssqQYKBBEBAs3o3v5hM/e8k86h/GL+cm53bhjelbq11YwuXCgQRIQte3OZ8lYKn2bs4OSYQDO6uLZqRhduFAgiYaykxPHq0o08PH8lxSWOv10Yz7hTYtWMLkwpEETC1NodB5iUmMLS9bs5rVtrHrq0Lx1bNvK7LPGRAkEkzBQVl/D8l+t4bMEq6tepxSOX9+OKgdFqOyEKBJFwkr5lPxMSl5O6eT8jerflvtF9iIxQMzoJUCCIhIH8omKe/CSTZz5bQ/NGdXn6NwM4v087HRXIjygQRGq4ZRt2MzExhcztB7h0QAfuuiCeFmpGJ2VQIIjUUAfzi3j0wwxmLF5P+2YNefmaQQzvEel3WRLEFAgiNdAXq3cweU4KWXtyGTesE7eP7EmT+vrvLsemV4hIDbLvUCH3v5fOm8uy6NKmMW/eMIxBsS39LktChAJBpIb4IHUrd81NZffBAv4wvCu3nB2nZnRyQhQIIiFue04eU+el8X7KVuKjInjpt4Po06GZ32VJCFIgiIQo5xyJ323mvnfTyS0s5vYRPRh/Rhc1o5Nyq9Arx8yam9lsM1tpZivMbJiZtTSzBWa22vveotT2k80s08wyzGxEqfGBZpbirXvC9OZokWPK2nOIcS99y1/fXE5cZBPev+V0bjqrm8JAKqSir57HgQ+ccz2B/sAKYBKw0DkXByz0bmNm8cBYoDcwEnjazA6f4HwGGA/EeV8jK1iXSI1UUuKY8fV6zntsEUnrd3PPxb2Zdf0wukU28bs0qQHKfcrIzCKAM4DfAjjnCoACMxsNDPc2mwF8BkwERgOvO+fygXVmlgkMNrP1QIRzbrF3vzOBS4D55a1NpCZas+MAE2cnk7RhD2d0b8ODY/oQ3ULN6KTyVOQaQhdgB/CSmfUHlgG3Am2dc9kAzrlsMzv8SZgOwDel9s/yxgq95SPHf8LMxhM4kiAmJqYCpYuEjsLiEqYvWsvjC1fTsG5t/n5Ffy4b0EFtJ6TSVSQQ6gADgD8655aY2eN4p4eOoqxXrzvG+E8HnZsOTAdISEgocxuRmiR18z4mzE4mPXs/o/q2Y+rFvYlsqmZ0UjUqEghZQJZzbol3ezaBQNhmZlHe0UEUsL3U9h1L7R8NbPHGo8sYFwlbeYXFPL5wNdMXraVFo3o8e+UARvaJ8rssqeHKfVHZObcV2GRmPbyhs4F0YB4wzhsbB8z1lucBY82svpl1JnDxeKl3einHzIZ67y66utQ+ImHn2/W7GfX4Fzzz2RouPbkDC/98psJAqkVFP4fwR+BVM6sHrAWuIRAys8zsOmAjcAWAcy7NzGYRCI0i4CbnXLF3PzcCLwMNCVxM1gVlCTsH8ot45IOVzFy8gegWDfnPdYM5Pa6N32VJGDHnQvNUfEJCgktKSvK7DJFK8fmqHdwxJ4Ut+3IZNyyW20f0oLGa0UkVMLNlzrmEstbpFSfio72HCrj33XTmfLeZrm0aM/uGYQzspGZ04g8FgogPnHPMT93K3+amsvdQITef1Y2bf9FNzejEVwoEkWq2fX8ed81N5cO0bfTpEMGMawfTu72a0Yn/FAgi1cQ5x5vLsrj/3XTyi0qYdH5PfndaZ+qo/5AECQWCSDXYtPsQk+ek8GXmTgbHtmTaZX3p0kb9hyS4KBBEqlBxiWPm4vU88kEGtQzuu6QPvxkcQ61aajshwUeBIFJFMrfnMGF2Mt9t3MvwHm14YExfOjRv6HdZIkelQBCpZIXFJTz72Rr+/UkmjerX5rFf9ueSk9SMToKfAkGkEqVk7eP22ctZuTWHC/tFMfXi3rRuUt/vskSOiwJBpBLkFRbz2MereG7RWlo3qc/0qwZyXu92fpclckIUCCIVtGTtLibNSWHdzoOMHdSRyaN60axhXb/LEjlhCgSRcsrJK+ThD1byyjcb6diyIa/+bgindmvtd1ki5aZAECmHT1du5463Uti6P4/rTuvMX87rTqN6+u8koU2vYJETsPtgAfe+k8bbP2whLrIJiTeewoCYFn6XJVIpFAgix8E5x7vJ2Uydl8a+3EJuOTuOm87qSv06akYnNYcCQeRnbNufx5S3Uvl4xTb6RTfj1d8PoWe7CL/LEql0CgSRo3DO8ca3m3jg/RUUFJUwZVQvrjk1Vs3opMZSIIiUYeOuQ0yak8zXa3YxpHNLHr6sH7GtG/tdlkiVUiCIlFJc4njpq3X8/aMM6tSqxYNj+jJ2UEc1o5OwoEAQ8WRszWFCYjLLN+3lFz0jeWBMH6KaqRmdhA8FgoS9gqISnv4sk6c+zaRpg7o8PvYkLu7fXs3oJOwoECSsLd+0lwmzk8nYlsPok9rztwvjaaVmdBKmFAgSlnILivnnggxe+HIdkU0b8PzVCZwT39bvskR8pUCQsPP1mp1MnpPChl2H+PWQGCad35OIBmpGJ6JAkLCxP6+Qh95fyWtLN9KpVSP++/shnNJVzehEDlMgSFj4OH0bU95OYUdOPuPP6MJt53SnYT21nRApTYEgNdquA/nc804685ZvoWe7pky/KoH+HZv7XZZIUFIgSI3knGPe8i1MnZfGgfwibjunOzcO70q9Omo7IXI0CgSpcbL35XLnW6ksXLmdkzo255HL+9G9bVO/yxIJehX+dcnMapvZ92b2rne7pZktMLPV3vcWpbadbGaZZpZhZiNKjQ80sxRv3ROmTwRJOZSUOF5dsoFz/7mIr9bs5M4LepF44ykKA5HjVBnHz7cCK0rdngQsdM7FAQu925hZPDAW6A2MBJ42s8NX9Z4BxgNx3tfISqhLwsi6nQf51XPfMOWtVPpFN+OjP53J707vQm31IBI5bhUKBDOLBi4Ani81PBqY4S3PAC4pNf66cy7fObcOyAQGm1kUEOGcW+ycc8DMUvuIHFNRcQnTF61h5L8WkZ69n4cv68urvxtCTKtGfpcmEnIqeg3hX8AEoPQxeVvnXDaAcy7bzCK98Q7AN6W2y/LGCr3lI8d/wszGEziSICYmpoKlS6hbkb2fiYnJJGft49z4ttx/SR/aRjTwuyyRkFXuQDCzC4HtzrllZjb8eHYpY8wdY/yng85NB6YDJCQklLmN1Hz5RcU89ekanv40k2YN6/Lkr0/mgr5RakYnUkEVOUI4FbjYzEYBDYAIM3sF2GZmUd7RQRSw3ds+C+hYav9oYIs3Hl3GuMhPfLdxDxNnJ7N6+wHGnNyBv10YT4vG9fwuS6RGKPc1BOfcZOdctHMulsDF4k+cc1cC84Bx3mbjgLne8jxgrJnVN7POBC4eL/VOL+WY2VDv3UVXl9pHBIBDBUXc+046lz3zNQfyi3jpt4N47JcnKQxEKlFVfA5hGjDLzK4DNgJXADjn0sxsFpAOFAE3OeeKvX1uBF4GGgLzvS8RAL7K3MmkOcls2p3LlUNjmDiyJ03VjE6k0lngjT2hJyEhwSUlJfldhlShfbmFPPjeCt5I2kTn1o2ZdmlfhnRp5XdZIiHNzJY55xLKWqdPKktQ+ihtK3e+ncqugwXccGZX/nROHA3qqhmdSFVSIEhQ2ZGTz9R30ngvOZteURG8MG4QfaOb+V2WSFhQIEhQcM7x1vebuffddA7lF/PX87pz/ZldqVtbzehEqosCQXy3eW8uU95K4bOMHQyICTSj6xap/kMi1U2BIL453Ixu2vyVlDi4+6J4rh4Wq/5DIj5RIIgv1u44wKTEFJau383pca15cExfOrZU/yERPykQpFoVFZfw3BfreOzjVTSoU4tHL+/H5QOj1XZCJAgoEKTapG/Zz4TE5aRu3s+I3m25b3QfItWMTiRoKBCkyuUVFvPkJ5k8+/kamjeqxzO/GcD5faP8LktEjqBAkCq1bMNuJsxOZs2Og1w2IJq7LuxF80bqPyQSjBQIUiUO5hfx6IcZzFi8nvbNGjLj2sGc2b2N32WJyDEoEKTSLVq1g8lzUtiyL5erh3bi9pE9aVJfLzWRYKf/pVJp9h0q5L730pm9LIsubRoz6/phDIpt6XdZInKcFAhSKT5IzeauuWnsPljAH4Z35Zaz1YxOJNQoEKRCtufkcffcNOanbiU+KoKXfjuIPh3UjE4kFCkQpFycc8xelsX9760gt7CY20f0YPwZXdSMTiSEKRDkhG3afYg73krhi9U7SejUgmmX9aNbZBO/yxKRClIgyHErKXHMXLyeRz7MwIB7R/fmyiGdqKVmdCI1ggJBjkvm9gNMSkwmacMezujehgfH9CG6hZrRidQkCgQ5psLiEqYvWsvjH6+mYb3a/OOK/lw6oIOa0YnUQAoEOarUzfuYMDuZ9Oz9jOrbjnsu7kObpvX9LktEqogCQX4ir7CYxxeuZvqitbRsXI9nrxzAyD5qRidS0ykQ5Ee+Xb+bibOTWbvzIP+TEM2UUfE0a1TX77JEpBooEASAA/lFPPLBSmYu3kB0i4a8ct0QTotr7XdZIlKNFAjCpxnbmTInhez9eVxzaix/Pa8HjdWMTiTs6H99GNtzsID73k1nzveb6RbZhNk3nMLATi38LktEfKJACEPOOd5P2crd81LZe6iQP/6iGzf/ohv166gZnUg4UyCEme3787jz7VQ+St9G3w7NmHntEOLbR/hdlogEAQVCmHDO8WZSFve9l05BUQmTz+/Jdad1po6a0YmIp9w/Dcyso5l9amYrzCzNzG71xlua2QIzW+19b1Fqn8lmlmlmGWY2otT4QDNL8dY9YfoYbKXatPsQV72wlAmJyfSKimD+radz/ZldFQYi8iMV+YlQBPzFOdcLGArcZGbxwCRgoXMuDljo3cZbNxboDYwEnjazwyetnwHGA3He18gK1CWe4hLHi1+u47zHFvHDpr3cf0kfXv/9ULq0UWdSEfmpcp8ycs5lA9neco6ZrQA6AKOB4d5mM4DPgIne+OvOuXxgnZllAoPNbD0Q4ZxbDGBmM4FLgPnlrU1g9bYcJiQm8/3GvQzv0YYHx/SlffOGfpclIkGsUq4hmFkscDKwBGjrhQXOuWwzi/Q26wB8U2q3LG+s0Fs+crysxxlP4EiCmJiYyii9xikoKuHZz9fw5CeZNK5fm3/98iRGn9RezehE5GdVOBDMrAmQCPzJObf/GD94ylrhjjH+00HnpgPTARISEsrcJpwlZ+1lwuxkVm7N4aL+7bn7onhaN1EzOhE5PhUKBDOrSyAMXnXOzfGGt5lZlHd0EAVs98azgI6ldo8Gtnjj0WWMy3HKKyzmsQWreO6LtbRpWp/nrk7g3Pi2fpclIiGmIu8yMuAFYIVz7p+lVs0DxnnL44C5pcbHmll9M+tM4OLxUu/0Uo6ZDfXu8+pS+8jP+GbtLkb+axH/u2gtvxzUkY9uO1NhICLlUpEjhFOBq4AUM/vBG7sDmAbMMrPrgI3AFQDOuTQzmwWkE3iH0k3OuWJvvxuBl4GGBC4m64Lyz8jJK2Ta/JW8umQjMS0b8d/fDeGUbmpGJyLlZ86F5qn4hIQEl5SU5HcZvvhk5TamvJXKtv15XHtqZ/58Xnca1dNnDEXk55nZMudcQlnr9FMkhOw+WMC976Tx9g9biItswtM3nsLJMWpGJyKVQ4EQApxzvJOczdR5aeTkFXLr2XH84ayuakYnIpVKgRDktu4LNKP7eMU2+kc34+HLh9CznZrRiUjlUyAEKeccr3+7iQffW0FhSQlTRvXi2tM6U7uWPmAmIlVDgRCENuw6yKTEFBav3cXQLi2Zdmk/Yls39rssEanhFAhBpLjE8dJX6/j7RxnUrVWLB8f0ZeygjtTSUYGIVAMFQpDI2BpoRrd8017O7hnJ/WP6ENVMzehEpPooEHxWUFTC059l8tSnmTRtUJcnfnUyF/WLUjM6Eal2CgQf/bBpLxNnJ5OxLYfRJ7Xn7ot607JxPb/LEpEwpUDwQW5BMf/4KIMXv1pHZNMGvDAugbN7qf+QiPhLgVDNvl6zk0mJKWzcfYhfD4lh0vk9iWhQ1++yREQUCNVlf14hD72/gteWbqJTq0a89vuhDOvayu+yRET+PwVCNfg4fRtT3k5hR04+48/owm3ndKdhPbWdEJHgokCoQrsO5DP1nXTeWb6Fnu2aMv2qBPp3bO53WSIiZVIgVAHnHHN/2MI976RxIL+IP5/bnRvO7Eq9OuX+e0QiIlVOgVDJtuzN5c63U/lk5XZO6ticRy7vR/e2Tf0uS0TkZykQKklJieO/Szcybf5Kikscd10Yz29PiVUzOhEJGQqESrBu50EmJSazZN1uTu3WiofG9COmVSO/yxIROSEKhAooKi7hhS/X8c8Fq6hXpxYPX9aX/0noqLYTIhKSFAjltCJ7PxMTk0nO2se58W25/5I+tI1o4HdZIiLlpkA4QflFxTz1SSZPf7aG5o3q8tSvBzCqbzsdFYhIyFMgnIBlG/YwMTGZzO0HuPTkDtx1YTwt1IxORGoIBcJxOFRQxKMfZvDy1+uJimjAS9cM4qwekX6XJSJSqRQIP+PL1TuZNCeZrD25XDW0ExNG9qCpmtGJSA2kQDiKfbmFPPBeOrOSsujcujFvjB/KkC5qRiciNZcCoQwfpm3lrrdT2XWwgBuHd+XWs+NoUFfN6ESkZlMglLIjJ5+p89J4LyWbXlERvDBuEH2jm/ldlohItVAgEGhGN+e7zdz7bjq5BcXcPqIH48/oQt3aakYnIuEj7ANh895c7piTwuerdjAgJtCMrlukmtGJSPgJmkAws5HA40Bt4Hnn3LSqfLySEscrSzbw8PyVOGDqRfFcNUzN6EQkfAVFIJhZbeAp4FwgC/jWzOY559Kr4vHW7DjApMRkvl2/h9PjWvPgmL50bKlmdCIS3oIiEIDBQKZzbi2Amb0OjAYqPRBmfbuJO+em0qBOLR69vB+XD4xW2wkREYInEDoAm0rdzgKGHLmRmY0HxgPExMSU64E6t2nM2T0juWd0byKbqhmdiMhhwRIIZf2K7n4y4Nx0YDpAQkLCT9Yfj0GxLRkU27I8u4qI1GjB8r7KLKBjqdvRwBafahERCUvBEgjfAnFm1tnM6gFjgXk+1yQiElaC4pSRc67IzG4GPiTwttMXnXNpPpclIhJWgiIQAJxz7wPv+12HiEi4CpZTRiIi4jMFgoiIAAoEERHxKBBERAQAc65cn+/ynZntADaUc/fWwM5KLCeYhdNcIbzmq7nWTFU9107OuTZlrQjZQKgIM0tyziX4XUd1CKe5QnjNV3Otmfycq04ZiYgIoEAQERFPuAbCdL8LqEbhNFcIr/lqrjWTb3MNy2sIIiLyU+F6hCAiIkdQIIiICBCGgWBmI80sw8wyzWyS3/UcLzNbb2YpZvaDmSV5Yy3NbIGZrfa+tyi1/WRvjhlmNqLU+EDvfjLN7Anz/n6omdU3sze88SVmFlvN83vRzLabWWqpsWqZn5mN8x5jtZmN82muU81ss/f8/mBmo0J9rmbW0cw+NbMVZpZmZrd64zX1eT3afEPnuXXOhc0Xgdbaa4AuQD1gORDvd13HWft6oPURY48Ak7zlScDD3nK8N7f6QGdvzrW9dUuBYQT+St184Hxv/A/As97yWOCNap7fGcAAILU65we0BNZ631t4yy18mOtU4K9lbBuycwWigAHeclNglTefmvq8Hm2+IfPchtsRwmAg0zm31jlXALwOjPa5pooYDczwlmcAl5Qaf905l++cWwdkAoPNLAqIcM4tdoFX0cwj9jl8X7OBsw//VlIdnHOLgN1HDFfH/EYAC5xzu51ze4AFwMjKnl9pR5nr0YTsXJ1z2c6577zlHGAFgb+fXlOf16PN92iCbr7hFggdgE2lbmdx7CcsmDjgIzNbZmbjvbG2zrlsCLwYgUhv/Gjz7OAtHzn+o32cc0XAPqBVFczjRFTH/ILpNXGzmSV7p5QOn0apEXP1Tm2cDCwhDJ7XI+YLIfLchlsglPUbb6i87/ZU59wA4HzgJjM74xjbHm2ex5p/KP3bVOb8gmXezwBdgZOAbOAf3njIz9XMmgCJwJ+cc/uPtWkZYyE1VyhzviHz3IZbIGQBHUvdjga2+FTLCXHObfG+bwfeInD6a5t3eIn3fbu3+dHmmeUtHzn+o33MrA7QjOM/rVFVqmN+QfGacM5tc84VO+dKgOcIPL8co76QmKuZ1SXww/FV59wcb7jGPq9lzTekntuqvMgSbF8E/mToWgIXcA5fVO7td13HUXdjoGmp5a8JnB98lB9fnHvEW+7Njy9WreX/LlZ9Cwzl/y5WjfLGb+LHF6tm+TDPWH58obXK50fgItw6AhfiWnjLLX2Ya1Sp5dsInFsO6bl6dc0E/nXEeI18Xo8x35B5bqv1P3wwfAGjCFz9XwNM8bue46y5i/fCWQ6kHa6bwLnDhcBq73vLUvtM8eaYgfcOBW88AUj11j3J/31avQHwJoELW0uBLtU8x9cIHE4XEvht57rqmh9wrTeeCVzj01z/A6QAycC8I36IhORcgdMInLZIBn7wvkbV4Of1aPMNmedWrStERAQIv2sIIiJyFAoEEREBFAgiIuJRIIiICKBAEBERjwJBREQABYKIiHj+H7Fc8jQ2Uq9RAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(t_vec)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# import validation data\n", "val_OL_M1_p = interpolate_validation_data(comp_val_t_vec_TB,comp_val_OL_M1_p,t_vec)\n", "val_OL_M1_LA = interpolate_validation_data(comp_val_t_vec_TB,comp_val_OL_M1_LA,t_vec)\n", "val_OL_M2_p = interpolate_validation_data(comp_val_t_vec_TB,comp_val_OL_M2_p,t_vec)\n", "val_OL_M2_LA = interpolate_validation_data(comp_val_t_vec_TB,comp_val_OL_M2_LA,t_vec)\n", "val_OL_level = interpolate_validation_data(comp_val_t_vec_UT,comp_val_OL_level,t_vec)\n", "val_UL_M1_p = interpolate_validation_data(comp_val_t_vec_UT,comp_val_UL_M1_p,t_vec)\n", "val_UL_M1_LA = interpolate_validation_data(comp_val_t_vec_UT,comp_val_UL_M1_LA,t_vec)\n", "val_UL_M2_p = interpolate_validation_data(comp_val_t_vec_UT,comp_val_UL_M2_p,t_vec)\n", "val_UL_M2_LA = interpolate_validation_data(comp_val_t_vec_UT,comp_val_UL_M2_LA,t_vec)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# create objects\n", "\n", "# downstream turbines\n", "UL_T1 = Francis_Turbine(UL_T1_Q_nenn,UL_T1_p_nenn,UL_T1_closingTime,Pip_dt,pUnit_conv)\n", "UL_T2 = Francis_Turbine(UL_T2_Q_nenn,UL_T2_p_nenn,UL_T2_closingTime,Pip_dt,pUnit_conv)\n", "\n", "KW_UL = Kraftwerk_class()\n", "KW_UL.add_turbine(UL_T1)\n", "KW_UL.add_turbine(UL_T2)\n", "\n", "KW_UL.set_steady_state_by_LA(LA_vec=[val_UL_M1_LA[0],val_UL_M2_LA[0]],ss_pressure=np.max([val_UL_M1_p,val_UL_M2_p])) \n", "flux_init = KW_UL.get_current_Q()\n", "\n", "# Upstream reservoir\n", "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,pUnit_conv,Res_level_crit_lo,Res_level_crit_hi,rho)\n", "reservoir.set_steady_state(flux_init,level_init)\n", "\n", "# pipeline\n", "pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_head,Pip_n_seg,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n", "pipe.set_steady_state(flux_init,reservoir.get_current_pressure())\n", "\n", "# influx setting turbines\n", "OL_T1 = Francis_Turbine(OL_T1_Q_nenn,OL_T1_p_nenn,OL_T1_closingTime,Pip_dt,pUnit_conv)\n", "OL_T2 = Francis_Turbine(OL_T2_Q_nenn,OL_T2_p_nenn,OL_T2_closingTime,Pip_dt,pUnit_conv)\n", "\n", "KW_OL = Kraftwerk_class()\n", "KW_OL.add_turbine(OL_T1)\n", "KW_OL.add_turbine(OL_T2)\n", "\n", "KW_OL.set_steady_state_by_LA(LA_vec=[val_OL_M1_LA[0],val_OL_M2_LA[0]],ss_pressure=np.max([val_OL_M1_p,val_OL_M2_p])) \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# initialization for Timeloop\n", "\n", "# pipeline\n", "v_old = pipe.get_current_velocity_distribution() # storing the velocity from the last timestep\n", "v_min = pipe.get_lowest_velocity_per_node() # storing minimal flux velocity at each node\n", "v_max = pipe.get_highest_velocity_per_node() # storing maximal flux velocity at each node\n", "Q_old = pipe.get_current_flux_distribution() # storing the flux from the last timestep\n", "Q_min = pipe.get_lowest_flux_per_node() # storing minimal flux at each node\n", "Q_max = pipe.get_highest_flux_per_node() # storing maximal flux at each node\n", "p_old = pipe.get_current_pressure_distribution() # storing the pressure from the last timestep\n", "p_min = pipe.get_lowest_pressure_per_node() # storing minimal pressure at each node\n", "p_max = pipe.get_highest_pressure_per_node() # storing maximal pressure at each node\n", "p_0 = pipe.get_initial_pressure_distribution() # storing initial pressure at each node\n", "\n", "v_boundary_res = np.zeros_like(t_vec) # storing the boundary velocity at the reservoir\n", "v_boundary_tur = np.zeros_like(t_vec) # storing the boundary velocity at the turbine\n", "Q_boundary_res = np.zeros_like(t_vec) # storing the boundary flux at the reservoir\n", "Q_boundary_tur = np.zeros_like(t_vec) # storing the boundary flux at the turbine\n", "p_boundary_res = np.zeros_like(t_vec) # storing the boundary pressure at the reservoir\n", "p_boundary_tur = np.zeros_like(t_vec) # storing the boundary pressure at the turbine\n", "\n", "v_boundary_res[0] = v_old[0] # storing the initial value for the boundary velocity at the reservoir\n", "v_boundary_tur[0] = v_old[-1] # storing the initial value for the boundary velocity at the turbine\n", "Q_boundary_res[0] = Q_old[0] # storing the initial value for the boundary flux at the reservoir\n", "Q_boundary_tur[0] = Q_old[-1] # storing the initial value for the boundary flux at the turbine\n", "p_boundary_res[0] = p_old[0] # storing the initial value for the boundary pressure at the reservoir\n", "p_boundary_tur[0] = p_old[-1] # storing the initial value for the boundary pressure at the turbine\n", "\n", "# reservoir\n", "Q_in_vec = np.zeros_like(t_vec) # storing the influx to the reservoir\n", "Q_in_vec[0] = flux_init # storing the initial influx to the reservoir\n", "# Outflux from reservoir is stored in Q_boundary_res\n", "level_vec = np.zeros_like(t_vec) # storing the level in the reservoir at the end of each pipeline timestep\n", "level_vec[0] = level_init # storing the initial level in the reservoir\n", "volume_vec = np.zeros_like(t_vec) # storing the volume in the reservoir at the end of each pipeline timestep\n", "volume_vec[0] = reservoir.get_current_volume() # storing the initial volume in the reservoir\n", "\n", "# OL KW\n", " # manual input to modulate influx\n", "OL_T1_LA_soll_vec = val_OL_M1_LA\n", "\n", "OL_T2_LA_soll_vec = val_OL_M2_LA # storing the target value for the guide van opening\n", "\n", "OL_T1_LA_ist_vec = np.zeros_like(t_vec) # storing the actual value of the guide vane opening\n", "OL_T1_LA_ist_vec[0] = OL_T1.get_current_LA() # storing the initial value of the guide vane opening\n", "\n", "OL_T2_LA_ist_vec = np.zeros_like(t_vec) # storing the actual value of the guide vane opening\n", "OL_T2_LA_ist_vec[0] = OL_T2.get_current_LA() # storing the initial value of the guide vane opening\n", "\n", "# UL KW\n", "UL_T1_LA_soll_vec = val_UL_M1_LA # storing the target value of the guide vane opening\n", "UL_T1_LA_soll_vec[0] = UL_T1.get_current_LA() # storing the initial value of the guide vane opening\n", "\n", "UL_T2_LA_soll_vec = val_UL_M2_LA # storing the target value of the guide vane opening\n", "\n", "UL_T1_LA_ist_vec = np.zeros_like(t_vec) # storing the actual value of the guide vane opening\n", "UL_T1_LA_ist_vec[0] = UL_T1.get_current_LA() # storing the initial value of the guide vane opening\n", "\n", "UL_T2_LA_ist_vec = np.zeros_like(t_vec) # storing the actual value of the guide vane opening\n", "UL_T2_LA_ist_vec[0] = UL_T2.get_current_LA() # storing the initial value of the guide vane opening\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib qt5\n", "\n", "\n", "fig1,axs1 = plt.subplots(3,1)\n", "fig1.suptitle(str(0) +' s / '+str(round(t_vec[-1],2)) + ' s' )\n", "axs1[0].set_title('Pressure distribution in pipeline')\n", "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[0].set_ylabel(r'$p$ ['+pUnit_conv+']')\n", "axs1[1].set_title('Pressure distribution in pipeline \\n Difference to t=0')\n", "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[1].set_ylabel(r'$p$ ['+pUnit_conv+']')\n", "axs1[2].set_title('Flux distribution in pipeline')\n", "axs1[2].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[2].set_ylabel(r'$Q$ [$\\mathrm{m}^3 / \\mathrm{s}$]')\n", "lo_0, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,pUnit_calc, pUnit_conv),marker='.')\n", "lo_0min, = axs1[0].plot(Pip_x_vec,pressure_conversion(pipe.get_lowest_pressure_per_node(),pUnit_calc,pUnit_conv),c='red')\n", "lo_0max, = axs1[0].plot(Pip_x_vec,pressure_conversion(pipe.get_highest_pressure_per_node(),pUnit_calc,pUnit_conv),c='red')\n", "lo_1, = axs1[1].plot(Pip_x_vec,pressure_conversion(p_old-p_0,pUnit_calc, pUnit_conv),marker='.')\n", "lo_1min, = axs1[1].plot(Pip_x_vec,pressure_conversion(pipe.get_lowest_pressure_per_node()-p_0,pUnit_calc,pUnit_conv),c='red')\n", "lo_1max, = axs1[1].plot(Pip_x_vec,pressure_conversion(pipe.get_highest_pressure_per_node()-p_0,pUnit_calc,pUnit_conv),c='red')\n", "lo_2, = axs1[1].plot(Pip_x_vec,Q_old,marker='.')\n", "lo_2min, = axs1[2].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", "lo_2max, = axs1[2].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", "\n", "# axs1[0].autoscale()\n", "# axs1[1].autoscale()\n", "\n", "fig1.tight_layout()\n", "fig1.show()\n", "plt.pause(1)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# needed for turbine convergence\n", "convergence_parameters = [p_old[-2],v_old[-2],Pip_dia,Pip_area,Pip_angle,Pip_f_D,Pip_pw_vel,rho,Pip_dt,p_old[-1]]\n", "\n", "# loop through time steps of the pipeline\n", "for it_pipe in range(1,nt+1):\n", "\n", " KW_OL.update_LAs([OL_T1_LA_soll_vec[it_pipe],OL_T2_LA_soll_vec[it_pipe]])\n", " KW_OL.set_pressure(OL_T1_p_nenn)\n", " Q_in_vec[it_pipe] = KW_OL.get_current_Q()\n", " reservoir.set_influx(Q_in_vec[it_pipe])\n", "\n", "# for each pipeline timestep, execute Res_nt timesteps of the reservoir code\n", " # set initial condition for the reservoir time evolution calculted with the timestep_reservoir_evolution() method\n", " reservoir.set_pressure(p_old[0],display_warning=False)\n", " reservoir.set_outflux(Q_old[0],display_warning=False)\n", " # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n", " for it_res in range(Res_nt):\n", " reservoir.timestep_reservoir_evolution() \n", " level_vec[it_pipe] = reservoir.get_current_level() \n", " volume_vec[it_pipe] = reservoir.get_current_volume() \n", " \n", " # change the guide vane opening based on the target value and closing time limitation\n", " KW_UL.update_LAs([UL_T1_LA_soll_vec[it_pipe],UL_T2_LA_soll_vec[it_pipe]])\n", " OL_T1_LA_ist_vec[it_pipe], OL_T2_LA_ist_vec[it_pipe] = KW_OL.get_current_LAs()\n", " UL_T1_LA_ist_vec[it_pipe], UL_T2_LA_ist_vec[it_pipe] = KW_UL.get_current_LAs()\n", "\n", " # set boundary condition for the next timestep of the characteristic method\n", " convergence_parameters[0] = p_old[-2]\n", " convergence_parameters[1] = v_old[-2]\n", " convergence_parameters[9] = p_old[-1]\n", " KW_UL.set_pressure(p_old[-1])\n", " KW_UL.converge(convergence_parameters)\n", " p_boundary_res[it_pipe] = reservoir.get_current_pressure()\n", " v_boundary_tur[it_pipe] = 1/Pip_area*KW_UL.get_current_Q()\n", " Q_boundary_tur[it_pipe] = KW_UL.get_current_Q()\n", "\n", " # the the boundary condition in the pipe.object and thereby calculate boundary pressure at turbine\n", " pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n", " # pipe.v[0] = (0.8*pipe.v[0]+0.2*reservoir.get_current_outflux()/Res_area_out) # unnecessary\n", " p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n", " v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n", " Q_boundary_res[it_pipe] = pipe.get_current_flux_distribution()[0]\n", "\n", " # perform the next timestep via the characteristic method\n", " pipe.timestep_characteristic_method_vectorized()\n", "\n", " # prepare for next loop\n", " p_old = pipe.get_current_pressure_distribution()\n", " v_old = pipe.get_current_velocity_distribution()\n", " Q_old = pipe.get_current_flux_distribution()\n", "\n", " # plot some stuff\n", " # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n", " if it_pipe%25 == 0:\n", " lo_0.remove()\n", " lo_0min.remove()\n", " lo_0max.remove()\n", " lo_1.remove()\n", " lo_1min.remove()\n", " lo_1max.remove()\n", " lo_2.remove()\n", " lo_2min.remove()\n", " lo_2max.remove()\n", " # plot new pressure and velocity distribution in the pipeline\n", " lo_0, = axs1[0].plot(Pip_x_vec,pressure_conversion(pipe.get_current_pressure_distribution(),pUnit_calc,pUnit_conv),marker='.',c='blue')\n", " lo_0min, = axs1[0].plot(Pip_x_vec,pressure_conversion(pipe.get_lowest_pressure_per_node(),pUnit_calc,pUnit_conv),c='red')\n", " lo_0max, = axs1[0].plot(Pip_x_vec,pressure_conversion(pipe.get_highest_pressure_per_node(),pUnit_calc,pUnit_conv),c='red') \n", " lo_1, = axs1[1].plot(Pip_x_vec,pressure_conversion(pipe.get_current_pressure_distribution()-p_0,pUnit_calc,pUnit_conv),marker='.',c='blue')\n", " lo_1min, = axs1[1].plot(Pip_x_vec,pressure_conversion(pipe.get_lowest_pressure_per_node()-p_0,pUnit_calc,pUnit_conv),c='red')\n", " lo_1max, = axs1[1].plot(Pip_x_vec,pressure_conversion(pipe.get_highest_pressure_per_node()-p_0,pUnit_calc,pUnit_conv),c='red')\n", " lo_2, = axs1[2].plot(Pip_x_vec,pipe.get_current_flux_distribution(),marker='.',c='blue')\n", " lo_2min, = axs1[2].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", " lo_2max, = axs1[2].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", " fig1.suptitle(str(round(t_vec[it_pipe],2))+ ' s / '+str(round(t_vec[-1],2)) + ' s' )\n", " fig1.canvas.draw()\n", " fig1.tight_layout()\n", " fig1.show()\n", " plt.pause(0.00000001) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig2,axs2 = plt.subplots(1,1)\n", "axs2.set_title('Level and Volume reservoir')\n", "axs2.plot(t_vec,level_vec,label='level')\n", "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2.set_ylabel(r'$h$ [m]')\n", "x_twin_00 = axs2.twinx()\n", "x_twin_00.set_ylabel(r'$V$ [$\\mathrm{m}^3$]')\n", "x_twin_00.plot(t_vec,volume_vec)\n", "axs2.legend()\n", "\n", "fig2,axs2 = plt.subplots(1,1)\n", "axs2.set_title('LA')\n", "axs2.plot(t_vec,100*OL_T1_LA_soll_vec,label='OL_T1 Target',c='b')\n", "axs2.scatter(t_vec[::200],100*OL_T1_LA_ist_vec[::200],label='OL_T1 Actual',c='b',marker='+')\n", "axs2.plot(t_vec,100*OL_T2_LA_soll_vec,label='OL_T2 Target',c='g')\n", "axs2.scatter(t_vec[::200],100*OL_T2_LA_ist_vec[::200],label='OL_T2 Actual',c='g',marker='+')\n", "axs2.plot(t_vec,100*UL_T1_LA_soll_vec,label='UL_T1 Target',c='r')\n", "axs2.scatter(t_vec[::200],100*UL_T1_LA_ist_vec[::200],label='UL_T1 Actual',c='r',marker='+')\n", "axs2.plot(t_vec,100*UL_T2_LA_soll_vec,label='UL_T2 Target',c='k')\n", "axs2.scatter(t_vec[::200],100*UL_T2_LA_ist_vec[::200],label='UL_T2 Actual',c='k',marker='+')\n", "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2.set_ylabel(r'$LA$ [%]')\n", "axs2.legend()\n", "\n", "fig2,axs2 = plt.subplots(1,1)\n", "axs2.set_title('Pressure change vs t=0 at reservoir and turbine')\n", "axs2.plot(t_vec,pressure_conversion(p_boundary_res-p_boundary_res[0],pUnit_calc, pUnit_conv),label='Reservoir')\n", "axs2.plot(t_vec,pressure_conversion(p_boundary_tur-p_boundary_tur[0],pUnit_calc, pUnit_conv),label='Turbine')\n", "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n", "axs2.legend()\n", "\n", "fig2,axs2 = plt.subplots(1,1)\n", "axs2.set_title('Fluxes')\n", "axs2.plot(t_vec,Q_in_vec,label='Influx')\n", "axs2.plot(t_vec,Q_boundary_res,label='Outflux')\n", "axs2.scatter(t_vec[::200],Q_boundary_tur[::200],label='Flux Turbine',c='g',marker='+')\n", "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n", "axs2.legend()\n", "\n", "# fig2,axs2 = plt.subplots(1,1)\n", "# axs2.set_title('Min and Max Pressure')\n", "# axs2.plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n", "# axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n", "# axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "# axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n", "\n", "# fig2,axs2 = plt.subplots(1,1)\n", "# axs2.set_title('Min and Max Fluxes')\n", "# axs2.plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", "# axs2.plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", "# axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "# axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n", "\n", "\n", "fig2.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig3,axs3 = plt.subplots(2,2)\n", "axs3[0,0].set_title('Level and Volume reservoir')\n", "axs3[0,0].plot(t_vec,level_vec,label='level')\n", "axs3[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs3[0,0].set_ylabel(r'$h$ [m]')\n", "x_twin_00 = axs3[0,0].twinx()\n", "x_twin_00.set_ylabel(r'$V$ [$\\mathrm{m}^3$]')\n", "x_twin_00.plot(t_vec,volume_vec)\n", "axs3[0,0].legend()\n", "\n", "axs3[0,1].set_title('LA')\n", "axs3[0,1].plot(t_vec,100*OL_T1_LA_soll_vec,label='OL_T1 Target',c='b')\n", "axs3[0,1].scatter(t_vec[::200],100*OL_T1_LA_ist_vec[::200],label='OL_T1 Actual',c='b',marker='+')\n", "axs3[0,1].plot(t_vec,100*OL_T2_LA_soll_vec,label='OL_T2 Target',c='g')\n", "axs3[0,1].scatter(t_vec[::200],100*OL_T2_LA_ist_vec[::200],label='OL_T2 Actual',c='g',marker='+')\n", "axs3[0,1].plot(t_vec,100*UL_T1_LA_soll_vec,label='UL_T1 Target',c='r')\n", "axs3[0,1].scatter(t_vec[::200],100*UL_T1_LA_ist_vec[::200],label='UL_T1 Actual',c='r',marker='+')\n", "axs3[0,1].plot(t_vec,100*UL_T2_LA_soll_vec,label='UL_T2 Target',c='k')\n", "axs3[0,1].scatter(t_vec[::200],100*UL_T2_LA_ist_vec[::200],label='UL_T2 Actual',c='k',marker='+')\n", "axs3[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs3[0,1].set_ylabel(r'$LA$ [%]')\n", "axs3[0,1].legend()\n", "\n", "axs3[1,0].set_title('Fluxes')\n", "axs3[1,0].plot(t_vec,Q_in_vec,label='Influx')\n", "axs3[1,0].plot(t_vec,Q_boundary_res,label='Outflux')\n", "axs3[1,0].scatter(t_vec[::200],Q_boundary_tur[::200],label='Flux Turbine',c='g',marker='+')\n", "axs3[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs3[1,0].set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n", "axs3[1,0].legend()\n", "\n", "axs3[1,1].set_title('Pressure change vs t=0 at reservoir and turbine')\n", "axs3[1,1].plot(t_vec,pressure_conversion(p_boundary_res-p_boundary_res[0],pUnit_calc, pUnit_conv),label='Reservoir')\n", "axs3[1,1].plot(t_vec,pressure_conversion(p_boundary_tur-p_boundary_tur[0],pUnit_calc, pUnit_conv),label='Turbine')\n", "axs3[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs3[1,1].set_ylabel(r'$p$ ['+pUnit_conv+']')\n", "axs3[1,1].legend()\n", "\n", "fig3.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "figure = plt.figure()\n", "plt.plot(t_vec,pressure_conversion(p_boundary_tur,'Pa','mWS'))\n", "plt.plot(t_vec,pressure_conversion(val_p_vec,'Pa','mWS'),marker='+')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.8.13 ('DT_Slot_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.13" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48" } } }, "nbformat": 4, "nbformat_minor": 2 }