# Preventive Control Strategy of Cascading Fault considering Safety and Economy

Authors: Hui-Qiong Deng; Jie Luo; Chen-Chen Li; Pei-Qiang Li; Rong-Jin Zheng
Journal: Journal of Electrical and Computer Engineering
 (2021)

---

## Abstract

The operation and structure of the power system are becoming increasingly complex, and the probability of cascading fault increases. To this end, this paper proposes a cascading fault preventive control strategy that considers safety and the economy. First is to give a mathematical form to discriminate the cascading fault according to the action characteristics of the current-type backup protection. Second, the safety and economy of the system are evaluated in terms of power grid safety margin and generation operation cost, respectively, the initial faults are selected based on the power grid vulnerability and safety margin, and a cascading fault preventive control model is constructed for different initial faults’ scenarios. The model is a two-layer optimization mathematical model, with the inner model being solved by particle swarm optimization to minimize the power grid safety margin. The outer model is solved by the multiobjective algorithm to minimize generation cost and maximizing power grid safety margin. Finally, the calculated Pareto set is evaluated using fuzzy set theory to determine the optimal generator output strategy. The feasibility of the proposed method is verified by conducting a simulation study with the IEEE39 node system as an example.


---

## Body

## 1. Introduction

The modern power grid has evolved toward a highly intelligent, complex, and high-capacity Internet trend, promoting rapid social and economic development along with a series of fault problems. A fault in a small area can quickly involve the entire power grid and then develop into a major blackout accident. Therefore, learning the lessons from accidents, studying the chain fault development process, and formulating sound cascading fault defense decisions have important application values to enhance the safety and reliability of power grid operation.At the present stage, the most notable studies on modeling aspects of cascading faults are related models based on self-organized criticality theory [1–7], such as the OPA model based on power system weak links [1, 2], CASCADE model based on component cascading failure [3, 4], and Manchester model based on AC-DC currents [5–7]. These models explore the intrinsic mechanism of large outage accidents triggered by cascading faults and reveal that large outage accidents have self-organized critical characteristics, which are of great reference value for preventing cascading faults.Numerous studies have shown that, in the early stages of cascading faults, the action of relay protection and power flow redistribution are the main causes of further deterioration of the grid operation. To prevent cascading faults, literature [8] proposed a vulnerability-based cascading fault preventive control model, which quantifies the risk of initial faults leading to subsequent faults and improves the safety of power system operation. Literature [9] establishes a Bayesian network-based chain fault model to minimize fault risk and achieve preventive control using machine and load cutting. Literature [10] proposes a multiobjective algorithm-based chain trip prevention strategy by screening sensitive control nodes with the goal of less equipment adjustment and less load cutting. In literature [11], multistage dynamic game theory is introduced to analyze the process of a complex system cascading faults and give a risk assessment method of cascading fault accident chain to determine the best preventive control scheme. Literature [12, 13] established a coordinated control model for cascading faults considering safety and economy and obtained generation output adjustment and load cutting strategies by solving them. Although the above studies have good reference value for cascading fault analysis and prevention, most of these studies are based on accident chain control methods, which are computationally intensive and time-consuming, can hardly reflect the actual development of cascading faults, and may not be used to guide cascading fault prevention in practice.Based on the above analysis, this paper mainly focuses on the trigger phase of cascading faults, reflects the safety of the power grid by the distance between the nodal injection power in the current power grid operation state and the nodal injection power in the critical state, reflects the economy of the power grid by the generation operation cost, and establishes a cascading fault preventive control model for different initial faults based on the selection of initial faults using power grid vulnerability and safety margin. For this model, the Pareto solution set under the model is calculated by combining particle swarm algorithm and multiobjective algorithm, and the Pareto solution set is evaluated using fuzzy set theory to determine the optimal generator output strategy. The simulation analysis of the IEEE39 node system is carried out to verify the effectiveness of the method proposed in this paper.


## 2. Mathematical Expression of Cascading Fault

In the event of a blackout, the cascading fault of the power grid is usually caused by the backup protection action caused by the power flow transfer, which usually occurs in the early stage of the blackout. To facilitate the analysis, this paper mainly studies a common cascading fault mode: after the initial fault branch is cut off, with the power flow transfer, the cascading fault occurs due to the corresponding electrical quantity of the remaining branches entering the action area of the backup protection (mainly current protection and distance protection (Section3)). In a cascading fault sequence, if the initial fault is regarded as the first-level fault, the cascading trip after the initial fault is cut off can be regarded as the second-level fault, the second-level fault may also cause the third level fault, and so on. This paper only studies the second-level fault and considers it together with the initial fault.In this paper, taking each branch configured with current-type backup protection as an example, assuming that, in a certain power grid, after the initial fault of branchLa is cut off by the backup protection device at a certain time, whether any other branch Lb will have cascading fault can be judged according to whether the current flow detected by the configured current protection enters the protection action area. Let the branch Lb be between node i and node j, and its i-side is configured with current-type backup protection; according to the protection calibration value and measured current, the equation is shown in formula (1):(1)Ii.distb=Iib.set−Iib.m,where Iib.m is the current measured by i-side backup protection of branch Lb, Iib.set is the setting value of i-side backup protection of branch Lb, and Ii.distb is the distance between Iib.set and Iib.m. It can be seen from the specific performance of branch cascading fault that when Ii.distb>0, branch Lb will not have cascading fault; when Ii.distb<0, branch Lb will have cascading fault; when Ii.distb=0, branch Lb is just at the boundary of cascading fault. Similarly, if the j-side of branch Lb is also equipped with current protection, just change the superscript from i to j.Let the number of branchesLb in the power grid be b, and there are l remaining branches except for the initial fault. Considering all the remaining branches, we can further give the expression shown in the following formula:(2)Ib=minIi.distb,Ij.distbI=minIbb=1,2,…,l.It can be seen from formula (2) that, from the power grid level, when I<0, the power grid is in the cascading fault state, when I>0, the power grid is in the safe state, and when I=0, the power grid is just in the critical state of cascading fault. In this way, according to formula (2), the mathematical expression for judging the cascading fault of the power grid is given.


## 3. The Basic Model of Preventive Control

Preventive control is usually a control measure before faults occur when the initial fault has not yet occurred and the power grid dispatcher has enough time to adjust the current power grid operation to a safer operating state. Therefore, for the power grid interlocking fault scenario, it will be the problem studied in this paper to establish a preventive control model to ensure the improvement of power grid safety while taking proper care of the power grid economy.


### 3.1. System Safety

In the triggering stage of cascading fault, the injected power of power grid nodes before and after initial fault removal is unchanged. Therefore, for a given power grid and initial fault, a method for analyzing the safety of the power grid for cascading fault can be given according to the nodal injection power [14].Assume that the node power combination consists of active power and reactive power and the nodal injection power combination for which no cascading fault occurs is called the safe combination; the nodal injection power combination for which a cascading fault can be triggered is called the dangerous combination. When the power system operates normally, the farther the nodal injection power combination of the current power grid is from the dangerous combination, the higher the safety level of the power grid for cascading faults after the initial fault occurs is. If the shortest distance between the nodal injection power combination and the dangerous combination of the current power grid can be analyzed before the initial fault occurs, the safety of the power grid and the allowable variation range of the nodal injection power of the power grid can be observed. For this purpose, Figure1 gives a further explanation.Figure 1
Safety level of the power grid for cascading failure.In Figure1, the nodal injection power combination (P, Q) is denoted by S. S0 is the nodal injection power combination corresponding to the current power grid operation state (safety combination), denoted by ·, while the nodal injection power combination that can trigger cascading fault (dangerous combination) is denoted by ∗; S1 is the nodal injection power combination closest to S0 in the dangerous combination. With S0 being the center of the sphere, the distance between S0 and S1 is used as the radius to form a circular sphere A, and the area outside the sphere is denoted by B. From Figure 1 and the previous analysis, it is clear that, for a certain initial fault, no cascading fault occurs when the power grid is operated at any position within sphere A. If we can manage to calculate the radius of sphere A, we can observe a plausible safety region. Therefore, the radius R of the sphere A can be used to measure the current safety level of the grid for interlocking faults, and R can be expressed as(3)R=S0−S12,where R is the distance between S0 and S1. Furthermore, if it is extended to a multinode power system, formula (3) can be written in the form shown in the following formula:(4)R=min∑d=1DS0.d−S1.d2,where D is the total number of nodes of the system; R can be used to characterize the safety margin of the power grid. The greater the R value is, the safer the power grid will be. Therefore, this paper takes it as the safety index of the system.


### 3.2. System Economy

In this paper, the power generation operation cost is defined as the economic index of the system, and the power generation operation costC can be expressed as [15](5)C=∑u=1GEauPGu2+buPGu+cu,where GE is the set of generators participating in the adjustment; PGu is the active output of the generator at node u; au, bu, and cu are the operating cost coefficients.


### 3.3. Mathematical Model of Preventive Control

In this paper, the power output of the generator set is adjusted as a preventive control measure, the control variable is the nodal injection power of the generator, and the interlocking fault prevention control model considering safety and economy is established. The model can be divided into two levels (called inner layer and outer layer, respectively). Among them, the inner model mainly calculates the power grid safety margin, and the outer model optimizes the output of the generator set under the consideration of power grid safety margin and power generation operation cost.


#### 3.3.1. Inner Layer Optimization Model

The objective function of inner layer optimization is to minimize the safety margin of the power grid, to characterize the shortest distance between the initial operational state of the power grid and the critical state of cascading fault. The objective functionF1 is as follows:(6)F1=minR.Constraints:(1)
Critical constraints of the power gridTo search the critical state closest to the current operation state of the power grid, the cascading fault criticality constraint needs to be introduced. According to the analysis of formulas (1) and (2), under the action of S, when the initial fault La is cut off, the conditions under which the power grid is just in the critical state of cascading fault are as follows:(7)I=minIb=0.(2)
Power flow constraintsThe power grid also needs to ensure that the power flow constraints are satisfied before and after the fault in the cascading fault critical state. Among them, the equation constraint is the system node power flow constraint, while the inequality constraint includes generator active and reactive power output constraint, node voltage constraint, and branch transmission power constraint. For the convenience of description, they can be abbreviated to the form shown in the following formula:(8)s.t.m0x=0m1x=0g0x≥0g1x≥0,where m0 and m1 are the reduced form of the power flow equation constraint corresponding to the power grid before and after the initial fault, respectively; g0 and g1 are the reduced form of the power flow inequality constraint corresponding to the power grid before and after the initial fault, respectively; x is the state variable.


#### 3.3.2. Outer Layer Optimization Model

The outer layer optimization is based on the inner layer model to solve the power grid safety margin. The objective function is to maximize the power grid safety margin and minimize the power generation operation cost. The expression of the objective function is as follows:(9)F2=maxF1,F3=minC,where F2 and F3 are the objective functions of the outer model; F1 is the power grid safety margin calculated by the inner model; C is the power generation operation cost.Constraints:(1)
Safety constraints of the power gridAccording to the previous analysis, under the action ofS, when the initial fault La is cut off, the power grid is in a safe state without cascading fault. The conditions are as follows:(10)I=minIb>0.(2)
Power flow constraintsIn the safe state, the power grid also needs to meet certain power flow constraints. The constraints similar to formula (8) are given, which can be expressed as(11)s.t.d0x=0d1x=0h0x≥0h1x≥0,whered0 and d1 are the reduced form of the power flow equation constraint corresponding to the power grid before and after the initial fault, respectively; h0 and h1 are the reduced form of the power flow inequality constraint corresponding to the power grid before and after the initial fault, respectively.


## 3.1. System Safety

In the triggering stage of cascading fault, the injected power of power grid nodes before and after initial fault removal is unchanged. Therefore, for a given power grid and initial fault, a method for analyzing the safety of the power grid for cascading fault can be given according to the nodal injection power [14].Assume that the node power combination consists of active power and reactive power and the nodal injection power combination for which no cascading fault occurs is called the safe combination; the nodal injection power combination for which a cascading fault can be triggered is called the dangerous combination. When the power system operates normally, the farther the nodal injection power combination of the current power grid is from the dangerous combination, the higher the safety level of the power grid for cascading faults after the initial fault occurs is. If the shortest distance between the nodal injection power combination and the dangerous combination of the current power grid can be analyzed before the initial fault occurs, the safety of the power grid and the allowable variation range of the nodal injection power of the power grid can be observed. For this purpose, Figure1 gives a further explanation.Figure 1
Safety level of the power grid for cascading failure.In Figure1, the nodal injection power combination (P, Q) is denoted by S. S0 is the nodal injection power combination corresponding to the current power grid operation state (safety combination), denoted by ·, while the nodal injection power combination that can trigger cascading fault (dangerous combination) is denoted by ∗; S1 is the nodal injection power combination closest to S0 in the dangerous combination. With S0 being the center of the sphere, the distance between S0 and S1 is used as the radius to form a circular sphere A, and the area outside the sphere is denoted by B. From Figure 1 and the previous analysis, it is clear that, for a certain initial fault, no cascading fault occurs when the power grid is operated at any position within sphere A. If we can manage to calculate the radius of sphere A, we can observe a plausible safety region. Therefore, the radius R of the sphere A can be used to measure the current safety level of the grid for interlocking faults, and R can be expressed as(3)R=S0−S12,where R is the distance between S0 and S1. Furthermore, if it is extended to a multinode power system, formula (3) can be written in the form shown in the following formula:(4)R=min∑d=1DS0.d−S1.d2,where D is the total number of nodes of the system; R can be used to characterize the safety margin of the power grid. The greater the R value is, the safer the power grid will be. Therefore, this paper takes it as the safety index of the system.


## 3.2. System Economy

In this paper, the power generation operation cost is defined as the economic index of the system, and the power generation operation costC can be expressed as [15](5)C=∑u=1GEauPGu2+buPGu+cu,where GE is the set of generators participating in the adjustment; PGu is the active output of the generator at node u; au, bu, and cu are the operating cost coefficients.


## 3.3. Mathematical Model of Preventive Control

In this paper, the power output of the generator set is adjusted as a preventive control measure, the control variable is the nodal injection power of the generator, and the interlocking fault prevention control model considering safety and economy is established. The model can be divided into two levels (called inner layer and outer layer, respectively). Among them, the inner model mainly calculates the power grid safety margin, and the outer model optimizes the output of the generator set under the consideration of power grid safety margin and power generation operation cost.


### 3.3.1. Inner Layer Optimization Model

The objective function of inner layer optimization is to minimize the safety margin of the power grid, to characterize the shortest distance between the initial operational state of the power grid and the critical state of cascading fault. The objective functionF1 is as follows:(6)F1=minR.Constraints:(1)
Critical constraints of the power gridTo search the critical state closest to the current operation state of the power grid, the cascading fault criticality constraint needs to be introduced. According to the analysis of formulas (1) and (2), under the action of S, when the initial fault La is cut off, the conditions under which the power grid is just in the critical state of cascading fault are as follows:(7)I=minIb=0.(2)
Power flow constraintsThe power grid also needs to ensure that the power flow constraints are satisfied before and after the fault in the cascading fault critical state. Among them, the equation constraint is the system node power flow constraint, while the inequality constraint includes generator active and reactive power output constraint, node voltage constraint, and branch transmission power constraint. For the convenience of description, they can be abbreviated to the form shown in the following formula:(8)s.t.m0x=0m1x=0g0x≥0g1x≥0,where m0 and m1 are the reduced form of the power flow equation constraint corresponding to the power grid before and after the initial fault, respectively; g0 and g1 are the reduced form of the power flow inequality constraint corresponding to the power grid before and after the initial fault, respectively; x is the state variable.


### 3.3.2. Outer Layer Optimization Model

The outer layer optimization is based on the inner layer model to solve the power grid safety margin. The objective function is to maximize the power grid safety margin and minimize the power generation operation cost. The expression of the objective function is as follows:(9)F2=maxF1,F3=minC,where F2 and F3 are the objective functions of the outer model; F1 is the power grid safety margin calculated by the inner model; C is the power generation operation cost.Constraints:(1)
Safety constraints of the power gridAccording to the previous analysis, under the action ofS, when the initial fault La is cut off, the power grid is in a safe state without cascading fault. The conditions are as follows:(10)I=minIb>0.(2)
Power flow constraintsIn the safe state, the power grid also needs to meet certain power flow constraints. The constraints similar to formula (8) are given, which can be expressed as(11)s.t.d0x=0d1x=0h0x≥0h1x≥0,whered0 and d1 are the reduced form of the power flow equation constraint corresponding to the power grid before and after the initial fault, respectively; h0 and h1 are the reduced form of the power flow inequality constraint corresponding to the power grid before and after the initial fault, respectively.


## 3.3.1. Inner Layer Optimization Model

The objective function of inner layer optimization is to minimize the safety margin of the power grid, to characterize the shortest distance between the initial operational state of the power grid and the critical state of cascading fault. The objective functionF1 is as follows:(6)F1=minR.Constraints:(1)
Critical constraints of the power gridTo search the critical state closest to the current operation state of the power grid, the cascading fault criticality constraint needs to be introduced. According to the analysis of formulas (1) and (2), under the action of S, when the initial fault La is cut off, the conditions under which the power grid is just in the critical state of cascading fault are as follows:(7)I=minIb=0.(2)
Power flow constraintsThe power grid also needs to ensure that the power flow constraints are satisfied before and after the fault in the cascading fault critical state. Among them, the equation constraint is the system node power flow constraint, while the inequality constraint includes generator active and reactive power output constraint, node voltage constraint, and branch transmission power constraint. For the convenience of description, they can be abbreviated to the form shown in the following formula:(8)s.t.m0x=0m1x=0g0x≥0g1x≥0,where m0 and m1 are the reduced form of the power flow equation constraint corresponding to the power grid before and after the initial fault, respectively; g0 and g1 are the reduced form of the power flow inequality constraint corresponding to the power grid before and after the initial fault, respectively; x is the state variable.


## 3.3.2. Outer Layer Optimization Model

The outer layer optimization is based on the inner layer model to solve the power grid safety margin. The objective function is to maximize the power grid safety margin and minimize the power generation operation cost. The expression of the objective function is as follows:(9)F2=maxF1,F3=minC,where F2 and F3 are the objective functions of the outer model; F1 is the power grid safety margin calculated by the inner model; C is the power generation operation cost.Constraints:(1)
Safety constraints of the power gridAccording to the previous analysis, under the action ofS, when the initial fault La is cut off, the power grid is in a safe state without cascading fault. The conditions are as follows:(10)I=minIb>0.(2)
Power flow constraintsIn the safe state, the power grid also needs to meet certain power flow constraints. The constraints similar to formula (8) are given, which can be expressed as(11)s.t.d0x=0d1x=0h0x≥0h1x≥0,whered0 and d1 are the reduced form of the power flow equation constraint corresponding to the power grid before and after the initial fault, respectively; h0 and h1 are the reduced form of the power flow inequality constraint corresponding to the power grid before and after the initial fault, respectively.


## 4. Preventive Control Strategies under Different Initial Fault Scenarios

Cascading faults often start with an initial fault and then trigger subsequent faults. Therefore, it is important to identify the vulnerable branches in the system and to preventive control. In this paper, based on the power network topology and operation state, we give the definition and evaluation methods of vulnerability indicators, further filter out the important initial faults according to the safety margin, and then propose the preventive control strategies for different initial fault scenarios.


### 4.1. Selection of Initial Fault

#### 4.1.1. Initial Fault Preselection Based on Vulnerability

(1) Severity Index. According to the previous analysis, without considering the uncertainty of impedance relay operation and other blocking conditions, the severity of the initial fault can be effectively evaluated from the branch level by the remaining system branches subjected to the initial fault removal. Let the number of the initial faults in the power grid be a. The severity index shown in formula (12) can be obtained after further processing of formula (1) as(12)ka=wamaxIib.setIib.m+wbmaxIjb.setIjb.m,where wa and wb are weights; the value of weights are set as ωa = wb = 1, which represents that the branch i-side and j-side are of equal importance [16]. k(a) integrates the distance protection on both sides, and the maximum disturbance degree on both sides of the remaining branch is defined as the severity index of the initial fault.(2) Weighted Power Flow Entropy Index. Starting from the level of branch anti-impact ability and system power flow transfer, the branch load rate is used to weight the power flow entropy index. The weighted power flow entropy index considers the influence of branch load rate and the total number of branches on the system operation state, which overcomes the problem that the traditional model only considers the balance of transfer power flow distribution. Define the weighted power flow entropy Va after branch La is cut off as [17](13)Va=∑b=1lαabβablnβablnl,where βab is the power flow impact rate on branch Lb after branch La is cut off; αab is the load rate of branch Lb after branch La is cut off.(3) Comprehensive Vulnerability Index. To more comprehensively and reasonably evaluate the vulnerability of initial fault, this paper further considers the system efficiency index S(a) [18], branch dielectric B(a) [19], and voltage offset index U(a) [20]. After normalizing and weighting each index, the comprehensive vulnerability index F(a) is defined as(14)Fa=w1ka′+w2Va′+w3Sa′+w4Ba′+w5Ua′,where ka′, Va′, Sa′, Ba′, and Ua′ are the normalized values of each index, respectively. w1,w2,w3,w4,andw5 are the weights of each index, and the weight is 1. The selection of each index weight in the practical application can be determined according to field experience. The comprehensive vulnerability index focuses on identifying the branch that is prone to overload in the process of fault propagation and has a serious impact on the topology and operation state of the power grid after removal. The greater the value of F(a), the more likely it is to cause cascading faults after fault and the greater the impact on the stable operation of the power grid.(4) Vulnerability Assessment Methods. According to the above analysis, this paper takes each vulnerability index as the character input and plans to use a fuzzy C-means algorithm (FCM) to effectively divide the branches. FCM algorithm is a clustering algorithm based on classification. It obtains the membership degree of each sample to all class centers by optimizing the objective function, determining the category of samples, to realize the automatic classification of sample data [21]. The objective function I and constraints are as follows:(15)Z=∑O=1we∑v=1veavtObtO−cv2,∑O=1we∑v=1veavtOb=1,where to is the characteristic input quantity composed of each evaluation index; cv is the cluster center of class v; avto is the membership function of the O-th sample data for class v; b is the exponential factor of membership avto; we is the number of samples; ve is the number of sample classifications. The value of the classification number is determined by referring to the classification entropy proposed in the document [22]. The classification entropy H is defined as(16)H=−1we∑O=1we∑v=1veavtOlnavtO,where the definition of each variable is the same as that in formula (15), and the number with the smallest classification entropy is the best classification number.The process of selecting the initial fault based on the FCM algorithm is as follows:(1)
Initialize the membership function and clustering center, and determine the algorithm parameters and feature input. The feature input includes the normalized values of the above six evaluation indexes.(2)
The Euclidean distance from each sample to the cluster centertO−cv and the classification entropy H of the sample are calculated.(3)
Determine the sample classification numberve and calculate the objective function Z. When ΔZ is less than the given threshold, the calculation terminates and the clustering result is output. Otherwise, return to step 2 to continue the calculation.In this way, the FCM algorithm can synthesize various vulnerability indexes, preliminarily screen out more important initial faults, and select these initial faults into the expected initial fault setG1.


#### 4.1.2. Initial Faults Are Further Selected according to the Safety Margin

According to the selected expected initial fault setG1, the safety margin of the power grid corresponding to each initial fault is calculated separately, whereby the safety level of the power grid for cascading faults after each initial fault is cut can be observed. Further, the initial fault with a low safety margin value is selected into the key initial fault set G2, which can be expressed as(17)G2=La|Ra<ε,where ε is the threshold value for screening and ε is a positive number, which can be given according to experience; Ra is the power grid safety margin corresponding to any initial fault branch La. Any initial fault branch can be regarded as a key initial fault as long as it satisfies formula (17).


### 4.2. Preventive Control Model for Different Initial Fault Scenarios

Due to the large scale of the actual power grid, a large number of branches, it leads to a significant reduction in the computational speed of the algorithm. In this paper, the key initial fault setG2 is selected based on the power grid vulnerability and safety margin, and preventive control is taken only for the initial faults in G2, which overcomes the problem of excessive computation due to previous control strategies that do not screen the initial faults.For the key initial fault setG2, this paper mainly considers the preventive control strategy in two cases: (1) the safety margin and economic cost after the cutoff of each key initial fault are calculated according to the basic model of preventive control respectively; (2) the second case is to transform the basic model of preventive control into an extended model, which applies to all key initial faults in G2, its objective function remains unchanged, and the corresponding constraints change from a single initial fault to multiple initial faults. The extended model can be shown as follows:(18)F1=minR1,R2,…,ReF2=maxF1,F3=minCme0x=0,me1x=0ge0x≥0,ge1x≥0de0x=0,de1x=0he0x≥0,he1x≥0I≥0,e=1,2,…,E,where e is the e-th initial fault; E is the number of initial faults; the corresponding constraints are expanded into E sets, and each set of constraints corresponds to one initial fault; the objective function of the inner model is the minimum value of the power grid safety margin F1 under multiple initial faults, while the objective function of the outer model is the same as the basic preventive control model. After the preventive control is carried out, it is necessary to ensure that the nodal injection power of the power grid simultaneously satisfies the condition of high safety margin and low economic cost in case of fault of any branch in the key initial fault set G2.


## 4.1. Selection of Initial Fault

### 4.1.1. Initial Fault Preselection Based on Vulnerability

(1) Severity Index. According to the previous analysis, without considering the uncertainty of impedance relay operation and other blocking conditions, the severity of the initial fault can be effectively evaluated from the branch level by the remaining system branches subjected to the initial fault removal. Let the number of the initial faults in the power grid be a. The severity index shown in formula (12) can be obtained after further processing of formula (1) as(12)ka=wamaxIib.setIib.m+wbmaxIjb.setIjb.m,where wa and wb are weights; the value of weights are set as ωa = wb = 1, which represents that the branch i-side and j-side are of equal importance [16]. k(a) integrates the distance protection on both sides, and the maximum disturbance degree on both sides of the remaining branch is defined as the severity index of the initial fault.(2) Weighted Power Flow Entropy Index. Starting from the level of branch anti-impact ability and system power flow transfer, the branch load rate is used to weight the power flow entropy index. The weighted power flow entropy index considers the influence of branch load rate and the total number of branches on the system operation state, which overcomes the problem that the traditional model only considers the balance of transfer power flow distribution. Define the weighted power flow entropy Va after branch La is cut off as [17](13)Va=∑b=1lαabβablnβablnl,where βab is the power flow impact rate on branch Lb after branch La is cut off; αab is the load rate of branch Lb after branch La is cut off.(3) Comprehensive Vulnerability Index. To more comprehensively and reasonably evaluate the vulnerability of initial fault, this paper further considers the system efficiency index S(a) [18], branch dielectric B(a) [19], and voltage offset index U(a) [20]. After normalizing and weighting each index, the comprehensive vulnerability index F(a) is defined as(14)Fa=w1ka′+w2Va′+w3Sa′+w4Ba′+w5Ua′,where ka′, Va′, Sa′, Ba′, and Ua′ are the normalized values of each index, respectively. w1,w2,w3,w4,andw5 are the weights of each index, and the weight is 1. The selection of each index weight in the practical application can be determined according to field experience. The comprehensive vulnerability index focuses on identifying the branch that is prone to overload in the process of fault propagation and has a serious impact on the topology and operation state of the power grid after removal. The greater the value of F(a), the more likely it is to cause cascading faults after fault and the greater the impact on the stable operation of the power grid.(4) Vulnerability Assessment Methods. According to the above analysis, this paper takes each vulnerability index as the character input and plans to use a fuzzy C-means algorithm (FCM) to effectively divide the branches. FCM algorithm is a clustering algorithm based on classification. It obtains the membership degree of each sample to all class centers by optimizing the objective function, determining the category of samples, to realize the automatic classification of sample data [21]. The objective function I and constraints are as follows:(15)Z=∑O=1we∑v=1veavtObtO−cv2,∑O=1we∑v=1veavtOb=1,where to is the characteristic input quantity composed of each evaluation index; cv is the cluster center of class v; avto is the membership function of the O-th sample data for class v; b is the exponential factor of membership avto; we is the number of samples; ve is the number of sample classifications. The value of the classification number is determined by referring to the classification entropy proposed in the document [22]. The classification entropy H is defined as(16)H=−1we∑O=1we∑v=1veavtOlnavtO,where the definition of each variable is the same as that in formula (15), and the number with the smallest classification entropy is the best classification number.The process of selecting the initial fault based on the FCM algorithm is as follows:(1)
Initialize the membership function and clustering center, and determine the algorithm parameters and feature input. The feature input includes the normalized values of the above six evaluation indexes.(2)
The Euclidean distance from each sample to the cluster centertO−cv and the classification entropy H of the sample are calculated.(3)
Determine the sample classification numberve and calculate the objective function Z. When ΔZ is less than the given threshold, the calculation terminates and the clustering result is output. Otherwise, return to step 2 to continue the calculation.In this way, the FCM algorithm can synthesize various vulnerability indexes, preliminarily screen out more important initial faults, and select these initial faults into the expected initial fault setG1.


### 4.1.2. Initial Faults Are Further Selected according to the Safety Margin

According to the selected expected initial fault setG1, the safety margin of the power grid corresponding to each initial fault is calculated separately, whereby the safety level of the power grid for cascading faults after each initial fault is cut can be observed. Further, the initial fault with a low safety margin value is selected into the key initial fault set G2, which can be expressed as(17)G2=La|Ra<ε,where ε is the threshold value for screening and ε is a positive number, which can be given according to experience; Ra is the power grid safety margin corresponding to any initial fault branch La. Any initial fault branch can be regarded as a key initial fault as long as it satisfies formula (17).


## 4.1.1. Initial Fault Preselection Based on Vulnerability

(1) Severity Index. According to the previous analysis, without considering the uncertainty of impedance relay operation and other blocking conditions, the severity of the initial fault can be effectively evaluated from the branch level by the remaining system branches subjected to the initial fault removal. Let the number of the initial faults in the power grid be a. The severity index shown in formula (12) can be obtained after further processing of formula (1) as(12)ka=wamaxIib.setIib.m+wbmaxIjb.setIjb.m,where wa and wb are weights; the value of weights are set as ωa = wb = 1, which represents that the branch i-side and j-side are of equal importance [16]. k(a) integrates the distance protection on both sides, and the maximum disturbance degree on both sides of the remaining branch is defined as the severity index of the initial fault.(2) Weighted Power Flow Entropy Index. Starting from the level of branch anti-impact ability and system power flow transfer, the branch load rate is used to weight the power flow entropy index. The weighted power flow entropy index considers the influence of branch load rate and the total number of branches on the system operation state, which overcomes the problem that the traditional model only considers the balance of transfer power flow distribution. Define the weighted power flow entropy Va after branch La is cut off as [17](13)Va=∑b=1lαabβablnβablnl,where βab is the power flow impact rate on branch Lb after branch La is cut off; αab is the load rate of branch Lb after branch La is cut off.(3) Comprehensive Vulnerability Index. To more comprehensively and reasonably evaluate the vulnerability of initial fault, this paper further considers the system efficiency index S(a) [18], branch dielectric B(a) [19], and voltage offset index U(a) [20]. After normalizing and weighting each index, the comprehensive vulnerability index F(a) is defined as(14)Fa=w1ka′+w2Va′+w3Sa′+w4Ba′+w5Ua′,where ka′, Va′, Sa′, Ba′, and Ua′ are the normalized values of each index, respectively. w1,w2,w3,w4,andw5 are the weights of each index, and the weight is 1. The selection of each index weight in the practical application can be determined according to field experience. The comprehensive vulnerability index focuses on identifying the branch that is prone to overload in the process of fault propagation and has a serious impact on the topology and operation state of the power grid after removal. The greater the value of F(a), the more likely it is to cause cascading faults after fault and the greater the impact on the stable operation of the power grid.(4) Vulnerability Assessment Methods. According to the above analysis, this paper takes each vulnerability index as the character input and plans to use a fuzzy C-means algorithm (FCM) to effectively divide the branches. FCM algorithm is a clustering algorithm based on classification. It obtains the membership degree of each sample to all class centers by optimizing the objective function, determining the category of samples, to realize the automatic classification of sample data [21]. The objective function I and constraints are as follows:(15)Z=∑O=1we∑v=1veavtObtO−cv2,∑O=1we∑v=1veavtOb=1,where to is the characteristic input quantity composed of each evaluation index; cv is the cluster center of class v; avto is the membership function of the O-th sample data for class v; b is the exponential factor of membership avto; we is the number of samples; ve is the number of sample classifications. The value of the classification number is determined by referring to the classification entropy proposed in the document [22]. The classification entropy H is defined as(16)H=−1we∑O=1we∑v=1veavtOlnavtO,where the definition of each variable is the same as that in formula (15), and the number with the smallest classification entropy is the best classification number.The process of selecting the initial fault based on the FCM algorithm is as follows:(1)
Initialize the membership function and clustering center, and determine the algorithm parameters and feature input. The feature input includes the normalized values of the above six evaluation indexes.(2)
The Euclidean distance from each sample to the cluster centertO−cv and the classification entropy H of the sample are calculated.(3)
Determine the sample classification numberve and calculate the objective function Z. When ΔZ is less than the given threshold, the calculation terminates and the clustering result is output. Otherwise, return to step 2 to continue the calculation.In this way, the FCM algorithm can synthesize various vulnerability indexes, preliminarily screen out more important initial faults, and select these initial faults into the expected initial fault setG1.


## 4.1.2. Initial Faults Are Further Selected according to the Safety Margin

According to the selected expected initial fault setG1, the safety margin of the power grid corresponding to each initial fault is calculated separately, whereby the safety level of the power grid for cascading faults after each initial fault is cut can be observed. Further, the initial fault with a low safety margin value is selected into the key initial fault set G2, which can be expressed as(17)G2=La|Ra<ε,where ε is the threshold value for screening and ε is a positive number, which can be given according to experience; Ra is the power grid safety margin corresponding to any initial fault branch La. Any initial fault branch can be regarded as a key initial fault as long as it satisfies formula (17).


## 4.2. Preventive Control Model for Different Initial Fault Scenarios

Due to the large scale of the actual power grid, a large number of branches, it leads to a significant reduction in the computational speed of the algorithm. In this paper, the key initial fault setG2 is selected based on the power grid vulnerability and safety margin, and preventive control is taken only for the initial faults in G2, which overcomes the problem of excessive computation due to previous control strategies that do not screen the initial faults.For the key initial fault setG2, this paper mainly considers the preventive control strategy in two cases: (1) the safety margin and economic cost after the cutoff of each key initial fault are calculated according to the basic model of preventive control respectively; (2) the second case is to transform the basic model of preventive control into an extended model, which applies to all key initial faults in G2, its objective function remains unchanged, and the corresponding constraints change from a single initial fault to multiple initial faults. The extended model can be shown as follows:(18)F1=minR1,R2,…,ReF2=maxF1,F3=minCme0x=0,me1x=0ge0x≥0,ge1x≥0de0x=0,de1x=0he0x≥0,he1x≥0I≥0,e=1,2,…,E,where e is the e-th initial fault; E is the number of initial faults; the corresponding constraints are expanded into E sets, and each set of constraints corresponds to one initial fault; the objective function of the inner model is the minimum value of the power grid safety margin F1 under multiple initial faults, while the objective function of the outer model is the same as the basic preventive control model. After the preventive control is carried out, it is necessary to ensure that the nodal injection power of the power grid simultaneously satisfies the condition of high safety margin and low economic cost in case of fault of any branch in the key initial fault set G2.


## 5. Solution of the Model

In this paper, we solve the preventive control model, the inner model uses particle swarm optimization (PSO) [23], and the outer model uses a multiobjective particle swarm optimization based on adaptive grid algorithm (AGA-MOPSO) [24]. In the solution process, the key variables and results between the inner and outer layers will be transferred to each other as needed, and the specific transfer method is as follows: the inner layer model passes the power grid safety margin to the outer layer model, and the outer layer model optimizes the generator nodal injection power taking into account the power grid safety margin and economic cost and returns the optimized nodal injection power to the inner layer model; then the inner layer model searches for the cascading fault again based on this result. Then, the inner model searches for the critical state of the cascading faults again based on the results to obtain the new power grid safety margin; finally, through the optimization iterations of the inner and outer layers, the optimal generator output is solved to coordinate safety and economy [25].


### 5.1. Particle Update Criterion

In PSO and AGA-MOPSO algorithms, the position and rate of each particle are updated according to the following formula [24]:(19)yit+1=yit+vit+1,vit+1=wvit+c1r1yp.i−yit+c2r2yg−yit,where yi(t+1) and vit+1 are the position and rate of particle i at the (t+1)st time, respectively, where the particle position yi of the inner model corresponds to the nodal injection power vector in the critical state of cascading fault. The particle position yi of the outer model corresponds to the nodal injection power vector of the current state of the power grid; yp.i is the best individual position of particle i; yg is the best position of the whole particle swarm; w is the inertia coefficient; c1 and c2 are acceleration factors; r1 and r2 are random numbers within [0,1].


### 5.2. Fitness Function

In this paper, the penalty function is introduced to further deal with the constraint problem in the model [26]. The specific process of this method is as follows: in the inner model, I=0 in formula (7) is written as zx=0, and all power flow equality constraints and inequality constraints in formula (8) are written as mx=0 and gx≥0, respectively. In the outer model, I>0 in formula (10) is written as px>0, and all power flow equality constraints and inequality constraints in formula (11) are written as dx=0 and hx≥0, respectively. Add the penalty factor, and use the objective function and constraint function to jointly construct the fitness function as follows:(20)fa=F1+M1∑emin0,gex2+M2∑emin0,mex2+M3min0,zx2fb=1F2+M1∑emin0,hex2+M2∑emin0,dex2+M3min0,px2fc=F3+M1∑emin0,hex2+M2∑emin0,dex2+M3min0,px2,where fa is the fitness function of the inner model; fb and fc are the fitness functions of the outer model; F1 is the objective function of the inner model; F2 and F3 are the objective functions of the outer model; M1, M2, and M3 are penalty factors; gex, me(x), de(x), and he(x) are some components of gx, m(x), d(x), and h(x), respectively.


### 5.3. Improvement of the Solution Algorithm

#### 5.3.1. Power Flow Transfer Method

When solving the power flow, the Newton Raphson method is proposed to solve the power flow equation in the form ofm1x=0, to obtain more accurate calculation results. For the power flow equation in the form of m0x=0, considering that each particle is generated based on S0, and the solution of the power flow equation in the form of m0x=0 under S0 is generally known, we should try to use the power flow solution corresponding to S0 to obtain the power flow equation in the form of m0x=0 contained by any particle, to improve the calculation efficiency. Therefore, the following calculation ideas are designed in this paper.For ease of explanation, the power flow equation in the form ofm0x=0 is defined as [27](21)f0x=S,where S is the nodal injection power of the power grid, f0 is the mapping relation corresponding to the power flow equation of the form m0x=0, and x is the state variable corresponding to S. Based on the previous description, it is clear that there exists f0x0=S0, where x0 is the state variable under the current initial operating conditions of the power grid, and setting S−S0=ΔS, x−x0=Δx, it can be defined as(22)ΔxA0=ΔS,where A0=f′x|x=x0 is the Jacobian matrix. It can be further defined as(23)x=x0+A0−1ΔS.Thus, for each particle, according to the value of its nodal injection powerS, the solution of the tidal equation of the form of m0x=0, for each particle in the iterative process, can be found using formula (23).


#### 5.3.2. Adaptive Mutation Operation

To prevent the algorithm from falling into local optimum, this paper draws on the idea of literature [28] and introduces an adaptive mutation operation to mutate the current particles with a certain probability to achieve the goal of breaking up the overly aggregated particles, to search a new direction again and increase the chance of finding the global optimal solution. The specific mutation method is as follows.Calculate the rate of mutationPa:(24)Pa=1−iter−1Maxt−11/mt,where iter is the number of current iterations; Maxt is the maximum number of iterations; mt is the coefficient of mutation, and its value is taken as 0.1.Then calculate the mutation interval:(25)minvmin,yit+1−Δy,maxvmax,yit+1+Δy,where vmin and vmax are particle velocity limits; Δy can be defined as(26)Δy=Pa×vmax−vmin.Finally, execute the mutation operation, and the new variableY can be defined as(27)Y=unifrndminvmin,yit+1−Δy,maxvmax,yit+1+Δy,where unifrnd is a random array function uniformly distributed on the mutation interval.


#### 5.3.3. Choosing the Compromise Optimal Solution

After obtaining a set of Pareto optimal front (POF), the decision-maker can choose the corresponding compromise optimal solution according to the demand. There are many methods to choose the compromise optimal solution, such as the decision-making method based on membership function [29], the decision-making method based on distance evaluation index [30], and the comprehensive decision-making method based on the combination of FCM clustering and GRP [31, 32]. Because the decision-makers’ judgment on the weight of each objective function is uncertain, they often select the compromise optimal solution that can meet each objective function to the greatest extent in the Pareto solution set according to the membership function. In this paper, the membership degree of the Pareto solution is calculated by using fuzzy set theory. The membership degree λrm of Pareto solution in a one-dimensional objective function can be expressed as(28)λrm=1,frm≤fm.minfm.max−frmfm.max−fm.min,fm.min<frm<fm.max0frm≥fm.max,where frm is the m-dimensional objective function value of the r-th optimal solution, m=1,2,…,N; fm.min and fm.max are the minimum and maximum values of the m-dimensional objective function value, respectively.The standardized membership degree calculation formula of each Pareto solution is(29)λr=∑m=1Nλrm∑r=1M∑m=1Nλrm,where M is the number of Pareto solutions, λr is the comprehensive membership of the r-th Pareto solution, and the solution with the largest λr is selected as the optimal solution.


### 5.4. Algorithm Flow

Based on the PSO algorithm and AGA-MOPSO algorithm, this paper solves the double-layer optimized preventive control model. The overall algorithm flow is shown in Figure2:Figure 2
Algorithm flowchart.


## 5.1. Particle Update Criterion

In PSO and AGA-MOPSO algorithms, the position and rate of each particle are updated according to the following formula [24]:(19)yit+1=yit+vit+1,vit+1=wvit+c1r1yp.i−yit+c2r2yg−yit,where yi(t+1) and vit+1 are the position and rate of particle i at the (t+1)st time, respectively, where the particle position yi of the inner model corresponds to the nodal injection power vector in the critical state of cascading fault. The particle position yi of the outer model corresponds to the nodal injection power vector of the current state of the power grid; yp.i is the best individual position of particle i; yg is the best position of the whole particle swarm; w is the inertia coefficient; c1 and c2 are acceleration factors; r1 and r2 are random numbers within [0,1].


## 5.2. Fitness Function

In this paper, the penalty function is introduced to further deal with the constraint problem in the model [26]. The specific process of this method is as follows: in the inner model, I=0 in formula (7) is written as zx=0, and all power flow equality constraints and inequality constraints in formula (8) are written as mx=0 and gx≥0, respectively. In the outer model, I>0 in formula (10) is written as px>0, and all power flow equality constraints and inequality constraints in formula (11) are written as dx=0 and hx≥0, respectively. Add the penalty factor, and use the objective function and constraint function to jointly construct the fitness function as follows:(20)fa=F1+M1∑emin0,gex2+M2∑emin0,mex2+M3min0,zx2fb=1F2+M1∑emin0,hex2+M2∑emin0,dex2+M3min0,px2fc=F3+M1∑emin0,hex2+M2∑emin0,dex2+M3min0,px2,where fa is the fitness function of the inner model; fb and fc are the fitness functions of the outer model; F1 is the objective function of the inner model; F2 and F3 are the objective functions of the outer model; M1, M2, and M3 are penalty factors; gex, me(x), de(x), and he(x) are some components of gx, m(x), d(x), and h(x), respectively.


## 5.3. Improvement of the Solution Algorithm

### 5.3.1. Power Flow Transfer Method

When solving the power flow, the Newton Raphson method is proposed to solve the power flow equation in the form ofm1x=0, to obtain more accurate calculation results. For the power flow equation in the form of m0x=0, considering that each particle is generated based on S0, and the solution of the power flow equation in the form of m0x=0 under S0 is generally known, we should try to use the power flow solution corresponding to S0 to obtain the power flow equation in the form of m0x=0 contained by any particle, to improve the calculation efficiency. Therefore, the following calculation ideas are designed in this paper.For ease of explanation, the power flow equation in the form ofm0x=0 is defined as [27](21)f0x=S,where S is the nodal injection power of the power grid, f0 is the mapping relation corresponding to the power flow equation of the form m0x=0, and x is the state variable corresponding to S. Based on the previous description, it is clear that there exists f0x0=S0, where x0 is the state variable under the current initial operating conditions of the power grid, and setting S−S0=ΔS, x−x0=Δx, it can be defined as(22)ΔxA0=ΔS,where A0=f′x|x=x0 is the Jacobian matrix. It can be further defined as(23)x=x0+A0−1ΔS.Thus, for each particle, according to the value of its nodal injection powerS, the solution of the tidal equation of the form of m0x=0, for each particle in the iterative process, can be found using formula (23).


### 5.3.2. Adaptive Mutation Operation

To prevent the algorithm from falling into local optimum, this paper draws on the idea of literature [28] and introduces an adaptive mutation operation to mutate the current particles with a certain probability to achieve the goal of breaking up the overly aggregated particles, to search a new direction again and increase the chance of finding the global optimal solution. The specific mutation method is as follows.Calculate the rate of mutationPa:(24)Pa=1−iter−1Maxt−11/mt,where iter is the number of current iterations; Maxt is the maximum number of iterations; mt is the coefficient of mutation, and its value is taken as 0.1.Then calculate the mutation interval:(25)minvmin,yit+1−Δy,maxvmax,yit+1+Δy,where vmin and vmax are particle velocity limits; Δy can be defined as(26)Δy=Pa×vmax−vmin.Finally, execute the mutation operation, and the new variableY can be defined as(27)Y=unifrndminvmin,yit+1−Δy,maxvmax,yit+1+Δy,where unifrnd is a random array function uniformly distributed on the mutation interval.


### 5.3.3. Choosing the Compromise Optimal Solution

After obtaining a set of Pareto optimal front (POF), the decision-maker can choose the corresponding compromise optimal solution according to the demand. There are many methods to choose the compromise optimal solution, such as the decision-making method based on membership function [29], the decision-making method based on distance evaluation index [30], and the comprehensive decision-making method based on the combination of FCM clustering and GRP [31, 32]. Because the decision-makers’ judgment on the weight of each objective function is uncertain, they often select the compromise optimal solution that can meet each objective function to the greatest extent in the Pareto solution set according to the membership function. In this paper, the membership degree of the Pareto solution is calculated by using fuzzy set theory. The membership degree λrm of Pareto solution in a one-dimensional objective function can be expressed as(28)λrm=1,frm≤fm.minfm.max−frmfm.max−fm.min,fm.min<frm<fm.max0frm≥fm.max,where frm is the m-dimensional objective function value of the r-th optimal solution, m=1,2,…,N; fm.min and fm.max are the minimum and maximum values of the m-dimensional objective function value, respectively.The standardized membership degree calculation formula of each Pareto solution is(29)λr=∑m=1Nλrm∑r=1M∑m=1Nλrm,where M is the number of Pareto solutions, λr is the comprehensive membership of the r-th Pareto solution, and the solution with the largest λr is selected as the optimal solution.


## 5.3.1. Power Flow Transfer Method

When solving the power flow, the Newton Raphson method is proposed to solve the power flow equation in the form ofm1x=0, to obtain more accurate calculation results. For the power flow equation in the form of m0x=0, considering that each particle is generated based on S0, and the solution of the power flow equation in the form of m0x=0 under S0 is generally known, we should try to use the power flow solution corresponding to S0 to obtain the power flow equation in the form of m0x=0 contained by any particle, to improve the calculation efficiency. Therefore, the following calculation ideas are designed in this paper.For ease of explanation, the power flow equation in the form ofm0x=0 is defined as [27](21)f0x=S,where S is the nodal injection power of the power grid, f0 is the mapping relation corresponding to the power flow equation of the form m0x=0, and x is the state variable corresponding to S. Based on the previous description, it is clear that there exists f0x0=S0, where x0 is the state variable under the current initial operating conditions of the power grid, and setting S−S0=ΔS, x−x0=Δx, it can be defined as(22)ΔxA0=ΔS,where A0=f′x|x=x0 is the Jacobian matrix. It can be further defined as(23)x=x0+A0−1ΔS.Thus, for each particle, according to the value of its nodal injection powerS, the solution of the tidal equation of the form of m0x=0, for each particle in the iterative process, can be found using formula (23).


## 5.3.2. Adaptive Mutation Operation

To prevent the algorithm from falling into local optimum, this paper draws on the idea of literature [28] and introduces an adaptive mutation operation to mutate the current particles with a certain probability to achieve the goal of breaking up the overly aggregated particles, to search a new direction again and increase the chance of finding the global optimal solution. The specific mutation method is as follows.Calculate the rate of mutationPa:(24)Pa=1−iter−1Maxt−11/mt,where iter is the number of current iterations; Maxt is the maximum number of iterations; mt is the coefficient of mutation, and its value is taken as 0.1.Then calculate the mutation interval:(25)minvmin,yit+1−Δy,maxvmax,yit+1+Δy,where vmin and vmax are particle velocity limits; Δy can be defined as(26)Δy=Pa×vmax−vmin.Finally, execute the mutation operation, and the new variableY can be defined as(27)Y=unifrndminvmin,yit+1−Δy,maxvmax,yit+1+Δy,where unifrnd is a random array function uniformly distributed on the mutation interval.


## 5.3.3. Choosing the Compromise Optimal Solution

After obtaining a set of Pareto optimal front (POF), the decision-maker can choose the corresponding compromise optimal solution according to the demand. There are many methods to choose the compromise optimal solution, such as the decision-making method based on membership function [29], the decision-making method based on distance evaluation index [30], and the comprehensive decision-making method based on the combination of FCM clustering and GRP [31, 32]. Because the decision-makers’ judgment on the weight of each objective function is uncertain, they often select the compromise optimal solution that can meet each objective function to the greatest extent in the Pareto solution set according to the membership function. In this paper, the membership degree of the Pareto solution is calculated by using fuzzy set theory. The membership degree λrm of Pareto solution in a one-dimensional objective function can be expressed as(28)λrm=1,frm≤fm.minfm.max−frmfm.max−fm.min,fm.min<frm<fm.max0frm≥fm.max,where frm is the m-dimensional objective function value of the r-th optimal solution, m=1,2,…,N; fm.min and fm.max are the minimum and maximum values of the m-dimensional objective function value, respectively.The standardized membership degree calculation formula of each Pareto solution is(29)λr=∑m=1Nλrm∑r=1M∑m=1Nλrm,where M is the number of Pareto solutions, λr is the comprehensive membership of the r-th Pareto solution, and the solution with the largest λr is selected as the optimal solution.


## 5.4. Algorithm Flow

Based on the PSO algorithm and AGA-MOPSO algorithm, this paper solves the double-layer optimized preventive control model. The overall algorithm flow is shown in Figure2:Figure 2
Algorithm flowchart.


## 6. Example Analysis

The test example takes the IEEE39 node system as an example, and the system wiring is shown in Figure3.Figure 3
Diagram of IEEE39 node system.According to the system wiring diagram shown in Figure3, the analysis program is written in the Matlab environment according to the above algorithm. In the following calculation examples, the cost coefficients of each generator of the system can be given concerning literature [15], and the nodal injection power S0 corresponding to the parameters of each component of the system and the current power grid operation state can be given concerning literature [33]. To facilitate the calculation, this paper uses the unit value for representation. For the specific protection configuration data on the IEEE39 node system, this paper uses virtual data and sets the current protection setting value of each branch as three times the current of the branch in the base-state power flow state.


### 6.1. Initial Fault Analysis

Since the power flow does not converge immediately after the initial fault branchL16-19 (the branch between node 16 and node 19, etc.) is cut off, it indicates that there may be voltage instability in the system after the fault, which should be studied in combination with the analysis and control of voltage stability, so branch L16-19 is not considered as the initial fault for the time being.For each index weight in formula (14), this paper determines the weight according to the entropy weight method [34], and the final calculation results are w1 = 0.3566, w2 = 0.0561, w3 = 0.1354, w4 = 0.2047, and w5 = 0.2472. Select other nongenerator outlet branches n-1 except branch L16-19 as the initial fault, normalize each index according to the FCM algorithm process, and then cluster it. Finally, the initial fault is divided into two categories. The clustering centers are shown in Table 1.Table 1
Cluster center of the example.
CategoryF(a)ka′Va′Sa′Ba′Ua′A0.510.390.630.620.580.53B0.160.060.370.290.170.16From Table1, it can be seen that the indexes of class A are larger, so the branch of class A is selected into the expected initial fault set G1. There are 7 branches in class A. The vulnerability indexes are sorted by size, and the results are shown in Table 2.Table 2
Vulnerable branches of IEEE39 node system.
BranchF(a)Rankingka′RankingVa′RankingSa′RankingBa′RankingUa′RankingL15-160.7510.4630.7330.7240.96211L26-270.6720.5810.5050.6450.3240.325L21-220.6530.395110.317110.843L2-30.5540.4720.7420.8130.7830.325L2-250.3750.4540.5740.8820.1950.097L14-150.3360.0170.2770.6060.0970.882L16-170.3270.1160.406110.1260.434As can be seen from Table2, there is a big difference between the ranking of each vulnerability index individually and the ranking of the integrated vulnerability index F(a), and the branch with the top ranking of the integrated vulnerability index often contains both strong and weak vulnerability indexes, which indicates that the use of a single index for branch vulnerability evaluation is not comprehensive enough and may ignore some weak links, while the integrated vulnerability index takes into account both the network topology and the branch after opening and breaking. The comprehensive vulnerability index takes into account both the network topology and the operation status of the branch after opening and closing, which has certain superiority. To verify the rationality and effectiveness of the initial fault screening, ranking results of the comprehensive vulnerability index of the branch are compared with literature [35, 36], and the comparison results show that more than half of the identified vulnerable branches are the same, which indicates that the expected initial fault set G1 obtained in this paper is more reasonable and comprehensive.To further select the initial faults, the safety margin after cutoff of each initial fault in the expected initial fault setG1 is calculated separately, and then the conditions are judged to be satisfied according to formula (17), where the value of ε is taken as 1. The final calculation results are shown in Table 3.Table 3
Safety margin for each initial failure.
Initial faultL14-15L16-17L2-3L26-27L2-25L15-16L21-22Safety margin3.352.631.741.391.290.270.19It can be seen from Table3 that the safety margins of the initial fault branches L15-16 and L21-22 after cutoff are less than 1, which meets the conditions for screening key initial faults. Therefore, branches L15-16 and L21-22 are selected into the key initial fault set G2, and preventive control will be taken later.


### 6.2. Analysis of Preventive Control Results

In this paper, we study preventive control, and we do not know which branch is going to have cascading faults, so we cannot know which initial fault to give the optimization result by. Therefore, it is necessary to make the preventive control method applicable to all key initial faults, and based on obtaining the key initial fault setG2, we consider the preventive control strategies under two different initial fault scenarios, solve the Pareto solutions for the initial fault branches L15-16 and L21-22 in G2 according to the basic preventive control model and the extended preventive control model, respectively, and record the generator node power corresponding to each Pareto optimal solution. Set PSO and AGA-MOPSO algorithm parameters: the maximum number of iterations is 200, the number of populations is 50, the acceleration factors c1 and c2 are 2, the inertia coefficient w is 0.5, the penalty factors M1 and M2 are 1000, and M3 is 10000. The Pareto optimal front (POF) obtained under different initial fault scenarios is as follows.Figures4 and 5 are the Pareto solutions calculated according to the basic preventive control model for the initial fault branches L15-16 and L21-22, respectively, and Figure 6 is the Pareto solutions obtained according to the extended preventive control model for the initial fault branches L15-16 and L21-22 at the same time. From the simulation results, it can be seen that the extended model established in this paper is suitable for all initial faults screened, can meet the conditions of high safety margin and low economic cost under any fault, and ensure that the controlled power grid is in a safer operation state. The Pareto solution obtained by the algorithm forms a POF curve. The POF curve has good uniform ductility. Each Pareto solution has a different economic cost and power grid safety margin. The POF curve is from left to right, and the safety margin and economic cost of the corresponding Pareto solution are from small to large. The membership degree of each Pareto solution is calculated through the membership function, and the respective compromise optimal solution is found according to the value of the comprehensive membership degree. The results of the optimal solution are shown in Table 4:Figure 4
POF from initial failureL15-16.Figure 5
POF from initial failureL21-22.Figure 6
POF from initial failuresL15-16 and L21-22 at the same time.Table 4
Preventive control optimization results.
Generator no.L15-16L21-22L15-16; L21-22G30881.93899.31477.33G32634.51581.78166.43G33168.36389.95639.32G34250.43558.72506.11G35714.05562.74503.09G36218.75268.34286.76G37550.46624.76962.73G38885.92631.39354.63G39860.50631.85832.60Safety margin9.178.738.55Generation cost1101.121085.341051.16The method given in this paper is the simulation and tested on several other node systems as well, and the results are satisfactory, which shows that the method in this paper is suitable for various systems. From the above calculation example, it is clear that the safety and economy of the power grid are significantly improved after preventive control for all possible key initial fault sets. At this point, the power grid enters a safer operating stateS0, and when the power grid nodal injection power in an arbitrary variation within a hypersphere with S0 as the center of the sphere and the safety margin value R as the radius, the power grid is safe for different key initial fault impacts for the triggering of cascading faults; that is, no cascading faults are triggered. In conclusion, compared with the existing methods, the method in this paper can facilitate the observation of the safety of the power grid and the reasonable variation range of the nodal injection power of the power grid and minimizes the generation cost during the power grid adjustment process, providing the operators with the optimal generator output strategy that balances safety and economy.


## 6.1. Initial Fault Analysis

Since the power flow does not converge immediately after the initial fault branchL16-19 (the branch between node 16 and node 19, etc.) is cut off, it indicates that there may be voltage instability in the system after the fault, which should be studied in combination with the analysis and control of voltage stability, so branch L16-19 is not considered as the initial fault for the time being.For each index weight in formula (14), this paper determines the weight according to the entropy weight method [34], and the final calculation results are w1 = 0.3566, w2 = 0.0561, w3 = 0.1354, w4 = 0.2047, and w5 = 0.2472. Select other nongenerator outlet branches n-1 except branch L16-19 as the initial fault, normalize each index according to the FCM algorithm process, and then cluster it. Finally, the initial fault is divided into two categories. The clustering centers are shown in Table 1.Table 1
Cluster center of the example.
CategoryF(a)ka′Va′Sa′Ba′Ua′A0.510.390.630.620.580.53B0.160.060.370.290.170.16From Table1, it can be seen that the indexes of class A are larger, so the branch of class A is selected into the expected initial fault set G1. There are 7 branches in class A. The vulnerability indexes are sorted by size, and the results are shown in Table 2.Table 2
Vulnerable branches of IEEE39 node system.
BranchF(a)Rankingka′RankingVa′RankingSa′RankingBa′RankingUa′RankingL15-160.7510.4630.7330.7240.96211L26-270.6720.5810.5050.6450.3240.325L21-220.6530.395110.317110.843L2-30.5540.4720.7420.8130.7830.325L2-250.3750.4540.5740.8820.1950.097L14-150.3360.0170.2770.6060.0970.882L16-170.3270.1160.406110.1260.434As can be seen from Table2, there is a big difference between the ranking of each vulnerability index individually and the ranking of the integrated vulnerability index F(a), and the branch with the top ranking of the integrated vulnerability index often contains both strong and weak vulnerability indexes, which indicates that the use of a single index for branch vulnerability evaluation is not comprehensive enough and may ignore some weak links, while the integrated vulnerability index takes into account both the network topology and the branch after opening and breaking. The comprehensive vulnerability index takes into account both the network topology and the operation status of the branch after opening and closing, which has certain superiority. To verify the rationality and effectiveness of the initial fault screening, ranking results of the comprehensive vulnerability index of the branch are compared with literature [35, 36], and the comparison results show that more than half of the identified vulnerable branches are the same, which indicates that the expected initial fault set G1 obtained in this paper is more reasonable and comprehensive.To further select the initial faults, the safety margin after cutoff of each initial fault in the expected initial fault setG1 is calculated separately, and then the conditions are judged to be satisfied according to formula (17), where the value of ε is taken as 1. The final calculation results are shown in Table 3.Table 3
Safety margin for each initial failure.
Initial faultL14-15L16-17L2-3L26-27L2-25L15-16L21-22Safety margin3.352.631.741.391.290.270.19It can be seen from Table3 that the safety margins of the initial fault branches L15-16 and L21-22 after cutoff are less than 1, which meets the conditions for screening key initial faults. Therefore, branches L15-16 and L21-22 are selected into the key initial fault set G2, and preventive control will be taken later.


## 6.2. Analysis of Preventive Control Results

In this paper, we study preventive control, and we do not know which branch is going to have cascading faults, so we cannot know which initial fault to give the optimization result by. Therefore, it is necessary to make the preventive control method applicable to all key initial faults, and based on obtaining the key initial fault setG2, we consider the preventive control strategies under two different initial fault scenarios, solve the Pareto solutions for the initial fault branches L15-16 and L21-22 in G2 according to the basic preventive control model and the extended preventive control model, respectively, and record the generator node power corresponding to each Pareto optimal solution. Set PSO and AGA-MOPSO algorithm parameters: the maximum number of iterations is 200, the number of populations is 50, the acceleration factors c1 and c2 are 2, the inertia coefficient w is 0.5, the penalty factors M1 and M2 are 1000, and M3 is 10000. The Pareto optimal front (POF) obtained under different initial fault scenarios is as follows.Figures4 and 5 are the Pareto solutions calculated according to the basic preventive control model for the initial fault branches L15-16 and L21-22, respectively, and Figure 6 is the Pareto solutions obtained according to the extended preventive control model for the initial fault branches L15-16 and L21-22 at the same time. From the simulation results, it can be seen that the extended model established in this paper is suitable for all initial faults screened, can meet the conditions of high safety margin and low economic cost under any fault, and ensure that the controlled power grid is in a safer operation state. The Pareto solution obtained by the algorithm forms a POF curve. The POF curve has good uniform ductility. Each Pareto solution has a different economic cost and power grid safety margin. The POF curve is from left to right, and the safety margin and economic cost of the corresponding Pareto solution are from small to large. The membership degree of each Pareto solution is calculated through the membership function, and the respective compromise optimal solution is found according to the value of the comprehensive membership degree. The results of the optimal solution are shown in Table 4:Figure 4
POF from initial failureL15-16.Figure 5
POF from initial failureL21-22.Figure 6
POF from initial failuresL15-16 and L21-22 at the same time.Table 4
Preventive control optimization results.
Generator no.L15-16L21-22L15-16; L21-22G30881.93899.31477.33G32634.51581.78166.43G33168.36389.95639.32G34250.43558.72506.11G35714.05562.74503.09G36218.75268.34286.76G37550.46624.76962.73G38885.92631.39354.63G39860.50631.85832.60Safety margin9.178.738.55Generation cost1101.121085.341051.16The method given in this paper is the simulation and tested on several other node systems as well, and the results are satisfactory, which shows that the method in this paper is suitable for various systems. From the above calculation example, it is clear that the safety and economy of the power grid are significantly improved after preventive control for all possible key initial fault sets. At this point, the power grid enters a safer operating stateS0, and when the power grid nodal injection power in an arbitrary variation within a hypersphere with S0 as the center of the sphere and the safety margin value R as the radius, the power grid is safe for different key initial fault impacts for the triggering of cascading faults; that is, no cascading faults are triggered. In conclusion, compared with the existing methods, the method in this paper can facilitate the observation of the safety of the power grid and the reasonable variation range of the nodal injection power of the power grid and minimizes the generation cost during the power grid adjustment process, providing the operators with the optimal generator output strategy that balances safety and economy.


## 7. Conclusion

Considering safety and economy, this paper presents a multiobjective preventive control model and specific solution process for preventing cascading faults of the power grid and verifies it with a numerical example. The main conclusions are as follows:(1)
The preventive control model established in this paper is essentially a double-layer optimization model with dynamic interaction and interaction between internal and external layers. The optimal generator output strategy is obtained by solving the double-layer model so that the power grid can take into account the requirements of safety and economy at the same time.(2)
In this paper, the safety margin of the power grid is mainly measured by the distance between the current operating state and the critical state of the power grid, which is in the form of the distance between the nodal injection power vectors of the power grid in the two states. This form of safety margin representation can visually quantify and compare the safety level of the power grid for cascading faults before and after taking preventive measures.(3)
In this paper, for the selection of initial faults, the selection is mainly made considering the vulnerability indexes of the network topology and operation state as well as the safety margin, ensuring that the important initial faults are selected.In conclusion, the method proposed in this paper can effectively select the initial faults that are more important to the power grid and preventive control according to the idea of improving the safety margin of the power grid and ensuring the economic cost, providing a reference for further research.



---

*Source: 4554236-2021-10-11.xml*