Topological Modeling Research on the Functional Vulnerability of Power Grid under Extreme Weather

: The large-scale interconnection of the power grid has brought great beneﬁts to social development, but simultaneously, the frequency of large-scale fault accidents caused by extreme weather is also rocketing. The power grid is regarded as a representative complex network in this paper to analyze its functional vulnerability. First, the actual power grid topology is modeled on the basis of the complex network theory, which is transformed into a directed-weighted topology model after introducing the node voltage together with line reactance. Then, the algorithm of weighted reactance betweenness is proposed by analyzing the characteristic parameters of the power grid topology model. The product of unit reliability and topology model’s characteristic parameters under extreme weather is used as the index to measure the functional vulnerability of the power grid, which considers the extreme weather of freezing and gale and quantiﬁes the functional vulnerability of lines under wind load, ice load, and their synergistic effects. Finally, a simulation using the IEEE-30 node system is implemented. The result shows that the proposed method can effectively measure the short-term vulnerability of power grid units under extreme weather. Meanwhile, the example analysis veriﬁes the different effects of normal and extreme weather on the power grid and identiﬁes the nodes and lines with high vulnerability under extreme weather, which provides theoretical support for preventing and reducing the impact of extreme weather on the power grid.


Introduction
The power grid is the foundation of the national power supply, and its safe and reliable operation exerts a profound impact on national security, stability, and economic development [1]. In recent years, the rapid development of China's economy and society has been witnessed, along with a decline in regional isolated grid operation. The completion and operation of a 500 kV substation in Jilong county indicate that the main power grid in Tibet has achieved complete regional coverage, further indicating the success of the national interconnection of the power grid in mainland China [2].
The trend for power system development worldwide lies in the large-scale interconnection of the power grid, which has brought great benefits to the development of the national economy as well as other fields. Power grid accidents are characterized by their rapid spread, wide influence, and serious accident consequences, and thus large-scale interconnection has left demanding requirements for the safe and reliable operation of the power grid. Power grid accidents have been occurring frequently recently, inflicting enormous losses on physical and economic security and seriously hindering social harmony, stability, and long-term development. The complex and changing climate, together with seasonal weather in China, increases the frequency of grid failures caused by extreme weather [3,4]. Table 1 lists the typical large-scale grid failure accidents caused by extreme conditions worldwide [5][6][7][8][9][10][11][12][13].
Most of the extreme conditions affecting the power grid are extreme weather. Therefore, extreme weather will accelerate the failure rate of power grid units and lead to large-scale power grid failures; it is thus necessary to carry out research on improving power grid reliability and reducing the occurrence of large-scale power grid faults under extreme weather. As an important indicator evaluating whether the power grid system is strong, power grid vulnerability is appealing in the field of power grid reliability. The power grid has the characteristics of a small-world network, and complex network modeling has been extensively applied into power grid vulnerability modeling and analysis, providing a developed method for power grid vulnerability assessment [14][15][16]. At present, the research on the power grid mainly focuses on reliability, vulnerability, and resilience; reliability and vulnerability measure the stability of the power system, and resilience considers its ability to maintain its function after damage [17]. Power grid vulnerability can be divided into functional vulnerability and structure vulnerability. Functional vulnerability is associated with power grid units, and structural vulnerability is associated with the whole power grid system [18]. Taking into account the functional features critical to largescale blackouts and exploring the failure law of power grid units, power grid functional vulnerability has attracted great research interest recently [19]. Therefore, the study of power grid functional vulnerability under extreme weather can provide theoretical support for the solution of large-scale power grid faults, yet recent research seldom considers the impact of extreme weather synergism.
Consequently, considering the deficiencies of current research on power grid functional vulnerability, this paper establishes a topology model of the power grid based on complex network modeling to study the functional vulnerability of power grid units, and an analysis method of power grid functional vulnerability under extreme conditions, including gale and freezing weather, is proposed to measure the short-term vulnerability of power grid units. The major contributions of this paper are the power network topology modeling based on complex network modeling and power system functional vulnerability analysis based on a reliability and example analysis. The key points of the proposed method are as follows.
(1) First, the topology of the power grid is modeled, analyzing the characteristics of the power grid from a macro perspective, in which the internal structure of the power grid is ignored. The grid model is transformed into a directed weighted topology model through grid nodes classification and line reactance weighting, of which the characteristic parameters are analyzed after improving the betweenness algorithm. (2) Second, extreme weather (gale weather and freezing weather) is introduced into the power grid functional vulnerability, and the power grid functional vulnerability model based on reliability under extreme weather is proposed. (3) Last, on the foundation of the functional vulnerability model under extreme weather based on reliability, the topology modeling and functional vulnerability of the power grid system are analyzed through the IEEE-30 (Institute of Electrical and Electronics Engineers) node system. The simulation result, outputting the nodes and edges with high vulnerability under extreme weather, verifies the rationality of the proposed power grid functional vulnerability model.
The results of the final example analysis indicate that our proposed model can serve as theoretical support for the research field of power grid large-scale faults under extreme conditions, improving the adaptability of traditional models to different situations so as to fit the real situation. The framework of the proposed method is shown in Figure 1. model through grid nodes classification and line reactance weighting, of which the characteristic parameters are analyzed after improving the betweenness algorithm.
(2) Second, extreme weather (gale weather and freezing weather) is introduced into the power grid functional vulnerability, and the power grid functional vulnerability model based on reliability under extreme weather is proposed. (3) Last, on the foundation of the functional vulnerability model under extreme weather based on reliability, the topology modeling and functional vulnerability of the power grid system are analyzed through the IEEE-30 (Institute of Electrical and Electronics Engineers) node system. The simulation result, outputting the nodes and edges with high vulnerability under extreme weather, verifies the rationality of the proposed power grid functional vulnerability model.
The results of the final example analysis indicate that our proposed model can serve as theoretical support for the research field of power grid large-scale faults under extreme conditions, improving the adaptability of traditional models to different situations so as to fit the real situation. The framework of the proposed method is shown in Figure 1.
The rest of this paper is constructed as follows. Section 2 presents the research status related to extreme weather and power grid functional vulnerability. Sections 3 and 4 introduce the power grid topology modeling method based on complex network modeling and power grid functional vulnerability based on reliability, respectively. Section 5 performs an example analysis using the IEEE-30 node system. The final analysis and conclusion are presented in Section 6. The rest of this paper is constructed as follows. Section 2 presents the research status related to extreme weather and power grid functional vulnerability. Sections 3 and 4 introduce the power grid topology modeling method based on complex network modeling and power grid functional vulnerability based on reliability, respectively. Section 5 performs an example analysis using the IEEE-30 node system. The final analysis and conclusion are presented in Section 6.

Related Work
There is extensive relevant research on the functional vulnerability of power grid units and relevant theoretical models under extreme weather. Extreme weather is the main influence on power grid failure. The general modeling methods of power grid reliability, vulnerability, and resilience under extreme weather are two-state models and three-state models according to the severity of extreme weather, performing a comparative analysis of reliability under different weather conditions [20].
The influences of extreme weather on the reliability and vulnerability of the power grid are mainly the destruction and increased failure rate of power grid units. Zhou T. et al. [21] utilized the detrended fluctuation method to analyze the impact of extreme weather conditions on grid unit failures. Billiton R. et al. [22] researched the normal case and the extreme weather of a gale disaster simultaneously, obtaining the grid reliability models under the two states through comparative analysis. R. Billinton et al. [23] established a three-state climate model to study the reliability of the power grid. Lightning weather [24,25], gale weather [26,27], and line icing [28][29][30] were separately researched to search the impact on power grid reliability. Yong Liu et al. [31] integrated the extreme conditions of gale disasters and lightning, studying the reliability models under different severities. Bhuiyan M. R. et al. [32] integrated the three-state climate model into the establishment of the failure rate model, in which the weather superposition state of transmission lines was considered at the same time, and the Monte Carlo method was utilized to sample and simulate the climate region and unit state. Kiel E. S. et al. [33] connected severe weather with power grid vulnerability, exploring the influence of severe weather on the power system reliability index through the failure rate model.
Power grid resilience involves the four stages of anticipation, absorption, recovery, and adaptability, and each stage is easily influenced by extreme weather. The energy management system, optimal power flow, and dynamic control are widely applied to the resilience modeling of the power grid under extreme weather [34]. In terms of frequently concerning typhoon weather, Nguyen H. T. et al. [35] utilized a stochastic optimization model based on resilience improvement, and Zhou J. et al. [36] proposed a coordinated optimization method based on the regional integrated energy system to realize resilience improvement. The mathematical programming model is an important optimization method to study the resilience of power grid; Hou H. et al. [37] used machine learning methods to analyze the impact of wind disasters on transmission lines, and Karangelos E. et al. [38] combined the failure rate of grid units sensitive to weather with an optimization model based on vulnerability identification through the Monte Carlo method. Under extreme weather, the power grid capabilities related to resilience include physical hardiness and operational capability, and each capability enhancement strategy improves the overall power grid system resilience [39]. With the rapid development of smart grid technology, the power grid resilience evaluation method based on modularity can effectively reduce the failure rate of power grid units under extreme weather [40].
In the above works researching the impact of extreme weather on the power grid, the extreme weather considered is limited, mainly focusing on lightning, icing, and gale disasters, and it lacks consideration of extreme conditions with the synergism of multi-disaster competition. In particular, most case studies of power grid resilience have been carried out under stormy weather. In the research on power grid reliability and vulnerability, although the failure rate of power grid units in different extreme weather conditions is taken into account, there are few specific physical models carried out in quantitative research.
Various research methods describing the different characteristics of power grid units have attracted more and more attention in terms of power grid functional vulnerability. Functional vulnerability analysis methods can be divided into complex network theory and power flow analysis. Nima Amja et al. [41] studied node voltage vulnerability based on energy function. Ubeda J. R. et al. [42] studied power grid functional vulnerability using an analytical method, yet this method failed to work when there were many network nodes with high complexity; this can be solved through the Monte Carlo method when the complexity of the power grid is high. Singh C. et al. [43] modeled reliability based on the uncertainty probability and guarantee time of hidden failure. Yu Xingbin et al. [44] defined grid vulnerability parameters, including the bus isolation rate, load loss rate, and load reduction expectation, and proposed a functional vulnerability evaluation index of the system level based on these parameters. By analyzing the cascading failure model of the power grid based on the Monte Carlo method, Su Sheng et al. [45] identified the key vulnerable nodes of the power grid system and optimized the distribution of system compensation devices. However, although the Monte Carlo method can deal with the high complexity of the power grid, it also brings a huge computational burden. Nan L. et al. [46] proposed a vulnerability analysis method considering both topological and functional characteristics of the power grid to identify vulnerable units. Liu Y. et al. [47] Energies 2021, 14, 5183 5 of 25 evaluated the dependence between grid units using a dynamic network flow redistribution model. Functional vulnerability deeply studies the dynamic process of cascading faults; Abedi A. et al. [48] established a power flow model by introducing a novel vulnerability assessment index to analyze cascading faults, and Fang J. et al. [49] designed a multiobjective optimization algorithm to evaluate power grid vulnerability with regard to cascading faults. Among the methods above, the functional vulnerability of the power grid is rarely analyzed from the physical and mechanical points of view.
The existing research rarely combines the power grid functional vulnerability with extreme weather, which makes the modeling deviate from reality and unable to effectively study the power grid cascading fault mechanism. Table 2 presents some details in the investigated literature on extreme weather and functional vulnerability. Identifying how to carry out topology modeling that fits the reality is the major challenge when analyzing vulnerability; this means it must fully consider related physical and mechanical characteristics. Therefore, extreme conditions with the synergism of multi-disaster competition need to be considered to analyze power grid functional vulnerability. Complex network modeling is an efficient modeling method for the power grid that does not ignore significant characteristic information of the power grid during simulation [50]. In conclusion, in this paper, we used the complex network topology theory to model the power grid according to the deficiencies of existing research. The extreme weather of gale and freezing was considered, and the feasibility of the proposed topology model was verified by the quantitative analysis of an example.

Topology Modeling Rules
In order to study the vulnerability of the power grid from a macro perspective, the first task is to transform the power grid structure into a topology model based on graph theory.
According to the structure and characteristics of different types of power grid units, power grid units are divided into two types represented by nodes and edges following the common modeling method, forming a powerless undirected topology model. First, the generator, transformer, and user load are taken as network nodes with a total number n while modeling, which are divided into three categories according to the energy flow direction, including the source node set S G , the sink node set S L , and the link node set S t . Then, the high voltage line L G and the transformer branch L T are regarded as the edges with a total number k. The nodes connecting each edge are denoted by i and j, and the n-order incidence matrix A is defined to represent the connection relationship between nodes. The element of the i-th row and j-th column in matrix A is the connection coefficient a ij .
When i = j or there is no connecting edge i and j, the value of a ij is 0. Otherwise, the value of a ij is 1. From a macro perspective, on the premise of ensuring that it is suitable for the actual case, the model is appropriately simplified through the following three steps: (1) the power grid structure is simplified, mainly considering the 220 kV and above high-voltage transmission network; (2) the parallel lines on the same tower are merged, and the shunt capacitor lines are ignored to eliminate the case of multiple edges and self-direction in the power grid; and (3) regardless of the internal physical properties of each unit, all grounding wires are ignored.
In the undirected and unweighted topology model, the attributes of all nodes are the same, the energy can be transmitted along the edges of the network in two directions, and the weight of each edge is set to 1 [51,52]. When there is a parallel connection, the same number of edges leads to the same weight for multi-path, which makes it difficult to determine the minimum path and has a great impact on the analysis of network parameters.
In order to better adapt to reality, the power grid node is divided into three types: when the output power of the node is greater than the load power, the node is the source node; when the output power of the node is smaller than the load power, the node is the sink node; and other nodes are link nodes. Energy can only be transferred from the node with higher potential to the node with lower potential, transforming the undirected topology model into the directed topology model. The reactance value x ij0 is equivalent to the edge weight w ij , the effect of reactance is similar to that of the resistance in a direct current. When the reactance value of the line is larger, the transmitted power is smaller, transforming the unweighted topology model is into a weighted topology model. The n-order weighted square matrix W is redefined to represent the connection relationship between nodes while weighting. The element of i-th row and j-th column in the matrix W is W ij : When i is equal to j, the value of W ij is 0, and when i = j and there is a direct connecting edge between i and j, the value of W ij is equal to x ij . Otherwise, the value of W ij is infinite.

Characteristic Parameters of Power Grid Topology Model
After network weighting, the distance d ij between node i and node j is the reactance sum of the edges contained in the shortest path between two nodes. The weights of edge ij and edge jk are w ij and w jk , respectively, and Equation (3) can be deduced to calculate the distance d ik between nodes i and k. The network average distance L is obtained from the average distance of all pairs of nodes.
The weighted betweenness, measuring the importance of nodes and edges in the whole network, is divided into node betweenness C Bv and edge betweenness C Be . Node betweenness C Bv is the proportion of the weighted shortest paths passing through a node in the network to all the weighted shortest paths in the network. Edge betweenness C Be is the proportion of the weighted shortest paths passing through an edge in the network to all the weighted shortest paths in the network.
where σ ij refers to the number of shortest paths between node i and node j, and σ ijv and σ ije are the number of shortest paths between node i and node j containing node v or edge e, respectively. When node v is the link node, the number of shortest paths to enter node v is equal to the number of shortest paths to leave node v, which can be directly calculated by Equation (6). If node v belongs to the source node, the path only leaves from this node, as the source node is the starting node of the shortest path, leaving the number of shortest paths entering the node less than that leaving the node. In the same way, if node v belongs to the sink node, the number of shortest paths entering the node is larger than that leaving the node. In this case, a virtual node is introduced to transform node v into a link node, as shown in Figure 2.
where ij  refers to the number of shortest paths between node i and node j, and v ij  and e ij  are the number of shortest paths between node i and node j containing node v or edge e, respectively.
When node v is the link node, the number of shortest paths to enter node v is equal to the number of shortest paths to leave node v, which can be directly calculated by Equation (6). If node v belongs to the source node, the path only leaves from this node, as the source node is the starting node of the shortest path, leaving the number of shortest paths entering the node less than that leaving the node. In the same way, if node v belongs to the sink node, the number of shortest paths entering the node is larger than that leaving the node. In this case, a virtual node is introduced to transform node v into a link node, as shown in Figure 2. When the node v is the source node, the shortest path from v to any sink node only leaves from the edge and does not enter from the edge. The betweenness CBG of source node v is as follows: When the node v is the source node, the shortest path from v to any sink node only leaves from the edge and does not enter from the edge. The betweenness C BG of source node v is as follows: where n L is the number of sink nodes. When node v is the sink node, the shortest path from any source node to v only enters from the edge and does not leave from the edge. The betweenness C BL of sink node v is as follows: where n G is the number of source nodes.
In the reactance weighted network, the node degree k i , which is defined as the number of edges connected to the node, remains unchanged. The average degree K of a network is the average degree of all nodes and edges k i .

Power Grid Topology Modeling and Characteristic Parameter Analysis Process
The process of power grid topology modeling and characteristic parameter analysis is shown in Figure 3, and the detailed steps are as follows: Step 1: The actual power grid model is transformed into the power grid topology model according to power grid topology modeling rules in Section 3.1. The abstract class element is established, and then the node subclass and link subclass are established to store the initial data of the node and line, based on which the weighted model is established. Table 3 shows the data structure information of the nodes and lines in the power grid topology model. Step 2: All nodes are traversed, and then the value of power_out and power_in are judged. Nodes are divided into the source node, sink node, and link node, and then the directed weighted model is obtained.
Step 3: The characteristic parameters of the power grid are analyzed. The nodes and the line in the loop are traversed. When the traversed node is the previous node, the leaving degree increases by one. Similarly, if the traversed node is the succeeding node, the entering degree increases by one. The node degree is the sum of the entering degree and the leaving degree.
Step 4: The Floyd algorithm is applied to calculate the minimum distance between nodes. In the line reactance matrix, the nodes are traversed from the top left to the bottom right. If there is a direct connection between the node pairs, the distance between the node pairs is the line reactance. Conversely, if there is no direct connection, the distance between the node pairs is determined by Equation (3), and the shortest distance between the node pairs is obtained by traversing all nodes.
Step 5: The number of link nodes appearing in the node pair is counted while traversing the nodes. If the node is a source node, the betweenness is half the value calculated by the number of link nodes plus the number of sink nodes, according to Equation (7). If the node is a sink node, the betweenness is the half the value calculated by the number of link nodes plus the number of source nodes, according to Equation (8).
Step 6: The characteristic parameters of the nodes and lines are ranked, the cumulative distribution diagram is drawn.
Energies 2021, 14, 5183 9 of 27 Step 3: The characteristic parameters of the power grid are analyzed. The nodes and the line in the loop are traversed. When the traversed node is the previous node, the leaving degree increases by one. Similarly, if the traversed node is the succeeding node, the entering degree increases by one. The node degree is the sum of the entering degree and the leaving degree.
Step 4: The Floyd algorithm is applied to calculate the minimum distance between nodes. In the line reactance matrix, the nodes are traversed from the top left to the bottom right. If there is a direct connection between the node pairs, the distance between the node pairs is the line reactance. Conversely, if there is no direct connection, the distance between the node pairs is determined by Equation (3), and the shortest distance between the node pairs is obtained by traversing all nodes.
Step 5: The number of link nodes appearing in the node pair is counted while traversing the nodes. If the node is a source node, the betweenness is half the value calculated by the number of link nodes plus the number of sink nodes, according to Equation (7). If the node is a sink node, the betweenness is the half the value calculated by the number of link nodes plus the number of source nodes, according to Equation (8).
Step 6: The characteristic parameters of the nodes and lines are ranked, the cumulative distribution diagram is drawn. In this paper, the power grid nodes were divided into the source node, sink node, and link node according to the relationship between the energy inflow and outflow, and the line reactance value was introduced to transform the model into a directed weighted model. The characteristic parameters of power grid units were analyzed by the calculation method of reset betweenness, laying the foundation for the calculation of power grid vulnerability.   In this paper, the power grid nodes were divided into the source node, sink node, and link node according to the relationship between the energy inflow and outflow, and the line reactance value was introduced to transform the model into a directed weighted model. The characteristic parameters of power grid units were analyzed by the calculation method of reset betweenness, laying the foundation for the calculation of power grid vulnerability. The function vulnerability of the power grid system is defined as the disturbance degree to the function of the power grid unit in the specified time conditions. It is crucial to select the appropriate topology weight of power grid units when measuring functional vulnerability. When an extreme condition occurs, the external environment conditions and the state of the power grid unit vary, and the frequency of this variation is faster compared with the normal case. When the original betweenness index is used to measure the functional vulnerability of the power grid, it is difficult to measure the disturbance of the external environment to the power grid units under extreme conditions. However, the node vulnerability obtained from the betweenness analysis in the normal case may not be consistent with the actual load level or the reliability under extreme conditions. In particular, the node vulnerability in the normal case is low, yet the actual load is very large and the reliability under extreme conditions is low, leading to high risk. Therefore, the power grid functional vulnerability model based on reliability was established according to reference [46], and the time parameters were selected to weigh the power grid functional vulnerability. Under extreme conditions, if there is no time to repair the power grid lines in a short period of time, the unit reliability R i can be calculated following Equation (10):

Power Grid Functional Vulnerability Based on Reliability
where the reliability R i is a function of time t, λ i0 is the unit failure rate (1/h). When the reliability of the element obeys exponential distribution, the unit failure rate λ i0 is a constant. The reliability R i can be obtained from Equation (11).
Energies 2021, 14, 5183 10 of 25 According to the above modeling, in the weighted network, the line functional vulnerability B Bei is defined as the product of unit unreliability F ei and betweenness B ei : where the unreliability F ei can be obtained from the reliability R ei . The unreliability F ei is selected to measure the power grid functional vulnerability, for which the unreliability F ei is directly proportional to the vulnerability. The larger F ei is, the more vulnerable the unit becomes. The reliability R ei is inversely proportional to the vulnerability.

Evaluation Index of Power Grid Node Functional Vulnerability
Under extreme conditions, power grid nodes will be affected by line load fluctuations in the operation of the power grid system, on which the influence of the neighborhood cannot be reflected through the original betweenness index measuring functional vulnerability [53]. Therefore, the functional vulnerability measured by the node degree B Dvi is defined as the product of the power grid line influence weight on the node F vi and the node degree k i .
The functional vulnerability under extreme conditions measured by the betweenness of power grid nodes is defined as the product of power grid line influence weight on the node F vi and the node betweenness B vi .
The influence weight of the power grid line on the nodes F vi is defined as follows. The lines flowing through the nodes are divided into the current inflow line set E in and the current outflow line set E out according to the sign of the current. In this paper, the reliability of the inflow line set E in and the outflow line set E out were calculated in order. Using E in as an example, E in is composed of M elements (e in1 , e in2 , ..., e inm ). The influence weight of the grid node inflow F in is defined as Equation (15).
where R ei (t) is the reliability of unit e i , and the calculation method of F out is as same as F in . Only when the two subsystems of inflow, line set E in and outflow line set E out , are both normal is the system in the normal state. Therefore, the system conforms to the series system model, and the influence weight F vi of power grid nodes is defined as Equation (16).
If the node v is the source node or the sink node, either F in or F out does not exist, and the influence weight of the source node F Gi is defined as Equation (17).
Similarly, the grid node influence weight of sink node F Li is defined as Equation (18). As a country with vast territory, a complex climate, and a large temperature difference between the South and the North, there are extensive differences in climate and natural disasters in China. In addition, there are many trans-regional and long-distance transmission facilities in the power grid system, most of which are constructed in open air and greatly affected by the climate, causing the occurrence of extreme weather disturbances to be more frequent and the grid functional vulnerability to be far higher than the normal case [54]. According to the standard of the IEEE 346, weather intensity is divided into the following three types: normal weather, severe weather, and extreme severe weather [55]. When the grid unit is exposed to different weather intensities, the corresponding failure rate is different. The types of severe weather mainly include freezing, gale, lightning, rainstorm, and drought, of which freezing and gale are the major causes of power grid unit failure [56,57].
The occurrence of severe weather has obvious regional and temporal characteristics. In China, gale disasters mainly occur in the coastal areas of the South and East, and freezing disasters mainly occur in the northeast, north, and central regions. Therefore, the ability of transmission lines to withstand severe weather in different regions is different because of the different geographical locations and meteorological conditions. The geographical and seasonal factors endow a certain time characteristic to severe weather, resulting in the concentrated occurrence of specific severe weather conditions. Ice disasters often occur in winter, and gale weather occurs in summer and winter. By contrast, the frequency of severe weather in the spring and autumn is lower, and the intensity is lower than that in the summer and winter.
In summary, the impact of extreme conditions caused by gale weather and freezing weather on the functional vulnerability of the power grid is the major research direction of this paper. According to the relationship among line operation conditions, the line load, and its own attributes, the calculation model of line functional vulnerability was established. A node functional vulnerability calculation model was established to analyze the vulnerability distribution of the power grid system under extreme weather based on the relationship between the current inflow and outflow lines, together with the reliability of the power grid nodes.
The ice storm weather load model was analyzed consistent with [58], utilizing the circle model with maximum central strength to express the load of a transmission line section. The bivariate function in the following formula has appropriate properties and can be regarded as the basic model to describe the wind load and freezing load: where A refers to weather severity (amplitude), µ x and µ y are the coordinates of weather center, and σ x , σ y are the load calculation parameters related to the influence radius of weather.

Vulnerability Model of Gale Weather
The main influencing factor in gale weather is the wind load in the power grid. The wind load refers to the instantaneous wind force in severe weather, also known as gust wind. Since the wind load is 0 at the absolute center of the storm, Equation (19) is not applicable. A more practical model is obtained by subtracting an extra function with a smaller amplitude A and load calculation parameter σ [59].
where w is the influence of wind direction on transmission line wind load, A 1 and A 2 are parameters of gale weather severity (A 1 > A 2 ), and σ x1 , σ y1 , σ x2 , σ y2 are load calculation parameters (σ x1 > σ x2 , σ y1 > σ y2 ). In Equation (20), the influence of wind direction on transmission line wind load w(t) considers the wind direction factor composed of the angle β between the wind direction and the transmission line. In the northern hemisphere, the wind direction rotates anticlockwise around the weather center, perpendicular to the wind load of the power grid line. Therefore, the worst weather case is represented by β = π/2. Figure 4 shows the relative position between the line and the weather center as well as the effect of the wind load.
In Equation (20), the influence of wind direction on transmission line wind load w(t) considers the wind direction factor composed of the angle β between the wind direction and the transmission line. In the northern hemisphere, the wind direction rotates anticlockwise around the weather center, perpendicular to the wind load of the power grid line. Therefore, the worst weather case is represented by β = 2 /  . Figure 4 shows the relative position between the line and the weather center as well as the effect of the wind load.  In Figure 4, α is the angle between vector r = [(x 1 , y 1 ), (µ x , µ y )] and vector u = [(x 1 , y 1 ), (x 2 , y 2 )]. The wind direction is always perpendicular to vector r, and the angle β between the wind direction and the transmission line is acute or right. α = arccos ru |r||u| (21) The influence of the wind direction on the wind load of the transmission line is as follows: When β is 0, the wind load influence coefficient is 0, for which the wind load parallel to the line will not affect the line in practice. When the line is vertically subjected to the wind load, β is π/2. The influence coefficient of the wind load is 1, indicating that the wind load has the greatest influence at this time.
The line failure rate under extreme weather λ W (times/h, KM) is related to the line wind load.
According to reference [60], the discrete relationship between the wind load and the line failure rate under extreme weather is calculated as shown in Table 4, in which dl W is the line design value of the wind load.

Vulnerability Model of Freezing Weather
The main influencing factor in freezing weather is the icing of power grid lines. Icing will lead to an increase in the line load, leading to a change in the conductor cross-section and even the collapse of the high-voltage line tower, broken lines, and relaxation of the vibration of the iced line [60,61]. According to Equation (19), the ice load of power grid lines should also be related to the degree of severe weather and the distance from the line to the weather center. However, the ice load caused by icing is a cumulative process and is related to the duration of freezing weather. As shown in Equation (25), it can be feasible to solve the above problems using integrals.
))du (25) where A 3 is a parameter of the freezing weather severity. If only considering the influence of ice load and ignoring the melting of ice in a short time under extreme weather, the ice load of the line increases with the duration of extreme weather until the freezing weather does not affect the line and the ice load reaches the maximum, becoming constant.
The line failure rate under extreme weather λ I (times/h, KM) is related to the line ice load. λ The discrete relationship between the ice load and the line failure rate under extreme weather can be obtained according to Table 5 [61], in which dl I is the line design value of the ice load.

Collaborative Vulnerability Model of Gale and Freezing Weather
Gale and freezing weather often occur together in winter, leaving power grid lines to bear both a wind load and an ice load. Moreover, the failure modes of the wind load and the ice load act on the power grid lines independently. Therefore, this collaborative vulnerability model conforms to the competitive fault model: where λ is the total failure rate, and a W , a I are the influence weights of wind load and ice load, respectively.

Power Grid Functional Vulnerability Analysis Process
The process of the power grid functional vulnerability analysis is shown in Figure 5, and the detailed steps are as follows: Energies 2021, 14, 5183 14 of 25 constructed in this paper. Focusing on the typical extreme weather of gale weather and freezing weather, the load of the line section in these two cases was constructed, which was converted into the failure rate according to a certain functional relationship. Then, the line reliability was obtained and multiplied by the weighted betweenness of the line to obtain the line functional vulnerability. The node influence weight was obtained using the line reliability and multiplied by the node degree and betweenness to obtain the node reliabilities measured by degree and betweenness, respectively.

Example Analysis
In this paper, the IEEE-30 bus system was applied to the example analysis. First, the topology of the power grid was modeled to obtain the characteristic parameters. Then, the functional vulnerability and structural vulnerability of the power grid were calculated to analyze the vulnerability of the power grid load imbalance under extreme weather. Last, the constructed model was verified through a simulation experiment. Step 1: Set the type of extreme situation. In this paper, the extreme case was freezing and gale weather, and the calculations of freezing weather and gale weather were carried out respectively. Set the extreme case parameters, initial position, and wind direction.
Step 2: Traverse the line and calculate the load of the line under the wind force. Set the duration of wind force, calculate the angle between the line and wind direction according to Equations (21)- (23), and then calculate the load according to Equation (20). Finally, record the maximum load of each line within the duration.
Step 3: Traverse the line and calculate the load of the line under the action of freezing. Set the freezing action time, accumulate the line load under freezing action according to Equation (25), and finally, record the maximum load of each line within the duration.
Step 4: Calculate the line failure rate according to Equation (27).
Step 5: Calculate the line reliability and the functional vulnerability measured by betweenness according to Equations (10)- (12).
Step 6: Traverse the nodes and calculate the node function vulnerability influence weight according to Equations (15)- (19). Then, calculate the degree and betweenness to measure the node functional vulnerability according to Equations (13) and (14).
Step 7: Traverse the nodes and edges, ranking the functional fragility. The calculation method of unit functional vulnerability under extreme weather was constructed in this paper. Focusing on the typical extreme weather of gale weather and freezing weather, the load of the line section in these two cases was constructed, which was converted into the failure rate according to a certain functional relationship. Then, the line reliability was obtained and multiplied by the weighted betweenness of the line to obtain the line functional vulnerability. The node influence weight was obtained using the line reliability and multiplied by the node degree and betweenness to obtain the node reliabilities measured by degree and betweenness, respectively.

Example Analysis
In this paper, the IEEE-30 bus system was applied to the example analysis. First, the topology of the power grid was modeled to obtain the characteristic parameters. Then, the functional vulnerability and structural vulnerability of the power grid were calculated to analyze the vulnerability of the power grid load imbalance under extreme weather. Last, the constructed model was verified through a simulation experiment. Figure 6 shows the wiring diagram of the IEEE-30 system, which is composed of 30 nodes and 41 lines. The nodes can be divided into 5 source nodes, 19 sink nodes, and 6 link nodes, according to their properties. Specifically, although node 5 is both an engine node and a load node, its load power is greater than its generator input power, which means that node 5 is a sink node. The specific type of each node is shown in Table 6.  Figure 6 shows the wiring diagram of the IEEE-30 system, which is composed of 30 nodes and 41 lines. The nodes can be divided into 5 source nodes, 19 sink nodes, and 6 link nodes, according to their properties. Specifically, although node 5 is both an engine node and a load node, its load power is greater than its generator input power, which means that node 5 is a sink node. The specific type of each node is shown in Table 6.   6 6, 9, 22, 25, 26, 27 The topology model of the power grid was established according to Section 3, transforming the power grid into a directed weighted topology model. Figure 7 shows the schematic diagram of the IEEE-30 system topology model, in which the coordinate axis is established, and the unit length is 50 km. The nodes and lines are labeled, and the coordinate positions were determined respectively, distinguished by different colors according to node types.   6 6, 9, 22, 25, 26, 27 The topology model of the power grid was established according to Section 3, transforming the power grid into a directed weighted topology model. Figure 7 shows the schematic diagram of the IEEE-30 system topology model, in which the coordinate axis is

Characteristic Parameter Analysis
All the lines of the power grid were traversed in turn, accumulating the occurrence times of each node in the previous node and subsequent nodes of the line.
If the occurrence time in the previous nodes of node v is m, the leaving degree of node v is m. Similarly, if the occurrence time in the subsequent nodes of node v is n, the entering degree of node v is n. The node degree is the sum of the entering degree and leaving degree. Table 7 indicates the ranking of the top seven nodes by degrees, and it is obvious that the nodes with larger degrees connect more lines, approach more neighboring nodes, and bear heavier loads. Table 7. Node degree ranking.

Ranking
Node ID Node Degree  1  6  7  2  10  6  3  12  5  4  4  4  4  2  4  4  15  4  4 27 4 The node degree distribution is shown in Figure 8, in which the abscissa and ordinate indicate the node degree and the probability accumulation of the degree, respectively.

Characteristic Parameter Analysis
All the lines of the power grid were traversed in turn, accumulating the occurrence times of each node in the previous node and subsequent nodes of the line.
If the occurrence time in the previous nodes of node v is m, the leaving degree of node v is m. Similarly, if the occurrence time in the subsequent nodes of node v is n, the entering degree of node v is n. The node degree is the sum of the entering degree and leaving degree. Table 7 indicates the ranking of the top seven nodes by degrees, and it is obvious that the nodes with larger degrees connect more lines, approach more neighboring nodes, and bear heavier loads. Table 7. Node degree ranking.

Ranking
Node ID Node Degree The node degree distribution is shown in Figure 8, in which the abscissa and ordinate indicate the node degree and the probability accumulation of the degree, respectively.
After traversing all nodes of the power grid in turn, the shortest path between all node pairs can be identified according to the Floyd algorithm described in Section 3. The weighted betweennesses of the node and the line are calculated according to Equations (5) and (6), respectively. Table 8 shows the ranking of the top 5 nodes according to betweenness, and Table 9 shows the ranking of the top 7 lines according to betweenness. After traversing all nodes of the power grid in turn, the shortest path between all node pairs can be identified according to the Floyd algorithm described in Section 3. The weighted betweennesses of the node and the line are calculated according to Equations (5) and (6), respectively. Table 8 shows the ranking of the top 5 nodes according to betweenness, and Table 9 shows the ranking of the top 7 lines according to betweenness.  Table 9. Line weighted betweenness ranking.

Ranking
Line ID Line Betweenness  1  14  38  2  11  36  3  4  31  4  41  31  5  6  30  5  7  30 The distribution of node betweenness is shown in Figure 9. The abscissa and the ordinate are the node betweenness and the probability accumulation of the betweenness, respectively.   Table 9. Line weighted betweenness ranking.

Ranking
Line ID Line Betweenness   1  14  38  2  11  36  3  4  31  4  41  31  5  6  30  5  7  30 The distribution of node betweenness is shown in Figure 9. The abscissa and the ordinate are the node betweenness and the probability accumulation of the betweenness, respectively. The betweenness distribution of the lines is shown in Figure 10. The abscissa and the ordinate are the line betweenness and the probability accumulation of the betweenness, respectively. The betweenness distribution of the lines is shown in Figure 10. The abscissa and the ordinate are the line betweenness and the probability accumulation of the betweenness, respectively. The betweenness distribution of the lines is shown in Figure 10. The abscissa and the ordinate are the line betweenness and the probability accumulation of the betweenness, respectively.

Simulation and Analysis of Power Grid Functional Vulnerability
The extreme case set in this paper was severe weather with a synergistic effect of gale and freezing. The radius of the gale action was 200 km, and the radius of the freezing action was 130 km. Under the action of gale, the load calculation parameters σx1, σy1, σx2, and σy2 satisfied σx1 = σy1 = 0.4 RW and σx2 = σy2 = 0.05 RW; the amplitudes of severe weather A1, A2, and A3 were 20 m/s, 4 m/s, and 30 mm/h, respectively; and the amplitudes of severe weather A1, A2, and A3 were 25 m/s, 5 m/s, and 50 mm/h, respectively. The influence weights of the wind load aW and the ice load aI were 0.1 and 0.9, respectively. The initial center position of gale weather and freezing weather was (6,10), and the center moving speed of gale weather and freezing weather was 100 km/h in the southern direction. Figure 11 shows the schematic diagram of the extreme weather scope.

Simulation and Analysis of Power Grid Functional Vulnerability
The extreme case set in this paper was severe weather with a synergistic effect of gale and freezing. The radius of the gale action was 200 km, and the radius of the freezing action was 130 km. Under the action of gale, the load calculation parameters σ x1 , σ y1 , σ x2 , and σ y2 satisfied σ x1 = σ y1 = 0.4 R W and σ x2 = σ y2 = 0.05 R W ; the amplitudes of severe weather A 1 , A 2 , and A 3 were 20 m/s, 4 m/s, and 30 mm/h, respectively; and the amplitudes of severe weather A 1 , A 2 , and A 3 were 25 m/s, 5 m/s, and 50 mm/h, respectively. The influence weights of the wind load a W and the ice load a I were 0.1 and 0.9, respectively. The initial center position of gale weather and freezing weather was (6,10), and the center moving speed of gale weather and freezing weather was 100 km/h in the southern direction. Figure 11 shows the schematic diagram of the extreme weather scope.

Load Calculation
According to the computing method proposed in Section 3, the wind load and ice load of each line in the power grid under extreme weather can be obtained. The four typical lines 6, 14, 21, and 24 were selected to perform the load analysis. Figure 12 shows the wind load on these four lines.
The following could be inferred from the wind load curve: (1) Because of the different distances and wind direction angles between the line and the gale weather center, the wind load curve of each line was a single peak or a double peak curve. (2) The arrival time of the wind load extreme value was related to the relative position between the line and the beginning of the weather center. The relative positions of the four lines were different from that of the weather center, and the arrival time of wind load extreme value of these four lines also varied. Line 6 was to the north of line 14, and the distance from the initial coordinate of the weather center was shorter than that of line 14; therefore line 6 was first to be affected by the wind and reach the extreme value. Line 14 and line 21 were located at the same latitude and were almost simultaneously affected by the wind.

Load Calculation
According to the computing method proposed in Section 3, the wind load and ice load of each line in the power grid under extreme weather can be obtained. The four typical lines 6, 14, 21, and 24 were selected to perform the load analysis. Figure 12 shows the wind load on these four lines. The following could be inferred from the wind load curve:

Load Calculation
According to the computing method proposed in Section 3, the wind load and ice load of each line in the power grid under extreme weather can be obtained. The four typical lines 6, 14, 21, and 24 were selected to perform the load analysis. Figure 12 shows the wind load on these four lines. The following could be inferred from the wind load curve:   Figure 13 shows the ice load on the four lines. The following were indicated from the analysis of the ice load curve: (1) Since the effect of ice load on the line was cumulative, the ice load curve of each line presents S-type growth, and a moment of ice load maximum value exists, after which the ice load of the line did not increase.

Functional Vulnerability Calculation
Based on the characteristic parameters of the power grid topology model in the previous section and the load calculated in this section, the design values dlW and dli of the wind and ice loads for all lines were set to 18.95 m/s and 50 mm, respectively. The line failure rate could be obtained according to the discrete relationship between the load and line failure rates in Tables 4 and 5. Then, the reliability of the transmission line in this extreme weather could be obtained from Equation (10). Finally, the betweenness of lines was calculated using Equation (12) to measure the functional vulnerability, and the ranking was carried out. Table 10 shows the top 10 lines in terms of function vulnerability measured by betweenness in severe weather and extreme severe weather.

Functional Vulnerability Calculation
Based on the characteristic parameters of the power grid topology model in the previous section and the load calculated in this section, the design values dl W and dl i of the wind and ice loads for all lines were set to 18.95 m/s and 50 mm, respectively. The line failure rate could be obtained according to the discrete relationship between the load and line failure rates in Tables 4 and 5. Then, the reliability of the transmission line in this extreme weather could be obtained from Equation (10). Finally, the betweenness of lines was calculated using Equation (12) to measure the functional vulnerability, and the ranking was carried out. Table 10 shows the top 10 lines in terms of function vulnerability measured by betweenness in severe weather and extreme severe weather.
According to the index definition of power grid node functional vulnerability in Section 3, the mean reliability value of the node inflow line set and node outflow line set was calculated by Equation (15), and the influence weight of power grid node vulnerability was calculated by Equations (16)- (18). Then, the degree of nodes obtained from Equation (13) was used to measure the functional vulnerability and was ranked. Table 11 shows the top 10 lines in terms of functional vulnerability measured by the node degree in severe weather and extreme severe weather. The node degrees obtained from Equation (14) were used to measure the functional vulnerability, and functional vulnerability was then ranked. Table 12 shows the top 10 lines of function vulnerability measured by node betweenness in severe weather and extreme severe weather. (1) Functional vulnerability mainly reflects the difficulty of power grid units exiting the system operation due to internal causes or external interference in a short time. It can be inferred from Tables 10-12 that with increased proximity to the weather center, the power grid unit will bear a greater weather disturbance the higher the functional vulnerability. However, the units far away from the weather center are generally less vulnerable. (2) Table 10 shows that the vulnerability of lines in extreme severe weather is more than four times that in severe weather, while the betweenness of the lines in two extreme conditions is the same. Hence, the lines are generally vulnerable to environmental interference in extreme severe weather.
(3) Figure 14 shows the distribution of functional vulnerability measured by degree and betweenness in extreme severe weather, of which the horizontal axis is the node degree measurement index, the vertical axis is the node betweenness measurement index, and the area of scattered points is the grid betweenness. As shown in Figure 14, there is a linear positive correlation between the node degree measuring index and the betweenness measuring index. The degree measuring index considers the influence of the node on its neighborhood, and the betweenness measuring index considers the overall influence of the node on the power grid system. (4) In the normal case, the nodes with high characteristic parameters in the topology model are not guaranteed to have a higher functional vulnerability and vice versa. For instance, nodes 10 and node 25 have lower betweenness in the normal case, but their functional vulnerability is high. Therefore, the nodes that need to be protected and monitored are those with high functional vulnerability, not those with high betweenness.
Energies 2021, 14, 5183 24 of 27 there is a linear positive correlation between the node degree measuring index and the betweenness measuring index. The degree measuring index considers the influence of the node on its neighborhood, and the betweenness measuring index considers the overall influence of the node on the power grid system. (4) In the normal case, the nodes with high characteristic parameters in the topology model are not guaranteed to have a higher functional vulnerability and vice versa. For instance, nodes 10 and node 25 have lower betweenness in the normal case, but their functional vulnerability is high. Therefore, the nodes that need to be protected and monitored are those with high functional vulnerability, not those with high betweenness. Figure 14. Functional vulnerability distribution of power grid nodes measured by degree and betweenness.

Conclusions
In this paper, the topology of the power grid was modeled based on the complex network theory, which considers the impact of extreme weather on power grid and researches the functional vulnerability of the power grid. The main contributions of this paper are as follows: a topology model of the power grid based on complex network modeling was constructed; a functional vulnerability model based on reliability was established; and the proposed model was verified through a simulation experiment.
The major conclusions are as follows: (1) A power grid topology model based on complex network modeling was constructed. According to the relationship between output power and load power, the nodes were divided into source nodes, sink nodes, and link nodes, and the undirected topology model was transformed into a directed topology model according to node voltage. By introducing the line reactance, the weighted model was transformed into the weighted model, which can better identify the vulnerability of different units. In addition, the node degree, node weighted betweenness, and line weighted betweenness were defined as the characteristic parameters of the power grid topology model to analyze its topology characteristics. (2) A functional vulnerability model based on reliability was established. The reliability

Conclusions
In this paper, the topology of the power grid was modeled based on the complex network theory, which considers the impact of extreme weather on power grid and researches the functional vulnerability of the power grid. The main contributions of this paper are as follows: a topology model of the power grid based on complex network modeling was constructed; a functional vulnerability model based on reliability was established; and the proposed model was verified through a simulation experiment.
The major conclusions are as follows: (1) A power grid topology model based on complex network modeling was constructed. According to the relationship between output power and load power, the nodes were divided into source nodes, sink nodes, and link nodes, and the undirected topology model was transformed into a directed topology model according to node voltage. By introducing the line reactance, the weighted model was transformed into the weighted model, which can better identify the vulnerability of different units. In addition, the node degree, node weighted betweenness, and line weighted betweenness were defined as the characteristic parameters of the power grid topology model to analyze its topology characteristics. (2) A functional vulnerability model based on reliability was established. The reliability of the line was obtained by calculating the load in gale weather, freezing weather, and freezing gale weather. The functional vulnerability of the line was calculated by weighting the line reliability with the line betweenness. The node functional vulnerability weight was defined by the line reliability of the node neighborhood, and then the node functional vulnerability model measured by degree and betweenness was obtained based on the above work.
(3) The model and algorithm proposed in this paper were verified by a simulation experiment. Taking the IEEE-30 bus system as an example, the topology model of power grid system was established to simulate the load and vulnerability of the power grid in freezing and gale weather. The power grid functional vulnerability level under different units and modes was analyzed to put forward suggestions on the maintenance of the power grid, which was verified through the simulation experiment. The simulation results show that the vulnerability of the power grid unit reaches the highest value under the synergistic effect of extreme weather.