01_Data_Labeling_and_Preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data Labeling\n",
    "===============================\n",
    "## Description\n",
    "In this notebook, the labels for the data are created so that supervised learning algorithms can be trained and tested. This is done by applying expert-domain rules based on the damping factor $\\zeta$. The first part is to import the required libraries and the raw eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run -i import_libraries.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing Raw Eigenvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We have 19815 number of scenarios, each with 49 eigenvalues\n",
      "Dataframe loaded!\n"
     ]
    }
   ],
   "source": [
    "# Loading simulation results\n",
    "eigs19branches = np.load(\"99_Data/19Branches_all_eigs.npy\")\n",
    "print(f\"We have {eigs19branches.shape[0]} number of scenarios, each with {eigs19branches.shape[1]} eigenvalues\")\n",
    "\n",
    "current_wd = os.getcwd()\n",
    "\n",
    "if not os.path.exists(os.path.join(current_wd, \"99_Data/df_eigenvalues.pkl\")):\n",
    "    \n",
    "    print(\"Creating dataframe...\")\n",
    "    # Creating column labels\n",
    "    col_labels = []\n",
    "    for i in range(1, eigs19branches.shape[0] + 1):\n",
    "        col_labels.append(\"Scenario {}\".format(i))\n",
    "\n",
    "    # Organizing the data in Pandas\n",
    "    df_eigenvalues = pd.DataFrame([], columns = col_labels)\n",
    "\n",
    "    for sc in range(1, eigs19branches.shape[0] + 1):\n",
    "        df_eigenvalues[\"Scenario {}\".format(sc)] = eigs19branches[(sc - 1)][:]\n",
    "\n",
    "    # The DataFrame has not been created in memory - Saving it as a .pkl file\n",
    "    with open(\"99_Data/df_eigenvalues.pkl\", 'wb') as f:\n",
    "        pickle.dump(df_eigenvalues, f, pickle.HIGHEST_PROTOCOL)\n",
    "    print(\"Dataframe created!\")\n",
    "    \n",
    "if os.path.exists(os.path.join(current_wd, \"99_Data/df_eigenvalues.pkl\")):\n",
    "    # DataFrame exists as a .pkl file - Loading it\n",
    "    with open(\"99_Data/df_eigenvalues.pkl\", 'rb') as f:\n",
    "        df_eigenvalues = pickle.load(f)\n",
    "    print(\"Dataframe loaded!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalizing Eigenvalues (Data Cleaning)\n",
    "\n",
    "Given the significative irregularity of the eigenvalues, it is necessary to perform some data cleaning to bring them to a form where the most important information is preserved. The relevant information about the power system operating condition specified by a given eigenvalue is contained in the _angle_ it forms with the positive real-axis. Hence, a magnitude normalization seems to be more convenient for this case. \n",
    "\n",
    "Despite this, if the real-part of the eigenvalue is positive, it is _unstable_. In other words, it describes an unstable operating condition of the power system. A standard normalization will map the real eigenvalues close to the imaginary axis (i.e., those having a small negative real-part) to $-1+0j$. This is not desirable since it hiddens the fact that the eigenvalue is close to the region of instability. For this reason, the normalization procedure carried in the cell below is an ad-hoc solution that normalizes only the eigenvalues outside of the unit circle in the complex plane: that is, only the eigenvalues whose magnitude is strictly larger or equal than one are normalized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normalized eigenvalues dataframe loaded!\n"
     ]
    }
   ],
   "source": [
    "current_wd = os.getcwd()\n",
    "\n",
    "if not os.path.exists(os.path.join(current_wd, \"99_Data/df_normalized_eigenvalues.pkl\")):\n",
    "\n",
    "    print(\"Normalizing eigenvalues...\")\n",
    "    # Creating new Pandas DataFrame for normalized eigenvalues\n",
    "    df_normalized_eigenvalues = pd.DataFrame([], columns = \\\n",
    "                                             list(df_eigenvalues.columns.values))\n",
    "\n",
    "    # Lambda expression to normalize only the eigenvalues outside the unit circle\n",
    "    norm_scalar = lambda eig : (eig / np.sqrt(np.power(eig.real, 2) + \\\n",
    "                                              np.power(eig.imag, 2))\n",
    "                        if np.sqrt(np.power(eig.real, 2) + np.power(eig.imag, 2)) >= 1\n",
    "                        else eig)\n",
    "\n",
    "    # Vectorizing the lambda expresion defined before\n",
    "    norm_vector = np.vectorize(norm_scalar)\n",
    "\n",
    "    # Normalizing eigenvalues\n",
    "    for sc in list(df_eigenvalues.columns.values):\n",
    "        df_normalized_eigenvalues[sc] = norm_vector(df_eigenvalues[sc][:])\n",
    "\n",
    "    # Saving normalized eigenvalues\n",
    "    with open(\"99_Data/df_normalized_eigenvalues.pkl\", 'wb') as f:\n",
    "        pickle.dump(df_normalized_eigenvalues, f, pickle.HIGHEST_PROTOCOL)\n",
    "    print(\"Normalized eigenvalues saved!\")\n",
    "    \n",
    "if os.path.exists(os.path.join(current_wd, \"99_Data/df_normalized_eigenvalues.pkl\")):\n",
    "    # DataFrame exists as a .pkl file - Loading it\n",
    "    with open(\"99_Data/df_normalized_eigenvalues.pkl\", 'rb') as f:\n",
    "        df_normalized_eigenvalues = pickle.load(f)\n",
    "    print(\"Normalized eigenvalues dataframe loaded!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating Ground Truth\n",
    "\n",
    "The labeling is carried out by identifying whether an eigenvalue is real- and complex-valued on the first place. Then, depending on whether the eigenvalue lies on the right-half or left-half plane (positive or negative real part), it is labeled as real/complex unstable or stable, respectively. All unstable eigenvalues are assigned a negative damping factor arbitrarily ($\\zeta = -1.1$). For stable complex-valued eigenvalues, we compute the damping factor using a closed-form expression given by\n",
    "\n",
    "$$\\begin{equation}\n",
    "    \\zeta_i = -\\frac{\\text{Re}\\{\\lambda_i\\})}{\\sqrt{(\\text{Re}\\{\\lambda_i\\})^2 + (\\text{Im}\\{\\lambda_i\\})^2}}\n",
    "\\end{equation}$$\n",
    "\n",
    "where $\\text{Re}\\{\\lambda_i\\}$ and $\\text{Im}\\{\\lambda_i\\}$ are the real and imaginary parts of the eigenvalue $\\lambda_i$. Stable real-valued eigenvalues have a $\\zeta$ larger than 1. Therefore, whenever an eigenvalue is identified as being stable and real, a damping factor $\\zeta = 1.1$ is assigned to it.\n",
    "\n",
    "Once the damping ratio is computed, we use it to create the labels of each of the categories, namely:\n",
    "* **Unstable** ($\\zeta < 0$): the eigenvalue lies on the right-half plane.\n",
    "* **Stable but critical ($0 \\leq \\zeta < 0.05$)**: the eigenvalue is stable. However, it exhibits light damping (since the damping ratio is small). Therefore, there could be significant oscillations in the system which is an undesired behavior and would require a corrective action.\n",
    "* **Acceptable operation ($0.05 \\leq \\zeta \\leq 0.1$)**: this operation corresponds to the scenario in which an eigenvalue is not as close to the imaginary axis and/or the oscillations caused by it do not compromise the stability of the system.\n",
    "* **Good operation ($0.1 \\leq \\zeta < 1$)**: the response caused by the eigenvalue does not compromise the current operating condition of the system.\n",
    "* **Satisfactory operation ($\\zeta \\geq 1$)**: the damping of the eigenvalue is large enough so that no oscillation will be caused on the system. \n",
    "\n",
    "The labeling results are stored in an array called `tag_label`. The function below employs the aforementioned damping ratio-based criterion to label the eigenvalues into the five pre-established categories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_ground_truth(df_eigenvalues, current_wd, df_name):\n",
    "    \n",
    "    # Path of the target file for storing the ground truth data\n",
    "    path_output = os.path.join(current_wd, \"99_Data/{}_ground_truth_data.pkl\".format(df_name))\n",
    "\n",
    "    if not os.path.exists(path_output):\n",
    "        \n",
    "        t_0 = time.time() # Initializing time measurement of algorithm execution\n",
    "        \n",
    "        # Number of scenarios\n",
    "        n_scenarios = df_eigenvalues.shape[1]\n",
    "        # Number of eigenvalues per scenario\n",
    "        n_eigs = df_eigenvalues.shape[0]\n",
    "\n",
    "        # Name of each scenario\n",
    "        scenarios = list(df_eigenvalues.columns)\n",
    "\n",
    "        # Creating label container\n",
    "        label = np.empty((n_eigs, n_scenarios), dtype = object)\n",
    "        # Creating damping ratio container\n",
    "        damping_ratio = np.zeros((n_eigs, n_scenarios))\n",
    "        # Creating tag (i.e., category) container\n",
    "        tag_label = np.zeros((n_eigs, n_scenarios))\n",
    "\n",
    "        for n_sc, sc in enumerate(scenarios):\n",
    "            for n_ev in range(0, n_eigs):\n",
    "                # Classification between real and complex conjugate eigenvalues\n",
    "                if (np.imag(df_eigenvalues[sc][n_ev]) == 0):\n",
    "                    label[n_ev][n_sc] = \"Real,\"\n",
    "                else:\n",
    "                    if (np.real(df_eigenvalues[sc][n_ev]) == 0):\n",
    "                        label[n_ev][n_sc] = \"Pure imaginary,\"\n",
    "                    else:\n",
    "                        label[n_ev][n_sc] = \"Complex conjugate,\" \n",
    "\n",
    "                # Classification between stable and unstable eigenvalues\n",
    "                if (np.real(df_eigenvalues[sc][n_ev]) > 0):\n",
    "                    label[n_ev][n_sc] = label[n_ev][n_sc] + \" unstable\"\n",
    "                else:\n",
    "                    label[n_ev][n_sc] = label[n_ev][n_sc] + \" stable\"\n",
    "\n",
    "                # Computation of damping ratio\n",
    "                if (label[n_ev][n_sc] == \"Complex conjugate, unstable\" or label[n_ev][n_sc] == \"Complex conjugate, stable\"):\n",
    "                    damping_ratio[n_ev][n_sc] = -np.real(df_eigenvalues[sc][n_ev]) / np.sqrt(np.square(np.real(df_eigenvalues[sc][n_ev])) + np.square(np.imag(df_eigenvalues[sc][n_ev])))\n",
    "                else:\n",
    "                    if (label[n_ev][n_sc] == \"Real, stable\"):\n",
    "                        damping_ratio[n_ev][n_sc] = 1.1\n",
    "                    else:\n",
    "                        damping_ratio[n_ev][n_sc] = -1.1\n",
    "\n",
    "                # Creating labels based on the damping ratio\n",
    "                if (damping_ratio[n_ev][n_sc] < 0):\n",
    "                    tag_label[n_ev][n_sc] = 1 # 1 - Unstable operation\n",
    "                elif (damping_ratio[n_ev][n_sc] < 0.05):\n",
    "                    tag_label[n_ev][n_sc] = 2 # 2 - Stable but critical operation\n",
    "                elif (damping_ratio[n_ev][n_sc] >= 0.05 and damping_ratio[n_ev][n_sc] < 0.1):\n",
    "                    tag_label[n_ev][n_sc] = 3 # 3 - Acceptable operation\n",
    "                elif (damping_ratio[n_ev][n_sc] >= 0.1 and damping_ratio[n_ev][n_sc] < 1.0):\n",
    "                    tag_label[n_ev][n_sc] = 4 # Good operation\n",
    "                elif (np.abs(np.real(df_eigenvalues[sc][n_ev])) < 0.5 and damping_ratio[n_ev][n_sc] >= 1.0):\n",
    "                    if np.abs(np.real(df_eigenvalues[sc][n_ev])) < 0.01:\n",
    "                        tag_label[n_ev][n_sc] = 6 # Not relevant\n",
    "                    else:\n",
    "                        tag_label[n_ev][n_sc] = 2 # Critical operation; real eigenvalue close to the imaginary axis\n",
    "                else:\n",
    "                    tag_label[n_ev][n_sc] = 5 # Satisfactory operation\n",
    "\n",
    "        # Saving ground truth data of eigenvalues    \n",
    "        ground_truth_data = {'text_label' : label,\n",
    "                          'damping_ratio' : damping_ratio,\n",
    "                          'tag_label' : tag_label}\n",
    "\n",
    "        with open(path_output, 'wb') as f:\n",
    "            pickle.dump(ground_truth_data, f, pickle.HIGHEST_PROTOCOL)\n",
    "        print(\"Ground truth data saved!\\n\")\n",
    "\n",
    "        t_f = time.time() # Finishing time measurement of algorithm execution\n",
    "        \n",
    "        print(\"Time elapsed: {} s\".format(t_f-t_0))\n",
    "        print(\"{} Eigenvalues (shape): {} (eigs) x {} (scenarios)\".format(df_name, df_eigenvalues.shape[0], df_eigenvalues.shape[1]))\n",
    "        print(\"{} Text labels (shape): {} (eigs) x {} (scenarios)\".format(df_name, label.shape[0], label.shape[1]))\n",
    "        print(\"{} Damping ratio (shape): {} (eigs) x {} (scenarios)\".format(df_name, damping_ratio.shape[0], damping_ratio.shape[1]))\n",
    "        print(\"{} Tag label (shape): {} (eigs) x {} (scenarios)\\n\".format(df_name, tag_label.shape[0], tag_label.shape[1]))\n",
    "    else:\n",
    "        print(\"Ground truth already exists!\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Improved Labeling Implementation\n",
    "This is a more efficient implementation in terms of computational efficiency for the hard-coded classifier. The time below is the one used for the score computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "import_libraries.py:1: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  # Mathematical tools\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time elapsed: 5.476852893829346 s\n"
     ]
    }
   ],
   "source": [
    "damping_scalar = lambda eig : (-np.real(eig) / np.sqrt(np.power(np.real(eig), 2) + np.power(np.imag(eig), 2)))\n",
    "\n",
    "damping_vector = np.vectorize(damping_scalar)\n",
    "t_0 = time.time()\n",
    "res = damping_vector(df_eigenvalues.values)\n",
    "t_f = time.time()\n",
    "\n",
    "print(\"Time elapsed: {} s\".format(t_f-t_0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ground truth data saved!\n",
      "\n",
      "Time elapsed: 46.803407430648804 s\n",
      "raw Eigenvalues (shape): 49 (eigs) x 19815 (scenarios)\n",
      "raw Text labels (shape): 49 (eigs) x 19815 (scenarios)\n",
      "raw Damping ratio (shape): 49 (eigs) x 19815 (scenarios)\n",
      "raw Tag label (shape): 49 (eigs) x 19815 (scenarios)\n",
      "\n",
      "Ground truth data saved!\n",
      "\n",
      "Time elapsed: 47.08603644371033 s\n",
      "normalized Eigenvalues (shape): 49 (eigs) x 19815 (scenarios)\n",
      "normalized Text labels (shape): 49 (eigs) x 19815 (scenarios)\n",
      "normalized Damping ratio (shape): 49 (eigs) x 19815 (scenarios)\n",
      "normalized Tag label (shape): 49 (eigs) x 19815 (scenarios)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "create_ground_truth(df_eigenvalues, current_wd, 'raw')\n",
    "create_ground_truth(df_normalized_eigenvalues, current_wd, 'normalized')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}