{ "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": [ "\$$\n", "p_u = \\sum_v L_{u,v} \\theta_v\n", "\$$" ] }, { "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 }