{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Energy System Modelling - Solutions to Tutorial II"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "**Imports**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "***\n",
    "## Problem II.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "**(a) Compile the nodes list and the edge list.**\n",
    "> **Remark:** While graph-theoretically both lists are unordered sets, let's agree on an ordering now which can serve as basis for the matrices in the following exercises: we sort everything in ascending numerical order, i.e.\\ node 1 before node 2 and edge (1,2) before (1,4) before (2,3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "nodes = [0, 1, 2, 3, 4, 5]\n",
    "edges = [(0, 1), (1, 2), (1, 4),\n",
    "         (1, 5), (2, 3), (2, 4),\n",
    "         (3, 4), (4, 5)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "**(b) Determine the order and the size of the network.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "N = len(nodes)\n",
    "E = len(edges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Order: 6\n",
      "Size: 8\n"
     ]
    }
   ],
   "source": [
    "print(\"Order: {}\\nSize: {}\".format(N, E))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "**(c) Compute the adjacency matrix $A$ and check that it is symmetric.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 1., 0., 0., 0., 0.],\n",
       "       [1., 0., 1., 0., 1., 1.],\n",
       "       [0., 1., 0., 1., 1., 0.],\n",
       "       [0., 0., 1., 0., 1., 0.],\n",
       "       [0., 1., 1., 1., 0., 1.],\n",
       "       [0., 1., 0., 0., 1., 0.]])"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.zeros((N, N))\n",
    "\n",
    "for u, v in edges:\n",
    "    A[u, v] += 1\n",
    "    A[v, u] += 1\n",
    "\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A == A.T).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "**(d) Find the $k_n$ of each node $n$ and compute the average degree of the network.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 4., 3., 2., 4., 2.])"
      ]
     },
     "execution_count": 147,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k = A.sum(axis=1)\n",
    "k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.6666666666666665"
      ]
     },
     "execution_count": 148,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "**(e) Determine the incidence matrix $K$ by assuming the links are always directed from smaller-numbered node to larger-numbered node, i.e.\\ from node 2 to node 3, instead of from 3 to 2.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [-1.,  1.,  1.,  1.,  0.,  0.,  0.,  0.],\n",
       "       [ 0., -1.,  0.,  0.,  1.,  1.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0., -1.,  0.,  1.,  0.],\n",
       "       [ 0.,  0., -1.,  0.,  0., -1., -1.,  1.],\n",
       "       [ 0.,  0.,  0., -1.,  0.,  0.,  0., -1.]])"
      ]
     },
     "execution_count": 150,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K = np.zeros((N,E))\n",
    "\n",
    "for i, (u, v) in enumerate(edges):\n",
    "    K[u,i] = 1\n",
    "    K[v,i] = -1\n",
    "    \n",
    "K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "**(f) Compute the Laplacian $L$ of the network using $k_n$ and $A$. Remember that the Laplacian can also be computed as $L=KK^T$ and check that the two definitions agree.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "D = np.diag(A.sum(axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1., -1.,  0.,  0.,  0.,  0.],\n",
       "       [-1.,  4., -1.,  0., -1., -1.],\n",
       "       [ 0., -1.,  3., -1., -1.,  0.],\n",
       "       [ 0.,  0., -1.,  2., -1.,  0.],\n",
       "       [ 0., -1., -1., -1.,  4., -1.],\n",
       "       [ 0., -1.,  0.,  0., -1.,  2.]])"
      ]
     },
     "execution_count": 152,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = D - A\n",
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 153,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array_equal(L, K.dot(K.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "**(g) Find the diameter of the network by looking at the graph.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "By looking: Between nodes $0$ and $3$, f.ex. $0 \\to 1 \\to 2 \\to 3$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "## Problem II.3\n",
    "\n",
    "If you map the nodes to `0=DK, 1=DE, 2=CH, 3=IT, 4=AT,5=CZ` the network represents a small part of the European electricity network (albeit very simplified). On the [course homepage](https://nworbmot.org/courses/complex_renewable_energy_networks/), you can find the power imbalance time series for the six countries for January 2017 in hourly MW in the file `imbalance.csv`. They have been derived from physical flows as published by [ENTSO-E](https://transparency.entsoe.eu/transmission-domain/physicalFlow/show)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The linear power flow is given by\n",
    "$$p_i = \\sum_j \\tilde{L}_{i,j}\\theta_j \\qquad \\text{and} \\qquad f_l = \\frac{1}{x_l} \\sum_i K_{i,l}\\theta_i, \\qquad \\text{where} \\qquad \\tilde{L}_{i,j}= \\sum_l = K_{i,l}\\frac{1}{x_l} K_{j,l}$$\n",
    "is the weighted Laplacian. For simplicity, we assume identity reactance on all links $x_l = 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "**Read data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>DK</th>\n",
       "      <th>DE</th>\n",
       "      <th>CH</th>\n",
       "      <th>IT</th>\n",
       "      <th>AT</th>\n",
       "      <th>CZ</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DateTime</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2017-01-01 00:00:00</th>\n",
       "      <td>-79.73</td>\n",
       "      <td>3486.25</td>\n",
       "      <td>-3403.37</td>\n",
       "      <td>-97.5</td>\n",
       "      <td>-2867.33</td>\n",
       "      <td>2961.68</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-01-01 01:00:00</th>\n",
       "      <td>-109.14</td>\n",
       "      <td>4248.17</td>\n",
       "      <td>-3829.59</td>\n",
       "      <td>-105.9</td>\n",
       "      <td>-3149.98</td>\n",
       "      <td>2946.44</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-01-01 02:00:00</th>\n",
       "      <td>-101.49</td>\n",
       "      <td>5314.56</td>\n",
       "      <td>-4131.66</td>\n",
       "      <td>-199.2</td>\n",
       "      <td>-3646.87</td>\n",
       "      <td>2764.66</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-01-01 03:00:00</th>\n",
       "      <td>-197.96</td>\n",
       "      <td>5450.71</td>\n",
       "      <td>-4009.61</td>\n",
       "      <td>-656.0</td>\n",
       "      <td>-3345.56</td>\n",
       "      <td>2758.42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-01-01 04:00:00</th>\n",
       "      <td>-526.15</td>\n",
       "      <td>5695.72</td>\n",
       "      <td>-3924.43</td>\n",
       "      <td>-646.2</td>\n",
       "      <td>-3393.12</td>\n",
       "      <td>2794.18</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         DK       DE       CH     IT       AT       CZ\n",
       "DateTime                                                              \n",
       "2017-01-01 00:00:00  -79.73  3486.25 -3403.37  -97.5 -2867.33  2961.68\n",
       "2017-01-01 01:00:00 -109.14  4248.17 -3829.59 -105.9 -3149.98  2946.44\n",
       "2017-01-01 02:00:00 -101.49  5314.56 -4131.66 -199.2 -3646.87  2764.66\n",
       "2017-01-01 03:00:00 -197.96  5450.71 -4009.61 -656.0 -3345.56  2758.42\n",
       "2017-01-01 04:00:00 -526.15  5695.72 -3924.43 -646.2 -3393.12  2794.18"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imbalance = pd.read_csv('imbalance.csv', index_col=0, parse_dates=True)\n",
    "imbalance.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "\\begin{equation}\n",
    "p_u = \\sum_v L_{u,v} \\theta_v\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "**(a) Compute the voltage angles $\\theta_j$ and flows $f_l$ for the first hour in the dataset with the convention of $\\theta_0 = 0$; i.e. the slack bus is at node 0.**\n",
    "> **Remark:** Linear equation systems are solved efficiently using `numpy.linalg.solve`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Calculate the *voltage angles* first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 3486.25, -3403.37,   -97.5 , -2867.33,  2961.68])"
      ]
     },
     "execution_count": 155,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imbalance.iloc[0].values[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([    0.        ,    79.73      , -2302.97857143, -1995.25809524,\n",
       "       -1590.03761905,   725.68619048])"
      ]
     },
     "execution_count": 156,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta = np.r_[0, np.linalg.solve(L[1:,1:], imbalance.iloc[0].values[1:])]\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Then, calculate the *flows*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  -79.73      ,  2382.70857143,  1669.76761905,  -645.95619048,\n",
       "        -307.72047619,  -712.94095238,  -405.22047619, -2315.72380952])"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flows = K.T.dot(theta)\n",
    "flows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "**(b) Determine the average flow on each link for January 2017 and draw it as a directed network**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "flows = K.T.dot(np.vstack([np.zeros((1, len(imbalance))), np.linalg.solve(L[1:,1:], imbalance.values[:,1:].T)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 744)"
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flows.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  451.02294355,  2402.42223502,  2124.94871608,   426.46578277,\n",
       "         421.1171115 ,  -277.47351895,  -698.59063044, -1698.48293331])"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "avg_flows = flows.mean(axis=1)\n",
    "avg_flows"
   ]
  }
 ],
 "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.6.4"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": false,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  },
  "toc_position": {
   "height": "47px",
   "left": "595px",
   "right": "20px",
   "top": "-14px",
   "width": "318px"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}