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

Faculty of Electrical Engineering Technology, Industrial University of Ho Chi Minh City, Ho Chi Minh City, Vietnam Faculty of Engineering and Technology, Quy Nhon University, Quy Nhon City, Binh Dinh, Vietnam Department for Management of Science and Technology Development, Ton Duc (ang University, Ho Chi Minh City, Vietnam Faculty of Electrical and Electronics Engineering, Ton Duc (ang University, Ho Chi Minh City, Vietnam


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. e 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]. e impact of DGs on the voltage stability in the long and short terms is examined in [5][6][7][8][9], respectively. In addition, the DGs' optimal allocation, considering uncertainty in load demand and growth, is studied in [10][11][12].
e 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]. e 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. e 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].

Literature Review.
e blackouts in the power system can cause great damage to socioeconomic and brought inconvenience to people. erefore, 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]. e 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]. e aforesaid methods based on the power flow as the modal analysis have been introduced in [23]. is 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. e 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]. ese 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. erefore, 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].

Contribution and Study Organization.
is 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. is 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. e 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.
e remaining parts of the paper can be divided into four sections as follows.
e voltage stability problem is formulated in Section 2 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. e effectiveness of the proposed method is verified in Section 4. Finally, the conclusions are given in Section 5.

Problem Formulation.
e 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. e power flow calculation aims to determine voltage magnitudes and angles for a given generation, load, and grid condition. e 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 the P 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 2 Complexity of unknowns which is 2(n − 1) − (g − 1) can be calculated by using the reactive and active power balance equations. erefore, the active and reactive power balance equations injected into the network at bus j can be defined as follows [30]: where P Generation,j and Q Generation,j are the generator's active and reactive powers at bus j, respectively; P Load,j and Q Load,j are the load's active and reactive powers at bus j, respectively. e complex power at bus j can be given by where U j is the voltage at bus j and I j is the current injected into bus j which can be written as follows: where U k is the voltage at bus k and equation (3) can be rewritten under the following relationship between the bus voltage and current as follows: Substituting equation (3) into equation (2), one gets where Y jk is mutual admittance between buses j and k. e active and reactive power balance equation at bus j can be written as follows: where G jk and B jk are, respectively, the real and imaginary part of the element in the bus admittance matrix corresponding to the jth row and kth column. e 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.

Analyses.
e 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. e CPF determines the consecutive load flow solutions under a load scenario using the prediction and correction steps. e prediction-correction procedure is illustrated in Figure 1. 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. is process will be finished when the critical point is reached; in other words, the critical point is one where tangent vector is zero.
e CPF is performed based on the predictor-corrector procedure via the following system of nonlinear equation [22]: where e nonlinear power flow equation in equation (8) is augmented by the loading factor λ as follows: where U is the vector of bus voltage magnitudes and δ is the vector of bus voltage angles. e dimension of F is 2n + g, as mentioned, n denotes the number of PQ buses, and g denotes the number of PV buses. e power flow equation in equation (1) is rewritten, comprising a loading factor λ as follows: where P Load o,j and Q Load o,j are, respectively, the original load demands at bus j, P Generation o,j is the original generation demands at bus j, k Generation,i is the generation factor at bus j, and c 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. is can be shown from the tracing of the PV curve based on the predictor-corrector scheme as shown in Figure 1.
e predictor-corrector continuation process procedure can be realized in two steps.

e First
Step.
is step realizes the predictor process to calculate the tangent vector and can be expressed in the following form: e tangent vector t with a corresponding dimension 2n + g + 1 is determined from the following form: where F x 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: where e i is a row vector having an appropriate dimension, in which the ith element is the only nonzero element. From equation (13), t i 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: 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.

e 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: where x i is the ith state variable that is selected as the continuation parameter and μ is the predicted value of x i . 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:

Proposed Method
As analyzed in Section 2, 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. en, 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. e executorial predictor-corrector procedure is illustrated in Figure 2, 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). e process of finding solutions to the proposed method is presented in the following.

Predictor.
e secant method and arc-length parameterization are used to realize this step. e function of the predictor is to find an approximate point for the next solution. Figure 2(a) shows that if the continuation process at ith step is the current corrector solution (x (i) c ; λ (i) c ), the predictor solution will be found in an approximate point and can be obtained as follows: 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 so- . However, the process of  Complexity finding solutions can quickly reach the voltage collapse point or not depending on the step length control σ (i) . e 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 . Figure 2(a), in order to find the next correct solution (

Corrector. From
, the corrector step is performed by adding a line equation that is perpendicular to line connected between the two previous correct solutions ( ). e added line equation is expressed as follows: ) on the bifurcation manifold can be calculated on the following set of equations (9) and (18) for x and λ as follows: 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). e flowchart of the continuation power flow is shown in Figure 3 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 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.

Voltage Stability Index.
e voltage stability condition in a power system can be represented by using the voltage stability index. is index can be calculated based on the tangent vector of the P-V curve at each operating state. e voltage stability degree at buses, experiencing when the system reached the state of the voltage collapse, is considered. us, 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. e 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]: Successful corrector Unsuccessful corrector where dP Load,j is the differential change in load jth, k Load,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 λ. e weakest bus hth is determined as follows: When the weakest bus hth reaches its steady-state voltage stability limit, the differential of changes dλ closes to zero, and the ratio of dV h /Cdλ will become infinite. is ratio is defined as the voltage stability index at the bus hth. e ratio of dV h /Cdλ in equation (21) shows the magnitude comparison of the elements dV j /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.

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. e detailed data of these systems have been introduced in [31][32][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. is 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]. e placed locations of DGs at buses of each tested system are listed in Table 1. e numerical simulation results are realized in the following case studies.

Test on IEEE 14-Bus System.
e 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. e 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. e P-V curve of load bus 14 is shown in Figure 4 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.
From the P-V curve in Figure 4, 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. e 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. e 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 Step length control σ = 1 Note: x = [δ, U] T Predictor At (i + 1) using equation (17) Corrector At (i + 1) using equation (19) Check if critical point is reached 6 Complexity addition, the stability margin is expanded with two cases when using MCPF. Figure 5(a) plots the voltage magnitude profile of the 14bus 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.
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. at is the reason why the proposed method is applied to consider its effectiveness based on the P-V curve. e voltage at buses 4, 5, 7, 9, 10, and 14 is plotted concerning the loading factor as shown in Figure 6. As the loading factor is increased, the voltage at load buses decreases. e most reduction in bus voltages occurs on bus 14.
e 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 Figure 7. 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;     Figures 5-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. e 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.

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. is 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. e total load demand is 1,250.8 MW and 336.4 MVAr. More  Figure 6: e 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. 8 Complexity 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 Figure 8 considering either installing three 20 MW DGs at buses 26, 39, and 54 or not. e numerical results are summarized in Table 3. e 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 Table 3; 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. e 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. e 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. e zoom-in for this case is shown in Figure 8. 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. Figure 9 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 10 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 11 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.

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. e 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. e load bus 44 is selected to  Figure 12 and the numerical results are summarized in Table 4.  From 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. e 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 predictorcorrector 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. e 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. e stability margin is expanded when using MCPF. e obtained results are also summarized in Table 4.    Complexity e 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. e stability margin is     Table 4. Figure 13 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 14 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 15 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. erefore, the simulated results in Figures 13-15 could conclude that buses 33 and 44 are identified as critical ones.

Conclusions
is 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.
is proposed method is developed based on the continuation power flow (CPF) method. e 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. e proposed method is realized on the predictor and corrector procedures to draw Voltage (p.u.)   the P-V curves at buses according to a specified generationload scenario. ree 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. e numerical results are simulated for all study cases based on the loadincreasing process and voltage stability index by using MATLAB software on a PC with Intel(R) Core processor (TM) i7, 3.2 GHz. e 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.

Data Availability
e data used to support the study are presented in [31][32][33].

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.