Nodal vulnerability to targeted attacks in power grids

Due to the open data policies, nowadays, some countries have their power grid data available online. This may bring a new concern to the power grid operators in terms of malicious threats. In this paper, we assess the vulnerability of power grids to targeted attacks based on network science. By employing two graph models for power grids as simple and weighted graphs, we first calculate the centrality metrics of each node in a power grid. Subsequently, we formulate different node-attack strategies based on those centrality metrics, and empirically analyse the impact of targeted attacks on the structural and the operational performance of power grids. We demonstrate our methodology in the high-voltage transmission networks of 5 European countries and in commonly used IEEE test power grids.


Introduction
The unavailability of electrical power can severely disrupt daily life and result in substantial economic and social costs (Baldick et al. 2008). This vital importance encourages a robust design and operation of power grids (Alvarado and Oren 2002). Robust power grids are able to anticipate, adapt to and/or rapidly recover from a disruptive event or a failure.
In current practice, flow-based simulations play an essential role in the security analysis of power grids. Given the generation and demand profiles, the steady-state analyses estimate the operation of power grids. Additionally, many countries require that the power grids should withstand the scheduled and unscheduled outages of its most critical lines or other components. In these contingency analyses, the component outages are also simulated to determine whether the power grids can still function properly under the failure and consequent loss of an element. However, as power grids grow in size and get more complex, the number of contingencies increases significantly, increasing the associated computational time. This motivates research towards alternative techniques Gutierrez et al. 2013;Hines et al. 2010a).
Disruptions in networks can be caused by unintentional failures or intentional attacks. Unintentional failures can include manufacturing defects, malfunction in network elements or human error. These kinds of failures can occur randomly throughout the grid and are characterized as random failures (Trajanovski et al. 2013). Intentional attacks or targeted attacks, on the other hand, are not random and are aimed at maximizing damage (Rueda et al. 2017). A major challenge in power grids is to evaluate the vulnerability of a power system to these intentional hazards, starting by quantifying the importance of electrical buses and the impact of the attacks on the network performance.
Power grids are amongst the largest and the most complex man-made systems and, like many other complex networks, it features a specific topology which characterizes its connectivity and influences the dynamics of processes executed on the network. The complex nature of power grids and its underlying structure make it possible to analyse power grids relying on network science (Strogatz 2001;Barabási and Albert 1999). The application of network science on power grids has shown the promising potential to capture the interdependencies between components and to understand the collective emergent behaviour of complex power grids (Crucitti et al. 2004;Chassin and Posse 2005).
Topological investigations of power grids have demonstrated that power grids have several components with significant importance compared to the rest of the network (Koç et al. 2014a). These components are crucial for the grid as their removal can significantly disrupt the operation of the power grids. Identifying these critical components in advance can enable power grid operators to improve system robustness by monitoring and protecting these components continuously. Koç et al. 2014a) Currently, many studies use a complex networks perspective in analysing power system vulnerabilities (Cuadra et al. 2015;Crucitti et al. 2004;Bompard et al. 2008). A significant part of these studies investigates the relationship between the topology and specific performance metrics in the underlying graph of power grids (Crucitti et al. 2004;Negeri et al. 2015;Rosas-Casals et al. 2007). Such studies focus on the basic structural properties of a graph (such as nodal degree, clustering coefficient (Hernández and Van Mieghem 2011)), which typically ignore the electrical properties, such as flow allocation according to Kirchhoff 's laws or the impedance values of transmissions lines in the grid. Mainly, two different aspects are important in the operation and consequent robustness of power grids: the topology of the network formed by electrical buses and their interconnections, and the operating conditions such as supply and demand distributions (Bompard et al. 2012;Cetinay et al. 2016). Consequently, these purely topological metrics could result in misleading research results, which may be far from real physical behaviours of power grids (Arianos et al. 2009;Hines et al. 2010b;Koç et al. 2014a).
To include the electrical properties of the grid in the analyses, several studies propose extended metrics (such as effective graph resistance (Koç et al. 2014b), the electrical centrality (Hines and Blumsack 2008) and the net-ability (Bompard et al. 2012)) by introducing a set of link weights (such as distance or resistance (Qi et al. 2015)) and node properties (such as the electrical demand and supply (Bompard et al. 2012)). Additionally, other studies have used topological and electrical metrics to rank the electrical buses and lines in power grids as a selective contingency analysis (Gorton et al. 2009;Nasiruzzaman et al. 2011).
Motivated by the increasing need of alternative studies to the flow-based analyses and the merits of network science on the investigation of power grids, in this paper, we combine both of the aforementioned approaches: First, we present two different graph representations for a power grid: a simple graph and an extended graph representation that takes the electrical properties of power grids into account. Next, we develop a methodology to identify the critical electrical buses (nodes) in power grids, and compare the impact of targeted node-attacks in detail for European high-voltage transmission networks (Cetinay 2018) and for the publicly available IEEE test power grids (IEEE 2018). Our contributions can be summarized as follows: (i) we consider two different graph models for power grids based on either purely topological information or by including the link weight information and the linearised DC power flow equations; (ii) we employ these two graph models to formulate the standard and the extended centrality metrics of nodes in power grids; (iii) we formulate 8 different attack scenarios according to these centrality metrics and empirically investigate the impact of targeted node-attacks on the structural and operational performance of power grids.

Power grids and network science
In this section, we provide details about power grids, the steady-state power flow equations and our models for power grids as simple and weighted graphs.

Power grids preliminaries
Power grids consist of nodes (electrical buses) and interconnecting links (transmission lines and transformers). The status of each node i is represented by its voltage v i = |v i |e iθ i in which |v i | is the voltage magnitude, θ i is the phase angle, and i denotes the imaginary unit. Each line l has a predetermined capacity C l that bounds its flow f l under a normal operation of the system. In the steady-state of a power grid with N nodes and L links, the injected apparent power p i + iq i at node i, where p i is the active power and q i is the reactive power, is calculated using the AC power flow equations (Grainger and Stevenson 1994): where θ ik = θ i − θ k and y (R) ik = Re(y ik ) and y (I) ik = Im(y ik ) are the real and the imaginary parts of the element y ik in the bus admittance matrix Y corresponding to the i th row and k th column, respectively.
Each node in a power grid contains a number of electrical devices and according to those, two basic types of nodes can be defined (Bergen and Vittal 1999): • Supply node: A supply node generates the active power p i and controls the voltage magnitude |v i | at its node i. • Demand node: At a demand node, it is possible to specify the extracted active p i and the reactive powers q i from the type of the electrical loads that are connected to that node. There are also nodes without a supply or a demand connected, which can be modelled as a demand node with no injected power, i.e., p i = 0 and q i = 0.
Due to the impedance of transmission elements, there are power losses during the operation in power grids. As the losses are dependent on the system state-the supply and demand dispatches-they cannot be calculated in advance. Therefore, a slack node among the supply nodes is assigned in power grids to compensate for the difference between the total supply and the total demand plus the losses.

DC power flow equations
The AC power flow Eqs. (1) and (2) are non-linear and the solution process is generally iterative. A linear set of equations is more desirable whenever fast and repetitive solutions are needed. Linearisation can be reasonably accurate when the following conditions are met Van Hertem et al. 2006): 1 The difference between the phase angles of neighbouring nodes is small such that sin θ ik ≈ θ ik and cos θ ik ≈ 1. 2 The active power losses are negligible, and therefore, the bus admittance matrix can be approximated as Y ≈ iY (I) where Y (I) is the imaginary part of the admittance matrix Y, calculated neglecting the line resistances. 3 The variations in the voltage magnitudes |v i | are small and, can be assumed as |v i | = 1 for all nodes.
If these conditions are approximately met, the AC power flow equations can be simplified to the so-called the DC power flow equations: Although the DC power flow solution is less accurate than the AC power flow solution, in practice, the differences in high-voltage transmission networks between the phase angles of neighbouring buses and the variations in voltage magnitudes are relatively small, thus the error is assumed to be negligible (Van Hertem et al. 2006).

Graph representations of power grids
This section presents our models for power grids as simple and weighted graphs.

Power grid as a simple graph
A simple graph is an unweighted, undirected graph containing no self-loops or multiple links. A power grid can be modelled as a graph G(N , L) where N denotes the set of N nodes and L denotes the set of L links in which multiple lines connecting the same pair of nodes are modelled as one link. The N × N adjacency matrix A specifies the interconnection pattern of the graph G(N , L): a ik = 1 only if the pair of nodes i and k are connected by a direct link; otherwise a ik = 0. The N × N Laplacian matrix Q is defined as

Power grid as a weighted graph
Alternatively, a power grid can be modelled as a weighted graph where each link is assigned a weight that is related to the admittance of the transmission line it represents.
We model a power grid as a weighted graph G(N , L) where N denotes the set of N nodes and L denotes the set of L links 1 . By writing the DC power flow equations in (3) in terms of the adjacency matrix A of G(N , L) we introduce the weighted adjacency matrixÃ, where each nonzero elementã ij = a ij y represents both the connectivity and the admittance between nodes i and j. Equation (4) can then be written as: Since (5) holds for every node i, the corresponding matrix representation is where P = [ p 1 . . . p N ] T is the vector of net active power injection at the nodes with a balanced supply and demand, i.e. u T P = 0 where u is all-one vector,˜ is the weighted diagonal degree matrix, and = [θ 1 . . . θ N ] T is the vector of phase angles at the nodes. Finally, introducing the weighted LaplacianQ =˜ −Ã into (6) yields P =Q (7) where the weighted LaplacianQ is a symmetric, positive semi-definite matrix that possesses non-negative eigenvalues apart from the smallest eigenvalue, which is zero (Van Mieghem 2010). The inversion of the active power -phase angle relation P =Q in (7) is not possible due to the fact that detQ = 0, which follows from the characteristic propertỹ Qu = 0 of the weighted Laplacian. Although the inverse of the weighted Laplacian matrix does not exist, the active power -phase angle relation inversion can be shown to be = Q † P + u T N u, where Q † is the pseudo-inverse of the weighted LaplacianQ, obeying QQ † = Q †Q = I − 1 N J with the identity I and all-one matrix J = uu T . By choosing the average phase angle in the graph θ av = u T N = 0 as the reference (Cetinay et al. 2016), the phase angle -active power relation takes the elegant form of While the weighted LaplacianQ and its pseudo-inverse Q † are derived here based on the linearised DC power flow equations in power grids, their applicability is far wider . A weighted LaplacianQ can describe many processes, that are linear in or proportional to the network topology such as electrical circuits, water flow networks, mechanical or thermal systems. The process equivalence between those systems are given in Table 1.

Targeted attacks on power grids
The threats for power grids can be classified by using multiple criteria considering the causes of the threat, their consequences or the preventive actions to manage the hazards (Ciapessoni et al. 2017). One example of such threats are targeted attacks on power grids, which involve intentional, criminal actions to destroy the network. In modelling these threats, we assume that the attacks are performed with the knowledge of power grid layout and with the intention to maximally disrupt the network performance while attacking as few nodes (electrical buses) as possible. Throughout this section, we describe how network science can be employed to formulate such attack strategies, where target nodes correspond to most critical or most vulnerable nodes whose removal significantly disrupts the network functioning. We first describe the standard centrality metrics, which are purely based on the underlying topology of power grids, and then we extend these metrics to include the information on the link weights, i.e. the admittances of the transmission lines, and the DC power flow equations in power grids.

Ranking nodes in the simple graph representation of a power grid
In this section, we review some of the existing topological centrality metrics in order to rank the importance and the centrality of nodes in the underlying simple graph of power grids.

Degree centrality
The degree d i of a node i in the graph G(N , L) is equal to the number of its neighbouring nodes (Freeman 1978). The degree d i can be calculated using the adjacency matrix A:

Eigenvector centrality
The eigenvector centrality of a node is a global centrality metric that depends not only on the number of its neighbouring nodes, but also on the number of 2-hop neighbouring nodes, 3-hop neighbouring nodes, and so on (Bonacich 1987;1991). The eigenvector centrality x i of node i is equal to the i th component of the eigenvector corresponding to the largest eigenvalue λ 1 of the adjacency matrix A. The principal eigenvector centralities thus follow from the linear equations:

Betweenness centrality
Another metric to assess node importance or centrality is the betweenness centrality (Freeman 1977). In calculating the betweenness centrality, it is assumed that information or services are transmitted over shortest paths between node pairs. Hence, if many shortest paths pass through a certain node, this node takes a central role in the network.
If |P s→t | is the number of all possible shortest paths from node s to node t, and |P s→t (i)| is the number of those paths that pass through node i, then the betweenness b i of node i is equal to In other words, the betweenness centrality of a node i shows the fraction of all shortest paths between any pair (s, t) of nodes, that pass through node i.

Closeness centrality
In calculating the closeness centrality, the hopcount H(P i→j ) that is the number of links in the shortest path P i→j between a pair of nodes i and j, is used. The closeness centrality c i of a node i is defined as (Freeman 1977): which is the reciprocal of the sum of the hopcounts of node i to all other nodes. A large closeness centrality value thus corresponds to a "central" node that is well-connected by a few hops to other nodes.

Ranking nodes in the weighted graph representation of a power grid
While the standard centrality metrics are based on purely topological information, it is possible to extend the definition of these metrics by including the link weight information and the power flow equations in power grids. Different definitions of extended centrality metrics (extended betweenness (Bompard et al. 2010), modified betweenness and closeness centrality (Gutierrez et al. 2013), electical degree (Hines and Blumsack 2008)) exist 2 and are evaluated by simulations via power flow solvers or by calculating power transfer distribution factors (PTDF) in power grids. Such simulation-based definitions are generally computationally expensive and formulations with the absence of slack node(s) may not fully explain the analogy between the extended centrality definitions and the weighted graph model for power grids. Extended metrics were also defined before (Ellens et al. 2011;Newman 2005) based on the voltage -current relation in electrical circuits. Since the phase angle -active power relations in (7) and (8) in power grids obey the same linear relation as those in electrical circuits (as described before in Table 1), these metrics can identify central nodes in power grids.
We take here a graph theoretical approach using the slack-node-independent weighted graph representation for power grids described in the previous section. This weighted graph model facilitates both the analogy between the standard and the extended centrality metrics, and the enhanced linear algebra to formulate the closed-form expression of centrality metrics via graph-related matrices.

Weighted degree centrality
Similar to the topological definition in (9), the weighted degree centrality is related to the number of neighbours of a node. However, rather than only considering the number of neighbours, the weighted degreed i also includes the information of the admittancesã ij of the transmission lines that link the nodes, which leads to the definition: A large value of the weighted degreed i corresponds to larger values of the admittance directly connected to that node, which indicate that node i is well connected to its neighbours.

Weighted eigenvector centrality
In analogy with the eigenvector centrality in (10), the weighted eigenvector centralitỹ x i not only captures the total admittance of all lines connected to node i, but is also influenced by the admittance of all lines connected to its neighbours, their neighbours and so on. The weighted eigenvector centralities correspond to the eigenvector of the highest eigenvalueλ 1 of the weighted adjacency matrixÃ. Thus, the principal weighted eigenvector centralityx i is given by the equation:

Flow betweenness centrality
While in the standard definition of the betweenness and the closeness centrality in (11) and (12), information exchange and other processes are assumed to travel over shortest paths, in the case of the DC power flow equations (or in the equivalent linear systems in Table 1), the flow distribution obeys Kirchhoff 's and Ohm's laws. Therefore, the standard betweenness and closeness centrality based on shortest paths may not fully capture the operation of power grids. Instead, the flow betweenness centralityb i of node i depends on the total flow running through that node, as proposed by Newman (2005): where B(i) denotes the direct neighbours of node i, and |f s→t (i, j)| is the magnitude of the power flow through the link between i and j when a unit active power is injected at node s and extracted from node t. In Appendix 2, we show how these flows can be calculated from the weighted graph representation of a power grid. Higher values of the flow betweenness centralityb i indicate the importance of a node with respect to the electrical power transmission in power grids.

Electrical closeness centrality
Similar to the definition of the closeness centrality in (12), the electrical closeness centrality of a node is an indicator of the average distance of that node to all other nodes. However, since the flow in a power grids obeys Kirchhoff 's laws, the effective resistance (Ellens et al. 2011;Cetinay et al. 2016) is a more appropriate distance metric between nodes than the hopcount of shortest-path. The effective resistance ij between a pair of nodes can be calculated from the pseudo-inverse Laplacian matrix as (Van Mieghem 2010): and captures the effect of the active power transfer p ij and the phase angle difference θ i −θ j between a pair of nodes, when active power is only injected at and extracted from nodes i or j: Since the effective resistance satisfies the properties 3 of a distance function (Klein and Randić 1993) and obeys the flow equations in power grids, it can be used to define a distance-based centrality metric. The electrical closeness centralityc i of a node equals the reciprocal of the total effective resistance of that node to all other nodes 4 : Compared to the shortest-path hopcount H(P i→j ), the effective resistance ij does not depend only on the shortest path, but also incorporates the information of all possible paths between node i and j, where the contribution of each possible path follows from the linear flow equations. In the case of the unweighted tree networks, the effective resistance ij equals the hopcount H(P i→j ) for all nodes. Thus, for tree-like power grids with equal admittances, the electrical closeness centrality closely resembles the topological closeness centrality, while for power grids with many loops (i.e. non-tree-like) both metrics could differ significantly.
Each of the centrality metrics we present captures a certain aspect of the structural and the operational centrality in the network, such as the strength of a direct connectivity (degree and eigenvector centrality), being a part of many important paths (betweenness centrality) or being close to other nodes (closeness centrality). In recent years, another conceptual definition of centrality has emerged. Based on optimal percolation theory (Morone and Makse 2015), which considers the problem of "finding the smallest set of nodes whose removal fragments the network in small disconnected pieces", a number of new metrics have been proposed (such as the collective influence (CI) (Morone and Makse 2015; Morone et al. 2016), belief propagation decimation (BPD) (Mugisha and Zhou 2016) and CoreHD (Zdeborová et al. 2016)). Such metrics reflect the importance of a node for the global structural coherence as well as their influence in spreading behaviour. However, to the best of our knowledge, extended metrics based on percolation theory have not been studied yet. Therefore, in this work, we focus on the generally accepted and adopted centrality metrics to the power grids.

Identifying the effect of node removals in power grids
In this section, we empirically compare the effects of the targeted node removals based on the centrality metrics presented in the previous section. To evaluate the change in the network functioning, we use two performance metrics that can quantify both the topological and the operational characteristics of the grid after targeted attacks. We consider the networks from 5 real-world power grids of European countries (Cetinay 2018) and 5 synthetic power grids from the IEEE test case database (IEEE 2018).

Performance metrics
In an ideal power grid which is robust to targeted attacks, the removal of nodes should not significantly alter the network functioning. In some cases, removing a node from the power grid can partition the network into several components, which are disconnected from each other. This is undesirable as this partition both adversely affects (i) the structure: as the size (i.e. the number of nodes) of the connected component of the network is decreased, and (ii) the operation: since the disconnected structure disrupts the service and the capacity of the network. In this work, we present two performance metrics in our case studies, the size and the capacity of the giant component, to assess both the topological and the operational performance aspect in the network.

The size of the giant component
The giant component (Molloy and Reed 1998) is the connected component of a graph that contains the largest fraction of the entire graph's nodes. The size of the giant component in the graph reflects the disruptive effect of node removals on the structure of the network.
We assume that the underlying graph of the initial network is connected, thus the initial size of the giant component is N. Then, we calculate the normalized size σ of the giant component after each node removal as the ratio between the size of the current giant component and the initial network size N, in other words where 1

The capacity of the giant component
Each transmission line in a power grid is associated with a maximum flow carrying capability. For the safe operation of a network, the flows through the network links should be below these capability. If the flow limits are exceeded, the situation is detected by protection relays, the circuit breakers are tripped, and the corresponding element is taken out of service. The possibility and the negative impact of cascading failures in power grids increases when the operating point of a power grid is close to the flow carrying capabilities of its links (Koç et al. 2013;Cetinay et al. 2017). Consequently, a network with a high flow carrying capability is desired. We calculate the total capacity of the network as the sum of the maximum flow carrying capabilities of links in the largest connected component of the graph. When multiple lines are connecting the same pair of nodes, we consider an equivalent capacity between those nodes. This equivalent capacity represents the maximum power that can be transferred between these nodes such that the resulting power flow through each single line is at most at its capacity. In Appendix 3, we describe how this equivalent capacity is calculated.
The capacity of the giant component depends on the number of links in the giant component as well as the flow carrying capability of links, which are closely related to the electrical demands and supplies at the neighbouring nodes 5 . We calculate the normalized capacity of the giant component γ after each node removal as the ratio between the total capacity of the current giant component and the total capacity of the initial network, in other words (18)

Properties of the networks used in simulations
We considered the high-voltage transmission networks of 5 European countries in our case studies: Austrian, Belgian, Dutch, French and German power grids. In addition, we included 5 widely used test power grids from IEEE database (IEEE 2018). In all networks, multiple lines connecting the same pair of nodes are represented as an equivalent single link using the equivalent admittance (23) and the equivalent maximum flow carrying capability (25). The degree distributions of the underlying graphs are shown Figs. 1 and 2.
Additionally, more details of the power grids in our case study are available on our GitHub page (Cetinay 2018).

The effects of targeted node removals in power grids
We apply both the standard and the extended centrality metrics as node-attack strategies in power grids. For each centrality metric, we start the attacks by removing the node (and all its links) with the highest value of the chosen centrality metric. After each node removal, we recalculate the values of the centrality metric, and continue by removing the node with the highest value of the centrality metric in the current giant component of the graph. Note that, during the successive node removals, we do not take the cascading dynamics (such as overloading of links or demand/supply redistribution due to cascading failures ) into account. In other words, we focus on the instant just after the removal of nodes to identify the effects on the structure and the operational performance indicators of the power grid. Figures 3 and 4 show the changes in the normalized size and the capacity of the giant component when we sequentially remove the nodes according to 8 different centrality metrics. We observe that the betweenness and the flow betweenness centrality are the best attack strategies as they can maximally disrupt the network functioning with fewer attacked nodes. On the other hand, the degree centrality may not always successfully assign an important node. Compared to the degree centrality, the betweenness and the flow betweenness centrality give more fine-grained centrality values for each node, whereas, multiple nodes with the same degree exist, as illustrated in Figs. 1 and 2, making them indistinguishable. In addition, we observe that the eigenvector and the weighted eigenvector centrality are the least effective attack strategies: Targeted attacks according to these centrality metrics destroy the network slower than other attack strategies. Fig. 1 The degree distribution of the simple graphs of 5 European power grids Fig. 2 The degree distribution of the simple graphs of 5 IEEE power grids Next, in order to compare the attack strategies and to further quantify the topological and operational changes in power grids, we calculate the average of the structural and the operational performance indicators (or the energy ε values of a graph (Trajanovski et al. 2013)) that are the normalized sums of the size and the capacity of the giant component over successive targeted attacks, respectively. Thus, the average valueσ of the structural performance indicator of the power grid over K successive node-attacks can be calculated asσ where σ (k) is the normalized size of the giant component after k successive attacks. Similarly, the average valueγ of the structural performance indicator of the power grid over K successive node-attacks is  (19) and (20) are evaluated on a score between 0 and 1: In an ideal power grid which is robust to targeted attacks, the node removals should have slight effects on the network performance. Thus, a performance indicator close to 1 is desirable by the network operators. On the other hand, a lower performance indicator over successive node-attacks indicates a powerful (destructive) attack-strategy in which few important nodes of the network are identified and removed, with negative operational and structural consequences. In Figs. 5 and 6, we present the average values of the performance indicators in European and IEEE test power grids (PG) after different attack strategies that remove 10% of the initial network nodes, respectively. Higher values in Figs. 5 and 6 represent higher robustness to the targeted attacks, whereas lower values indicate vulnerability or a severe disrupt in network functioning. We observe that targeted attacks based on the flow betweenness and the betweenness centrality followed by the closeness and the electrical closeness centrality are the best attack strategies to decrease the structural and the operational performances of the power grids. As an example, the targeted node attacks

Main lessons learned from the analyses
In this section, we summarize the insights obtained in the previous sections. The main lessons learned from the analyses of the targeted node-attacks based on the different centrality metrics from the simple and the weighted graph representations of power grids in the tested networks are as follows: • The degree centrality (9) only provides information on the local structure around a node. Similarly, the weighted degree centrality (13) reflects local connectivity information. Thus, a node that is connected to many other nodes (with high admittance) is not necessarily a central node for the whole network. Therefore, as illustrated by the targeted attack simulations, the degree and the weighted degree centralities cannot always indicate the important nodes.
• The betweenness centrality (11) incorporates information about the global network structure, and in the analyses of the test networks, high betweenness centrality values were found to efficiently indicate the nodes whose removal would significantly disrupt the network performance. While successfully indicating vulnerable nodes, the betweenness centrality (11) is based on the shortest paths only. This means that the betweenness centrality does not discriminate nodes that are positioned "close" to many shortest paths (and would be considered central), and peripheral nodes. This limitation is partly addressed by the flow betweenness centrality (15), in which the flows through the network links are distributed throughout the network according to the Kirchhoff's laws. In the analyses of the test networks, removing nodes with a high flow betweenness usually resulted in the most destructive effects on the network.
• The closeness centrality (12) reflects the average shortest path distance from a node to all other nodes in the network. Higher closeness centrality values thus indicate nodes which can easily reach the other nodes in the network. Similarly, higher values of electrical closeness centrality (16) show a node that is on average close to the other nodes in the network, based on the operationally inspired effective resistance distance instead of the shortest-path distance. In the analyses of targeted attacks, the performance of the closeness and the electrical closeness centrality in identifying the important nodes in the tested power grids are found to be similar.
• The eigenvector centrality (10) can rarely identify the critical nodes, and thus, the targeted attacks based on the eigenvector centrality are generally the worst destructive strategy among the traditional centrality metrics in the tested networks. Similarly, the weighted eigenvector centrality (14) seems not to successfully indicate important nodes.
The analyses of the targeted node-attacks show that centrality metrics, in particular the (flow) betweenness and (electrical) closeness, are very successful in indicating the critical nodes whose removals sharply decrease the selected performance indicators (the size and the capacity of the giant component) of power grids. Identifying these critical components in advance can enable power grid operators to improve system robustness by monitoring and protecting these components continuously. Additionally, although the effect of targeted attacks are more significant when the centrality information is updated after each node removal, the information based on the initial calculation of the centrality metrics is also fairly successful in finding the important nodes. The degree centrality is a good indicator to fragment the network to decrease the structural and operational performance indicators of power grids (See Appendix 1).

Conclusion
In this paper, we took a network science approach to investigate the vulnerability of power grids to malicious targeted attacks. First, we presented two different graph models for power grids: simple and weighted graphs. Subsequently, using these graph models, we ranked the importance of each node according to the standard and the extended centrality metrics that take into account the electrical properties of the grids such as the admittance of the transmission lines and the flow allocation according to the DC power flow equations. Via case studies in both real-world and test power grids, we show that the power grids are highly vulnerable to targeted attacks: sequentially removing the nodes with the highest centrality is a good strategy to fragment the power grids, and to maximally decrease its operational performance. In almost all power grids in our case study, removing approximately 15% of the nodes according to the flow betweenness centrality destroys the network almost completely. Grid operators can use the proposed methodology to analyse the current vulnerability of their network to targeted attacks and to take necessary measures by protecting the important nodes in their networks.  Knowing the phase angle at each node, it is then possible to calculate the flow f s→t (i, j) through the link between nodes i and j as The flow betweenness centralityb i of a node i is the sum of the absolute flows that pass through that node i, over all possible pairs of source and target nodes 6 : b i = s,t∈N \{i} j∈B(i) ã ij (e i − e j ) T Q † (e s − e t ) .
Introducing the phase angle difference from Eq. (24) then leads to f l = f ij y (I) l y (I) ij for the flow f l through line l. Since the maximum flow through each line is constrained by its flow capacity: f l ≤ C l , we find that the total flow f ij between node i and j is constrained by an equivalent capacity C ij equal to:

Funding
This work was supported in part by Alliander N.V.

Availability of data and materials
The data set used in this article is available in the cited references (Cetinay 2018;IEEE 2018).