Correlated Failures of Power Systems: Analysis of the Nordic Grid

(1)

Correlated Failures of Power Systems:

Analysis of the Nordic Grid

MARTIN ANDREASSON

Master's Degree Project

Stockholm, Sweden February 2011

(2)

Abstract

The emphasis of this master's thesis is modeling and simulation of failures in large-scale power grids. The linear DC-model governing the active power ows is derived and dis-cussed, and the optimal load shedding problem is introduced. Using Monte Carlo simulations, we have analyzed the eects of correlations between failures of power lines on the total system load shed. Correlations are introduced by a Bernoulli failure model with its rst two ordinary moments given explicitly. The total system load shed is determined by solving the optimal load shedding problem in a MATLAB environment using YALMIP and the GLPK solver. We have introduced a Monte Carlo simula-tion framework for sampling the statistics of the system load shed as a funcsimula-tion of stochastic network parameters, and provide explicit guarantees on the sampling accu-racy. This framework has been applied to a 470 bus model of the Nordic power grid. It has been found that increased correlations between Bernoulli failures of power lines can dramatically increase the expected value as well as the variance of the system load shed.

(3)

Introduction

1.1 Previous work

Power systems are among the largest and most complex systems created by mankind. The complexity of power systems keeps increasing as power grids are expanded and new

functionalities are added, as with the development of the SmartGrid [Massoud and Wollenberg, 2005]. Because many vital parts of today's society require reliable supply of electricity, the

reliable and secure operation of power systems is essential [Amanullah et al., 2005, Amin, 2002]. We give a brief overview of the research in the two distinctive areas of security against adversarial attacks and reliability of power systems.

In the area of security, research on characterizing optimal attack and defense strate-gies has gained momentum over the past years. [Salmeron et al., 2004] consider the optimization problem of maximizing the power outage, for a given number of power transmission lines that an adversary is capable of disconnecting. The system opera-tor is assumed to take the best action to minimize the damage in form of compulsory load shedding. The problem is by nature game theoretic and gives rise to a maximin optimization problem, where the outer maximization seeks the most disruptive attack for a given budget of the adversary, and the inner minimization solves the optimal load shedding problem which minimizes the consequences of an attack. Because of the non-convexity and the existence of integer variables, the problem is inherently hard to solve for large systems. [Pinar et al., 2010] approximates the nonlinear mixed integer bi-level program by a mixed integer linear program, and derive an upper bound on the severity of adversarial attacks. [Arroyo and Galiana, 2005] considers a fairly more general formulation of the terrorist threat problem, where the terrorist's and system operator's objectives are not necessarily antagonistic. This is a natural extension of the maximin optimization formulation, and the formulation leads to a bi-level optimization program. The authors have considered the linear DC model, and thus the bi-level opti-mization problem is linear. By replacing the inner optiopti-mization by it's KKT conditions

(7)

and rewriting them as linear inequalities, the problem is transformed into a single-level mixed integer linear optimization program, which is solved by standard commercial software.

Traditionally, the reliability of power systems has often been characterized by de-terministic means, such as the widely uses N − k criterion [Billinton and Allan, 1995]. A power system satisfying the N − k criterion is able to withstand any contingency consisting of k outages. The advantage with deterministic reliability criteria is that rubustness is guaranteed explicitely, as with the N − k criterion. The main drawback of deterministic reliability criteria however is that they do not take into account the probabilities of contingencies. Furthermore the number of events which have to be con-sidered grows exponentially in k, making the reliability evaluation computationally in-tractable. More recent research on the reliability of power systems has emphasized that many events governing the reliability of power systems are by nature stochastic, e.g. demands and generation capacities. Various statistical and sampling based methods for evaluating the reliability of power systems have been developed to analyze stochas-tic phenomena in power systems [Roy Billinton, 1994, Allan and Billinton, 2000]. In [Billinton and Wang, 1999] a two state Markov model for the failure of various power system elements is considered, and the statistics of the power system are calculated using Monte Carlo techniques. [Chertkov et al., 2010] considers the power demand as a random variable, to determine the likelihood of the most common causes of power shortage in a power system by a heuristic optimization algorithm. The pa-per focuses on the linear DC model and assumes Gaussian demand distributions. [Billinton and Khan, 1992] studies the reliability of a 5-bus power system by consider-ing nearly all combinations of contconsider-ingencies, and calculatconsider-ing the load shed needed to restore the system to a secure operational state. The power system reliability is dened by the likelihood of the discrete operational states of normal, restorative, alert, emer-gency and extreme emeremer-gency respectively, as dened in the paper. The amount of load shedding is determined by a linear program, by using the linearized DC-model. One drawback with this approach is the dramatically increased complexity as the network size increases. In [da Silva et al., 2004] the authors use a Monte Carlo simulation based approach to quantify the reliability of a power system. Reliability here is quantied by the states healthy, marginal and at risk. The proposed algorithm samples states based on a probability distribution, after which optimal power dispatch is performed. In [Anghel et al., 2007], the authors use Monte Carlo methods to model cascading fail-ures and major blackouts of power grids. Their model covers random removal of system components, as well as random load uctuations. The load uctuations are translated into heat variations of the power transmission lines by the use of a heat diusion model, which causes random failures of the power transmission lines. To capture the dynam-ics of cascading failures, the system is simulated over a time period, and recoveries as

(8)

well as protective actions by the transmission operator (TO) are modeled. The algo-rithms are applied to a 100-node power transmission system and examples of predicted blackouts are demonstrated. Other than failure correlations due to cascading failures [Kinney et al., 2005], we found no model which attempts to model correlated failures in power systems.

The Monte Carlo techniques in this work rely on the deterministic operation of op-timal load shedding in power systems. Opop-timal load shedding for system protection is a form of an optimal power ow problem and is a classical problem in the power systems community. For a survey of studied problems see e.g. [Huneault and Galiana, 1991]. The optimal load shedding problem is a system protection measure to, in some mean-ingful sense, optimally dispatch power in the network to avoid system faults. The classical formulations of the optimal load shedding problem consider a static network, where the decision variables are the phase angles of the nodes in the power network, and the objective is a (weighted) sum of power decits in the nodes. [Fisher et al., 2008] consider a more general formulation of the optimal load shedding problem for system protection, where they take into account the switching of power transmission lines, so called optimal transmission switching (OTS). The binary variables associated with the on-o state of each power line together with the linearized power ow model, give rise to a mixed integer linear program, which may be solved by standard commercial software.

The OTS can be generalized by considering the use of exible AC transmission sys-tems (FACTS). FACTS is a collection of technologies used to control the AC power ows in transmission lines. The power ow through the power network is determined solely by the phase angles and voltages of the buses, and the power transmission line admittances. See e.g. [Abur and Exposito, 2004] for a review of the Kircho voltage law (KVL) and the power ow equations. Due to the KVL there are no direct ways of routing electricity, as with e.g. trac ows, and hence the control of power transmission networks is fairly inexible. FACTS technologies however, open the possibility to control the line admit-tances in close to real time, allowing to control the power ow from the nodes, without changing the phase angles of any other nodes. The line admittances in FACTS are con-trolled by changing adjustable inductors and capacitors installed in the transmission lines. The FACTS and applications to power network control technologies are described more in detail in [Hingorani, 1993, Gotham and Heydt, 1998, Xiao et al., 2002]. At-tempts to incorporate FACTS devices in existing optimal power dispatch models have been made in [Preedavichit and Srivastava, 1998,Lu and Abur, 2002].

(9)

1.2 Overview of work

This work aims at studying the eects of increasing correlations between failures of power system components. In contrast to previous papers on the subject, we introduce a failure model which explicitly takes into account correlated failures. In particular, we consider a Bernoulli model of correlated power line failures. We measure the impact of correlations of failures by the covariances between the failures.

Many of the SCADA systems used in controlling power systems are running general-purpose operating systems, e.g. Microsoft Windows based operating systems and general-purpose software. This makes SCADA systems vulnerable to software bugs within the operating system and control applications. When identical software is de-ployed in several system components, software failures are likely to be more correlated between these components. We feel that the eects of correlated failures of power sys-tems have not been adequately studied, although there is reason to believe that many failures in power systems could be correlated. In this work, we measure the impact of a system failure by the minimum system load shed required to restore the system to a safe state. This formulation gives rise to an optimization problem, which under some conditions can be made linear. For dierent values of the correlations of the failure distribution, we compute the sampled statistics of the total system load shed by Monte Carlo techniques, and provide guarantees on the convergence rate of the sampled statis-tics. In particular, we use a weighted sum of the mean and variance of the total system load shed as a risk measure of the failure statistics. To obtain statistical data from a realistic power system, we apply our techniques to a 470 bus model of the Nordic power system, acquired from publicly available sources. Certain information of the Nordic power grid are available and can be put together to build a full model, whereas other information is unavailable and needs to be estimated. Applying the Monte Carlo simulation techniques to the Nordic grid model, we have found that increasing corre-lations between Bernoulli failures of power lines lead to increased expected value and variance of the system load shed.

1.3 Main contributions

This work presents a novel model of the Nordic power grid, used for studying correlated failures of power lines. While examples from the Nordic power grid have been used widely in power systems research, we know of no full-scale publicly available model of the Nordic power system for simulating static power ows.

Using the novel model of the Nordic grid, we have simulated correlated Bernoulli failures of power lines. Besides the novelty of the model and problem formulation, we have found that increased correlations lead to increased costs measured in the mean and

(10)

the variance of the system load shed. We furthermore show sucient conditions, for when increasing correlations between Bernoulli failures of power lines lead to increased expected system load shed and increased variance of the system load shed.

(11)

Chapter 2

Power ow analysis

2.1 Graph models for power networks

Denition 2.1. A weighted, directed graph is a set G = (V, E)where E ∈ V × V

v ∈ V are called the vertices of, and e ∈ E are called the edges of the graph. The number of vertices V and edges E are denoted |V | and |E|.

Denition 2.2. An edge eij going from vertex vito vertex vj with weight wij is denoted

eij = vi, vj, wij

Denition 2.3. An graph is said to be unweighted if wij = 1 ∀ i, j : eij ∈ E, and

undirected if eij = eji ∀ i, j : eij ∈ E

Denition 2.4. The degree, deg(·) of a vertex is the number of edges connected to the vertex.

Denition 2.5. The vertex-edge incidence matrix for a directed graph is dened as: A = [Aij] where Aij =      1 if ei= (vj, u) ∈ E −1 if ei= (u, vj) ∈ E 0 otherwise for some node u ∈ V .

Denition 2.6. The weighted Laplacian is dened as

Lw = [lij] where lij =            P

jwij for i where vj is adjacent to vi if i = j

−wij if eij ∈ E

−wji if eji ∈ E

(12)

Theorem 1.

Lw = ATDwA (2.1)

where Dw= diag([wij]).

For more details on graphs and for proofs, the reader is referred to AppendixC.1. Denition 2.7. For power systems, the power buses, also referred to as nodes, are modeled as vertices. The power transmission lines are modeled as edges.

2.2 Power ow equations

Power systems can be described using graph models, by letting the buses of the power system be modeled as vertices, and the power transmission lines be modeled as edges in the graph. We use the following notation for the power ow model:

Pij : Real power ow from bus i to j

Qij : Reactive power ow from bus i to j

P: Vector of total active power ow injected to the buses Vi : Voltage of bus i

gij+ jbij : Admittance of the series branch between bus i and j

gsj+ jbsj : Admittance of the shunt branch at bus i

θi : Phase angle of bus i

θ : Vector of phase angles, θ = [θi]

θij = θi− θj : Phase angle dierence between bus i and j

From [Abur and Exposito, 2004] we obtain the nonlinear model for the real and reactive power ow from bus i to bus j:

Pij = Vi2(gsi+ gij) − ViVj(gijcos θij− bijsin θij) (2.2)

Qij = −Vi2(bsi+ bij) − ViVj(gijsin θij − bijcos θij) (2.3)

By adding the incoming and subtracting the outgoing power ows into each node, we obtain the equations for the real and reactive power injections for a node i.

Pi = X j:(ej,ei)∈E V2 i (gsi+ gij) − ViVj(gijcos θij − bijsin θij) − X j:(ei,ej)∈E V2 i (gsi+ gij) − ViVj(gijcos θij − bijsin θij) (2.4) Qi = X j:(ej,ei)∈E −V2 i (bsi+ bij) − ViVj(gijsin θij − bijcos θij − X j:(ei,ej)∈E −V_i2(bsi+ bij) − ViVj(gijsin θij − bijcos θij) (2.5)

(13)

2.2.1 Nonlinear DC-model

Following [Pinar et al., 2010], we assume the resistance-to-reactance ratio to be suf-ciently small. Furthermore we assume that the admittance in the shunt branch is negligible. In (2.2)-(2.3) this translates to: gij = gsj = bsj = 0. Imposing these

approximations, we get the following power ow equations:

Pij = ViVjbijsin θij (2.6)

Qij = ViVjbij(cos θij − 1) (2.7)

Equations (2.6) - (2.7) are often referred to as the nonlinear DC-model. We may write these equations in terms of power ow into the nodes, by using a graph model of the power transmission network. Let the topology of the power transmission network be given by the vertex-edge incidence matrix A. The line admittances are represented as weights of the power transmission lines as B = diag(bij), and the node voltages as

Vline _{= diag V} iVj

. The phase angles are represented as as Θ = θij

. It is easily seen that Θ = Aθ, where θ is a vector of the absolute phase angles in the nodes and A is the node art incidence matrix. Hence we can write:

Pline _{= V}line_{B sin (Aθ)} _(2.8)

Qline _{= V}line_{B cos (Aθ) − 1}

(2.9) where Pline ₌_P

ij , Qline =Qij

and sin(x) = sin(xi) , cos(x) = cos(xi)

. Note that Vline _{is easily calculated once A is known, since the rows of A contain the}

infor-mation of which nodes are connected. 2.2.2 Linear DC-model

For linearizing equations (2.6)- (2.7), we dene the following functions:

FP(θij|V, bij) = Vlinebijsin θij (2.10)

FQ(θij|V, bij) = Vlinebij(cos θij − 1) (2.11)

Linearizing (2.6)- (2.7) around θij = 0 ∀i, j, we obtain the following equations, valid

for suciently small values of θij:

Pij = FP(0|V, bij) + θij ∂FP(θij|V, bij) ∂θij θij=0 = Vlinebijθij (2.12) Qij = FQ(0|V, bij) + θij ∂FQ(θij|V, bij) ∂θij θij=0 = 0 (2.13)

(14)

Equation (2.12) is commonly referred to as the linear DC model. As in the nonlinear case, we may write equation (2.12) in vector form:

Pline= VlineBAθ (2.14)

By adding the power ows to each node, we may write the power injections in the nodes as:

P = ATVlineBAθ = LBθ (2.15)

where V = [Vi]Note that LB = ATVlineBAis analogue to the weighted Laplacian of

the power transmission network.

2.2.3 Transmission constraints in power systems

All power systems are endowed with hard physical transmission constraints, restricting the active power ows. For simplicity but without loss of generality, we assume that the nodes in the power system are partitioned into generator nodes, Pg ₀_{, and load}

nodes, Pl ₀_{, i.e.} P = " Pg Pl # (2.16) Line capacity constraints

Typically, there are physical capacity bounds on the power lines that take the form:

− P_maxline Pline P_maxline (2.17)

Generation constraints

Each generator i has a maximum capacity Pg

max,i, limiting the maximum power

pro-duction. Furthermore we assume that the minimum production for each generator is 0. Thus the constraints on the productions take the form

0 Pg P_maxg (2.18)

where 0 = [0, . . . , 0]T

Load constraints

We assume that each load node i has a predened demand, Pl

d,i. Furthermore we

assume that load nodes cannot be supplied with more power than the demand, and that load nodes, per denition, cannot produce any power. Hence the load constraints take the form

Pl

d Pl 0 (2.19)

(15)

DC-model constraints

For the DC-model to be valid, the phase angle dierences must be suciently small. By using the vertex-edge incidence matrix, this can be formulated as a single vector constraint

− ∆θmax· 1 Aθ ∆θmax· 1 (2.20)

where 1 = [1, . . . , 1]T_{, and 0 < ∆θ}

max ≤ π₂ is a suciently small real number.

2.2.4 Transmission constraints under the linear DC-model

By using the linearized power ow equations (2.14)-(2.15), we can express Pline_{, P}l

and Pg _{as linear functions of the phase angles θ. Furthermore, since all constraints}

mentioned in section 2.2.3 are linear inequalities, we may write the intersection of all constraints as a set of linear inequalities. Dene npas the number of power transmission

lines, ngand nlas the number of generator and load nodes respectively. Let n = ng+nl

be the total number of nodes. The feasibility constraints of equations (2.17)-(2.20) can be written as Cθ d (2.21) where C =           Vline_BA −Vline_BA LB −LB A −A           d =                Pline max Pline max Pmaxg 0nl×1 0ng×1 −Pl d ∆θmax· 1np×1 ∆θmax· 1np×1                (2.22)

Denition 2.8. The feasible region of a power network is dened as F =_{θ ∈ R}n_{|Cθ d}

Denition 2.9. A state θ ∈ Rn _{is called feasible i}

θ ∈ F If a state θ is not feasible, it is called infeasible.

The power ow in a power transmission network is dependent on the topology of the power network and load and generation data. We dene the conguration C of a power system.

(16)

Denition 2.10. The conguration C of a power system is the set C =nV, A, B, P_maxline, P_dl, P_maxg o

Depending on the conguration C, it is either possible or not possible to satisfy all power demands. We make the following denition to formalize the intuition.

Denition 2.11. A conguration C is called satisable i ∃θ ∈ F : HlLBθ = Pdl

A conguration which is not satisable is called unsatisable.

2.3 Linear optimal load shedding

Recall that a conguration is satisable i there exists a feasible state such that all load demands are satised. If a conguration is unsatisable, not all demands can be met. In such cases, it is up to the transmission system operator (TSO) to cut the loads for some load nodes. This procedure is called load shedding. When load shedding is necessary, the TSO will try to minimize the total amount of load shedding needed to restore the system to a feasible state. We dene the associated optimization problem, called the optimal load shedding problem. It is easily veried that the corresponding nonlinear optimal load shedding problem is in general non-convex, and hence inherently hard to solve. For this reason, we will only consider the linear optimal load shedding problem henceforth.

Denition 2.12. The linear optimal load shedding problem is dened as the linear program min θ c T_θ _(2.23) s.t. Cθ d where c =h01×ng ₁1×nl i LB (2.24)

and c, C and d are as dened in equation (2.22). For a review of linear programming, the reader is referred to appendix C.4.

The objective function ensures that the load of the demand nodes is minimized, since the demand loads by denition are negative. The linear constraints represent the constraints of the power system, discussed in 2.2.3. Given a disturbance of the

(17)

power system, which changes the conguration of the system, the resulting change to the solution to (2.23) can be interpreted as a measure of how severe the disturbance is. When no load shedding is necessary, this corresponds to normal operation of the power system, while 100% load shed corresponds to a complete blackout where no power is transmitted.

Example 2.3.1. Consider a 6-node and 7-line power transmission network, where the nodes are partitioned into generator nodes (g1, g2), and load nodes (l1, l2, l3, l4) For simplicity, we assume unit voltage V = 1 for all nodes. Given a power demand

g1 l1 l2 l3 l4 g2 l1 l2 l3 l4 l5 l6 l7

Figure 2.1: An example of a power network where the maximum line capacities are: Pline

max,1 = 1, Pmax,2line = 1, Pmax,3line = 1, Pmax,4line = 2, Pmax,5line = 1, Pmax,6line = 3, Pmax,7line = 3.

The line admittances are: b1 = 1, b2 = 1, b3 = 1, b4 = 1, b5 = 2, b6 = 2, b7 = 2. We

assume no generation capacity bounds −Pl

d= [2, 2, 1, 1]T, the optimal power ow policy is given by the phase angles:

θ =           2.6436 1.6782 1.2055 0.5891 −0.2055 0.0000           (2.25)

(18)

and the power generations and ows are: Pg = " 2.8764 3.1236 # Pline=             0.6164 0.7945 −0.2055 −1.2055 −0.9455 2.1782 2.8764             (2.26)

Since all demands are met, the system is feasible, i.e. load shedding is not necessary. The value of the objective function is hence cT_{θ =} _{P P}l

d. The optimal power ow is

illustrated in gure 2.2. g1 l1 l2 l3 l4 g2 0.6164 0.7945 0.2055 1.2055 0.9455 2.1782 2.8764

Figure 2.2: The optimal power ows shown graphically

2.4 Switched optimal linear load shedding

(19)

g1

l1 l2

l3

Figure 2.3: A 3 node, 3 line power transmission network.

The generation capacity of g1 is 2, and the demands of l1 and l2 are both 2. The

transmission line capacities for the lines l1, l2 and l3 are 2, 2, 0.1 respectively, and

the line admittances are 1, 2, 1 respectively. The solution to the linear optimal load shedding problem is illustrated in gure 2.4

g1

l1 l2

0.65 1.1

0.1

Figure 2.4: The optimal power ows shown graphically

Since node l1 only receives 0.75 but its demand is 1, the total load shed is 0.25.

Example 2.4.2. Now consider the same power transmission network as in example 2.4.1, but with line l3 cut o.

(20)

g1

l1 l2

Figure 2.5: A 3 node, 2 line power transmission network.

The solution to the linear optimal load shedding problem is illustrated in gure 2.6 g1

l1 l2

1 1

Figure 2.6: The optimal power ows shown graphically Since both demands are satised, the total load shed is 0.

These preceding examples show that it may be benecial for the network operator to switch o power transmission lines in order to reduce the system load shed. We dene the associated optimization problem where it is possible to switch power transmission lines on or o.

Denition 2.13. The switched optimal load shedding problem is dened as the mixed integer program

min

θ,Y c

T_θ _(2.27)

(21)

where c and d are as dened in equation (2.22), and Cs=           Vline_{BY A} −Vline_{BY A} LS B −LS B A −A           (2.28) where LS_B= ATVlineBY A (2.29) and

Y = diag([y1, . . . , ynp]) where yi∈ {0, 1} ∀i ∈ {1, . . . , np} (2.30)

Lemma 2. The optimal value of the optimization program (2.27) is less than or equal to the optimal value of the optimization program (2.23).

Proof. Assume, for the sake of contradiction, that the optimal solution of (2.23) is strictly lower than the solution of (2.27). Clearly the optimizer for (2.23) is also a feasible solution to (2.27), with Y = I. Thus, the optimal solution for (2.27) is at most equal to the optimal solution of (2.23), contradicting the previous assumption.

Even though power line switching may be practically feasible for power transmission networks, the switched optimal linear load shedding problem becomes computationally hard as the number of transmission lines grows. Also, since the placement of new power transmission lines is carefully considered before construction, it is unlikely that the disconnection of some transmission line(s) will reduce the system load shed in practice. Nevertheless, the switched optimal linear load shedding problem is a generalization of the linear optimal load shedding problem and hence worth mentioning.

2.5 Power planning with uncertain demand and generation

Typically in most power grids, both demand and generation capacities vary over time, and are not completely predictable. Demands usually vary depending on the time of the day, the weather, season and other factors. Generation capacity is approximately constant for traditional generators, such as gas turbines and nuclear power plants, whereas the availability of many renewable generators such as water and wind power exhibit major variations over time. By using their balancing controls, TSOs can take real time actions to reschedule power ows when demand and generation uctuations

(22)

occur. These uctuations require sucient balancing power capacities. We will show that the necessary balancing power can be approximated by using the framework of chance constrained optimization. For simplicity we only assume variations of the de-mands, but our methods can be extended to the case where some generation capacities also are random. Suppose the hard constraint Pl

d Plis removed from equation (2.22),

and instead we add probabilistic constraints of the form: P rhPl

d,i ≥ Pil

i

≥ 1 − i ∀i (2.31)

This new probabilistic constraint states that the demand in each node should be at least satised with equality, with probability of at least i. The problem with this

new constraint can be formulated as a chance constrained linear optimization program where the objective is to minimize the necessary overcapacity.

min

θ c T

sθ (2.32)

s.t. Cdθ dd (2.33)

P rCs,iθ ≤ ds,i ≥ 1 − i ∀i (2.34)

where Cd=           Vline_BA −Vline_BA LB −LB A −A           dd=           Pline max Pline max 0ng×1 Pmaxg 0nl×1 ∆θmax· 1np×1 ∆θmax· 1np×1          

Cs,i= −(HlLSB)i ds,i= −Pd,il

(2.35) cs= − h 01×ng ₁1×nl i LS_B (2.36)

Note that the generation capacity bounds have been removed, because we wish to determine the needed balancing generation capacities. We assume that the demands are Gaussian random variables with mean ¯Pl

d,i and variance Σld,i, although this assumption

can be relaxed. Dening:

ξi = ¯ξi+ σiu (2.37) ¯ ξi = (Pdl− HlLBθ)i (2.38) σ2 i = Σld,i (2.39) u ∼ N (0, 1) (2.40)

(23)

the probabilistic constraints in equation (2.34) can be written as P r [ξi≤ 0] ≥ 1 − i ⇔ (2.41) P r " u ≤ −ξ¯i σi # ≥ 1 − i ⇔ (2.42) Erf −ξ¯i σi ! ≤ i ⇔ (2.43) Erf    −(Pl d,i+ HlLSBθ)i q Σl d,i   ≤ i ⇔ (2.44) (HlLSBθ)i ≤ Pd,il + Erfinv(i) q Σl d,i (2.45) where Erf(t) = √1 2π Z t −∞ exp −x 2 2 ! dx (2.46) Erfinv(t) = Erf−1(t) (2.47)

Thus, the chance constrained LP in equation (2.32)-(2.34) becomes a standard LP: min θ c T sθ (2.48) s.t. Ccθ dc (2.49) where C =                     Vline_BA −Vline_BA HgLSB −HgLS_B HlLSB −(HlLSB)1 ... −(HlLSB)nl A −A                     d =                       Pline max Pline max Pmaxg 0ng×1 0nl×1 −Pl d,1− Erfinv(1) q Σl d,1 ... −Pl d,nl− Erfinv(nl) q Σl d,nl ∆θmax· 1np×1 ∆θmax· 1np×1                       (2.50) and Hg = h Ing×ng ₀ng×nl i Hl= h 0nl×ng _Inl×nl i (2.51)

(24)

Note that if i= 1₂ ∀i, i.e. we tolerate load shedding with probability 1₂, this constraint

is equivalent to the demand constraints in equation (2.35). If however i ≤ 1₂, the

chance constraint will force the available load for the demand nodes to be higher than the actual demand, to account for demand uctuations. Even though we have assumed independence between the demands, it is not a particularly severe assumption. Since the probability of load shedding should be less than i for all demand nodes i, it also

handles the case where all demands are very close to the feasibility boundary. Still, a solution to the optimal power ow problem as solved in section2.3is not guaranteed to exist for all demands within the safe bounds, since the power ows are governed by the Kircho voltage law (KVL). One would in theory have to solve a LP for every demand to check feasibility for that demand. However, if the demands are perfectly correlated, i.e.:

P_d,il = ¯P_dl+ Σl_d,iu ∀i (2.52)

u ∼ N (0, 1) (2.53)

the convexity of the feasible polyhedron implies that all demands up to the critical demand will allow for a solution of the optimal power ow problem.

Example 2.5.1. Consider the same power network as in example 2.3.1, with the only dierence that the demands are now a random variable Pl

d ∼ N ( ¯Pdl, Σld), with − ¯Pdl = [1.2, 2, 0.5, 1]T _and Σl_d=       0.12 0 0 0 0 0.23 0 0 0 0 0.08 0 0 0 0 0.15       (2.54) By applying the chance constrained optimization framework and letting = 0.05. The power demands are

− Pl =       1.346 2.203 0.6196 1.164       (2.55)

By solving (2.32), the total needed generation capacity is found to be Pg= " 2.5028 2.8296 # (2.56) The ratio between total generation capacity needed to meet unforeseen demand uctua-tions and total power demand is 1.13.

(25)

Chapter 3

Reliability of power systems

3.1 Deterministic reliability measures

Traditionally the reliability of power systems has been evaluated by deterministic ro-bustness measures. We will here describe the most widely used deterministic reliability criterion, the N − k criterion.

3.1.1 N − k criterion

A widely adopted deterministic reliability measure in power transmission systems is the N − k criterion.

Denition 3.1. Given a power transmission network and a discrete set N of com-ponents with |N | = N, the power transmission network is said to satisfy the N − k criterion if the conguration after the disconnection of any set of components Nl⊆ N

with |Nl| = l is satisable for l = 0, 1, . . . , k.

A commonly used special case of the N − k criterion is the case k = 1, in which case the power transmission network should be satisable for all congurations where a single component is disconnected. An advantage of the N − k criterion is that it guarantees robustness against any number of k contingencies. However, it does not take the probabilities of these contingencies into account, making it a rather conservative risk measure. Furthermore the complexity of computing all k contingencies is exponential in k, and quickly becomes infeasible for k > 1.

3.2 A probabilistic reliability measure

We here introduce a probabilistic risk measure of power systems based on the statis-tical properties of the state of the power system. The risk measure will be based on

(26)

the amount of compulsory load shedding necessary to bring the system to a feasible conguration. Recall the optimization problem given by (2.23). We dene the optimal amount of load load shed as a function of the matrices C and d, characteristic for a power system. S∗(C, d) = min θ n cTθ|Cθ do− 11×nl_{· P}l d (3.1)

By introducing a probability measure on the topology of the power system, i.e. on the matrices C and d, this corresponds to introducing a probability measure on the system load shed. Indeed, the topology of the power system can be formalized by viewing the matrices C and d as random variables endowed with the product probability measures µC _{and µ}d_{, dened over the Lebesgue measurable spaces Y}

C and Yd respectively.

Lemma 3. S∗_{(C, d)} _{is continuous in C and d.}

Proof. See e.g. [Bereanu, 1976] for a proof. Lemma 4. S∗_{(C, d)} _{is a random variable}

Proof. Since S∗_{(C, d)} _{is continuous in C and d, it is Borel measurable. Hence S}∗_{(C, d)}

is also a random variable endowed with a probability measure µS= µC × µd.

Thus, having shown that also S∗_{(C, d)}_{is a random variable, we can dene operations}

on the probability measure of S∗_{(C, d). Indeed, since 0 ≤ S}∗_{(C, d) ≤ −P}l

d, the mean

and the variance of S∗_{(C, d)} _{are bounded, and we can dene them as:}

Denition 3.2. Given the random variables C and d with the product measure µS =

µC_{× µ}d _{, the expected load shed ¯}_S∗ _{is dened as}

¯ S∗ =

Z

Yc×Yd

S∗(C, d) dµs (3.2)

Denition 3.3. Given the random variables C and d with the product measure µS =

µS = µC× µd, and the expected load shed ¯S∗, the variance of the load shed ¯S∗ is dened

as σ2_S∗ = Z Yc×Yd S∗(C, d) − ¯S∗2 dµs (3.3)

One can use the the rst two moments of the distribution of S∗_{(S, d)} _{as a vector}

valued risk measure. To be able to compare power systems of dierent size, one would like to dene a risk measure that does not scale with the size of the power system, like ¯S∗ _{and σ}2

S∗ generally will. This can be achieved by normalizing ¯S∗ and σ_S2∗ by

a scaling factor. Also, for comparisons between dierent distributions S∗_{(C, d), one}

would rather have a scalar risk measure than a vector valued in some cases. Motivated by this discussion, we here dene three dierent risk measures for power systems.

(27)

Denition 3.4. The absolute risk double Rabs is dened as

Rabs =

¯ S∗, σS∗

Denition 3.5. The relative risk double Rabs is dened as

Rrel= " ¯ S∗ 11×nl· Pl d , σS∗ 11×nl· Pl d #

Denition 3.6. The weighted mean and standard deviation, WMSD, is dened as ¯

S∗+ α · σS∗

where α ∈ R+ _{is a scaling factor of the standard deviation.}

We will show that ¯S∗_+α·σ

S∗is closely related to the value at risk (VaR) [Due and Pan, 1997].

Denition 3.7. VaR(S∗) is dened as

VaR(S∗) = inf{l ∈ R : Pr(S∗ > l) ≤ 1 − }

The intuitive meaning of VaR(X) is that the probability of a loss greater than

VaR(X)is less than or equal . Given the denition of VaR, we can state the following

proposition relating VaRα(X)to ¯X and σ(X).

Proposition 1.

VaRα(X) ≤ ¯X(X) +

1 √

α · σ(X) To prove the proposition, we need Chebyshev's inequality.

Lemma 5. (Chebyshev's inequality) Given any random variable X with mean ¯X and standard deviation σ

Pr(|X − ¯X| ≥ α · σ) ≤ 1 α2

The proof now follows as a corollary of Chabyshev's inequality. Proof. (of proposition 1) Note that

VaRα(X) ≤ ¯X + 1 √ ασ ⇔ Pr X < ¯X +√1 α · σ ≥ 1 − α

(28)

which is easily shown using Chebyshev's inequality Pr X < ¯X +√1 α · σ ≥ Pr X − ¯X < 1 √ α · σ = 1 − Pr X − ¯X ≥ 1 √ α · σ ≥ 1 − α

Thus the weighted sum of the mean and standard deviation of a random variable X is an upper bound of the value at risk of the random variable X.

3.3 Monte Carlo methods for probabilistic reliability

mea-sures

To be be able to calculate the expected value and the variance of S∗_{(C, d)} _{over the}

probability measure µs, we will use a Monte Carlo sampling method. Formally, we

have the following problem.

S∗(C, d) = min θ {c T_{θ|Cθ d}} _(3.4) ¯ S∗(µs) = Z Yc×Yd S∗(C, d) dµs (3.5) σ_S2∗(µ_s) = Z Yc×Yd S∗(C, d) − ¯S∗2 dµs (3.6)

In most applications, the event spaces are either continuous (in the case of continuously varying demands), or very large (for example the total number of power line failures is 2N_{, where N is the number of power lines). Since computing the expected loss ¯}_S

∗

requires integrating over a continuous or at least intractably large set of solutions of a LP, analytically computing the statistics of S∗_{(C, d)} _{is practically infeasible. Instead,}

we will use Monte Carlo methods to approximate equation (3.5)-(3.6) with the sampled mean of S∗(C, d), built from a nite number of samples Ci and di of C and d. S∗(c, D)

and ¯σS∗(µ_s) can be approximated by a nite number of samples by

¯ S∗(µs) ≈ ˆS∗,N(µs) = 1 N N X i=1 S∗(Ci, di) (3.7) σ2_S∗(µ_s) ≈ ˆσ2 S∗,N(µs) = 1 N N X i=1 S∗(Ci, di) − ˆS∗,N(µs) 2 (3.8)

(29)

By (2.18) we have the following bounds on S∗_{(C, d)} _{which will be useful for bounding}

the variance of S∗_{(C, d)}_:

0 ≤ S∗(C, d) ≤ ˆS = −X

i

P_d,il

The following lemma will be needed to prove the convergence of the sampled mean ˆ

S∗(µs) to the mean to S∗(C, d).

Lemma 6. The variance of a random variable X with compact support [a, b] is bounded by:

Var[X] ≤ (b − a)

2

4

Proof. Consider the random variable dened by Y = X − a+b

2 . We have

Var[X] = Var[Y ] = E[Y2_{] − E[Y ]}2

≤ b − a 2 2 = (b − a) 2 4

To see that the bound is actually tight, consider the probability distribution given by

pY(y) =

(

−b−a₂ with probability 1₂

b−a 2 with probability 1 2 (3.9) with E[Y ] = −1 2· b − a 2 + 1 2· b − a 2 = 0 (3.10) Var[Y ] = 1 2 · −b − a 2 − 0 2 +1 2 · b − a 2 − 0 2 = (b − a) 2 4 (3.11)

To be able to prove the main theorems of this section, we rst need some denitions. Denition 3.8. The sampled mean of the random variable X from N samples {X1, . . . , XN},

is given by ˆ XN = 1 N N X i=1 Xi

Denition 3.9. The sampled variance of the random variable X from N samples {X1, . . . , XN}, is given by ˆ σXN = 1 N N X i=1 Xi− ˆXN 2

(30)

The following lemmas state that the error of the sampled mean and variance of bounded random variables can be easily bounded.

Lemma 7. Given > 0, δ > 0, we have for a random variable X with compact support on [a, b] Prh ˆ XN − ¯X ≥ i ≤ δ for N ≥ & (b − a)2 4δ2 '

Proof. By standard probability theory we have, due to the assumption of iid random variables with compact support, we have

Var[ ˆXN] = Var   1 N N X i=1 Xi  = 1 N2Var   N X i=1 Xi  = N σ2_(X) N2 ≤ (b − a)2 4N By Chebyshev's inequality we have

Prn XˆN − ¯X ≥ o ≤ Var[ ˆXN] 2 ≤ (b − a)2 4N 2 ≤ δ

Lemma 8. Given > 0, δ > 0, we have for a random variable X with compact support on [a, b] Prh ˆσ 2 XN − σ 2 X ≥ i ≤ δ for N ≥ & (b − a)4 8δ2 '

Proof. By e.g. [Bar-Shalom et al., 2002], the variance of the sampled variance is given by

Varhσˆ_X2_Ni= 2σ

4 X

N By lemma6, the variance Var[X] = σ2

X is bounded by

σ2_X ≤ (b − a)

2

4 Thus, by Chebyshev's inequality

Prh ˆσ 2 XN− σ 2 X ≥ i ≤ Var[σ 2 XN] 2 = 2σ4 XN N 2 ≤ (b − a)4 8N 2 ≤ δ

(31)

Corollary 9. Given > 0, δ > 0, we have for a random variable X with compact support on [a, b] Pr ˆσXN − σX ≥ ≤ δ for N ≥ & (b − a)4 8δ4 '

Proof. By concavity of √·, Chebyshev's inequality and lemma8 Pr σˆ_X_N− σ_X ≥ ≤ Pr h ˆσ 2 XN − σ 2 X ≥ 2i_≤ Var[σ 2 XN] 4 = 2σ4 XN N 4 ≤ (b − a)4 8N 4 ≤ δ

By applying lemma7and corollary 9, we can state and prove the main theorem of this section.

Theorem 10. Given > 0, δ > 0, the number of samples N1 and N2 which assure

Prh Sˆ ∗,N1_(µ s) − ¯S∗(µs) ≥ i ≤ δ Pr ˆσ_S∗,N2(µs) − σS∗(µ_s)≥ ≤ δ are N1 ≥ & ˆ S2 4δ2 ' N2≥ & ˆ S4 8δ4 '

Proof. Follows by combining lemma7and corollary9wih the fact that 0 ≤ S∗_{(C, d) ≤}

−P

iPd,il .

Corollary 11. The total running time for obtaining estimates for ¯S∗(µs) and σS∗(µ_s)

with accuracy determined by and δ is polynomial in 1 ,

1

δ and the number of

transmis-sion lines np for connected power systems.

Proof. Clearly, both N1 and N2 are polynomial in both 1 and 1_δ. For each sample, a

LP is solved with the inequality constraint Cθ ≤ d. Since n ≤ np + 1 for connected

graphs, the encoding size of C (and hence d) is polynomially bounded by np. Since

polynomial time algorithms for solving LPs exist, the running time for each sampling is polynomially bounded by np.

With theorem10 we have guaranteed bounds on both the expected value and vari-ance of the load shed. The theorem can of course also be used in the reverse direction. Given numbers N1 and N2, we can obtain bounds on the numbers δ and . These

the-orems are the main theoretical foundation for the deployment of Monte Carlo methods to estimate the mean and variance of the load shed.

(32)

3.3.1 Sampling Bernoulli line failures

In particular, we consider Bernoulli failures of power lines henceforth. We model the disconnection of power line i as a binary random variable Xi ∈ {0, 1} where Xi = 0

corresponds to line i being fully functional with all parameters set to default, and Xi = 1 corresponds to line i being disconnected, i.e. the ith row of A being 0. Thus,

the failure statistics of the whole power system are given by Pr(X1 = Y1, . . . , Xnp = Ynp) ∀ Yi ∈ {0, 1}

Example 3.3.1. Consider the same power network and numerical data as in exam-ple 2.5.1, where the available power production is the one determined by the chance constrained optimization problem in example 2.5.1, i.e.

P_maxg = " 2.5028 2.8296 # (3.12) Furthermore, we introduce failure statistics on the power lines. We let E[Xi] = 0.01 ∀i.

By the Monte Carlo sampling algorithm of section 3.3, the sampled probability distri-bution of the load shed S∗ is determined. By and xing δ = 0.05, the mean and the

variance of S∗ are bounded by:

Prh Sˆ∗,N(µs) − ¯S∗(µs) ≥ 0.03 · X −Pd i ≤ 0.05 (3.13) Prh ˆσS∗,N(µs) − ¯σS∗(µs) ≥ 0.15 · X −Pd i ≤ 0.05 (3.14)

The empirical relative risk double is

Rr_rel= [0.0274, 0.0088] (3.15)

3.4 Sampling correlated Bernoulli failures

When identical software and hardware components are deployed in several system com-ponents, failures are likely to be correlated between those components. For example a software bug is likely to aect multiple components simultaneously running the same software.

To be able to handle correlations of binary distributions, we use the sampling algo-rithm for correlated Bernoulli random variables, described in [Macke et al., 2009]. The algorithm assumes the rst two moments of the Bernoulli distribution to be known, i.e. the mean and covariance matrix of the joint Bernoulli distribution. The algorithm uses a dichotomized multivariate Gaussian distribution to sample the multivariate Bernoulli

(33)

distribution. Given a Gaussian random variable U ∼ N(γ, Λ), the dichotomized ran-dom variable X is dened as

Xi= 0 if Ui≤ 0

Xi= 1 if Ui> 0

Assuming, wlog, unit variance for U, the mean (r) and variance (Σ) for X are given by:

ri= Φ(γi) (3.16)

Σii= Φ(γi)Φ(−γi) (3.17)

Σij = Ψ(γi, γj, Λij) i 6= j (3.18)

where Ψ(x, y, Λ) = Φ2(x, y, Λ) − Φ(x)Φ(y). Here Φ(·) is the cumulative distribution

function of a univariate standard Gaussian distribution, and Φ2(x, y, λ)is the bivariate

standard Gaussian distribution with correlation λ. For a given mean r and variance Σ, one must solve equations (3.16)-(3.18) The point of this algorithm is that sampling from multivariate Gaussian is computationally attractive, and can be done using standard functions in MATLAB [MATLAB, 2010]. Using this algorithm, the number of variables that need to be stored is np+n2p, compared to 2npfor storing the complete joint Bernoulli

distribution. Here np is the dimension of X. One has to choose the covariance matrix

Σ carefully. Obviously one must require Σ ≥ 0, but not all Σ ≥ 0 can be associated with a valid Bernoulli distribution. No general condition on r and Σ is known that guarantees a valid Bernoulli distribution, hence the sampling algorithm is in general not guaranteed to produce a set of correlated Bernoulli samples.

Example 3.4.1. We demonstrate the sampling algorithm with a small example. Con-sider the same power network as in example 3.3.1, with failure statistics of the power lines given by:

r =             0.05 0.05 0.05 0.05 0.05 0.05 0.05             Σ =             0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02            

The failure statistics are sampled using 5000 Monte Carlo samples. The empirical relative risk double is found to be

Rs

(34)

Note that both the mean and the variance of S∗ _{are greater than in example} _3.3.1_,

although both distributions of the line failures have the same mean.

Example 3.4.2. Consider again the same power network as in example 3.3.1, but with constant demands. In this example we will show how the reliability of the system de-creases as the correlations between the failures inde-creases. Consider correlated Bernoulli failures of power transmission lines, with mean µ = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]T_{. Let}

Σij = Σs. For Σs = 0, 0.003, . . . , 0.087, 0.09, sample 1000 line congurations and

cal-culate the empirical risk double for each value of Σs. The resulting means and standard

deviations of S∗ are shown in gure 3.1 for dierent Σs. While the expected value

of S∗ increases for small Σs, it remains nearly constant for bigger Σs. The standard

deviation however increases with increasing Σs for all Σs. Thus an increasing

correla-tion between the failures, results in greater variance of the losses, and thus more severe extreme losses and greater risk.

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 Σs µ σ

Figure 3.1: Mean ( ¯S∗_{) and standard deviation (σ}

S∗) of the load shed for correlated

(35)

Chapter 4

The Nordic Power grid

4.1 Overview of the Nordic grid

The power grid of Norway, Sweden, Finland and Denmark are highly interconnected, with signicant power exchange. For example in Sweden the electric power import accounted for more than 18 % of the total consumption in 2009. Due to the high degree of interconnection, the trading of electricity is partly handled by the common Nordpool electric market, where consumers and producers can trade power over the whole Nordic power grid. The common market and high interconnection requires a tight cooperation between the system operators, in order to ensure a well functioning electricity market. The responsibilities of the individual system operators are regulated in a system operation agreement [Nordel, 2006]. The trend over the past decades has been going towards a more unied operation of the whole Nordic power grid. The following is a citation from the system operation agreement between the Nordic transmission system operators (TSOs).

The interconnection of the individual subsystems into a common system means in-creased security and lower costs. The delivery capacity of the system as a whole is higher than the sum of the individual delivery capacities of the subsystems. As a result of the expansion of transmission capacity between the subsystems, the interconnected Nordic electric power system operates increasingly as a single entity. [Nordel, 2006] 4.1.1 Market structure in the Nordic power grid

The Nordic power market, Nord Pool, was when established in 1993 only foreseen for the Norwegian power grid [Flatabo et al., 2003]. As the trading was extended to include Sweden in 1996, the worlds rst truly international electricity market was formed. Subsequently Finland joined Nord Pool in 1998 and Denmark in 2000.

(36)

company. The trading on the Nord Pool electricity market is optional, but however in 2009 72 % of the electricity produced in the Nordic countries was traded within Nord Pool [Nordpool, 2009]. Moreover, since Nordpool is the only common market between all Nordic countries, its prices are to a high degree indicative for other trading markets and bilateral contracts. The actual trading of electric power is carried out on the Nord Pool Spot market. Here, producers and consumers bid for day-ahead hourly electricity contracts. Based on the bids from the Spot market, and without taking the actual power ow into account, the system price and electric power ows are calculated. In case of capacity limitations, the power grid is divided into predened subareas, for which the system price of electricity is adjusted in such a way that the bottlenecks vanish. These predened areas are somewhat arbitrary, and consist of individual countries and regions with historical power decit or surplus. The areas of the Nordic power grid are illustrated in gure B.2. Historically the areas have consisted of the Nordic countries, but gradually the areas are being subdivided into smaller areas. In 2010 Denmark led a complaint against Sweden for protecting its own market by not subdividing its own power system. The complaint forced Sweden to also divide its power system into four sub areas with dierent spot prices [Fredriksson and Johansson, 2010].

After the Nord Pool Spot market has closed, power trades take place on the El-bas market, allowing to adjust the Spot prices. Two hours before delivery, also the Elbas market is closed for trade, and the responsibility of balancing electric power transmissions is now transfered to the TSOs.

4.1.2 Operation of the Nordic power network

Although the Nordic power market operates as one unit, the transmission system oper-ation is for historical reasons handled individually by the Nordic countries. The TSOs are state owned but independent companies controlling the whole national power grid. However, since the Nordic countries trade on a common market and a considerable amount of electricity passes the borders, this requires a high degree of cooperation between the TSOs. Driven by the homogenization of the electricity market, the TSOs have formed a TSO association, Nordel. The operations of the TSOs are regulated by a system operation agreement between within Nordel [Nordel, 2006].

The TSOs have the responsibility of balancing power in the operational phase, when the actual power delivery is taking place. At this point, there is no option for the market to adjust power production and consumption. Still, there are unforeseen uctuations in power demand and production due to the non-deterministic nature of consumers and certain electric power producers such as wind power plants. When power production and demand deviate from the volumes determined by the Spot and Elbas markets, the TSOs are responsible of balancing the power supply. The power balancing is primarily carried out through bids on the regulating power market a

(37)

spe-cial balancing market, which represent producers willingness to increase or decrease production, and consumers willingness to cut consumption. The latter mainly applies to power intensive industry, like e.g. smelting plants. When an unbalance between production and consumption loads occurs, the TSOs match the available bids with the current power decit or surplus to restore power balance. In extreme situations, as a last resort, TSOs may stipulate the electricity generation and use load shedding to decrease demand.

4.1.3 Modernizations of the Nordic power grid

Motivated by the prospect of more ecient utilization of the current power network, eorts have been made to increase the deployment the use of modern technology in the Nordic power grid. For example, the use of demand response technology is anticipated to reduce peak energy consumption, as well as peak electricity prices [Albadi and El-Saadany, 2007]. Wide Area Measurement System (WAMS) technolo-gies increases the availability of real time accurate and time-stamped voltage, pha-sor and frequency measurements across the power system. These real time measure-ments help improve current state estimation algorithms, and improve stability con-trol of power systems [Bertsch et al., 2005]. There are eorts to increase power sys-tem capacities with the deployment of WAMS technologies [SvK, 2009]. PMUs, an integer part of WAMS, are currently being deployed in power networks around the world, most notably in the US power grid. In the Nordic power grid, the TSOs are collaborating with ABB in the deployment of WAMS technologies, including place-ment of PMUs [Leirbukt et al., 2008]. The TSO have not published the exact loca-tions of all PMUs, however most PMU localoca-tions are known and available from other sources. In eastern Denmark there are currently 4 PMUs installed, in the substa-tions of Hovegï¾1 2rd 400 kV, Hovegï¾ 1 2rd 132 kV, Asnï¾ 1 2svï¾ 1

2rket 400 kV and Radsted

132 kV [Xu et al., 2008]. In Norway there are 4 PMUs installed in the substations Hasle 400kV, Kristiansand 400 kV, Fardal 300 kV and Nedre Rï¾1

2ssï¾ 1

2ga 400 kV

[Uhlen et al., 2008]. In southern Sweden there is a PMU installed adjacent to the Konti-Skan 400 kV HVDC power line connecting eastern Denmark and southern Swe-den [Daniel Karlsson and Johannesson, 2008]. Approximate locations of PMU instal-lations in Finland, Sweden and western Denmark are provided in [Saarinen, 2010]. In gure 4.1the current 23 PMU locations mentioned above are shown.

(38)

Node PMU installation 400 kV power line 300 kV power line 220 kV power line 132 kV power line

Figure 4.1: PMU installations (black circles) in the Nordic power grid

4.2 Reliability of the Nordic power grid

Since the Nordic TSOs have agreed upon the N −1 criterion, we believe that the Nordic power grid is resilient to independent failures. Indeed, the probability of two ore more

(39)

failures occurring are second and higher order events, whose probability is small if the failure probabilities of all components are suciently small.

However we believe that the N − 1 criterion is a insucient reliability criterion when failures may be correlated. In the presence of correlated failures, the need for simulation based tools for determining system reliability increases. We will use the previously described Monte Carlo sampling technique to study the eects of correlated system failures in the Nordic grid. To perform the simulations however, we rst need an accurate reliable model of the Nordic power grid.

4.3 Model of the Nordic power grid

In order to evaluate the reliability of the Nordic power transmission system, an adequate model of the system is needed. However, the TSOs do not publish any of their internal system models, mainly due to security concerns. Due to the restrictions on the internal data of the TSOs, we will use an alternative system model. Our objective is to build an approximate model of the Nordic power transmission system, with as high accuracy as possible, using publicly available data. While there have been eorts to model other interconnected power transmission systems, such as the main European power grid [Zhou and Bialek, 2005], there are no known models of the Nordel power grid that are publicly available. The data needed to build a system model is

• Network topology

• Transmission line parameters. • Power generation data. • Power demand data.

Of this data needed, the only available directly from the TSOs is the network topology. 4.3.1 Collecting data of network topology

Since the TSOs provide no other data of their power transmission lines, the only pos-sibility was to manually enter data of node and power line locations from the maps in gure B.2and B.3. Data collected from the gures was:

• Relative map coordinates of the nodes.

• Connectivity of the nodes through power transmission lines. • Voltage of the power lines.

(40)

Afterwards, the coordinates of the Danish nodes were transformed into coordinates of the other Nordic networks, and the relative coordinates of the image were transformed into distance coordinates. The complete model has a total of 470 nodes and 717 power transmission lines.

4.3.2 Collecting power generation data

Neither the Nordel TSO association, nor the individual TSOs publish any complete data of power generation capacities in their networks, other than the total power generation capacities. The power generation in the Nordic countries is very heterogeneous in both electricity sources and geographic distribution of generation. Norway relies to 99 % on hydro power for it's electricity generation. The largest hydro power plant has a capacity of 1480 MW, but most have a generation capacity smaller than 200 MW. Most hydro power plants are located in the southern parts of Norway. Sweden has a mix of hydro power, nuclear power and other thermal power plants. Most hydro power plants are located in the middle and in the north of Sweden, whereas the four nuclear power plants are situated close to the largest metropolitan areas of Sweden. Finland also has hydro power plants in the north, and two nuclear power plants close to metropolitan areas. Finland however is relying on other thermal power plants for its electricity generation to a larger extent than Sweden, and has to import fossil fuels for it's electricity generation. Denmark has a large share of wind power for its electricity generation, but is mainly relying on thermal power plants for its electricity generation. The thermal power plants are mostly close to metropolitan areas, whereas the main wind farms are found oshore. Details of all power plants with at least 100MW capacity were found online. Since only major power plants with more than 100MW capacity are reported, minor power plants are not included in this data. It turns out that only data about smaller thermal and wind power plants was missing. We have made the assumption that this remaining thermal power generation is located in populated areas, and hence proportional to the population in each node. For the remaining wind power capacity we have assumed that the wind power generation is uniformly distributed over the land surface, and hence over the nodes.

4.3.3 Estimating power demand data

Since there is no electricity demand data available, other than cumulative data for the sub areas of the Nordic power grid, the power demands have to be estimated by other means. Following [Zhou and Bialek, 2005] we used population census data to estimate the power demand. This methodology relies on the following assumptions:

• Household power demand is proportional to the population connected to a sub-station.

(41)

• Industry power demand is also proportional to the population connected to a substation, because the workforce will settle relatively close to the industries. This may however not necessarily be the case for certain energy-intensive industries which are usually co-located with energy sources, nor for certain location-specic in-dustries such as forestry or the oil industry. Having made these assumptions, we collect population statistics from the Bureau of Statistics of the respective countries. We have collected population statistics for the major administrative regions of each country, and assumed that the demands are distributed uniformly over the demand nodes within each region. The number of administrative regions in each country used was between 12 and 21. Using smaller regions would introduce diculties in assigning the right population to each substation, since the actual connections are not known for this level. Since total power consumption data for all Nordic countries is publicly available, we have an approximation of the power demand of the nodes. Both the yearly average and yearly maximum of the daily maximum power consumption were used to estimate the power demands, yielding two dierent models for the power demands. The approximative model of the Nordic power grid is illustrated in gure 4.2, where all transmission lines are drawn, and power generation and demands are illustrated.

(42)

Production node Demand node 400 kV power line 300 kV power line 220 kV power line 132 kV power line

Figure 4.2: Approximative model of the Nordel high voltage power transmission net-work. The area of the production and demand nodes are proportional to the power production and demand respectively.

(43)

Production node Demand node 400 kV power line 300 kV power line 220 kV power line 132 kV power line

Figure 4.3: Approximative model of the Nordel high voltage power transmission net-work in Denmark.

4.3.4 Estimating line admittances

The two line parameters needed to determine the optimal power dispatch under the DC-model are the line admittances, and line capacities, none which are publicly available for the power transmission lines of the Nordic power grid. However, the admittances of power lines can be estimated by the length of the power line. Typically the reactance of high voltage power transmission lines is approximately 0.20 Ω/km [Brugg, 2006]. The lengths of a power line from a node with coordinates x to a node with coordinates y is estimated by the euclidean 2-norm as

l = dist(x, y) = kx − yk₂ =p(x1− y1)2+ (x2− y2)2 (4.1)

which is an underestimate of the actual line lengths. The total admittance Y of a power line with length l [km] is given by:

Y = 5

(44)

4.3.5 Estimating line capacities

Cross-border transmission line capacity constraints are available from Nordel [ENTSOE, 2007]. Transmission constraints are also available for the cross-sections of the sub areas within the countries The transmission capacities at 0◦ _C_{for the cross-border transmission lines}

are given in the following table.

From To Voltage Capacity

Ofoten (N) Rietsem (S) 400 kV 1350 MW N.Rï¾1 2ssï¾ 1 2ga (N) Gejmï¾ 1 2n-Ajaure (S) 220 kV 415 MW Nea (N) Jï¾1 2rpstrï¾ 1 2mmen (S) 300 kV 650 MW Hasle (N) Borgvik (S) 400 kV 2200 MW Halden (N) Skogssï¾1 2ter (S) 400 kV -Sildvik (N) Tornehamn (S) 132 kV 120 MW Eidskog (N) Charlottenberg (S) 132 kV 100 MW Lutufallet (N) Hï¾1 2ljes (S) 132 kV 40 MW Kalix (S) Ossauskoski (SF) 220 kV 1600 MW Letsi (S) Petï¾1 2jï¾ 1 2skoski (SF) 400 kV -Svartbyn (S) Keminmaa (SF) 400 kV -Teglstrupgï¾1 2rd 1 (DK) Mï¾ 1 2rarp 1 (S) 132 kV 175 MW Teglstrupgï¾1 2rd 2 (DK) Mï¾ 1 2rarp 2 (S) 132 kV 175 MW Gï¾1 2rlï¾ 1 2segï¾ 1 2rd (DK) Sï¾ 1 2derï¾ 1 2sen (S) 400 kV 800 MW Hovegï¾1 2rd (DK) Sï¾ 1 2derï¾ 1 2sen (S) 400 kV 800 MW Varangerboth (N) Ivalo (SF) 220 kV 120 MW Tjele (DK) Kristiansand (N) 250 kV DC 1040 MW Tjele (DK) Kristiansand (N) 350 kV DC -Vester Hassing (DK) Gï¾1 2teborg (S) 250 kV DC 290 MW

Vester Hassing (DK) Lindome (S) 280 kV DC 360 MW

Forsmark (S) Raumo (SF) 400 kV DC 550 MW

Table 4.1: Cross-border power transmission line capacity constraints

Some transmission capacities are given only implicitly by the capacity of a cross-section consisting of several lines. Therefore we assume that the 220 kV line from Kalix to Ossauskoski has a capacity of 200 MW, which means the 400 kV lines Letsi-Petï¾1

2jï¾ 1

2skoski and Svartbyn-Keminmaa have a cumulated capacity of 1400 MW. The

capacity of the both lines are assumed to be equal, i.e. 700 MW each. We also assume that the two power transmission lines connecting Sweden and southern Norway have the same capacity, 1100 MW each. We assume that these transmission capacities are valid for the whole cross-border links, and not just for the substations closest to the

(45)

borders. As for the power lines within the countries, the transmission capacities are only given for whole cross-sections, which makes it hard to estimate the individual line capacities.

Cross section 400 kV 300 kV 220 kV Capacity

Sweden 1-2 4 0 0 3300 MW Sweden 2-3 8 0 4 7300 MW Sweden 3-4 5 0 0 5300 MW Ofoten-Balsfjord (N) 1 0 0 500 MW Tunnsjï¾1 2dal-Rï¾ 1 2ssï¾ 1 2ga (N) 0 2 0 1000 MW NO 1-3 0 1 0 500 MW NO 5-1,2 1 3 0 1800 MW NO 1-2 2 3 0 2700 MW

Table 4.2: Transmission capacity constraints of sub area cross-sections, and number of transmission lines by capacity of each cross-section.

Instead of estimating the line capacities from the cross-section data, we estimate the line capacities using the more detailed cross-border data from gure 4.1. Under the assumption that the transmission capacity of each type of power line is the average capacity of the cross-border lines, the individual transmission capacities can be esti-mated. Here, the 220 kV line connecting Varangerboth (N) with Ivalo (SF) and the 132 kV line connecting Lutufallet (N) with Hï¾1

2ljes (S) are excluded due to their low

capacities. Voltage Capacity 400 kV 1030 MW 300 kV 650 MW 220 kV 415 MW 132 kV 143 MW

Table 4.3: Estimated transmission line capacities. 4.3.6 Evaluating the model

The linear optimal load shedding problem was applied to the above derived model of the Nordic power transmission grid. By using the YALMIP [Löfberg, 2004] interface with the GLPK [Makhorin, 2006] solver, the associated linear program was solved. With the yearly daily maximum load parameters, the total system load shed is found to be only 2 % of the total power demand, and with the yearly average of the daily

Correlated Failures of Power Systems: Analysis of the Nordic Grid