Problem 4.1

Integrated Energy Grids

Problem 4.1

This is a continuation of Problem 3.2 from Lecture 3: Let us assume now that we are in an hour with an excess of wind generation in Denmark and a deficit in other countries so that the power injection \(p_i\) of the different countries is as follows: Germany= - 2 MW, DK1=5 MW, DK2=6 MW, Norway = - 8 MW, Sweden = -1 MW.

Determine the voltage angles \(\theta_i\) and the flows \(p_l\) in the lines of the network. Assume that \(\theta_0\)=0; i.e. the reference bus is at node 0 (Germany); and the reactance in all links are \(x_l\)=1.

We will use the package numpy to operate with matrices.

import numpy as np
import numpy.linalg

Calculate list of nodes, links, degree, adjacency and Laplacian matrix. (this was already implemented in Problem 3.2)

nodes=[0,1,2,3,4]
links=[(0,1), (1,2), (1,3), (1,4), (2,4)]

D = np.zeros((len(nodes), len(nodes)))

for node in nodes:
    D[node, node] = sum([1 if node in link else 0 for link in links])

A = np.zeros((len(nodes), len(nodes)))

for node_a, node_b in links:
    A[node_a, node_b] = 1
    A[node_b, node_a] = 1

L = D - A
L

array([[ 1., -1.,  0.,  0.,  0.],
       [-1.,  4., -1., -1., -1.],
       [ 0., -1.,  2.,  0., -1.],
       [ 0., -1.,  0.,  1.,  0.],
       [ 0., -1., -1.,  0.,  2.]])

K = np.zeros((len(nodes),len(links)))

for i, (node_a, node_b) in enumerate(links):
    K[node_a,i] = 1
    K[node_b,i] = -1
    
K

array([[ 1.,  0.,  0.,  0.,  0.],
       [-1.,  1.,  1.,  1.,  0.],
       [ 0., -1.,  0.,  0.,  1.],
       [ 0.,  0., -1.,  0.,  0.],
       [ 0.,  0.,  0., -1., -1.]])

Power flow analysis

We know the power injection pattern for the nodes \(p_i\) and we want to determine the power flows \(p_l\) in the lines of the network

The reactance in every link is \(x_l\)=1 so \(L_{ij}=\sum_{l} K_{il} \frac{1}{x_l}K_{lj}=\sum_{l} K_{il} K_{lj}\) or \(L_{ij}=D_{ij}-A_{i,j}\)

First, we calculate the voltage angles, by solving the set of linear equations. \(p_i=\sum L_{i,j} \theta_j\)

p_i = [-2, 5, 6, -8, -1]

The node 0 (Germany) is the slack bus, so \(\theta_0\)=0.

theta = np.r_[0, np.linalg.solve(L[1:,1:], p_i[1:])]
theta

array([ 0.        ,  2.        ,  5.66666667, -6.        ,  3.33333333])

Second, we use the voltage angles to compute the flows.

\(p_l=\frac{1}{x_l}\sum_{j}K_{lj}\theta_j\)

p_l= K.T.dot(theta)
p_l

array([-2.        , -3.66666667,  8.        , -1.33333333,  2.33333333])

Additional discussion on the inversion of the Laplacian matrix

The Laplacian matrix is not invertible.

We can only invert a matrix if it is non-singular, i.e., the determinant should not be zero.

You can also check that the rank of the matrix is 4 but the matrix has 5 columns, so the determinant is zero and we can not find the inverse

np.linalg.matrix_rank(L)

np.int64(4)

Instead of selecting Germany as the slack bus, we could have used the Moore Penrose pseudo-inverse to invert the Laplacian.

In this case, we don’t get bus zero as slack and the theta values in the solution are shifted with respect to the previous solution, but we get the same results for the power flows.

theta_ = np.linalg.pinv(L).dot(p_i)
theta_

array([-1.        ,  1.        ,  4.66666667, -7.        ,  2.33333333])

P_l= K.T.dot(theta_)
P_l

array([-2.        , -3.66666667,  8.        , -1.33333333,  2.33333333])