gen_load_pv.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'PF_ALG': 1, 'PF_TOL': 1e-08, 'PF_MAX_IT': 10, 'PF_MAX_IT_FD': 30, 'PF_MAX_IT_GS': 1000, 'ENFORCE_Q_LIMS': False, 'PF_DC': False, 'CPF_PARAMETERIZATION': 3, 'CPF_STOP_AT': 'NOSE', 'CPF_STEP': 0.05, 'CPF_ADAPT_STEP': False, 'CPF_ERROR_TOL': 0.001, 'CPF_STEP_MIN': 0.0001, 'CPF_STEP_MAX': 0.2, 'CPF_PLOT_LEVEL': 0, 'CPF_PLOT_BUS': '', 'CPF_USER_CALLBACK': '', 'CPF_USER_CALLBACK_ARGS': '', 'OPF_ALG': 0, 'OPF_VIOLATION': 0.0005, 'OPF_FLOW_LIM': 1, 'OPF_IGNORE_ANG_LIM': False, 'OPF_ALG_DC': 0, 'VERBOSE': 0, 'OUT_ALL': 0, 'OUT_SYS_SUM': True, 'OUT_AREA_SUM': False, 'OUT_BUS': True, 'OUT_BRANCH': True, 'OUT_GEN': False, 'OUT_ALL_LIM': -1, 'OUT_V_LIM': 1, 'OUT_LINE_LIM': 1, 'OUT_PG_LIM': 1, 'OUT_QG_LIM': 1, 'RETURN_RAW_DER': 0, 'PDIPM_FEASTOL': 0, 'PDIPM_GRADTOL': 1e-06, 'PDIPM_COMPTOL': 1e-06, 'PDIPM_COSTTOL': 1e-06, 'PDIPM_MAX_IT': 150, 'SCPDIPM_RED_IT': 20}\n",
      "{'case_name': 'case14', 'load_resolution': '5min', 'fpr': 0.02, 'noise_ratio_power': 0.02, 'noise_ratio_voltage': 0.001, 'pv_bus': array([ 3,  4, 10, 12]), 'measure_type': 'HALF_RTU'}\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from gen_data.gen_data import improve_resolution, add_cloud, gen_case, gen_load, gen_pv, gen_measure\n",
    "from configs.config import sys_config, opt\n",
    "from pypower.api import runopf\n",
    "from utils.class_se import SE\n",
    "from configs.config_mea_idx import define_mea_idx_noise\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(opt)\n",
    "print(sys_config)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load raw data\n",
    "raw_load_dir = 'gen_data/raw_data/load.csv'\n",
    "raw_pv_dir = 'gen_data/raw_data/pv.csv'\n",
    "raw_load = pd.read_csv(raw_load_dir)\n",
    "raw_pv = pd.read_csv(raw_pv_dir)\n",
    "raw_load = raw_load[(raw_load['DateTime'] <= '2012-09-30 23:45:00') & (raw_load['DateTime'] >= '2012-06-00 00:00:00')]\n",
    "raw_pv = raw_pv[(raw_pv['DateTime'] <= '2019-09-30 23:45:00') & (raw_pv['DateTime'] >= '2019-06-00 00:00:00')]\n",
    "raw_load.reset_index(inplace=True, drop=True)\n",
    "raw_pv.reset_index(inplace=True, drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Improve data resolution\n",
    "high_load_dir = 'gen_data/raw_data/load_high.csv'\n",
    "high_pv_dir = 'gen_data/raw_data/pv_high.csv'\n",
    "load_high, pv_high = improve_resolution(load_raw = raw_load, pv_raw = raw_pv, res = '5min')\n",
    "\n",
    "# Add cloud (large deviations) to the pv dat\n",
    "pv_new = add_cloud(pv_high, unchange_rate=0.5, max_reduce=0.8)\n",
    "\n",
    "# Save\n",
    "#load_high.to_csv(high_load_dir, index=False)\n",
    "#pv_high.to_csv(high_pv_dir, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Default active load: [  0.    54.25 235.5  119.5   19.    28.    75.    75.    73.75  22.5\n",
      "   8.75  15.25  33.75  37.25]\n",
      "Default reactive load: [ 0.  12.7 19.  -3.9  1.6  7.5  0.   0.  16.6  5.8  1.8  1.6  5.8  5. ]\n",
      "Default total active load: 797.5\n",
      "Default generation limit: 772.4\n",
      "Success: False\n",
      "Voltage phase angle: [  0.          -5.65621323 -15.70237604 -14.21631462 -11.56732536\n",
      " -15.78823337 -21.03030736 -18.45556237 -22.75527478 -22.22793426\n",
      " -19.33417167 -18.08111832 -18.71131473 -23.3454873 ]\n"
     ]
    }
   ],
   "source": [
    "# Modify case\n",
    "case = gen_case(sys_config['case_name'])\n",
    "print(f'Default active load: {case[\"bus\"][:,2]}')\n",
    "print(f'Default reactive load: {case[\"bus\"][:,3]}')\n",
    "print(f'Default total active load: {case[\"bus\"][:,2].sum()}')\n",
    "print(f'Default generation limit: {case[\"gen\"][:,8].sum()}')\n",
    "# Test OPF on it\n",
    "result = runopf(case, opt)\n",
    "print(f'Success: {result[\"success\"]}')\n",
    "print(f'Voltage phase angle: {result[\"bus\"][:,8]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate load for case\n",
    "high_load_dir = 'gen_data/raw_data/load_high.csv'\n",
    "high_pv_dir = 'gen_data/raw_data/pv_high.csv'\n",
    "\n",
    "load_high = pd.read_csv(high_load_dir)\n",
    "pv_high = pd.read_csv(high_pv_dir)\n",
    "\n",
    "load_active_dir = f'gen_data\\{sys_config[\"case_name\"]}\\load_active.npy'\n",
    "load_reactive_dir = f'gen_data\\{sys_config[\"case_name\"]}\\load_reactive.npy'\n",
    "load_active, load_reactive = gen_load(case = case, load_raw = load_high)\n",
    "#np.save(load_active_dir, load_active, allow_pickle=True)\n",
    "#np.save(load_reactive_dir, load_reactive, allow_pickle=True)\n",
    "\n",
    "pv_active_dir = f'gen_data\\{sys_config[\"case_name\"]}\\pv_active.npy'\n",
    "pv_reactive_dir = f'gen_data\\{sys_config[\"case_name\"]}\\pv_reactive.npy'\n",
    "pv_active, pv_reactive = gen_pv(sys_config['pv_bus'], pv_high, load_active, penetration_ratio = 0.3)\n",
    "#np.save(pv_active_dir, pv_active, allow_pickle=True)\n",
    "#np.save(pv_reactive_dir, pv_reactive, allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max active load: 582.6279883399718\n",
      "Max pv load: 174.7883965019915\n",
      "PV penetration: 0.29999999999999993\n",
      "z shape: (68, 1)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f4e04cce7f0>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD6CAYAAABJTke4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5WklEQVR4nO3deXyU5bnw8d81M5nsZCMhQICEfRXQuK+VasFTpVas2s2+ta96emw9x/ac4unbnh6PXTw9R7vZxVarpYu2aCtVhKpYNxQJCMgWCBAgQCAkJASyzGTmfv94ngmTyUwyEybLTK7v58PHmed5MrmD4Zp7rvu6r0eMMSillBpeHIM9AKWUUgNPg79SSg1DGvyVUmoY0uCvlFLDkAZ/pZQahjT4K6XUMBRV8BeRhSJSKSJVIrI0zPlUEXnGPr9ORErt458SkU1Bf/wiMi++P4JSSqlYSW91/iLiBHYB1wA1wHrgNmPM9qBrvgicY4y5W0RuBW40xtwS8jpzgL8YYyb19P1GjhxpSktL+/KzKKXUsLVhw4bjxpjCaK93RXHNBUCVMWYvgIg8DSwGtgddsxj4lv14OfATERHT9Z3lNuDp3r5ZaWkpFRUVUQxLKaVUgIjsj+X6aNI+Y4GDQc9r7GNhrzHGdABNQEHINbcAf4hlcEoppfrHgCz4isiFQIsxZmuE83eKSIWIVNTV1Q3EkJRSaliLJvgfAsYFPS+xj4W9RkRcQA5QH3T+VnqY9RtjHjPGlBtjygsLo05ZKaWU6qNogv96YIqIlImIGyuQrwi5ZgVwu/14CbAmkO8XEQfwCaLI9yullBoYvS74GmM6ROQeYDXgBJ4wxmwTkQeACmPMCuBxYJmIVAENWG8QAVcABwMLxkoppQZfr6WeA628vNxotY9SSsVGRDYYY8qjvV53+Cql1DCkwV/FrMXTwbMbahhqnxqVUtHT4K9i9vL2o3zlT5uprm8Z7KEopfpIg7+KWZvXB1ifAJRSiUmDv4qZx2ele9o7/IM8EqVUX2nwVzHz2EG/3avBX6lEpcFfxczrs4N/h2+QR6KU6isN/ipm3sDMX9M+SiUsDf4qZh6fBn+lEp0GfxWzzuDv1bSPUolKg7+KmUfTPkolPA3+KmZeTfsolfA0+KuYeTsCdf6a9lEqUWnwVzE7k/PXmb9SiUqDv4qZVvsolfg0+KuYnVnw1bSPUolKg7+KmS74KpX4NPirmHk1569UwtPgr2KmaR+lEp8GfxUzbemsVOLT4K9ipjt8lUp8UQV/EVkoIpUiUiUiS8OcTxWRZ+zz60SkNOjcOSLyjohsE5EPRCQtjuNXg8CrvX2USni9Bn8RcQKPAouAmcBtIjIz5LI7gBPGmMnAI8BD9te6gN8CdxtjZgFXAd64jV4NCq32USrxRTPzvwCoMsbsNcZ4gKeBxSHXLAaesh8vBxaIiADXAluMMZsBjDH1xhidLiY4TfsolfiiCf5jgYNBz2vsY2GvMcZ0AE1AATAVMCKyWkQ2isi/hfsGInKniFSISEVdXV2sP4MaYHonL6USX38v+LqAy4BP2f+9UUQWhF5kjHnMGFNujCkvLCzs5yGps9Wu9/BVKuFFE/wPAeOCnpfYx8JeY+f5c4B6rE8JbxhjjhtjWoCVwLlnO2g1uDTnr1Tiiyb4rwemiEiZiLiBW4EVIdesAG63Hy8B1hhjDLAamCMiGfabwpXA9vgMXQ0W3eSlVOJz9XaBMaZDRO7BCuRO4AljzDYReQCoMMasAB4HlolIFdCA9QaBMeaEiDyM9QZigJXGmBf76WdRA8DnN/itPV4681cqgfUa/AGMMSuxUjbBx74Z9LgNuDnC1/4Wq9xTJYFAysftcuDp8GOMwSrsUkolEt3hq2ISmO2PSHN1ea6USiwa/FVMAjP/rFQN/kolMg3+KiaBxd6szpm/LvoqlYg0+KuYBGb+2akpgNb6K5WoNPirmHSmfTTnr1RC0+CvYhII9tma9lEqoWnwVzHx2jdyydYFX6USmgZ/FZNuC76a81cqIWnwVzHpXPBNsxd8Ne2jVELS4K9i4vGF5vx15q9UItLgr2Li6Qid+WvwVyoRafBXMTlT5x/I+WvaR6lEpMFfxaT7Dl+d+SuViDT4q5h4NeevVFLQ4K9i4rHr/M80dtO0j1KJSIO/iknngq/29lEqoWnwVzEJvpmL2+XQtI9SCUqDv4pJYOaf4hRSXQ5N+yiVoDT4q5h4fX4cAi6ng1SXU2f+SiUoDf4qJh6fnxSn9WuT6nJozl+pBKXBX8XE0+HH7bKDf4qmfZRKVFEFfxFZKCKVIlIlIkvDnE8VkWfs8+tEpNQ+XioirSKyyf7z8ziPXw0wr8+Pu3Pmr2kfpRKVq7cLRMQJPApcA9QA60VkhTFme9BldwAnjDGTReRW4CHgFvvcHmPMvPgOWw0WT0dI2keDv1IJKZqZ/wVAlTFmrzHGAzwNLA65ZjHwlP14ObBARCR+w1RDhddnzqR9XA7t7aNUgoom+I8FDgY9r7GPhb3GGNMBNAEF9rkyEXlfRF4XkcvDfQMRuVNEKkSkoq6uLqYfQA0sa8HXel9PTdG0j1KJqr8XfI8A440x84H7gN+LyIjQi4wxjxljyo0x5YWFhf08JHU2rAVfJ6BpH6USWTTB/xAwLuh5iX0s7DUi4gJygHpjTLsxph7AGLMB2ANMPdtBq8FjLfjaM3/d5KVUwoom+K8HpohImYi4gVuBFSHXrAButx8vAdYYY4yIFNoLxojIRGAKsDc+Q1eDoeuCr1Pr/JVKUL1W+xhjOkTkHmA14ASeMMZsE5EHgApjzArgcWCZiFQBDVhvEABXAA+IiBfwA3cbYxr64wdRA8PrC63z1+CvVCLqNfgDGGNWAitDjn0z6HEbcHOYr3sWePYsx6iGEI/PkOEOLvXUtI9SiUh3+KqYdNnhq5u8lEpYGvxVTLru8HXg6fBjjBnkUSmlYqXBX8XEWvAN1Plbvz46+1cq8WjwVzHpsuBr1/tr8Fcq8WjwVzHxhrR0Br2Pr1KJSIO/ikl7lwVfO/hrrb9SCUeDv4pJlwXfFE37KJWoNPirmIS2dAZN+yiViDT4q6j5/Aa/oXvaR2f+SiUcDf4qal6fFeSDe/uA5vyVSkQa/FXUAjP84N4+1nFN+yiVaDT4q6gFZv7BLZ1B0z5KJSIN/ipqno4IaR8N/kolHA3+KmqdM/9udf6a9lEq0WjwV1HrtuCrvX2USlga/FXUui34atpHqYSVNMG/4bSH5zcd4tjJtsEeStLy+qzWzW7d5KVUwkua4F9zooV7n97ElpqmwR5K0uq+4Gv9t03r/JVKOEkT/PMz3YD1CUD1j9AFXxHBrbdyVCohJU3wL8hMBeD46fZBHkny8nQu+ErnsVSXQ3f4KpWAogr+IrJQRCpFpEpEloY5nyoiz9jn14lIacj58SJySkS+Gqdxd5PudpLhdtJwSmf+/cUTsuALeh9fpRJVr8FfRJzAo8AiYCZwm4jMDLnsDuCEMWYy8AjwUMj5h4GXzn64PcvPdGvapx+d2eEbHPw17aNUIopm5n8BUGWM2WuM8QBPA4tDrlkMPGU/Xg4sEBEBEJGPAfuAbXEZcQ8KMt0c1+Dfb0IXfMGq9deZv1KJJ5rgPxY4GPS8xj4W9hpjTAfQBBSISBbwNeA/z36ovSvISqVBc/79JnTBF+y0j+b8lUo4/b3g+y3gEWPMqZ4uEpE7RaRCRCrq6ur6/M3yM92a8+9HHrvOP0XTPkolPFcU1xwCxgU9L7GPhbumRkRcQA5QD1wILBGR/wZyAb+ItBljfhL8xcaYx4DHAMrLy00ffg7gTNrHGIOddVJxFH7BV9M+SiWiaIL/emCKiJRhBflbgU+GXLMCuB14B1gCrDHGGODywAUi8i3gVGjgj6f8TDeeDj+nPT6yUqP50VQswi74pjhpavUO1pCUUn3Ua4Q0xnSIyD3AasAJPGGM2SYiDwAVxpgVwOPAMhGpAhqw3iAGXEGWVevfcMqjwb8fnFnwDa3z17SPUokmqghpjFkJrAw59s2gx23Azb28xrf6ML6YFNi7fOtPtzO+IKO/v92w4/X5cQi4QnL+Hk37KJVwkmaHL5xp8VCvi779wuPzd1nsBd3kpVSiSqrgX5Cl/X36k6fD32WxFwJ1/pr2USrRJFfwt/v71Gvw7xden7/LYi9obx+lElVSBf90t5P0FCf1p3SjV3/wdGjaR6lkkVTBH6zUj6Z9+ofXZ7qnfVwOPD4/fn+ft2copQZB8gX/TLemffqJteDbdfNc4D6+gXbPSqnEkHTBPz/TTb329+kX1oKvs8uxzvv4at5fqYSShME/Vfv79BNrwTdk5q/38VUqISVd8B+ZZaV9rO4SKp7CL/gGgr/O/JVKJEkX/PMz3bR3+Gnx6Ew03ry+cHX+dtpHZ/5KJZSkDP6gu3z7g8dnIs782zTnr1RCSbrgPzIrsNFLF33jLewOX037KJWQki74B2b+Wusff+F3+GraR6lElLTBX9M+8Wct+Iav89eZv1KJJemCf6C5m270ir+wC76BtI/m/JVKKEkX/DPcLtJTnHoj937gjdDSGTTto1SiSbrgD/YuX037xF27LvgqlTSSMvgXZGl/n/4QbsE3rbPOX4O/UokkOYN/pnb27A9hd/gGFnz1Pr5KJZSkDP75mana0z/OfH6D36BpH6WSRFIG/wLt7xN3Xrtlc+jMP5AG0uCvVGKJKviLyEIRqRSRKhFZGuZ8qog8Y59fJyKl9vELRGST/WeziNwY5/GHVaD9feIuENxDZ/4iYt3KUat9lEoovQZ/EXECjwKLgJnAbSIyM+SyO4ATxpjJwCPAQ/bxrUC5MWYesBD4hYi44jT2iHSXb/wFZv6hLZ1B7+OrVCKKZuZ/AVBljNlrjPEATwOLQ65ZDDxlP14OLBARMca0GGM67ONpwIDkYQIbvY5r3j9uPB3h0z5gdfbUtI9SiSWa4D8WOBj0vMY+FvYaO9g3AQUAInKhiGwDPgDuDnoz6DcFmVZzN535x0/nzN8VJvhr2kephNPvC77GmHXGmFnA+cD9IpIWeo2I3CkiFSJSUVdXd9bfs7O/jwb/uIm04AuB4K8zf6USSTTB/xAwLuh5iX0s7DV2Tj8HqA++wBizAzgFzA79BsaYx4wx5caY8sLCwuhHH0Fnfx/d5Rs3kRZ8wWrxoDl/pRJLNMF/PTBFRMpExA3cCqwIuWYFcLv9eAmwxhhj7K9xAYjIBGA6UB2Xkfcgw+0iLcWh/X3iyOuzlmtCd/iCtdFL0z5KJZZeK2+MMR0icg+wGnACTxhjtonIA0CFMWYF8DiwTESqgAasNwiAy4ClIuIF/MAXjTHH++MHCVWQmappnzjqccFX0z5KJZyoyi6NMSuBlSHHvhn0uA24OczXLQOWneUY+6QgS5u7xVPPC75OGlu9Az0kpdRZSModvmAt+mq1T/ycmflHqvPXtI9SiSRpg39BZqoG/zjy9DTzT3F2vjkopRJD8gb/LDf1p9u1v0+cnNnhG13Ov+ZECx//6dtUHWsekPEppWKTtME/P9NNm1f7+8RL7wu+Xf+e39p9nI0HGvnX5Vvw+fUNWKmhJmmDf4H294mr3hZ8Q+v8d9ZaM/73DzTy5Nrqfh+fUio2yRv89UbucdVzb5/uaZ+dtSeZNy6XD00r5H9WV3KgvmVAxqmUik7SBv98u7+P3tQlPjyBTV4Revt4fH78dnrHGMPO2mZmjM7m2zfOwekQlj63RddflBpCkjb4F/Rzf58tNY3DKqXU84KvdR/fQEXQ0ZPtNLZ4mV48gjG56dx/3XTW7qnnmfUHu32tUmpwJG/wz+rfnP+nfrmOH6/Z3S+vPRT1VucPdOb9d9aeBGBacTYAt50/ngvL8vn2izuobWobiOEqpXqRtME/0N+nP9I+bV4fze0d7Kk7HffXHqq8Pj8OAVeEnD/QWfETWOydbgd/h0N46KZz8Pr9fO+lHQM0YqVUT5I2+EP/9fdpslsZVB8fPsHf0+EPu9gLZ9I+gUXfytpmRuekkZvh7rymdGQmV0wpZPuRk/0/WKVUr5I7+Gf1T4uHxhYr+NecaBk2O1s9Pn/YxV4ISvvYM/8dR052pnyC5We6OdGiPYCUGgqSOvjnZ7rZX99Ca5w3ejW2WG8ofmO9AQwHXp8/7GIvnAn+bV4/Xp+fPXWnmF48ott1eZluGls8WvWj1BCQ1MH/H+aMprr+NDf+9O24pmiagjpYVtcPj9RPj2mflDNpn711p/H6TGe+P1heRgpen+FUe7/fyVMp1YukDv43l4/j1587n9qTbVz/47dYva02Lq8b3L64+vhwmfmbqNI+gUqf6aPDBX9rDaBRUz9KDbqkDv4AV00r4oUvXUZZYSZ3LdvAd1/a0Vmz3ldNdvByOx3DbObfvcwTgoO/n521zbgcwsSRWd2uCwT/4bQ/QqmhKumDP0BJXgZ/uvtiPn3ReH7x+l6u//FbbDxwos+v19TqxekQpozKonqYtC2wFnydYc91Vvt4/ew8cpLJRVlhPyXk2RvvTrRo8FdqsA2L4A9WgHrwY3P4xWfOo6nVy00/W8vX//xB5yw+Fo2tHnLSUygbmTlsyj2tBd8IM/+gOv/K2uaw+X6wcv6gwV+poWDYBP+Aj8wq5uX7ruTzl5bxh/cOsODhv8e8FtDY4iXXDv7Dpdyz5zp/63hdczuHm9qYFqbSB6zqK4ATpzXnr9RgG3bBHyAr1cU3PjqTFfdcRmF2Gvc+/X63fvQ9aWr1MiI9hQkFmcOm3NPbY52/lfbZXNMEhF/sBRiRloJDdOav1FAwLIN/wOyxOdx3zVTavH7eP9AY9dc1tXrJzUihbGQGMDzKPXsu9bSOb6lpBIiY9nE4hNwMtwZ/pYaAqIK/iCwUkUoRqRKRpWHOp4rIM/b5dSJSah+/RkQ2iMgH9n+vjvP4z9oFZfk4BNbuqY/6awJpnwkFmcDwKPf0RFHqub++hZz0FIpHpEV8ndyMFE37KDUE9Br8RcQJPAosAmYCt4nIzJDL7gBOGGMmA48AD9nHjwPXG2PmALcDy+I18HjJSU9hztgc3tlzPOqvaWyxFnwLMt1kp7qGxcy/px2+wcenFWcjEn5hGKxyT535KzX4opn5XwBUGWP2GmM8wNPA4pBrFgNP2Y+XAwtERIwx7xtjDtvHtwHpIpIaj4HH08WTRvL+gUZaPL3vPPX5Dc3tHeRkuBERJozMGBblnj3V+YtI5+x/RoSUT0BeRv/0W1JKxSaa4D8WCL4LR419LOw1xpgOoAkoCLnmJmCjMaZbj2URuVNEKkSkoq6uLtqxx80lkwro8BvWV/de+9/c5sUYyE23yhZLC4ZHuWdPC75wJvUzfXT4Sp+AvIwU3eGr1BAwIAu+IjILKxV0V7jzxpjHjDHlxpjywsLCgRhSF+eX5pPiFNZGkfoJBK7cjDPBfziUe/a04Atn+vuE6+YZzOrsqc3dlBps0QT/Q8C4oOcl9rGw14iIC8gB6u3nJcCfgc8aY/ac7YD7Q7rbyfzxebwTxaJvoKlbTmDmP3J4lHv21NIZzsz8p43qOfjnZrhp7/DT6o1vp1WlVGyiCf7rgSkiUiYibuBWYEXINSuwFnQBlgBrjDFGRHKBF4Glxpi34zTmfnHJpAK2HmrqdcdvoKnbmZm/Ve65P8nz/j0t+IIV/MfnZ5CZ6urxdfIzA7t8NfWj1GDqNfjbOfx7gNXADuCPxphtIvKAiNxgX/Y4UCAiVcB9QKAc9B5gMvBNEdlk/ymK+08RB5dMGonfwLp9Pc/+A738c9Kt3aqlI61yz31JnvfvLe0zNi+DC8rye32dwN29Tuiir1KDqudpms0YsxJYGXLsm0GP24Cbw3zdg8CDZznGATFvXC5pKQ7W7qnn2lnFEa87GZL2GQ7lnj6/wW/oMe3zy8+ehxC5xDMgX5u7KTUkRBX8hwO3y8H5pfm95v0DC76B4D8cyj0Di9k9LvhG6PgZKtDcTcs9lRpcw7q9Q6hLJo2k8mgzdc3dqlE7NbZ6yXQ7u8yCk73c02Pf/6CnmX+09IYuSg0NGvyDXDLJ2prw7t7Is/+mVm/nrD8g2cs9Aze/idTSORaBvzud+Ss1uDT4B5k1ZgTZaa4e+/w0tnjJsWevAcle7hlN2idaLqeDEWmuzoVzpdTg0OAfxOV0cGFZQY99fppaPZ27ewOSvdzTG8e0D1iLvg0R0j5tXh8HG1rw+XUTmFL9SRd8Q1wyqYBXdhzlUGMrY3PTu51vavV2uz9tcLnnhwZklAMrnjN/sMo9I838H3hhO79fd4BUl4PJRVlMKcpiWvEIPnnh+G7pNqVU3+nMP8Qlk628f6Sqn8YWb+cGr4CCTDdZqS72J2m5ZzwXfMGe+UfI+e88cpKJhZl85qIJFGSl8t6+Bh5atZM/VRwMez2AMYaHX97F9sMn4zI+pYYDnfmHmFqUjdvpYPex5rDnG1u95IQEfxGhdGQG+5I27WOlYHra4RuL3IwUKmvD//0eaGhlwfQi/t9Hz3QNL3/wZXYdDX89wLHmdn706m6a27z8x5hZcRmjUslOZ/4hHA5hTG4ah060djvX5vXh6fCHTT9MSOJyz3inffIj9PRv8XRw/FQ74+01lIApRdnsOnoq4usF3kiSfZe1UvGkwT+MkrwMasIE/86OnunubufKkrjcM94LvnmZblo8PtpCmrsdbLD+zsflhwT/UVlUHTsVsRNo4FOBBn+loqfBP4yxuekcagwT/Fut2Wpozh9gclEWfgN76iLPUBPVmZn/2df5Q+SNXoE1kwndgn82p9o7ONLUFvb1AsH/YENyvvkq1R80+IcxNi+duub2bjPTppDWDsHOKckBztzEPJnEe8E3UouHAw3Wmsn4kOA/tciqroqU96+0U0J+c+Y1lFI90+AfRqDE83DI7L+xNXLwLy3IZESai00Hm/p/gAPszA7f+KV9gG7lngcbWshOdXX7ZDXVvkfA7jB5f7/fsPtoM+eOzwU09aNUtDT4h1GSZwX/0NRPU8hdvII5HMLccblsPtjY7+MbaPFe8A2kfRpaus/8x+VndLsBfF6mm5FZ7rAVWIcaW2nx+Fg42+rEujcJ025K9QcN/mGMDQT/kEXfps4buXRf8AWYW5JL5dFmWj3JdZequC/4ZoS/ocuBhpZuKZ+ASBU/gUqf8ybkMzLLrTN/paKkwT+M4hFpOB3Sbebf2OrB6RAy3eHbF88dl4vPb9h2OLlSP/2xwxe63tDF7zccPNHarcwzYGqEip9Kex1g6qgsykZmsleDv1JR0eAfhsvpoHhEWrdyz8YWL7npKd3SEgFz7UXfTUmW+vEENnnFaebvdjnISnV1qfU/1tyOp8MfeeZvV/wcDqn42XW0mbG56WSnpVA2MlNn/kpFSYN/BGPz0sOmfUJ39wYrGpHG6Jw0ttQk18w/3gu+AHmZKV1m/oEyz8hpH6viZ3dIxU9lbTNTR1nnykZmUdfcTnOb3itAqd5o8I+gJEytf7he/qHmluSyOc7lnvvrT0dshzAQ4l3nD9aib3DOP1KZZ0C4ih+vz8/eutNMLbbOldkN9qqPa7mnUr3R4B/B2Lx0jjS1ds564Uzapydzx+Wyv74lrjcoX/rsB3zm8XWD1ubY6/PjECsdFi95IS0eDja04BAYE6aTKgQqflK71Prvrz+Nx+dnmv3GMLHQCv57j2vFj1K90eAfwdjcdPwGaoNyzE2t3oiVPgFzx1l5/3jN/ts7fGw8cIJjze28XRX5PgP9ydPhj9tib0BeRkqX4H+goYXROek9ritMHZXFrmNnAntl7Sn7uBX8x+dnIAJ76zTvr1RvovoXLSILRaRSRKpEZGmY86ki8ox9fp2IlNrHC0TkNRE5JSI/ifPY+1VJnpV+CE79NLZ4ek37zBmbgwhsjtNmr62Hmmi30y5/fv9QXF4zVh6fP26LvQF5mW4aT3dN+0RK+QRMKcqi6mhzZ8VP5dFmHGK11gBIS3FSkpeui75KRaHXf9Ei4gQeBRYBM4HbRGRmyGV3ACeMMZOBR4CH7ONtwDeAr8ZtxAMkUOsfqPjx+Q0n2zp6Df7ZaSlMLsyK28x/ffUJAD4yaxSrttZyur0jLq8bC6/PH9fFXrDSPs3tHZ3rCQcaWnsP/qOyOe3xdVb87KptprQgk7SUM6W3ZSOzNPgrFYVo/kVfAFQZY/YaYzzA08DikGsWA0/Zj5cDC0REjDGnjTFvYb0JJJTROWnAmY1egQqScLt7Q80dl8uWmsaIXShjsX5fAxMLM/nC5RNp9fpYtbX2rF8zVv2S9gm0eGj1RGzlHCqQ3gnk/Xcdbe48FjDRLveMx9+9Usksmn/RY4Hg2yjV2MfCXmOM6QCagIJoByEid4pIhYhU1NXVRftl/SotxUlRdiqHGq3KkcYemrqFmjsul+OnPGE7g8bC7zdU7D/BBaX5lE/IY1x++qCkfrw+E/+0T2CX72lvr5U+AcHlnm1eH9X1Zyp9AspGZnKqvYO6U+1xHa9SyWZILPgaYx4zxpQbY8oLCwsHezidxuadKfdsbI1h5m9v9jrbvP/uY6doavVyfmk+IsKN80t4e89xjjSd3ZtKrKyZf/zKPOFMf58TLR4O1EcX/AMVP7uPnqLq2Cn8hs5Kn4BAuacu+ianxhYPv3pzLx0+bd19tqIJ/oeAcUHPS+xjYa8REReQA4S/CW4CGZub3pnzb+rs6NlztQ/A9OIRuJ2OiHl/r89PZW0zf3n/EN99aQev7jga9rr3qhsAOL80H4Ab54/FGHh+0+FYf5SzYi34hm9p0Vd5QS0eop35w5mKn0DqZ1pxVpfzgeCvef/k9NTa/Tz44g5e/ODIYA8l4UVzD9/1wBQRKcMK8rcCnwy5ZgVwO/AOsARYY5Ig6To2L52/bTuK32862w9Hk/ZxuxzMHDOiW5uHtXuO872XdrLzSHNnj3yA5zYeYu3Swm559YrqBkaNSGVcvrX4XDYyk/njc3luYw13XTExYpuJeLMWfOM8888809wtUivncKaOyuZPFQeprG3G7XQwoSCzy/kxuVa5qAb/5LRqm7Xm9cRb+7hh7pgB+zeQjHqd+ds5/HuA1cAO4I/GmG0i8oCI3GBf9jhQICJVwH1AZzmoiFQDDwOfE5GaMJVCQ1ZJXgYen5+6U+1BHT17D1AA88bl8kFNU+fH09+vO8BnH3+Pk61e/s+lpfzglnms+ufL+fmnz6WuuT3s7H/9vgbK7ZRPwMfPLWHX0VNsO3wyDj9hdPqnzj8o7ROhlXM4U0Zlcdrj47XKY0wszOw2LqdDKCvIjEvaZ2ftSe774yYu/M4rSdesLxHtrz/NjiMnmV6czeaaJjbsPzHYQ0po0cz8McasBFaGHPtm0OM24OYIX1t6FuMbVCW5gXLPlh7v4hXO3HE5PLm2msqjzSzfUMOv367myqmF/PiT8xmRduY1JhdmUTwijd+/d5CFs0d3Hq850cLhpjbuslM+AR+dM5oH/rqNP79/iNljc872R4yKtx/q/NNSnKSnODvTPqFVO5FMKQpU/Jxi8bwxYa8pG5kZtvd/NIwxvLO3nsfe2MvfK+tIT3Hicgr//twHPPfFS3E6dKY5WFbbs/4f3TafJT9byxNv76M85N+Hit6QWPAdqoJr/RtbvWS6nVHPgOeW5ALw+SfX8+u3q/n8pWU8fnt5l8APVsuET5w/jjd313Ew6BaEFXZ9//khv9x5mW6unl7E85sO09TipbK2mTU7j7LsnWpeqzzW1x+1R/0x8wfIz3TTcNpjtXKOIt8PdDZxsx6Hf8MoK8zkQENLnxYFv/SH9/nkL9ex9VATX7lmKmuXXs23b5zD5pomfvvu/phfT8XPqq21zB47gqmjsrntwvGs2lrb5d+Mio0G/x4Ebud4qLHV6uvTS2uHYIHbOtaf8vCdG+fwzetnRuyNc8v54xDgmfVnKmrfq24gO9XFtOLuAe7G+SUcP9XO3Af+xkd+8Aaff7KCbzy/jbuXbaDFE/9NYJ5+KPUEK4W2s7YZT4efcVEG/9wMN4XZqUD3Sp+AspGZeH2mW0vu3uysPckLW47wuUtKeetrV/OlBVPIy3Rz/TmjuWJqId9fXdml3YcaOEdPtrHxQCMLZ1l3bLv94lJEhKfWVg/uwBKYBv8eZKa6yMtIoeZEa1QdPYM5HMLPPn0ey//xEj554fgerx2bm85V04r4Y8XBzkZyFdUNnFeaFzbN8OEZRXz12qksXTSdH982n+e+eAk/+9S5tHf4+Xtl9PskjDE8tGonG/Y39Hhdf+zwBWvmH7gZS7Qzfzgz+w/3xgjWRi+IveLnt+/ux+1ycO+CKV12DYsIDy6ejdfn51srtsX0mio+/rbdWhP7iB38x+Smc92c0Tyz/iCnBmHXezLQ4N+LQF//ptbe+/qEunTySOaNy43q2tsuGM+x5nbW7DzGidMedh091S3lE+ByOrjn6incfeUkrp87hnPH53HtrGIKMt2sjKEEbseRZn729z38ZE1Vj9f1R50/WLP4QKfSWIL/nLG55Ge6Oz+Zheqs9Y8h+J9q7+DPGw9x/TljOncfBxtfkMG9H57Cqm21vLw9fGmu6j+rt9YysTCzs48TwB2XldHc3sGfKg728JUqEg3+vRhr9/W30j6xBf9YfGhaIaNGpPKH9w5QsT98vr8nTodw7axRvLbzGG3e6O4h/Nct1n6Bt6qOc7KHG6C0en39kvbJt/8+e2rlHM69C6bw4pcvwxFh8TU/082INBf7Ymjt/Of3D3Ha4+PTF0X+lPZ/L5/ItFHZ/MfzWwelx9Jw1dji4Z299SycVdylImzeuFzOHZ/Lr9+uHrR254lMg38vxuZmcMhe8O3P4O9yOrilfByv76rjL5sO4XY6OKcktmqeRbNHc9rj483dvbd+Nsbw182HGZOThtdnWLMj/GLx9sMnqWtuZ9aY+FcWBdZQemvlHCrd7WR0TuQ3CxFhYmH0Dd6MMfz2nf3MHjuix09qKU4H3/n4bA43tfHIy7uiHq86O6/sOIbPb1g4u7jbuTsum8iBhhZeibBRUkWmwb8XJXnptHp91DW3R7W792x84nxrI/WLW45wTklOl7xzNC6eVEBOegovbe099fP+wUZqTrTyL9dMpSg7NWLDuGc31pDiFG6YG76s8mwE+vtM6KWhW19MHJnJvihr/ddXn6DyaDOfuWhCr3sNzpuQz63nj+PJtdW6kWyArNpay5icNOaEKW3+yKxRjM1N5zfvVMf9+7Z3+Fj67Bbe3D00+o3Fmwb/XgTKPSH6Gv++KsnL4MqpVm+j88tir19OcTr48IxRvLL9aGer5Ej+uvkwbpeDhbOL+cisYv6+61i3SiGvz8/zmw5x9fSisHnwsxV4zVjy/dGaVJTF4aa2qAL0snf3k53m4oa5of0Kw7vv2qm4XQ7+e9XOsx3mkLdi82GeeGvfoH3/0+0dvLG7jmtDUj4BLqeDm8tLWLunnpoT8S37fL2yjqfXH+T//Ho9z26oietrDwUa/HsRvKjYn2mfgE9fOAGAiydG3RS1i0WziznZ1sHaPZFTPz6/4YUtR7h6WhHZaSksml1Mm9fP6yGVQm/squP4KQ83nVvSp7H0JrDLN9oyz1jcXF5CptvJQy/1HKDrmttZtfUIS84rId0d3Setouw07r5yEi9traWiuudKqUR27GQbX1u+hQde2D5oi6p/r6zD0+EPm/IJuOncEoyBP2+Mb8fbVVtryUlP4aKJBXzlT5t59LWqpGoVrsG/FyVBM//e7t8bDx+eOYqX7r2cy6eM7NPXXzZlJFmprh77/q/bV09dczvX26mcC8ryyctI4aWQr3l2Yw0FmW4+NL2oT2PpTbF9z4RJhZm9XBm7QIBeta2W9T0EaKu81vDpiybE9PpfuLyMUSNSefDFHUkVEII98spuOvx+5o3L5et/2cqWON2gKBarttVSkOnusfhhXH4GF08sYPnGmrj9v/B0+Hl5x1GumTmKJz53PovnjeH7qyv55vPbkmZxWYN/L3LSU8hKdXU+HggzRo/oc8OqtBQnV08v4m/bj0bc4frXzUfIcFvXgfXR+dqZxazZeYz2DqtSqLHFwyvbj3HDvDH9srsXrB26z9x5EdfMjDyrOxtfuHxijwHa5zf87t39XDq5gEmFWWFeIbIMt4uvXDuNTQcbk7LDZNWxZp5Zf4BPXzSBJz53PoVZqdy9bAPHB/A+CU0tXl7ZfpRrZ43qta3GkvNK2F/f0nnnu7P19p7jNLd1cN2cYtwuB498Yh53XTGRZe/u58tPv48/Cd4ANPj3QkQ6Uz85A5D2iYdFs4tpOO3pbAkdzNPh56WtR7hm5qguaY6Fs4s51d7ReZP4v245gsfn77eUT8CFEwv6rV9OutvJV6+dxuaDjbywpXuA/vXb+zjc1NaZaovVTeeWML04m4dW7ex80zwbp9s7+OLvNvCxR9/m1sfe4fYn3uOuZRX88JXdA/7p4nsvVZLpdvGlq6eQn+nmF585j/rTHu75/cYB66X/+/cO0Or18ZmLSnu9dtGcYjLdTpZviE966qUPjpCd6uLSydYncIdDuP+6Gfzbwmm8uOUITybBzmIN/lEILPrG0t5hMF01rYj0FGfY1M/bVcdpbPF2q965ZHIB2akuXvrA+ppnN9QwvTibWWNGDMiY+8vHzy1hxugRPLRqZ+f+B5/f8F8vbOfBF3dw9fQirpk5qk+v7XQI/37dDA42tLLsnbPv+/MfK7bx0tZastNc+I316Wv30VM88squAb2Hw3v7Gnhlx1HuvmoS+fai/OyxOXznxjm8u7eB7/WyjhIPng4/T67dx6WTC5gZxe9ghtvFP5wzmhe3HDnrFiden5+/bT/KghlFpIbcx+Ifr5zEh2cU8b1VO9lZO3CddfuDBv8oBPL+A5X2OVvpbidXTStk1dbabh9P/7r5MCPSXFw+pesd01JdThbMKOLlHUeprG1m08FGbjq3JOH7pTsdwtevm0HNiVZ+8041p9s7uGtZBY+/tY/PXVLKY585L2LPpWhcMbWQK6YW8qNXd3PitKfPr7Ni82GWb6jhng9NZtkdF/LHuy7m+Xsu45X7rmRuSQ7fWbljQNoYGGP4zsodFI9I4/OXlnU5d9N5Jdx+8QR+9dY+7nhyPau31Xa2I4m3Fz84zNGT7Xzh8olRf82S88Zx2nP297let7eBxhZvly67ASLC9246hxFpLv756U1Rb6gcijT4R2HBjFEsnGV9rEwUC2cXc6y5nSfe3ke1fUPzNq+P1dtqWTR7dNhNVQtnF9PY4mXpc1twOoTF8+Nf2z8YLpsykqumFfLjNVXc/PN3WLPzGP95wyy+dcOsswr8Af9+3XRaPD6WPrelT+mZgw0tfP25D5g/PpcvL5jS5ZzDIXzrhlkca27vtQ1HQKvHx/+sruTFMKmu3ry0tZZNBxu575qpYauf/t9HZ/LlqyfzwaEm7lq2gYu/+yrfWbmDrYea4rYQaozhV2/uY3JRFleGTFJ6cn5pHuPzM1h+lmWZL2211sSumhb+e4/MSuX7S+ays7aZ76+uPKvvNZii6uc/3F05tbCz/j5RLJgxivH5GTz44g4efHEHuRkpjM/P4LTHxw0R+uBfOdVKF71/oJGrpxdRlJ02wKPuP/cvmsGiH77B/vrTPH77+XGtYJpePIKvLZzOt1fuYNm7+/nsxaVRf22Hz8+9T78PwI9unR92cX3++DyWnFfC42/t5RPlJUzsYXH6g5om7n3mffbWncYh4HSc12OZZLD2Dh//vWonU0dlcdN54dd6UpwO7rt2Gl9eMIXXd9XxzPqDPPHWPh57Yy+ZbifnlORy7oRc5o/L47IpI2PeqAjw7t4Gth0+yXc/PidiC49wRIQl55Xw8Mu7OGjfIChWPr9h9bZaPjS9qMexf2h6EZ+9eAKPv7WPq6YVdvsknQh05p+kslJdrPnKlaz88uV89+NzWDirGE+Hn7njcrkwwgay9KDZTn8v9A60acXZPPX5C1jxpcv6pXT1jsvKuGpaIQ++sCOmu379aE0VGw808uCNs3sMVl9bOJ00l5P/emF72PM+v+HR16q48adv09Lu4/Hby5k3Lpcv/+F91lb13u5jf/1plvzsHarrW7j/uhm9LsK7nA4WzBjFY58t5537F/CDW+Zx03klNLd7+fnre/nCbyq44r9f45dv7I25D9Kv3txLQaabG+dHt+ku2MfPHYuIdWvUvqiobuD4KQ+LonjD/PfrZjC5KIuv/HEz9TFWQbV5fby7t57K2maaWryDUi4sQ61Guby83FRUVAz2MIatdXvr+cUbe/npp87t06xtOKs/1c6iH75JVqqLv37pMjJTz3yw3rC/gafW7sfT4Scj1UmG24nL4eA371TzsfljefgT83p9/V+9uZcHX9zB47eXs2CGtUjt9xvWVzfwP3+rZH31Cf7hnNF8+2Ozyc1w09ji4ZZfvEvNiRb+cOdFnGPfYCjUC1sOs/TZD3AIfP/muZ1tk/uq1ePj3X31/PKNvazdU09uRgqfu6SUz11S2mvRxJ66Uyz439e5d8EU/uWaqX36/p/61bscaGjh9a9+KOwnh4MNLfz071W8tLWWJeeW8K8Lp3Uu7H5rxTb+8N4BNn7jmi7//yLZdriJGx9dS0aqk3+8chKfvbi0x82Cfr/hufcP8b9/q+RI0L0h0lIcjBqRxqLZo1m6aHoffmoQkQ3GmPKor9fgr1T8rN1znE/9ah032gH9g5om/vflSv5eWUdeRgojs1Jp8fho9fpo8XQwrXgEv/vChZ17SXri6fCz6Idv0OE3PPyJuaz8oJYXtxyh9mQb2aku/nPxLG6cP7bLIv3Rk20s+flaTrf7+ONdF3dpidzm9fFfL2znd+sOMH98Lj++bT4lefHdbb3xwAl++loVr+w4RobbyS3nj+Pzl5ZF/JTz9T9/wJ821LB26dWMzErt0/d8ftMh7n16E2Ny0lgwYxQLZhRx8aQCjja185PXdvPcxkM4RDi/LI+3q+qZXpzND2+dz5SiLC753hrOKcnhsc9GHUPZeqiJ76+u5PVddRRlp/KlBVO4pXxct3W1N3fX8Z2VO9lx5CRzxubwxasm0eE3HD3ZxtGTbdSebGfO2BHcecWkPv3cGvyVGmQPv7yLH726m/Mm5LFh/wlyM1K4+8pJfPbiCWS4z26Z7c3ddXzm8fcASHEKV04t4vq5o/nwjFERZ6rVx0+z5Ofv0OH3U5Sdyul2H81tXk57fPj8hruunMhXr53Wb5v5AHYcOckv39jLis2H8RvDdXNG838vn8jUUdm0en20eX00nPaw5Odr+di8sXzvpnP6/L2MMTy38RCrttXy5u462rx+MtxO2jv8OB3CJy8Yz91XTqI4J401O4/yb8u3cLKtg0+Ul/Dbdw/wg1vm8bE+pJzW7a3n+6srqdh/ghFpLrLTUkh1OUhxOvAZQ9WxU5TkpfOvH5nG9eeMiWk9Ixr9EvxFZCHwQ8AJ/MoY872Q86nAb4DzgHrgFmNMtX3ufuAOwAd82RizuqfvpcFfJboOn59PP76ObYdO8oXLJ/L5y0rJTotfmfBv392P2+ngI7OKo954uLP2JI+8vAtByEx1kZXqJDPVxWWTR3LJ5L61EumLI02tPPl2Nb9fd4DmCGsBf/uXKyLenzlWbV4fa/ccZ83OY2S6XXz+sjJGjehayHD8VDtfW76FV3ceI8UpbPjGNd3utR0tYwx/r6zj5R1Haff68fj8eDp8eH2GSyYV8JmLJ3TbOxAvcQ/+IuIEdgHXADXAeuA2Y8z2oGu+CJxjjLlbRG4FbjTG3CIiM4E/ABcAY4BXgKnGmIjFsRr8VTLwdPjxG6PrJhE0t3lZsfkwzW0dpLkcpKU4SUtxUpKXTnkMNzGKF2MMyzfU4PWZXm+7OlTFGvyj+Qx6AVBljNlrf4OngcVAcNnBYuBb9uPlwE/ESjwuBp42xrQD+0Skyn69d6IdoFKJqD/ufJZMstNS+FQf22r0BxHh5vJxgz2MARXNb+hYILhhRo19LOw1xpgOoAkoiPJrlVJKDbAhMT0RkTtFpEJEKurqkvOuOUopNZREE/wPAcGfh0rsY2GvEREXkIO18BvN12KMecwYU26MKS8sTLydckoplWiiCf7rgSkiUiYibuBWYEXINSuA2+3HS4A1xlpJXgHcKiKpIlIGTAHei8/QlVJK9VWvC77GmA4RuQdYjVXq+YQxZpuIPABUGGNWAI8Dy+wF3QasNwjs6/6ItTjcAfxTT5U+SimlBoZu8lJKqSQQa6nnkFjwVUopNbA0+Cul1DA05NI+IlIHnM098UYCvfewHVp0zAMnEcediGOGxBx3Io4ZrHFnGmOiLpcccsH/bIlIRSx5r6FAxzxwEnHciThmSMxxJ+KYoW/j1rSPUkoNQxr8lVJqGErG4P/YYA+gD3TMAycRx52IY4bEHHcijhn6MO6ky/krpZTqXTLO/JVSSvUiaYK/iCwUkUoRqRKRpYM9nkhE5AkROSYiW4OO5YvIyyKy2/5v3mCOMZSIjBOR10Rku4hsE5F77eNDdtwikiYi74nIZnvM/2kfLxORdfbvyTN2v6ohRUScIvK+iLxgP0+EMVeLyAcisklEKuxjQ/b3I0BEckVkuYjsFJEdInLxUB63iEyz/44Df06KyD/3ZcxJEfztu409CiwCZgK32XcRG4qeBBaGHFsKvGqMmQK8aj8fSjqArxhjZgIXAf9k//0O5XG3A1cbY+YC84CFInIR8BDwiDFmMnAC6xajQ829wI6g54kwZoAPGWPmBZUcDuXfj4AfAquMMdOBuVh/70N23MaYSvvveB7WbXNbgD/TlzEbYxL+D3AxsDro+f3A/YM9rh7GWwpsDXpeCYy2H48GKgd7jL2M/3ms23omxLiBDGAjcCHWBh5XuN+bofAHq+35q8DVwAuADPUx2+OqBkaGHBvSvx9Yref3Ya99Jsq4g8Z5LfB2X8ecFDN/Ev+OYaOMMUfsx7XAqMEcTE9EpBSYD6xjiI/bTp9sAo4BLwN7gEZj3W0OhubvyQ+AfwP89vMChv6YAQzwNxHZICJ32seG9O8HUAbUAb+202y/EpFMhv64A27Fukc69GHMyRL8k4ax3rqHZAmWiGQBzwL/bIw5GXxuKI7bGOMz1sfjEqx7R08f3BH1TEQ+ChwzxmwY7LH0wWXGmHOxUq//JCJXBJ8cir8fWC3tzwV+ZoyZD5wmJF0yRMeNve5zA/Cn0HPRjjlZgn9Udwwbwo6KyGgA+7/HBnk83YhIClbg/50x5jn78JAfN4AxphF4DStlkmvfbQ6G3u/JpcANIlINPI2V+vkhQ3vMABhjDtn/PYaVg76Aof/7UQPUGGPW2c+XY70ZDPVxg/Umu9EYc9R+HvOYkyX4R3O3saEs+E5ot2Pl1IcMERGsG/bsMMY8HHRqyI5bRApFJNd+nI61RrED601giX3ZkBqzMeZ+Y0yJMaYU63d4jTHmUwzhMQOISKaIZAceY+WitzKEfz8AjDG1wEERmWYfWoB146khPW7bbZxJ+UBfxjzYixZxXPy4DtiFldf9+mCPp4dx/gE4AnixZh53YOV1XwV2A68A+YM9zpAxX4b1MXILsMn+c91QHjdwDvC+PeatwDft4xOxbiVahfWROXWwxxph/FcBLyTCmO3xbbb/bAv8+xvKvx9BY58HVNi/J38B8ob6uIFMrHuk5wQdi3nMusNXKaWGoWRJ+yillIqBBn+llBqGNPgrpdQwpMFfKaWGIQ3+Sik1DGnwV0qpYUiDv1JKDUMa/JVSahj6/9gIl5SKAgpfAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Test if the maximum loads can converge in OPF\n",
    "case = gen_case(sys_config['case_name'])\n",
    "load_active_dir = f'gen_data/{sys_config[\"case_name\"]}/load_active.npy'\n",
    "load_reactive_dir = f'gen_data/{sys_config[\"case_name\"]}/load_reactive.npy'\n",
    "\n",
    "pv_active_dir = f'gen_data/{sys_config[\"case_name\"]}/pv_active.npy'\n",
    "pv_reactive_dir = f'gen_data/{sys_config[\"case_name\"]}/pv_reactive.npy'\n",
    "\n",
    "noise_sigma_dir = f'gen_data/{sys_config[\"case_name\"]}/noise_sigma.npy'\n",
    "\n",
    "load_active = np.load(load_active_dir)\n",
    "load_reactive = np.load(load_reactive_dir)\n",
    "\n",
    "pv_active = np.load(pv_active_dir)\n",
    "pv_reactive = np.load(pv_reactive_dir)\n",
    "\n",
    "# Analysis\n",
    "load_max = np.max(np.sum(load_active, axis = 1))\n",
    "print(f'Max active load: {load_max}')\n",
    "pv_max = np.max(np.sum(pv_active, axis = 1))\n",
    "print(f'Max pv load: {pv_max}')\n",
    "print(f'PV penetration: {pv_max/load_max}')\n",
    "\n",
    "# Generate noise sigma\n",
    "idx, no_mea, noise_sigma = define_mea_idx_noise(case, choice = sys_config['measure_type'])\n",
    "case_class = SE(case, noise_sigma, idx, fpr = sys_config['fpr'])\n",
    "result = case_class.run_opf(verbose=False)\n",
    "z, z_noise, vang_ref, vmag_ref = case_class.construct_mea(result)\n",
    "print(f'z shape: {z_noise.shape}')\n",
    "noise_sigma = np.abs(z * sys_config['noise_ratio_power']).flatten() + 1e-3\n",
    "#np.save(noise_sigma_dir, noise_sigma, allow_pickle=True)\n",
    "plt.plot(noise_sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "Voltage angle: [ 0.         -0.0745777  -0.18425626 -0.15516804 -0.12961976 -0.20633625\n",
      " -0.19371142 -0.13625724 -0.22392069 -0.22510798 -0.21872675 -0.22216505\n",
      " -0.22471583 -0.24458601]\n",
      "total load: 575.4147047709449\n",
      "True\n",
      "total load: 570.2323340055104\n",
      "True\n",
      "total load: 568.9799803256682\n",
      "True\n",
      "total load: 568.5910048448958\n",
      "True\n",
      "total load: 567.6804255549789\n",
      "True\n",
      "total load: 566.0975120449448\n",
      "True\n",
      "total load: 565.4247274294111\n",
      "True\n",
      "total load: 565.2204259385719\n",
      "True\n",
      "total load: 565.0539753302744\n",
      "True\n",
      "total load: 564.4193616259566\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "noise_sigma_dir = f'gen_data/{sys_config[\"case_name\"]}/noise_sigma.npy'\n",
    "case = gen_case(sys_config['case_name'])\n",
    "idx, no_mea, noise_sigma = define_mea_idx_noise(case, choice = sys_config['measure_type'])\n",
    "noise_sigma = np.load(noise_sigma_dir)\n",
    "case_class = SE(case, noise_sigma, idx, fpr = sys_config['fpr'])\n",
    "\n",
    "load_active_dir = f'gen_data/{sys_config[\"case_name\"]}/load_active.npy'\n",
    "load_reactive_dir = f'gen_data/{sys_config[\"case_name\"]}/load_reactive.npy'\n",
    "pv_active_dir = f'gen_data/{sys_config[\"case_name\"]}/pv_active.npy'\n",
    "pv_reactive_dir = f'gen_data/{sys_config[\"case_name\"]}/pv_reactive.npy'\n",
    "\n",
    "load_active = np.load(load_active_dir)\n",
    "load_reactive = np.load(load_reactive_dir)\n",
    "pv_active = np.load(pv_active_dir)\n",
    "pv_reactive = np.load(pv_reactive_dir)\n",
    "\n",
    "pv_active_ = np.zeros((load_active.shape[0], load_reactive.shape[1]))\n",
    "pv_reactive_ = np.zeros((load_reactive.shape[0], load_reactive.shape[1]))\n",
    "pv_active_[:,sys_config['pv_bus']] = pv_active\n",
    "pv_reactive_[:,sys_config['pv_bus']] = pv_reactive\n",
    "\n",
    "largest_load_idx = np.flip(np.argsort(np.sum(load_active - pv_active_, axis = 1)))[:10]\n",
    "\n",
    "for i in range(10):\n",
    "    result = case_class.run_opf(verbose=False, load_active = load_active[i] - pv_active_[i], load_reactive = load_reactive[i])\n",
    "    print(result['success'])\n",
    "    print(f'Voltage angle: {result[\"bus\"][:,8]*np.pi/180}')\n",
    "    break\n",
    "\n",
    "\"\"\"\n",
    "To make sure the largest load case can converge\n",
    "\"\"\"\n",
    "\n",
    "for i in largest_load_idx:\n",
    "    result = case_class.run_opf(verbose=False, load_active = load_active[i] - pv_active_[i], load_reactive = load_reactive[i])\n",
    "    print(f'total load: {np.sum(load_active[i] - pv_active_[i])}')\n",
    "    print(result['success'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('event_trigger_mtd')",
   "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.9.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "7dac1cb1ce86535b23df7051fbd88cb254e53683203cb38e019b53c60e510804"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}