A Novel Self-Healing Strategy for Distribution Network with Distributed Generators Considering Uncertain Power-Quality Constraints

: Self-healing of smart distribution networks with distributed generators, which are usually operated as independent islands after fault, can improve power-supply reliability. As a hot research topic, a self-healing scheme is usually treated as the output of a nonlinear optimizuoation model. However, existing strategies have two main shortcomings. The first, high-optimization dimension, results in low-optimization efficiency. The second, the effects of power-quality issues, which are more serious on islands and may further threaten the safe operation of islands, is usually neglected. To quickly obtain a reliable self-healing scheme, a novel strategy is proposed. As the first step, the distribution network after a fault occurrence can be divided into several trouble-free subnets via the connectivity analysis; each subnet is called an initial island. Further, for each initial island, a two-step optimization model of self-healing, which consists of load-shedding optimization and network reconfiguration optimization, is proposed to obtain the self-healing strategy with lower searching space as well as higher solving efficiency. In detail, in load-shedding optimization, by means of heuristic differential evolution algorithm, larger total recovery capacity is achieved by considering the droop characteristic of distributed generators (DGs) within the permissible change in frequency. In network-reconfiguration optimization, based on the improved hybrid particle swarm optimization algorithm, a comprehensive set of power-quality constraints, including constraint of change in frequency, uncertain constraints of node voltage total harmonic distortion (THD), and negative sequence components of DGs, is developed to guarantee the reliability of each island. To evaluate whether the constraints are satisfied during the optimization procedure, an improved flexible power-flow algorithm is developed to calculate the power flow of each island under change in frequency. Further, 2m+1-point estimate method is employed for uncertainty analyses of the distributions of harmonic and negative sequence components caused by the uncertainty of corresponding sources. Finally, via a 94-node practical distribution network, the effectiveness and advantages of the proposed strategy in safety, recovery capacity, and optimization efficiency are verified.


Introduction
As a core feature of smart distribution network (SDN), fast power-supply restoration after faults, which can reduce the adverse impacts of outage and abrupt faults, has important implications for both energy-supply companies and consumers. In recent years, technological advances and wide applications of distributed generators (DGs) provided powerful technical support for self-healing recovery of SDN.
1) Total capacity of recovered loads [1-9]: As the basic goal of self-healing recovery, de-energized loads must be restored as quickly as possible.
3) Number of switching operations [6,7,9]: As switching operations are time consuming and have operational cost, the number of switching operations shall be minimized to realize faster recovery and lower cost. 4) Power loss [8,9]: For economical operation, the smaller the power loss of a self-healing scheme, the better.
Good results have been achieved in References [1][2][3][4][5][6][7][8][9]. However, existing strategies have two main shortcomings. The first, high-optimization dimension, results in low-optimization efficiency. The second, the effects of power quality issues, which are more serious in islanding and may further threaten the safe operation of islands, are usually neglected.
With wide application of power electronic devices and access to various nonlinear loads, a series of power-quality issues, such as harmonic pollution and negative sequence component injection, widely exist. As the capacity of each operating island is smaller compared with normally operated SDN, these issues will be more serious and may affect the safe operation of islands. For synchronous generator-based DG, serious negative sequence currents may result in being out of service [10,11]. For inverter-interfaced DG (IIDG), serious harmonic or negative sequence components may result in output oscillation or even in becoming out of service [12][13][14]. Serious harmonic component can lead to misoperation or could even cause harmonic sensitive load and equipment to stop working [15,16]. Therefore, power-quality issues may result in power unbalance of island and further lead to recovery failure. Thus, to prevent potential recovery failures caused by power-quality issues, it is necessary to uniformly consider the influences of harmonic and negative sequence components in the analysis of self-healing of SDN. Furthermore, as the injections of harmonic and negative sequence sources are uncertain and fluctuant, it is necessary to employ the uncertainty analysis method to obtain an optimal self-healing scheme in accordance with the distribution network actual operation.
In this paper, to realize faster, stronger, and safer recovery, a new self-healing recovery strategy is proposed. To improve optimization efficiency, a decoupled two-step optimization model of load-shedding and network-reconfiguration is established. To enlarge recovery capacity, the droop characteristic of DGs is considered in load-shedding within the permissible change in frequency. To prevent potential recovery failure caused by power-quality issues, a comprehensive set of power-quality restraints, including frequency constraint, uncertain constraints of harmonic, and negative sequence components, is established in the network reconfiguration optimization model. This paper is organized as follows: Section 2 puts forward the framework of the proposed self-healing recovery strategy. Detailed instructions of the load-shedding optimization model are illustrated in Section 3. The network reconfiguration optimization model, considering uncertain power quality constraints, is minutely explained in Section 4. In Section 5, the solved algorithms of the proposed strategy are illustrated. Effectiveness and superiorities of the proposed strategy are further verified in Section 6. Finally, conclusions are drawn in Section 7.

Framework of the Proposed Self-Healing Recovery Strategy
In this section, the framework of the self-healing recovery strategy proposed in this paper, which is as shown in Figure 1, is illustrated. Once a fault of SDN is detected by distribution automation system (DAS), the self-healing recovery strategy is put forward as follows.
1) The distribution network after a fault occurrence can be divided into several trouble-free subnets, each of which is called an initial island. At the first step, the initial island division is carried out based on the connectivity analysis of the distribution network and each initial island is the minimum application object of optimization model of self-healing, which consists of load-shedding optimization and network reconfiguration optimization; the following two points are important for island division.
 Keep the normal operation of widely distributed high-priority loads as normal as possible.
 The island with larger capacity and more adjustable equipment can achieve better power quality indexes. Based on the above understandings, the following initial island division strategy is adopted with the consideration of tie-switches. a) If the network is connected after the fault, the network is operated as a whole island. b) If the network is no longer connected after the fault but the connectivity can be recovered by the tie-switches, the network is still operated as a whole island. c) If the network is no longer connected after the fault and the connectivity cannot be recovered by the tie-switches, the network shall be divided into independent small islands.
In the proposed strategy, the method in Reference [17] is used to judge the network connectivity.
The algorithm to judge the network connectivity is explained as follows. c) Referring to the short connection in the circuit, the network diagram can be simplified in a similar way. If two nodes are directly connected, they can be merged into one node to form a new network graph.
For example, if node i is connected to node j, which means Z[i,j] = 1, the simplification of the network diagram can be operated as follows. In detail, for a row operation, apply the "OR" operation to the corresponding element of row i and row j, and the result is placed in row j. For a column operation, apply the "OR" operation to the corresponding element of column i and column j, and the result can be placed in column j. After the operations, the diagonal element can be set to zero. Since node i has been merged into node j, row i and column i are removed. The above simplification decreases the dimension of adjacency matrix Z from n to n-1. d) Repeat the simplification. The network diagram is not connected if there is a row in which all the elements are zeros during the process. Otherwise, the network diagram is connected, and when the dimension decreases to 2, all nondiagonal elements of adjacency matrix Z are 1.
Through the above simplification, each connected part of the network diagram is considered as an initial island.
After the initial island division, the proposed load-shedding optimization and network reconstruction optimization are applied for each initial island. In practical application, this procedure is in parallel for each island, thus the increase of the number of initial islands may lead to the growth of the total amount of computational burden, but it will not lead to the delay of strategy implementation for each initial island by taking advantage of parallel computation. In consequence, there is no specific restrictions on the number of islands.
2) For each initial island, the load-shedding optimization and network reconfiguration optimization are carried out as follows, respectively. a) In the load-shedding optimization model, on one hand, the droop characteristic of DGs is considered to enlarge the sum recovery capacity within the allowable change in frequency. On the other hand, load priority is introduced to guarantee the normal operation of the loads that are more important. Concrete illustrations of the load-shedding optimization model are proposed in Section 3. b) In the network reconfiguration optimization model, the comprehensive set of power-quality constraints is considered so the obtained recovery scheme of the island can achieve good power-quality indexes. By this way, the threats caused by power-quality issues, such as high node voltage total harmonic distortion and excessive negative sequence current, can be prevented to guarantee the safe operation of the obtained island. Concrete illustrations of the load-shedding optimization model are proposed in Section 4.

Optimization Variables
In the load-shedding optimization model, the variables are the shedding state of each load, which is treated as 0-1 type variable. Concretely, xik = 1 means that the load of node i on the kth island is shed, while xik = 0 means the load is maintained.

Objective Function
In this section, the aim is to minimize the sum capacity of shed loads with the consideration of load priority. Thus, the objective function of load-shedding is as follows.
where Ls, m, k, PLi, and λi are the minimum of load-shedding of all the nodes, the total number of loads, the number of the island to be studied, active power, and weighting coefficient of the load of node i on island k, respectively.
In detail, λi is set as follows in this paper according to a lot of tests. 1) For important loads (L2), λi = 100.

Constraints
1) Active power balance constraint considering the droop characteristic of DGs.
The active power shortage is supposed to be smaller than the result of DGs' generated active power minus PLos0.
where PLto_re represents the total active power capacity of the recovered loads, PLos_0 is the initial active power loss of the initial island to be studied, and PDGj_f is the output active power of jth DG under frequency f if the DG can participate frequency modulation, which can be calculated according to the droop characteristic as follows [18].
where PDGj_r is the rated output active power of jth DG under rated frequency fr and Kpfj is its active power-frequency adjustment coefficient. According to Reference [19], fr and the minimum value of f are set as 50 Hz and 49.5 Hz, respectively. Thus, to maximize the total recovery capacity, PDGj_max, which is the maximum value of PDGj_f, can be calculated by setting the value of f as 49.5 Hz as follows.  (2) can be modified as the equation below.
The reactive power shortage is supposed to be smaller than the result of DGs' generated reactive power plus Qcmax.
where xik is defined in Section 3.1, QLi is the reactive power of the ith load, and Qcmax represents the maximum total reactive compensation capacity of the initial island. QDGj_max is the maximum allowable reactive power of the jth DG, which is decided by the apparent power constraint as follow: where SDGj_max and PDGj_max are the maximum apparent power and active power of the jth DG.

Network-Reconfiguration Optimization Model
As the second step, network reconfiguration is based on the optimal load-shedding scheme obtained in Section 3.

Optimization Variables
In network-reconfiguration optimization, the state of each switch and the specific compensation capacity of each reactive power compensation device are treated as 0-1-type and integer-type variables, respectively.

Objective Function
Considering rapidity, reliability, and economy, network-reconfiguration optimization aims at minimizing the number of switching operations. Concrete objective function can be expressed as follows: min F s t abs s s abs t t where k1 and k2 are the amount of section switches and tie-switches of the island, respectively; sa and sa_ini are the optimized and initial states of the ath section switch, respectively; and tb and tb_ini are the optimized and initial states of the bth tie-switch.

Constraints of Network-Reconfiguration
1) Regular constraints a) DG apparent power constraint.
where SDGj is the apparent power of the jth DG. b) Branch distribution power constraint.
where Sc_max is the maximum allowable value of Sc, which is the distribution power of the cth branch. c) Node voltage amplitude constraint.
where Vi_min and Vi_max represent the lower and upper limits of Vi, which is the voltage amplitude of node i.

2) Comprehensive power quality constraints set
a) Uncertain node voltage total harmonic distortion (THD) constraint.
where THD-uli is the upper limit of the distribution range of the THD index of node i. THDmax is the permissible maximum value of THD, which is defined as 4% for the distribution network of 10 kV according to Reference [20]. b) Uncertain negative sequence current constraint for synchronous generator-based DG where j j j

Cune-ul Ine-ul In =
where Ine-ulj is the upper bound of distribution range of the actual negative sequence current of the jth synchronous generator-based DG while Inj is the rated positive sequence current. Cumax is the permissible maximum value, which is set as 10% referring to Reference [21]. c) Uncertain negative sequence voltage constraint for IIDG where Vone-ulj is the upper bound of the distribution range of the actual negative sequence voltage of the kth IIDG while Vk is the actual positive sequence voltage of the kth IIDG. Vomax is the permissible maximum value, which is set as 2% referring to Reference [22]. d) Constraint of island's change in frequency According to Reference [19], Δfmax, which is the maximum allowable value of change in frequency Δf, is 0.5Hz.

Solving Algorithms
In this paper, the heuristic differential evolution algorithm and improved hybrid particle swarm optimization algorithm [23] are used to seek the optimal load-shedding scheme and network-reconfiguration scheme, respectively. The optimal solutions of the two models constitute the whole recovery scheme.
To evaluate whether the constraints are satisfied during the optimization procedure, an improved flexible power-flow algorithm for distribution network including second-order items is developed to calculate the power flow of each island under change in frequency. Also, 2m+1 point estimate methods (PEM), which has higher computational efficiency while the accuracy is almost the same as the Monte Carlo method, is chosen to carry out uncertainty analysis.
Detailed illustrations of the two algorithms are as follows.

Improved Flexible Power-Flow Algorithm for Distribution Network including Second-Order Items
In Reference [24], in order to obtain the generic initial value of the distribution network power flow, the equations are transformed to quadratic matrix equations in which the second-order items and droop characteristic of DGs and load are considered. It is expressed as the following equation.
In the equation B is the imaginary component of the network node admittance matrix and G is the real component. Concrete derivations and statements of the above equation can be found in Reference [24]. For Equation (14), its solving efficiency is restricted by the fsolve function of MATLAB.
To overcome the shortage, the solution issue of Equation (14) is proposed to be converted into a nonlinear constrained optimization problem and solved via the sequential quadratic programming (SQP) method in this paper.
Actually, Equation (14) contains 2n equations. If one moves SL to the left side of the equation, Equation (14) can be written as follows.
Generally, Equation (15) can be written as the uniform expression below.
As SQP has the characteristics of good global convergence as well as excellent local super-linear convergence [25]; better performances of convergence and computing efficiency can be achieved by the improved solving method.

2m+1 Point Estimate Method
The 2m+1 PEM, in which m represents the number of input random variables, can be used to calculate the statistical moments of a random quantity, which is a function of one or several random variables [26]. Compared with the two-point estimate method, the accuracy of 2m+1 PEM is higher. To get the statistical moments of output random variables, the function has to be calculated 2m+1 times [27].
Let Xi denote a random variable and Y = f(X1, X2, …, Xi) be the function of Xi. The procedures for finding the moments of the output Y by 2m+1 PEM are as follows.
1) Find probability concentrations of each input random variable.
, , where xi,k is the concentration of Xi; μi and σi are the mean and standard deviation of Xi; and θi,k is the location-factor, which can be obtained as follow: where αi,r is the rth standardized central moment of Xi and can be obtained by the following Equation.
, , where E[(xi−μi) r ] is the rth central moment of Xi.
2) Determine the 2m+1 estimation points. According to Equations (17)- (19), each random variable can be estimated by 2 points; the 2m estimation points can be determined as follows: The special addition of the 2m estimation points is the point of xμ, which is constituted by the means of all the input random variables: 3) According to Equations (19)- (20), the weighting factors are as follows: ( 1) 1, 2 ( )

Uncertainty Analyses of Harmonic and Negative Sequence Components Based on the 2m+1 PEM
The load-shedding optimization model and network reconfiguration optimization model make strategies according to the current electrical quantity and network topology collected by the distribution network automation system. According to the description in chapter 2 of the original manuscript, this information is deterministic; thus, the optimization models also adopt deterministic models and power-flow algorithm for the fundamental component.
In the network reconfiguration model, power-quality indexes are required to meet the constraints not only at the current moment but also in the following period of time. This is also one of the main contributions of this paper. In the following period of time, the power-quality index is uncertain, so the harmonic power flow and the unbalanced power flow of the fundamental component are obtained based on the 2m+1 PEM method, in which the estimation accuracy is similar to the classical Monte Carlo method but computational burden is much less to accommodate the need for online application.
In this paper, uncertainty analyses of harmonic and negative sequence components are independent of each other. The harmonic power-flow calculation adopts the method used in References [10,28,29]. The negative sequence component calculation method used in References [30][31][32] is introduced to analyze the ranges of system's negative sequence currents.
Due to the similarity of the two analyses, harmonic uncertainty analysis is chosen as the example to give detailed descriptions. Specific steps are as follows.
1) The sequences of the measured harmonic current values of harmonic sources are presented to calculate the mean and variance of each variable.

4)
Multiply the results of step 3 and their corresponding weighting factors together to get the expectation μ(Yk) and standard deviation σ(Yk) of each output variable Yk according to Equations (24) and (25).

5) Obtain the distribution range of Yk as [ ( ) 2 ( ), ( ) 2 ( )]
as the distribution range of k Y are as follows.  For normal distributions, the probability of the observed value within The event with probability less than 5% generally can be treated as a small probability event.
If Yk is node voltage total harmonic distortion, negative sequence current of synchronous generator-based DG, and negative sequence voltage of IIDG, the upper bound ( ) 2 ( ) be THD-ul, Cune-ul, and Vone-ul, respectively.

Introductions of a Practical Complex Distribution Network and Scenarios
To fully demonstrate the effectiveness and advantages of the proposed strategy, lots of tests verified by simulations via ETAP power system simulation software have been done via several examples with different operation status. Due to limited space, a practical complex distribution network in China, which is shown as Figure 2, is used to give detailed descriptions. The practical complex distribution network contains 94 nodes, 93 branches, 4 tie lines (10-87, 25-37, 54-55, and 76-85), and 5 DGs. In the case of normal power supply, total capacity of all the loads is 13,379.29 kW + 4519.12 kvar. Two reactive power compensations are installed on node 11 and node 59, respectively. A fluctuant harmonic source and a negative sequence current source, which fit normal distribution, are set at node 51 and node 61, respectively. Limited by space, more detailed instructions of the distribution network are listed in Table A1 and Table A2 of Appendix A.
For the distribution network, the following benchmark scenario is selected as the typical situation to give detailed descriptions: Benchmark scenario: disconnected with network and DG1, DG2, and DG3 are out of the running.
It is worth mentioning that the methods proposed in this paper can be applied to the distribution network benchmarks with different scales and topologies as well as to different scenarios of the same benchmarks under electrical faults at different locations.

Initial Results of Island Division
Based on the island division strategy proposed in Section 1, under the above scenario, the distribution network can be divided into islands.
For the benchmark scenario, the connectivity of the test system is not changed; therefore, the distribution network does not need to be divided into multiple islands.

Comparisons of the Proposed Strategy and the Conventional Strategy
To illustrate the superiorities of the proposed self-healing strategy, the benchmark scenario is studied to provide detailed comparison. Scheme 1 is the self-healing scheme obtained by the proposed strategy. Scheme 2 is the self-healing scheme obtained by the conventional strategy, in which the power quality constraint set and droop characteristic of DGs are not considered. The values of power-quality indexes of Scheme 2 are calculated. Concretely, the two schemes are minutely listed in the following Table 1. From Table 1, the following advantages of the proposed strategy can be summarized. 1) Total recovery capacity is larger. Specifically, total load removal of Scheme 1 is 9830.21 kW, while the value of Scheme 2 is 10,429.45 kW. That is, the total recovery capacity has risen nearly 20.31% from 2949.84 kW to 3549.08 kW, which means the proposed strategy can realize stronger recovery.
2) The obtained recovery scheme has better power-quality indexes and higher reliability. working. That is, the reliability of Scheme 2 will not be guarantee.
To sum up, the proposed strategy can realize the larger total recovery capacity and higher reliability by considering the DG frequency adjustment characteristic and the comprehensive set of power-quality constraints.

Necessity of Harmonic and Negative Uncertainty Analyses
To illustrate the necessity of uncertainty analysis, Scheme 2 is further analyzed by contrasting the THD-ul and THD-me indexes and the Cune-ul and Cune-me indexes which are minutely listed in Table 2 and Table 3, respectively. THD-me and Cune-me are the midpoints of the distribution ranges of THD index and negative sequence current, respectively.

Cune-ul
Cune-me DG5 10.56% 8.72% Through the above two tables, it can be illustrated that the values of THD-me and Cu-me are more conservative. Thus, recovery failure may be caused, as the scheme with actually exceeded power-quality indexes is treated as a feasible scheme. Reasons are as follows.
1) As all the nodes in THD-me are smaller than 4%, Scheme 2 is a feasible scheme in the view of THD-me. However, meanwhile, the values of THD-ul index of all the nodes are within (4.41%, 4.54%) which are larger than 4%.
2) Cune-me index of DG5 is within 10%, so Scheme 2 can be treated as a feasible one from the perspective of Cune-me. However, the Cune-ul index of DG5 is clearly exceeded.
Thus, in accordance with operation practices of power grids, it is necessary to adopt the uncertainty analysis method to deal with the uncertainty of harmonic and negative sequence components.

Solution Efficiency Comparison
Moreover, a comparison is conducted to illustrate the advantage of the proposed two-step optimization model. The computing platform is a PC with Intel i7-4790@3.60GHZ CPU and 8G RAM (Intel, Santa Clara, CA, USA), and the software environment is Windows 10 Professional and MATLAB R2010b (MathWorks, Natick, MA, USA).
Concretely, the proposed two-step optimization model and the conventional one-layer optimization model are used to find the recovery for benchmark scenario. Specific information of the two recovery schemes is shown in Table 4 below. Scheme 3 is obtained by the conventional one-layer optimization model. From Table 4, the following can be seen: 1) The power-quality indexes of the two schemes satisfy the constraints, and the numbers of switching actions are both 2.
2) Theoretically, the conventional one-layer optimization model can obtain better optimal recovery scheme with enough time. However, as the searching space of the one-layer optimization model is too large, the optimal solution is hard to be obtained. Until the set maximum iteration value, the total recovery capacity of Scheme 3 is 3469.96 kW, which is smaller than 3549.08 kW.
3) Elapsed time for solving two-step optimization model is 0.52 minute, which is merely 5.22% of 9.96 minute. Therefore, via the proposed two-step optimization model, the time consumed can be reduced significantly, which can achieve the aim of rapid obtainment of a good recovery scheme.

Conclusions
In this paper, a novel self-healing strategy for distribution network with distributed generators considering uncertain power-quality constraints is proposed. It can realize faster, stronger, and safer recovery of an SDN. Based on case studies via a practical large-scale distribution network, the following characteristics of the proposed strategy are verified.
1) Higher optimization efficiency is realized by establishing the two-step optimization model, which consists of load-shedding optimization and network reconfiguration optimization.
2) Larger recovery capacity achieved by fully using the droop characteristic of DGs within the permissible change in frequency in load-shedding optimization.
3) Safer operation of each island is achieved by improving power-quality indexes in network reconfiguration optimization to prevent recovery failure caused by power-quality issues. 4) Moreover, uncertainty analysis of power-quality indexes based on 2m+1 PEM is more in accordance with operation practices of SDN.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Concrete information including the installation location, rated active power, active power-frequency adjustment coefficient, and reserve capacity of the 5 DGs (DG1~DG5) are shown in the following Table A1.
Among them, DG1 and DG4 are micro-turbine power stations with droop control; DG2 and DG3 are doubly fed induction generator and photovoltaic power station with maximum power point tracking, respectively; DG5 is small synchronous generator; and DG1, DG2, DG3, and DG4 are IIDGs.