{ "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": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DKDECHITATCZ
DateTime
2017-01-01 00:00:00-79.733486.25-3403.37-97.5-2867.332961.68
2017-01-01 01:00:00-109.144248.17-3829.59-105.9-3149.982946.44
2017-01-01 02:00:00-101.495314.56-4131.66-199.2-3646.872764.66
2017-01-01 03:00:00-197.965450.71-4009.61-656.0-3345.562758.42
2017-01-01 04:00:00-526.155695.72-3924.43-646.2-3393.122794.18
\n", "
" ], "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 }