# An Innovatory Method Based on Continuation Power Flow to Analyze Power System Voltage Stability with Distributed Generation Penetration

Authors: Le Van Dai; Ngo Minh Khoa; Le Cao Quyen
Journal: Complexity
 (2020)

---

## Abstract

With the penetration of distributed generation (DG) units, the power systems will face insecurity problems and voltage stability issues. This paper proposes an innovatory method by modifying the conventional continuation power flow (CCPF) method. The proposed method is realized on two prediction and correction steps to find successive load flow solutions according to a specific load scenario. Firstly, the tangent predictor is proposed to estimate the next predicted solution from two previous corrected solutions. And then, the corrector step is proposed to determine the next corrected solution on the exact solution. This corrected solution is constrained to lie in the hyperplane running through the predicted solution orthogonal to the line from the two previous corrected solutions. Besides, once the convergence criterion is reached, the procedure for cutting the step length control down to a smaller one is proposed to be implemented. The effectiveness of the proposed method is verified via numerical simulations on three standard test systems, namely, IEEE 14-bus, 57-bus, and 118-bus, and compared to the CCPF method.


---

## Body

## 1. Introduction

### 1.1. Motivation

In recent years, distributed generations (DGs) as renewable energy sources connected to power system networks have been growing rapidly. Generally, the units of DG have a small size compared to existing central power plants; they have little effect if their penetration level is low at about 1% to 5%. However, if the penetration level is about 20% to 30%, the impact of DG units will be deep [1]. With a rapid increase in the DGs penetration level, the electric power system stability can change and it is affecting the voltage profile, reliability, protection, stability, power quality, and power flow [2, 3]; among them, the voltage stability problem is an important aspect.The voltage stability is an important problem that has been reported for several years, and a lot of papers and reports of high quality on this subject have been published by well-known publishers. According to the CIGRE Study Committee and IEEE Power System Dynamic Performance Committee that have set up a joint task force, the topic of voltage stability is classified into two types that are the large and small disturbances considering the short and long terms [4]. The impact of DGs on the voltage stability in the long and short terms is examined in [5–9], respectively. In addition, the DGs’ optimal allocation, considering uncertainty in load demand and growth, is studied in [10–12].The voltage stability is the efficiency of an electric system that voltage at all buses can maintain steady acceptable values under normal operating conditions for all the power system depending on a disturbance [13]. The voltage instability occurs when load dynamics struggle to restore power consumption beyond the possibilities of the generation system. So it could lead to the trip of loads and/or transmission lines and the loss of synchronism of some generators [14].During the last periods, numerous blackouts of the power system in the world relating to the voltage instability problems have been occurring in various countries such as Japan, India, the USA, France, Sweden, and Germany. It is worth noting that the blackouts have occurred in the Greek island of Kefalonia on January 24, 2006 [15]; the aftereffect is the cut of electricity energy to be approximately 3,000 MWh. The 61,800 MW capacity lost lasted around 4 to 7 days due to occurring blackouts throughout Ohio, Michigan, Pennsylvania, New York, Vermont, Massachusetts, Connecticut, New Jersey, and Ontario, Canada [16], affecting more than 50 million people. A blackout, affecting the normal life and work of 670 million people, happened on the power system of northern India on July 30 and 31, 2012 [17]. Another blackout, affecting more than 10 million customers, happened on the Brazilian power system, leading to heavy economic losses [18, 19].


### 1.2. Literature Review

The blackouts in the power system can cause great damage to socioeconomic and brought inconvenience to people. Therefore, the need for a trustworthy method to assess voltage stability in a power system is developed. To achieve this, a number of methods have been used in the literature for voltage stability analysis. For example, the active power-voltage (P-V) curve developed in [20] and the reactive power-voltage (Q-V) curve introduced in [13]. The shortcomings of these two methods are no information about key contributing factors for voltage stability problems and high computational time for a large power system network. To overcome these shortcomings, the continuation power flow (CPF) method has been proposed in [21]. Actually, CPF is also time-consuming for a large power system network, even though it has a lot of strong points based on a prediction-correction technique to find the continuous power flow solutions that start from the base-load condition to the limit of steady-state voltage stability [22].The aforesaid methods based on the power flow as the modal analysis have been introduced in [23]. This proposed method depends on the power flow Jacobian matrix, in which the active power is kept constant and the reactive power margin and voltage instability contribute as factors, resulting in the reduced Jacobian matrix of the system for computation. So that it is more suitable for the large power system. The methods have used voltage stability indices including the fast voltage stability index [24], line stability index [25], voltage reactive power index [26], and L index [27]. These methods are helpful to determine the proximity of a given operating point to a voltage collapse point, but the computational steps and time are heavy. Therefore, their main disadvantage is time-consuming to calculate for a large power system network and they can only predict the voltage collapse point on the power network [28, 29].


### 1.3. Contribution and Study Organization

This paper proposes a feasible method, namely, the modified continuation power flow (MCPF) method, to analyze the voltage stability problem to the electric power system considering the penetration of DG units. This proposed method is realized on the base of the conventional continuation power flow (CCPF) method considering the prediction and correction steps to find successive load flow solutions according to the load scenarios. A next real solution is determined based on two previously defined adjacent solutions using the predictor technique based on the arc-length parameterization and the correction technique based on the local parameterization. The process of next solutions on the whole P-V curve will continue until determining the critical solution corresponding to the load parameter equal to the critical. Besides, in exceptional cases, if the real solution is not identified by the correction step, the previous prediction step should be adjusted until the real solution is determined.To verify the effectiveness of the proposed method, three standard test systems, namely, IEEE 14-bus, 57-bus, and 118-bus, are proposed to perform the different test cases. Based on the obtained results, this paper has the main contributions as follows:(i)
Can exactly determine the voltage collapse point for analyzing voltage stability for large power systems and especially for the large power systems with DGs penetration.(ii)
Cut the number of predictor-corrector steps and reduce computational time.The remaining parts of the paper can be divided into four sections as follows. The voltage stability problem is formulated in Section2 to provide inputs to develop the analysis method. Section 3 recalls the CCPF method and develops a novel method for analyzing the voltage stability issue. The effectiveness of the proposed method is verified in Section 4. Finally, the conclusions are given in Section 5.


## 1.1. Motivation

In recent years, distributed generations (DGs) as renewable energy sources connected to power system networks have been growing rapidly. Generally, the units of DG have a small size compared to existing central power plants; they have little effect if their penetration level is low at about 1% to 5%. However, if the penetration level is about 20% to 30%, the impact of DG units will be deep [1]. With a rapid increase in the DGs penetration level, the electric power system stability can change and it is affecting the voltage profile, reliability, protection, stability, power quality, and power flow [2, 3]; among them, the voltage stability problem is an important aspect.The voltage stability is an important problem that has been reported for several years, and a lot of papers and reports of high quality on this subject have been published by well-known publishers. According to the CIGRE Study Committee and IEEE Power System Dynamic Performance Committee that have set up a joint task force, the topic of voltage stability is classified into two types that are the large and small disturbances considering the short and long terms [4]. The impact of DGs on the voltage stability in the long and short terms is examined in [5–9], respectively. In addition, the DGs’ optimal allocation, considering uncertainty in load demand and growth, is studied in [10–12].The voltage stability is the efficiency of an electric system that voltage at all buses can maintain steady acceptable values under normal operating conditions for all the power system depending on a disturbance [13]. The voltage instability occurs when load dynamics struggle to restore power consumption beyond the possibilities of the generation system. So it could lead to the trip of loads and/or transmission lines and the loss of synchronism of some generators [14].During the last periods, numerous blackouts of the power system in the world relating to the voltage instability problems have been occurring in various countries such as Japan, India, the USA, France, Sweden, and Germany. It is worth noting that the blackouts have occurred in the Greek island of Kefalonia on January 24, 2006 [15]; the aftereffect is the cut of electricity energy to be approximately 3,000 MWh. The 61,800 MW capacity lost lasted around 4 to 7 days due to occurring blackouts throughout Ohio, Michigan, Pennsylvania, New York, Vermont, Massachusetts, Connecticut, New Jersey, and Ontario, Canada [16], affecting more than 50 million people. A blackout, affecting the normal life and work of 670 million people, happened on the power system of northern India on July 30 and 31, 2012 [17]. Another blackout, affecting more than 10 million customers, happened on the Brazilian power system, leading to heavy economic losses [18, 19].


## 1.2. Literature Review

The blackouts in the power system can cause great damage to socioeconomic and brought inconvenience to people. Therefore, the need for a trustworthy method to assess voltage stability in a power system is developed. To achieve this, a number of methods have been used in the literature for voltage stability analysis. For example, the active power-voltage (P-V) curve developed in [20] and the reactive power-voltage (Q-V) curve introduced in [13]. The shortcomings of these two methods are no information about key contributing factors for voltage stability problems and high computational time for a large power system network. To overcome these shortcomings, the continuation power flow (CPF) method has been proposed in [21]. Actually, CPF is also time-consuming for a large power system network, even though it has a lot of strong points based on a prediction-correction technique to find the continuous power flow solutions that start from the base-load condition to the limit of steady-state voltage stability [22].The aforesaid methods based on the power flow as the modal analysis have been introduced in [23]. This proposed method depends on the power flow Jacobian matrix, in which the active power is kept constant and the reactive power margin and voltage instability contribute as factors, resulting in the reduced Jacobian matrix of the system for computation. So that it is more suitable for the large power system. The methods have used voltage stability indices including the fast voltage stability index [24], line stability index [25], voltage reactive power index [26], and L index [27]. These methods are helpful to determine the proximity of a given operating point to a voltage collapse point, but the computational steps and time are heavy. Therefore, their main disadvantage is time-consuming to calculate for a large power system network and they can only predict the voltage collapse point on the power network [28, 29].


## 1.3. Contribution and Study Organization

This paper proposes a feasible method, namely, the modified continuation power flow (MCPF) method, to analyze the voltage stability problem to the electric power system considering the penetration of DG units. This proposed method is realized on the base of the conventional continuation power flow (CCPF) method considering the prediction and correction steps to find successive load flow solutions according to the load scenarios. A next real solution is determined based on two previously defined adjacent solutions using the predictor technique based on the arc-length parameterization and the correction technique based on the local parameterization. The process of next solutions on the whole P-V curve will continue until determining the critical solution corresponding to the load parameter equal to the critical. Besides, in exceptional cases, if the real solution is not identified by the correction step, the previous prediction step should be adjusted until the real solution is determined.To verify the effectiveness of the proposed method, three standard test systems, namely, IEEE 14-bus, 57-bus, and 118-bus, are proposed to perform the different test cases. Based on the obtained results, this paper has the main contributions as follows:(i)
Can exactly determine the voltage collapse point for analyzing voltage stability for large power systems and especially for the large power systems with DGs penetration.(ii)
Cut the number of predictor-corrector steps and reduce computational time.The remaining parts of the paper can be divided into four sections as follows. The voltage stability problem is formulated in Section2 to provide inputs to develop the analysis method. Section 3 recalls the CCPF method and develops a novel method for analyzing the voltage stability issue. The effectiveness of the proposed method is verified in Section 4. Finally, the conclusions are given in Section 5.


## 2. Voltage Stability Problem

### 2.1. Problem Formulation

The power flow is an important problem in the field of power system engineering, where the voltage magnitudes and angles at buses are desired and the voltage magnitudes and power levels at other buses are known considering the available mode of the network configuration. The power flow calculation aims to determine voltage magnitudes and angles for a given generation, load, and grid condition. The starting step of resolving the power flow problem is to determine the unknown and known variables in the power system. Based on these variables, as known in the power systems, there exist three bus types that are the load bus (PQ bus), in which theP and Q variables are detailed and U and δ variables are necessary to be solved; the generation bus (PV bus), in which the P and U variables are detailed and the U and δ variables are necessary to be solved; and the slack bus (swing bus), in which the U and δ variables are detailed and the P and Q variables are necessary to be solved. It can be observed that, for the power system having n buses and g generators (including DG units), the number of unknowns which is 2n−1−g−1 can be calculated by using the reactive and active power balance equations. Therefore, the active and reactive power balance equations injected into the network at bus j can be defined as follows [30]:(1)Pj=PGeneration,j−PLoad,j,Qj=QGeneration,j−QLoad,j,where PGeneration,j and QGeneration,j are the generator’s active and reactive powers at bus j, respectively; PLoad,j and QLoad,j are the load’s active and reactive powers at bus j, respectively.The complex power at busj can be given by(2)Sj=Pj+jQj=UjIj∗,where Uj is the voltage at bus j and Ij is the current injected into bus j which can be written as follows:(3)Ij=∑k=1nYjkUk,whereUk is the voltage at bus k and equation (3) can be rewritten under the following relationship between the bus voltage and current as follows:(4)I1I2⋮In=Y11Y12⋯Y1nY21Y22⋯Y2n⋮⋮⋱⋮Yn1Yn2⋯YnnU1U2⋮Un.Substituting equation (3) into equation (2), one gets(5)Sj=Pj+jQj=Uj∑k=1nYjkUk∗,where Yjk is mutual admittance between buses j and k. The active and reactive power balance equation at bus j can be written as follows:(6)Pj=∑k=1nUjUkGjkcosδjk+Bjksinδjk,Qj=∑k=1nUjUkGjksinδjk−Bjkcosδjk,whereGjk and Bjk are, respectively, the real and imaginary part of the element in the bus admittance matrix corresponding to the jth row and kth column. The voltage phase angle between buses j and k is δjk=δj−δk, in which δj and δk are the voltage phase angles at buses j and k, respectively.


### 2.2. Analyses

The CPF method is a numerical technique used to overcome the problems of the Newton-Raphson technique such that it is more reliable for obtaining the solution curve, specifically to the ill-conditioned power flow equations, and faster via the effective predictor and corrector scheme, step length control, and parameterization.The CPF determines the consecutive load flow solutions under a load scenario using the prediction and correction steps. The prediction-correction procedure is illustrated in Figure1. It can be summarized as, from a previous known solution, the tangent predictor step is used to determine a next predictor solution under considering the increase of the loading factor and then the predictor step is applied to find the exact solution. And after that, based on the new tangent vector, a new prediction for a load increase is built. This process will be finished when the critical point is reached; in other words, the critical point is one where tangent vector is zero.Figure 1
The predictor-corrector procedure of continuation power flow.The CPF is performed based on the predictor-corrector procedure via the following system of nonlinear equation [22]:(7)Fx=0,where(8)x=δ,UT.The nonlinear power flow equation in equation (8) is augmented by the loading factor λ as follows:(9)Fx,λ=0,0≤λ≤λcritical,where U is the vector of bus voltage magnitudes and δ is the vector of bus voltage angles. The dimension of F is 2n+g, as mentioned, n denotes the number of PQ buses, and g denotes the number of PV buses.The power flow equation in equation (1) is rewritten, comprising a loading factor λ as follows:(10)Pj−λ+γkGeneration,iPGeneration_o,j+λPLoad_o,j=0,Qj−QGeneration,j+λQLoad_o,j=0,where PLoad_o,j and QLoad_o,j are, respectively, the original load demands at bus j, PGeneration_o,j is the original generation demands at bus j, kGeneration,i is the generation factor at bus j, and γ is the generator participation coefficient.From equation (9), the base solution corresponding to the loading factor equal to zero started to determine the next solutions corresponding to different load levels and the critical solution corresponding to the load parameter equal to the critical. This can be shown from the tracing of the PV curve based on the predictor-corrector scheme as shown in Figure 1. The predictor-corrector continuation process procedure can be realized in two steps.


#### 2.2.1. The First Step

This step realizes the predictor process to calculate the tangent vector and can be expressed in the following form:(11)t=dx,dλT.The tangent vectort with a corresponding dimension 2n+g+1 is determined from the following form:(12)Fx,Fλt=0,where Fx and Fλ are the partial derivative of equation (9) with respect to x and λ, respectively, and T denotes the transpose operation.From equation (12), the left side is a partial derivative matrix that this matrix is nothing but the original Jacobian augmented by one column [Fλ] of an unknown variable of load λ and a vector of differentials t is the differentials representing the tangent vector. Since an unknown variable of the load is added to the nonlinear system (9), this will lead the augmented Jacobian to be singular at the point of maximum possible system load. In order to solve this problem, the tangent vector should be determined based on the satisfied augment Jacobian as follows:(13)FxFλeit=0ti,where ei is a row vector having an appropriate dimension, in which the ith element is the only nonzero element.From equation (13), ti is a nonzero element of the tangent vector; this parameter can be either +1 or −1 depending on how the ith state variable is changing as a solution that is being traced. An estimate of the next solution can be determined by solving system (13) as follows:(14)xeλe=xiλi+σdxdλ,where the subscripts “i” and “e” denote the current and estimated solutions of the next step, respectively, and σ is a scalar defining the predictor step length control.


#### 2.2.2. The Second Step

After the prediction is finished, if the predictor solution may not be exactly on a desired solution curve, this step is realized to correct the predicted solution that is determined in equation (14) using a technique, namely, the local parameterization. Now, for the corrector step, system (9) is introduced by one equation unit with the purpose of determining the value of one of the bus voltage magnitude, bus voltage angle, and loading factor. Hence, the new system is expressed as follows:(15)Fx,λxi−μ=00,where xi is the ith state variable that is selected as the continuation parameter and μ is the predicted value of xi.Equation (15) shows that it is involved due to one additional equation and a state variable. In order to solve this difficulty, this paper reproposes Newton’s method for this second step as follows:(16)FxFλeiΔxΔλ=Fx,λxi−η.


## 2.1. Problem Formulation

The power flow is an important problem in the field of power system engineering, where the voltage magnitudes and angles at buses are desired and the voltage magnitudes and power levels at other buses are known considering the available mode of the network configuration. The power flow calculation aims to determine voltage magnitudes and angles for a given generation, load, and grid condition. The starting step of resolving the power flow problem is to determine the unknown and known variables in the power system. Based on these variables, as known in the power systems, there exist three bus types that are the load bus (PQ bus), in which theP and Q variables are detailed and U and δ variables are necessary to be solved; the generation bus (PV bus), in which the P and U variables are detailed and the U and δ variables are necessary to be solved; and the slack bus (swing bus), in which the U and δ variables are detailed and the P and Q variables are necessary to be solved. It can be observed that, for the power system having n buses and g generators (including DG units), the number of unknowns which is 2n−1−g−1 can be calculated by using the reactive and active power balance equations. Therefore, the active and reactive power balance equations injected into the network at bus j can be defined as follows [30]:(1)Pj=PGeneration,j−PLoad,j,Qj=QGeneration,j−QLoad,j,where PGeneration,j and QGeneration,j are the generator’s active and reactive powers at bus j, respectively; PLoad,j and QLoad,j are the load’s active and reactive powers at bus j, respectively.The complex power at busj can be given by(2)Sj=Pj+jQj=UjIj∗,where Uj is the voltage at bus j and Ij is the current injected into bus j which can be written as follows:(3)Ij=∑k=1nYjkUk,whereUk is the voltage at bus k and equation (3) can be rewritten under the following relationship between the bus voltage and current as follows:(4)I1I2⋮In=Y11Y12⋯Y1nY21Y22⋯Y2n⋮⋮⋱⋮Yn1Yn2⋯YnnU1U2⋮Un.Substituting equation (3) into equation (2), one gets(5)Sj=Pj+jQj=Uj∑k=1nYjkUk∗,where Yjk is mutual admittance between buses j and k. The active and reactive power balance equation at bus j can be written as follows:(6)Pj=∑k=1nUjUkGjkcosδjk+Bjksinδjk,Qj=∑k=1nUjUkGjksinδjk−Bjkcosδjk,whereGjk and Bjk are, respectively, the real and imaginary part of the element in the bus admittance matrix corresponding to the jth row and kth column. The voltage phase angle between buses j and k is δjk=δj−δk, in which δj and δk are the voltage phase angles at buses j and k, respectively.


## 2.2. Analyses

The CPF method is a numerical technique used to overcome the problems of the Newton-Raphson technique such that it is more reliable for obtaining the solution curve, specifically to the ill-conditioned power flow equations, and faster via the effective predictor and corrector scheme, step length control, and parameterization.The CPF determines the consecutive load flow solutions under a load scenario using the prediction and correction steps. The prediction-correction procedure is illustrated in Figure1. It can be summarized as, from a previous known solution, the tangent predictor step is used to determine a next predictor solution under considering the increase of the loading factor and then the predictor step is applied to find the exact solution. And after that, based on the new tangent vector, a new prediction for a load increase is built. This process will be finished when the critical point is reached; in other words, the critical point is one where tangent vector is zero.Figure 1
The predictor-corrector procedure of continuation power flow.The CPF is performed based on the predictor-corrector procedure via the following system of nonlinear equation [22]:(7)Fx=0,where(8)x=δ,UT.The nonlinear power flow equation in equation (8) is augmented by the loading factor λ as follows:(9)Fx,λ=0,0≤λ≤λcritical,where U is the vector of bus voltage magnitudes and δ is the vector of bus voltage angles. The dimension of F is 2n+g, as mentioned, n denotes the number of PQ buses, and g denotes the number of PV buses.The power flow equation in equation (1) is rewritten, comprising a loading factor λ as follows:(10)Pj−λ+γkGeneration,iPGeneration_o,j+λPLoad_o,j=0,Qj−QGeneration,j+λQLoad_o,j=0,where PLoad_o,j and QLoad_o,j are, respectively, the original load demands at bus j, PGeneration_o,j is the original generation demands at bus j, kGeneration,i is the generation factor at bus j, and γ is the generator participation coefficient.From equation (9), the base solution corresponding to the loading factor equal to zero started to determine the next solutions corresponding to different load levels and the critical solution corresponding to the load parameter equal to the critical. This can be shown from the tracing of the PV curve based on the predictor-corrector scheme as shown in Figure 1. The predictor-corrector continuation process procedure can be realized in two steps.


### 2.2.1. The First Step

This step realizes the predictor process to calculate the tangent vector and can be expressed in the following form:(11)t=dx,dλT.The tangent vectort with a corresponding dimension 2n+g+1 is determined from the following form:(12)Fx,Fλt=0,where Fx and Fλ are the partial derivative of equation (9) with respect to x and λ, respectively, and T denotes the transpose operation.From equation (12), the left side is a partial derivative matrix that this matrix is nothing but the original Jacobian augmented by one column [Fλ] of an unknown variable of load λ and a vector of differentials t is the differentials representing the tangent vector. Since an unknown variable of the load is added to the nonlinear system (9), this will lead the augmented Jacobian to be singular at the point of maximum possible system load. In order to solve this problem, the tangent vector should be determined based on the satisfied augment Jacobian as follows:(13)FxFλeit=0ti,where ei is a row vector having an appropriate dimension, in which the ith element is the only nonzero element.From equation (13), ti is a nonzero element of the tangent vector; this parameter can be either +1 or −1 depending on how the ith state variable is changing as a solution that is being traced. An estimate of the next solution can be determined by solving system (13) as follows:(14)xeλe=xiλi+σdxdλ,where the subscripts “i” and “e” denote the current and estimated solutions of the next step, respectively, and σ is a scalar defining the predictor step length control.


### 2.2.2. The Second Step

After the prediction is finished, if the predictor solution may not be exactly on a desired solution curve, this step is realized to correct the predicted solution that is determined in equation (14) using a technique, namely, the local parameterization. Now, for the corrector step, system (9) is introduced by one equation unit with the purpose of determining the value of one of the bus voltage magnitude, bus voltage angle, and loading factor. Hence, the new system is expressed as follows:(15)Fx,λxi−μ=00,where xi is the ith state variable that is selected as the continuation parameter and μ is the predicted value of xi.Equation (15) shows that it is involved due to one additional equation and a state variable. In order to solve this difficulty, this paper reproposes Newton’s method for this second step as follows:(16)FxFλeiΔxΔλ=Fx,λxi−η.


## 2.2.1. The First Step

This step realizes the predictor process to calculate the tangent vector and can be expressed in the following form:(11)t=dx,dλT.The tangent vectort with a corresponding dimension 2n+g+1 is determined from the following form:(12)Fx,Fλt=0,where Fx and Fλ are the partial derivative of equation (9) with respect to x and λ, respectively, and T denotes the transpose operation.From equation (12), the left side is a partial derivative matrix that this matrix is nothing but the original Jacobian augmented by one column [Fλ] of an unknown variable of load λ and a vector of differentials t is the differentials representing the tangent vector. Since an unknown variable of the load is added to the nonlinear system (9), this will lead the augmented Jacobian to be singular at the point of maximum possible system load. In order to solve this problem, the tangent vector should be determined based on the satisfied augment Jacobian as follows:(13)FxFλeit=0ti,where ei is a row vector having an appropriate dimension, in which the ith element is the only nonzero element.From equation (13), ti is a nonzero element of the tangent vector; this parameter can be either +1 or −1 depending on how the ith state variable is changing as a solution that is being traced. An estimate of the next solution can be determined by solving system (13) as follows:(14)xeλe=xiλi+σdxdλ,where the subscripts “i” and “e” denote the current and estimated solutions of the next step, respectively, and σ is a scalar defining the predictor step length control.


## 2.2.2. The Second Step

After the prediction is finished, if the predictor solution may not be exactly on a desired solution curve, this step is realized to correct the predicted solution that is determined in equation (14) using a technique, namely, the local parameterization. Now, for the corrector step, system (9) is introduced by one equation unit with the purpose of determining the value of one of the bus voltage magnitude, bus voltage angle, and loading factor. Hence, the new system is expressed as follows:(15)Fx,λxi−μ=00,where xi is the ith state variable that is selected as the continuation parameter and μ is the predicted value of xi.Equation (15) shows that it is involved due to one additional equation and a state variable. In order to solve this difficulty, this paper reproposes Newton’s method for this second step as follows:(16)FxFλeiΔxΔλ=Fx,λxi−η.


## 3. Proposed Method

As analyzed in Section2, when using the CCPF, firstly, the tangential predictor step with the P-V curve at the previous solution with a constant step length control σ as seen in equation (14) is applied. Then, the corrector step based on the parameterization method is used by solving equation (16). Obviously, using the tangential predictive method with a constant step length control value may increase the number of predictor and corrector steps and lead to nonconvergence when solving equation (16).In order to overcome these disadvantages, this paper proposed a method by using the scant method to perform the predictive solution. The executorial predictor-corrector procedure is illustrated in Figure2, in which it will help the process of finding solutions to the critical point to be faster as shown in Figure 2(a). If in case that has not found the next solution from the previous solution, the step length control will be cut down and resumed using the corrector step to the next corrected solution as shown in Figure 2(b). The process of finding solutions to the proposed method is presented in the following.Figure 2
The secant predictor-corrector procedure of the proposed method: (a) smooth function and (b) flat corner.
(a)(b)


### 3.1. Predictor

The secant method and arc-length parameterization are used to realize this step. The function of the predictor is to find an approximate point for the next solution. Figure2(a) shows that if the continuation process at ith step is the current corrector solution (xci; λci), the predictor solution will be found in an approximate point (xpi+1; λpi+1) for the next corrector solution (xci−1; λci−1) and can be obtained as follows:(17)xpi+1λpi+1=xciλci+σixciλci−xci−1λci−1.From equation (17), the number of iterations is used to perform the prediction process for finding solutions on the P-V curve depending on the distance between the two solutions (xci; λci) and (xci−1; λci−1). However, the process of finding solutions can quickly reach the voltage collapse point or not depending on the step length control σi. The process of corrector may be ineffective; if the step length control σi is too large, the process of the previous solution will be kept. So that the step length control needs to be adjusted in order to find a new solution around the critical point λcritical.


### 3.2. Corrector

From Figure2(a), in order to find the next correct solution (xci+1; λci+1), the corrector step is performed by adding a line equation that is perpendicular to line connected between the two previous correct solutions (xci−1; λci−1) and (xci; λci) at the prediction solution (xpi+1; λpi+1). The added line equation is expressed as follows:(18)ρx,λ=x−xpi+1λ−λpi+1Txciλci−xpi+1λpi+1=0.The actual correct solution (xci+1; λci+1) on the bifurcation manifold can be calculated on the following set of equations (9) and (18) for x and λ as follows:(19)Fx,λ=0,0≤λ≤λcritical,ρx,λ=0,where ρ. is the additional equation that is expressed in equation (18).Equation (19) can also be solved by using a slightly modified Newton-Raphson power flow method. However, this corrector step cannot converge if the step length control is too large. To overcome this limit, this paper proposed a procedure for cutting the step length control down to a smaller one until the convergence criteria are reached as shown in Figure 2(b).The flowchart of the continuation power flow is shown in Figure3 and can be explained in detail step by step as follows:Step 1: Read all input data including line data and bus data of the power system.Step 2: Set the initial parameters including the initial solutionδ;U;λ0=0;1;0 and the initial step length control σ=1; note that x=δ,UT.Step 3: Calculate the initial power flow with equation (3) using the Newton-Raphson method and update the initial solution δ;U;λ0; note that x=δ,UT.Step 4: Predict: find the (i+1)th predicted solution with equation (17).Step 5: Correct: find the (i+1)th corrected solution by solving equation (19).Step 6: If the convergence criterion is reached, the step length control will be cut using equationσ=σ/2 and go to Step 4. If the convergence criterion is not reached, it will go to Step 7.Step 7: If the critical point is not reached, it will go to Step 4.Step 8: If the critical point is reached, it stops and puts the P-V curves.Figure 3
The flowchart of voltage stability analyses of the proposed method.


### 3.3. Voltage Stability Index

The voltage stability condition in a power system can be represented by using the voltage stability index. This index can be calculated based on the tangent vector of the P-V curve at each operating state. The voltage stability degree at buses, experiencing when the system reached the state of the voltage collapse, is considered. Thus, we have to find out the weakest bus with respect to the voltage stability, so that we can find suitable solutions to improve voltage stability. The weak bus is one that has the greatest ratio of differential change in voltage to differential change in active load for the whole system. Using the reformulated power flow equations, the differential change in active system load is given as follows [19]:(20)dPtotal=∑j=1ndPLoad,j=SΔbase∑j=1nkLoad,jcosφjdλ=Cdλ,where dPLoad,j is the differential change in load jth, kLoad,j is the multiplier to designate the rate of load change at bus jth as λ changes, φj is the power factor angle of load change at bus jth, and SΔbase is the apparent power selected to provide appropriate scaling of λ. The weakest bus hth is determined as follows:(21)dVhdPtotal=dVhCdλ=maxdV1Cdλ,dV2Cdλ,…,dVnCdλ.When the weakest bushth reaches its steady-state voltage stability limit, the differential of changes dλ closes to zero, and the ratio of dVh/Cdλ will become infinite. This ratio is defined as the voltage stability index at the bus hth. The ratio of dVh/Cdλ in equation (21) shows the magnitude comparison of the elements dVj/Cdλ of the tangent vector of the P-V curve at each operating state. If the ratio of differential change in voltage at a bus to differential change in the connected load is the largest one, this bus indicates the weakest one.


## 3.1. Predictor

The secant method and arc-length parameterization are used to realize this step. The function of the predictor is to find an approximate point for the next solution. Figure2(a) shows that if the continuation process at ith step is the current corrector solution (xci; λci), the predictor solution will be found in an approximate point (xpi+1; λpi+1) for the next corrector solution (xci−1; λci−1) and can be obtained as follows:(17)xpi+1λpi+1=xciλci+σixciλci−xci−1λci−1.From equation (17), the number of iterations is used to perform the prediction process for finding solutions on the P-V curve depending on the distance between the two solutions (xci; λci) and (xci−1; λci−1). However, the process of finding solutions can quickly reach the voltage collapse point or not depending on the step length control σi. The process of corrector may be ineffective; if the step length control σi is too large, the process of the previous solution will be kept. So that the step length control needs to be adjusted in order to find a new solution around the critical point λcritical.


## 3.2. Corrector

From Figure2(a), in order to find the next correct solution (xci+1; λci+1), the corrector step is performed by adding a line equation that is perpendicular to line connected between the two previous correct solutions (xci−1; λci−1) and (xci; λci) at the prediction solution (xpi+1; λpi+1). The added line equation is expressed as follows:(18)ρx,λ=x−xpi+1λ−λpi+1Txciλci−xpi+1λpi+1=0.The actual correct solution (xci+1; λci+1) on the bifurcation manifold can be calculated on the following set of equations (9) and (18) for x and λ as follows:(19)Fx,λ=0,0≤λ≤λcritical,ρx,λ=0,where ρ. is the additional equation that is expressed in equation (18).Equation (19) can also be solved by using a slightly modified Newton-Raphson power flow method. However, this corrector step cannot converge if the step length control is too large. To overcome this limit, this paper proposed a procedure for cutting the step length control down to a smaller one until the convergence criteria are reached as shown in Figure 2(b).The flowchart of the continuation power flow is shown in Figure3 and can be explained in detail step by step as follows:Step 1: Read all input data including line data and bus data of the power system.Step 2: Set the initial parameters including the initial solutionδ;U;λ0=0;1;0 and the initial step length control σ=1; note that x=δ,UT.Step 3: Calculate the initial power flow with equation (3) using the Newton-Raphson method and update the initial solution δ;U;λ0; note that x=δ,UT.Step 4: Predict: find the (i+1)th predicted solution with equation (17).Step 5: Correct: find the (i+1)th corrected solution by solving equation (19).Step 6: If the convergence criterion is reached, the step length control will be cut using equationσ=σ/2 and go to Step 4. If the convergence criterion is not reached, it will go to Step 7.Step 7: If the critical point is not reached, it will go to Step 4.Step 8: If the critical point is reached, it stops and puts the P-V curves.Figure 3
The flowchart of voltage stability analyses of the proposed method.


## 3.3. Voltage Stability Index

The voltage stability condition in a power system can be represented by using the voltage stability index. This index can be calculated based on the tangent vector of the P-V curve at each operating state. The voltage stability degree at buses, experiencing when the system reached the state of the voltage collapse, is considered. Thus, we have to find out the weakest bus with respect to the voltage stability, so that we can find suitable solutions to improve voltage stability. The weak bus is one that has the greatest ratio of differential change in voltage to differential change in active load for the whole system. Using the reformulated power flow equations, the differential change in active system load is given as follows [19]:(20)dPtotal=∑j=1ndPLoad,j=SΔbase∑j=1nkLoad,jcosφjdλ=Cdλ,where dPLoad,j is the differential change in load jth, kLoad,j is the multiplier to designate the rate of load change at bus jth as λ changes, φj is the power factor angle of load change at bus jth, and SΔbase is the apparent power selected to provide appropriate scaling of λ. The weakest bus hth is determined as follows:(21)dVhdPtotal=dVhCdλ=maxdV1Cdλ,dV2Cdλ,…,dVnCdλ.When the weakest bushth reaches its steady-state voltage stability limit, the differential of changes dλ closes to zero, and the ratio of dVh/Cdλ will become infinite. This ratio is defined as the voltage stability index at the bus hth. The ratio of dVh/Cdλ in equation (21) shows the magnitude comparison of the elements dVj/Cdλ of the tangent vector of the P-V curve at each operating state. If the ratio of differential change in voltage at a bus to differential change in the connected load is the largest one, this bus indicates the weakest one.


## 4. Study and Simulation Results

To evaluate the effectiveness of the proposed method, the numerical case studies were examined through using three IEEE 14-bus, 57-bus, and 118-bus systems which are considered as a quite simple and small-scale network case, a medium-scale one, and a complex and large-scale system, respectively. The detailed data of these systems have been introduced in [31–33]. All dynamic models such as generators, excitation systems, transmission lines, and load are modeled based on [13]. For each system, the case studies are examined considering conditions as with or without the penetration of DGs, in which the DG is used to be the doubly fed induction generator-based variable-speed wind turbines (DFIG-VSWT) with total installed capacity being 20 MW at each bus. This bus is chosen based on the condition of the placed location to be far in comparison with the existing generation sources. Note also that, in this paper, the DG is suggested to be a pure distributed source with the constant power model and its model and structure are introduced in [34]. The placed locations of DGs at buses of each tested system are listed in Table 1. The numerical simulation results are realized in the following case studies.Table 1
The modification of the test systems.
Test systemDG placed at busesIEEE 14-bus5IEEE 57-bus26, 39, 54IEEE 118-bus2, 21, 41, 57, 101


### 4.1. Test on IEEE 14-Bus System

The IEEE 14-bus system, which is considered as the small-scale network case, has two generation buses, three synchronous compensators at buses 3, 6, and 8 used only for supporting reactive powers, eleven load buses, three static compensators, and three transformers. The total load demand is 259 MW and 73.5 MVAr. More details on this system can be also found in [33].In this case study, we choose the load bus 5 to install the 20 MW DG. The P-V curve of load bus 14 is shown in Figure4 using CCPF and MCPF methods to analyze voltage stability considering either penetrating the installed the 20 MW DG at bus 5 or not and the numerical results are summarized in Table 2.Figure 4
The P-V curve at bus 14 during the load-increasing process: (a) without DG; (b) with DG.
(a)(b)Table 2
Comparisons between CCPF and MCPF of the IEEE 14-bus test system.
Output parameterWithout DGsWith DGsCCPFMCPFCCPFMCPFPredictor-corrector steps41323725Computational time (ms)811.62730.61868.93750.92Maximum loading factorλcritical2.74372.81783.63973.7412Stability margin (MW)451.62470.81683.68709.97Weakest bus14141414From the P-V curve in Figure4, for the case of being without DG, using the CCPF method, it takes too many predictor-corrector steps to pass the bifurcation point, from the numerical results in Table 2; it takes 41 predictor-corrector steps totally and spends 0.81162 seconds finishing the critical point calculation. Meanwhile, using MCPF, it takes 32 predictor-corrector steps and spends 0.73061 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 2.7437 which is small compared to that using MCPF which is 2.8178. In case of being with DG, it takes 37 predictor-corrector steps totally and spends 0.86893 seconds finishing the critical point calculation. Meanwhile, using MCPF, it takes 25 predictor-corrector steps and spends 0.75092 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 3.6397 which is small compared to that using MCPF which is 3.7412. In addition, the stability margin is expanded with two cases when using MCPF.Figure5(a) plots the voltage magnitude profile of the 14-bus system without DG when the loading point has reached a maximum value, which is obtained using the CCPF and MCPF methods. From this figure, when using MCPF, it can detect the unstable voltage corresponding to the maximum loading point at a voltage collapse point to be more sensitive when using CCPF. In case of being with DG, the simulation results are shown in Figures 4(b) and 5(b) and the rest of the numerical results are shown in Table 2; using the CCPF method, it takes too many continuation steps to pass the bifurcation point, the stability margin is limited, and the detection of the voltage stability properties is low compared to the proposed method.Figure 5
Voltage profile of 14-bus system: (a) without DG; (b) with DG.
(a)(b)In the analysis of voltage stability, the relation between the power transfer to the load and voltage of the load bus is not weak. As known, the variation in power transfer from one bus to the other ones will affect the bus voltages. That is the reason why the proposed method is applied to consider its effectiveness based on the P-V curve. The voltage at buses 4, 5, 7, 9, 10, and 14 is plotted concerning the loading factor as shown in Figure6. As the loading factor is increased, the voltage at load buses decreases. The most reduction in bus voltages occurs on bus 14.Figure 6
The P-V curve of buses 4, 5, 7, 9, 10, and 14 during the load-increasing process when using MCPF: (a) without DG; (b) with DG.
(a)(b)The buses have a large voltage stability index to the change in the total load active power in the system to be buses numbers 4, 5, 7, 9, 10, and 14, in which buses 5 and 14 are the weakest in voltage stability aspect under cases of being with and without DG, respectively, as shown in Figure7. It can be concluded that buses 5 and 14 are identified as a critical bus for cases of being without and with DG in the system, respectively. After installing DG at bus 5, bus 5 has become the strongest as shown in Figure 7; naturally, it is an expected result since there are huge generating units.Figure 7
The voltage stability index of the 14-bus system when using MCPF without and with DG at bus 5.From the simulated results in Figures5–7, it can be concluded that bus 5 is the weakest for the case of being without DG and meanwhile the bus 14 is the weakest for two cases. Buses 1, 2, 3, 6, and 8 are applied for the voltage control, in which bus 1 is the slack bus, so the voltage profile is constant despite increasing the loading factor and it is pretty obvious bus 1 is the strongest. The voltage of the remaining buses decreases since they are load buses. Besides, when DG is connected at bus 5, the stability margin and loading factor are expanded as shown in the fourth row of Table 2.


### 4.2. Test on IEEE 57-Bus System

In this section, in order to verify the proposed method, the IEEE 57-bus system can be chosen as a medium-scale system. This system has 7 synchronous machines with IEEE type-1 exciters, three of which are synchronous compensators, 64 buses, 65 transmission lines, 22 transformers, and 42 constant impedance loads. The total load demand is 1,250.8 MW and 336.4 MVAr. More details on the IEEE 57-bus test system can be also found in [32].For this test case, we propose three 20 MW DG units, which are installed at buses 26, 39, and 54, respectively, and the load bus 31 is selected to verify the proposed method. In order to compare easily, the P-V curve of load bus 31 is plotted out in Figure8 considering either installing three 20 MW DGs at buses 26, 39, and 54 or not. The numerical results are summarized in Table 3.Figure 8
The P-V curve at bus 31 during the load-increasing process: (a) without DG; (b) with DG.
(a)(b)Table 3
Comparisons between CCPF and MCPF of the IEEE 57-bus test system.
Output parameterWithout DGsWith DGsCCPFMCPFCCPFMCPFPredictor-corrector steps50382720Computational time (ms)934.64850.21745.66710.25Maximum loading factorλcritical1.7741.84342.79362.8915Stability margin (MW)968.121,054.932,243.442,365.89Weakest bus31313131The obtained results from the case of being without the integration of DGs show that, using the CCPF method, it takes too many predictor-corrector steps to pass the bifurcation point, from the numerical results in Table3; it takes 50 predictor-corrector steps totally and spends 0.93464 seconds finishing the critical point calculation. Meanwhile, using MCPF, it takes 38 predictor-corrector steps and spends 0.85021 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 1.774 which is small compared to that using MCPF which is 1.8434. In case of being with DGs, it takes 27 predictor-corrector steps totally and spends 0.74566 seconds finishing the critical point calculation. Meanwhile, using MCPF, it takes 20 predictor-corrector steps and spends 0.71025 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 2.7936 which is small compared to that using MCPF which is 2.8915. In addition, the stability margin is expanded with two cases when using MCPF as shown in the fourth row of Table 3.The zoom-in for this case is shown in Figure8. From this figure, the effect of step length control reduction is clearly shown, and exactly, in this case, the maximum loading factors reach 1.9241 without DGs and 2.9152 with DGs corresponding to the step length control in equation (17) which is cut down to values 1/2 and 1/4, respectively.Figure9 plots the voltage magnitude profile of the system without and with the integration of DGs when the loading point has reached a maximum value, which is obtained using the CCPF and MCPF methods. From this figure, when using MCPF, it can detect the unstable voltage corresponding to the maximum loading point at a voltage collapse point to be more sensitive when using CCPF.Figure 9
Voltage profile of the 57-bus system buses: (a) without DG at buses 26, 39, and 54; (b) with DG at buses 26, 39, and 54.
(a)(b)Figure10 shows the voltage profiles of buses 25 and 30–34 when using the proposed method. It is seen from this figure that, after the load point with maximum loading factors λcritical of 1.8434 and 2.8915 for the cases of being without and with the integration of DGs, respectively, the bus voltages start to decrease because of the deficient power generation. At these points, the system is defined as the operating conditions and after these points, the system enters into an unstable condition which can cause the phenomena of voltage collapse. When comparing the voltage profiles of buses, bus 34 seems to be the strongest, and bus 31 seems to be the weakest bus under voltage stability facet.Figure 10
Voltage profile of buses from 25 to 30–34 during the load-increasing process when using MCPF: (a) without DG at buses 26, 39, and 54; (b) with DG at buses 26, 39, and 54.
(a)(b)Figure11 plots the voltage stability index of the whole system. From the obtained result, buses 25 and 31 have the large voltage stability index to the change in the total load active power, where bus 31 is the weakest in the voltage stability aspect under cases of being with and without the integration of DGs in the system. From the obtained results in Figures 9–11, it can be concluded that bus 31 is identified as a critical bus.Figure 11
The voltage stability index of 57-bus system when using MCPF without and with DG at buses 26, 39, and 54.


### 4.3. Test on IEEE 118-Bus System

Lastly, to add practicality and to ensure the efficiency of the proposed method, this test case is considered on the IEEE 118-bus system presumed as a complicated system, having 19 generators, 35 synchronous condensers, 177 lines, 9 transformers, and 91 loads. The total load demand is 4,242 MW and 1,438 MVar. Many details can be also found in [31]. Besides the differences in results discussed in the two previous test cases, for this test case, the paper proposes five 20 MW DG units installed at buses 2, 21, 41, 57, and 101, respectively. The load bus 44 is selected to verify and ensure the efficiencies of the proposed method. The results of the P-V curve of load bus 44 using CCPF and MCPF methods considering either penetrating the installed five 20 MW DG units at buses 2, 21, 41, 57, and 101 or not are shown in Figure 12 and the numerical results are summarized in Table 4.Figure 12
The P-V curve at the load bus 44 during the load-increasing process: (a) without DGs and (b) with DGs.
(a)(b)Table 4
Comparisons between CCPF and MCPF of the IEEE 118-bus test system.
Output parameterWithout DGsWith DGsCCPFMCPFCCPFMCPFPredictor-corrector steps177162215198Computational time (ms)3,238.53,184.63,814.23,725.1Maximum loading factorλcritical10.582010.909120.258321.1519Stability margin (MW)40,646.8442,034.481,693.7185,484.36Weakest bus33334444From the obtained results, the use of the proposed method to analyze the voltage stability has the advantage compared to the use of the CCPF method. For example, the predictor-corrector steps are less, the processing time is quicker, and the effect of adding the new generation units to the system is also observed clearly. Without the integration of DGs, it takes 177 predictor-corrector steps totally and spends 3.2385 seconds finishing the critical point calculation when using the CCPF method. Meanwhile, using MCPF, it takes 162 predictor-corrector steps and spends 3.1846 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 10.5820 which is small compared to that using MCPF which is 10.9091. It is the same. For the integration of DGs, it takes 215 predictor-corrector steps totally and spends 3.8142 seconds finishing the critical point calculation when using the CCPF method. Meanwhile, using MCPF, it takes 198 predictor-corrector steps and spends 3.7251 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 20.2583 which is small compared to that using MCPF which is 21.1519. The stability margin is expanded when using MCPF. The obtained results are also summarized in Table 4. The zoom-in for this case is shown in Figure 12. From this figure, the effect of step length control reduction is clearly shown, and exactly, in this case, the maximum loading factors reach 10.9091 without DGs and 21.1519 with DGs corresponding to the step length control in equation (17) which is cut down to value 1/2. The stability margin is expanded when using MCPF and can be shown in the fourth row of Table 4.Figure13 plots the voltage magnitude profile of the system without and with the integration of DGs when the loading point has reached a maximum value, which is obtained using the CCPF and MCPF methods. From this figure, when using MCPF, it can detect the unstable voltage corresponding to the maximum loading point at a voltage collapse point to be more sensitive when using CCPF.Figure 13
Voltage profile of 118 buses: (a) without the integration of DGs and (b) with the integration of DGs.
(a)(b)Figure14 shows the voltage profiles of buses 33 and 43–47 when using the proposed method. It is seen from this figure that, after the load point with maximum loading factors λcritical of 10.9091 and 21.1519 for the cases of being without and with the integration of DGs, respectively, the bus voltages start to decrease because of the deficient power generation. At these points, the system is defined as the operating conditions and after these points, the system enters into an unstable condition which can cause the phenomena of voltage collapse. When comparing the voltage profiles of buses, bus 47 seems to be the strongest, and buses 33 and 44 seem to be the weakest ones under the voltage stability facet for the cases of being without and with the integration of DGs, respectively.Figure 14
Voltage profile of buses during the load-increasing process when using MCPF: (a) without the integration of DGs and (b) with the integration of DGs.
(a)(b)Figure15 plotted the voltage stability index of the whole system. From this figure, buses 33 and 44 have a large voltage stability index to the change in the total load active power in the system. For this case, buses 33 and 44 are the weakest buses in the voltage stability aspect under cases of being with and without the integration of DGs, respectively. Therefore, the simulated results in Figures 13–15 could conclude that buses 33 and 44 are identified as critical ones.Figure 15
The voltage stability index of 118 buses when using MCPF without and with the integration of DGs.


## 4.1. Test on IEEE 14-Bus System

The IEEE 14-bus system, which is considered as the small-scale network case, has two generation buses, three synchronous compensators at buses 3, 6, and 8 used only for supporting reactive powers, eleven load buses, three static compensators, and three transformers. The total load demand is 259 MW and 73.5 MVAr. More details on this system can be also found in [33].In this case study, we choose the load bus 5 to install the 20 MW DG. The P-V curve of load bus 14 is shown in Figure4 using CCPF and MCPF methods to analyze voltage stability considering either penetrating the installed the 20 MW DG at bus 5 or not and the numerical results are summarized in Table 2.Figure 4
The P-V curve at bus 14 during the load-increasing process: (a) without DG; (b) with DG.
(a)(b)Table 2
Comparisons between CCPF and MCPF of the IEEE 14-bus test system.
Output parameterWithout DGsWith DGsCCPFMCPFCCPFMCPFPredictor-corrector steps41323725Computational time (ms)811.62730.61868.93750.92Maximum loading factorλcritical2.74372.81783.63973.7412Stability margin (MW)451.62470.81683.68709.97Weakest bus14141414From the P-V curve in Figure4, for the case of being without DG, using the CCPF method, it takes too many predictor-corrector steps to pass the bifurcation point, from the numerical results in Table 2; it takes 41 predictor-corrector steps totally and spends 0.81162 seconds finishing the critical point calculation. Meanwhile, using MCPF, it takes 32 predictor-corrector steps and spends 0.73061 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 2.7437 which is small compared to that using MCPF which is 2.8178. In case of being with DG, it takes 37 predictor-corrector steps totally and spends 0.86893 seconds finishing the critical point calculation. Meanwhile, using MCPF, it takes 25 predictor-corrector steps and spends 0.75092 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 3.6397 which is small compared to that using MCPF which is 3.7412. In addition, the stability margin is expanded with two cases when using MCPF.Figure5(a) plots the voltage magnitude profile of the 14-bus system without DG when the loading point has reached a maximum value, which is obtained using the CCPF and MCPF methods. From this figure, when using MCPF, it can detect the unstable voltage corresponding to the maximum loading point at a voltage collapse point to be more sensitive when using CCPF. In case of being with DG, the simulation results are shown in Figures 4(b) and 5(b) and the rest of the numerical results are shown in Table 2; using the CCPF method, it takes too many continuation steps to pass the bifurcation point, the stability margin is limited, and the detection of the voltage stability properties is low compared to the proposed method.Figure 5
Voltage profile of 14-bus system: (a) without DG; (b) with DG.
(a)(b)In the analysis of voltage stability, the relation between the power transfer to the load and voltage of the load bus is not weak. As known, the variation in power transfer from one bus to the other ones will affect the bus voltages. That is the reason why the proposed method is applied to consider its effectiveness based on the P-V curve. The voltage at buses 4, 5, 7, 9, 10, and 14 is plotted concerning the loading factor as shown in Figure6. As the loading factor is increased, the voltage at load buses decreases. The most reduction in bus voltages occurs on bus 14.Figure 6
The P-V curve of buses 4, 5, 7, 9, 10, and 14 during the load-increasing process when using MCPF: (a) without DG; (b) with DG.
(a)(b)The buses have a large voltage stability index to the change in the total load active power in the system to be buses numbers 4, 5, 7, 9, 10, and 14, in which buses 5 and 14 are the weakest in voltage stability aspect under cases of being with and without DG, respectively, as shown in Figure7. It can be concluded that buses 5 and 14 are identified as a critical bus for cases of being without and with DG in the system, respectively. After installing DG at bus 5, bus 5 has become the strongest as shown in Figure 7; naturally, it is an expected result since there are huge generating units.Figure 7
The voltage stability index of the 14-bus system when using MCPF without and with DG at bus 5.From the simulated results in Figures5–7, it can be concluded that bus 5 is the weakest for the case of being without DG and meanwhile the bus 14 is the weakest for two cases. Buses 1, 2, 3, 6, and 8 are applied for the voltage control, in which bus 1 is the slack bus, so the voltage profile is constant despite increasing the loading factor and it is pretty obvious bus 1 is the strongest. The voltage of the remaining buses decreases since they are load buses. Besides, when DG is connected at bus 5, the stability margin and loading factor are expanded as shown in the fourth row of Table 2.


## 4.2. Test on IEEE 57-Bus System

In this section, in order to verify the proposed method, the IEEE 57-bus system can be chosen as a medium-scale system. This system has 7 synchronous machines with IEEE type-1 exciters, three of which are synchronous compensators, 64 buses, 65 transmission lines, 22 transformers, and 42 constant impedance loads. The total load demand is 1,250.8 MW and 336.4 MVAr. More details on the IEEE 57-bus test system can be also found in [32].For this test case, we propose three 20 MW DG units, which are installed at buses 26, 39, and 54, respectively, and the load bus 31 is selected to verify the proposed method. In order to compare easily, the P-V curve of load bus 31 is plotted out in Figure8 considering either installing three 20 MW DGs at buses 26, 39, and 54 or not. The numerical results are summarized in Table 3.Figure 8
The P-V curve at bus 31 during the load-increasing process: (a) without DG; (b) with DG.
(a)(b)Table 3
Comparisons between CCPF and MCPF of the IEEE 57-bus test system.
Output parameterWithout DGsWith DGsCCPFMCPFCCPFMCPFPredictor-corrector steps50382720Computational time (ms)934.64850.21745.66710.25Maximum loading factorλcritical1.7741.84342.79362.8915Stability margin (MW)968.121,054.932,243.442,365.89Weakest bus31313131The obtained results from the case of being without the integration of DGs show that, using the CCPF method, it takes too many predictor-corrector steps to pass the bifurcation point, from the numerical results in Table3; it takes 50 predictor-corrector steps totally and spends 0.93464 seconds finishing the critical point calculation. Meanwhile, using MCPF, it takes 38 predictor-corrector steps and spends 0.85021 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 1.774 which is small compared to that using MCPF which is 1.8434. In case of being with DGs, it takes 27 predictor-corrector steps totally and spends 0.74566 seconds finishing the critical point calculation. Meanwhile, using MCPF, it takes 20 predictor-corrector steps and spends 0.71025 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 2.7936 which is small compared to that using MCPF which is 2.8915. In addition, the stability margin is expanded with two cases when using MCPF as shown in the fourth row of Table 3.The zoom-in for this case is shown in Figure8. From this figure, the effect of step length control reduction is clearly shown, and exactly, in this case, the maximum loading factors reach 1.9241 without DGs and 2.9152 with DGs corresponding to the step length control in equation (17) which is cut down to values 1/2 and 1/4, respectively.Figure9 plots the voltage magnitude profile of the system without and with the integration of DGs when the loading point has reached a maximum value, which is obtained using the CCPF and MCPF methods. From this figure, when using MCPF, it can detect the unstable voltage corresponding to the maximum loading point at a voltage collapse point to be more sensitive when using CCPF.Figure 9
Voltage profile of the 57-bus system buses: (a) without DG at buses 26, 39, and 54; (b) with DG at buses 26, 39, and 54.
(a)(b)Figure10 shows the voltage profiles of buses 25 and 30–34 when using the proposed method. It is seen from this figure that, after the load point with maximum loading factors λcritical of 1.8434 and 2.8915 for the cases of being without and with the integration of DGs, respectively, the bus voltages start to decrease because of the deficient power generation. At these points, the system is defined as the operating conditions and after these points, the system enters into an unstable condition which can cause the phenomena of voltage collapse. When comparing the voltage profiles of buses, bus 34 seems to be the strongest, and bus 31 seems to be the weakest bus under voltage stability facet.Figure 10
Voltage profile of buses from 25 to 30–34 during the load-increasing process when using MCPF: (a) without DG at buses 26, 39, and 54; (b) with DG at buses 26, 39, and 54.
(a)(b)Figure11 plots the voltage stability index of the whole system. From the obtained result, buses 25 and 31 have the large voltage stability index to the change in the total load active power, where bus 31 is the weakest in the voltage stability aspect under cases of being with and without the integration of DGs in the system. From the obtained results in Figures 9–11, it can be concluded that bus 31 is identified as a critical bus.Figure 11
The voltage stability index of 57-bus system when using MCPF without and with DG at buses 26, 39, and 54.


## 4.3. Test on IEEE 118-Bus System

Lastly, to add practicality and to ensure the efficiency of the proposed method, this test case is considered on the IEEE 118-bus system presumed as a complicated system, having 19 generators, 35 synchronous condensers, 177 lines, 9 transformers, and 91 loads. The total load demand is 4,242 MW and 1,438 MVar. Many details can be also found in [31]. Besides the differences in results discussed in the two previous test cases, for this test case, the paper proposes five 20 MW DG units installed at buses 2, 21, 41, 57, and 101, respectively. The load bus 44 is selected to verify and ensure the efficiencies of the proposed method. The results of the P-V curve of load bus 44 using CCPF and MCPF methods considering either penetrating the installed five 20 MW DG units at buses 2, 21, 41, 57, and 101 or not are shown in Figure 12 and the numerical results are summarized in Table 4.Figure 12
The P-V curve at the load bus 44 during the load-increasing process: (a) without DGs and (b) with DGs.
(a)(b)Table 4
Comparisons between CCPF and MCPF of the IEEE 118-bus test system.
Output parameterWithout DGsWith DGsCCPFMCPFCCPFMCPFPredictor-corrector steps177162215198Computational time (ms)3,238.53,184.63,814.23,725.1Maximum loading factorλcritical10.582010.909120.258321.1519Stability margin (MW)40,646.8442,034.481,693.7185,484.36Weakest bus33334444From the obtained results, the use of the proposed method to analyze the voltage stability has the advantage compared to the use of the CCPF method. For example, the predictor-corrector steps are less, the processing time is quicker, and the effect of adding the new generation units to the system is also observed clearly. Without the integration of DGs, it takes 177 predictor-corrector steps totally and spends 3.2385 seconds finishing the critical point calculation when using the CCPF method. Meanwhile, using MCPF, it takes 162 predictor-corrector steps and spends 3.1846 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 10.5820 which is small compared to that using MCPF which is 10.9091. It is the same. For the integration of DGs, it takes 215 predictor-corrector steps totally and spends 3.8142 seconds finishing the critical point calculation when using the CCPF method. Meanwhile, using MCPF, it takes 198 predictor-corrector steps and spends 3.7251 seconds. The maximum loading factor (λcritical) at the critical point when using CCPF is 20.2583 which is small compared to that using MCPF which is 21.1519. The stability margin is expanded when using MCPF. The obtained results are also summarized in Table 4. The zoom-in for this case is shown in Figure 12. From this figure, the effect of step length control reduction is clearly shown, and exactly, in this case, the maximum loading factors reach 10.9091 without DGs and 21.1519 with DGs corresponding to the step length control in equation (17) which is cut down to value 1/2. The stability margin is expanded when using MCPF and can be shown in the fourth row of Table 4.Figure13 plots the voltage magnitude profile of the system without and with the integration of DGs when the loading point has reached a maximum value, which is obtained using the CCPF and MCPF methods. From this figure, when using MCPF, it can detect the unstable voltage corresponding to the maximum loading point at a voltage collapse point to be more sensitive when using CCPF.Figure 13
Voltage profile of 118 buses: (a) without the integration of DGs and (b) with the integration of DGs.
(a)(b)Figure14 shows the voltage profiles of buses 33 and 43–47 when using the proposed method. It is seen from this figure that, after the load point with maximum loading factors λcritical of 10.9091 and 21.1519 for the cases of being without and with the integration of DGs, respectively, the bus voltages start to decrease because of the deficient power generation. At these points, the system is defined as the operating conditions and after these points, the system enters into an unstable condition which can cause the phenomena of voltage collapse. When comparing the voltage profiles of buses, bus 47 seems to be the strongest, and buses 33 and 44 seem to be the weakest ones under the voltage stability facet for the cases of being without and with the integration of DGs, respectively.Figure 14
Voltage profile of buses during the load-increasing process when using MCPF: (a) without the integration of DGs and (b) with the integration of DGs.
(a)(b)Figure15 plotted the voltage stability index of the whole system. From this figure, buses 33 and 44 have a large voltage stability index to the change in the total load active power in the system. For this case, buses 33 and 44 are the weakest buses in the voltage stability aspect under cases of being with and without the integration of DGs, respectively. Therefore, the simulated results in Figures 13–15 could conclude that buses 33 and 44 are identified as critical ones.Figure 15
The voltage stability index of 118 buses when using MCPF without and with the integration of DGs.


## 5. Conclusions

This paper proposes an innovatory method for analyzing the voltage stability and specifically monitoring the bus voltage when the system is operating at a load point near the critical one. This proposed method is developed based on the continuation power flow (CPF) method. The voltage stability problem of the power system has been analyzed to establish the proposed method; the CPF method based on the tangent and local parameterization methods is recalled to be compared with the proposed method. The proposed method is realized on the predictor and corrector procedures to draw the P-V curves at buses according to a specified generation-load scenario.Three IEEE 14-bus, 57-bus, and 118-bus test systems are considered as a quite simple and small-scale network case, a medium-scale one, and a complex and large-scale system, respectively, and the distributed generation (DG) is used to be the doubly fed induction generator-based variable-speed wind turbines (DFIG-VSWT) with the constant power to verify the efficiency of the proposed method. The numerical results are simulated for all study cases based on the load-increasing process and voltage stability index by using MATLAB software on a PC with Intel(R) Core processor (TM) i7, 3.2 GHz. The obtained results show that using the proposed method to analyze the voltage stability has the advantages compared to the CCPF method; namely, the predictor-corrector steps are less, the processing time is quicker, the effect of adding the new generation units into the system is also observed clearly during the load-increasing process, and the stability margin and loading factor are expanded. In addition, the proposed method was shown to be effective on a large power system with the integration of many DG units and the voltage stability index indicated closer proximity to voltage collapse when the system is operating at a load point near the critical one.



---

*Source: 8037837-2020-09-07.xml*