#### 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.

#### 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 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. 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 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 generators (including DG units), the number of unknowns which is 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]:where and are the generator’s active and reactive powers at bus *j*, respectively; and are the load’s active and reactive powers at bus *j*, respectively.

The complex power at bus *j* can be given bywhere is the voltage at bus *j* and is the current injected into bus *j* which can be written as follows:

where 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 getswhere is mutual admittance between buses *j* and *k*. The active and reactive power balance equation at bus *j* can be written as follows:

where and are, respectively, the real and imaginary part of the element in the bus admittance matrix corresponding to the *j*th row and *k*th column. The voltage phase angle between buses *j* and *k* is , in which and 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 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. This process will be finished when the critical point is reached; in other words, the critical point is one where tangent vector is zero.

The CPF is performed based on the predictor-corrector procedure via the following system of nonlinear equation [22]:where

The 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. The dimension of *F* is , as mentioned, *n* denotes the number of PQ buses, and denotes the number of PV buses.

The power flow equation in equation (1) is rewritten, comprising a loading factor as follows:where and are, respectively, the original load demands at bus *j*, is the original generation demands at bus *j*, 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:

The tangent vector *t* with a corresponding dimension is determined from the following form:where and 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 [] 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 is a row vector having an appropriate dimension, in which the *i*th element is the only nonzero element.

From equation (13), is a nonzero element of the tangent vector; this parameter can be either +1 or −1 depending on how the *i*th 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.

###### 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:where is the *i*th state variable that is selected as the continuation parameter and is the predicted value of .

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:

#### 3. 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. 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 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). The process of finding solutions to the proposed method is presented in the following.

**(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. Figure 2(a) shows that if the continuation process at *i*th step is the current corrector solution (; ), the predictor solution will be found in an approximate point (; ) for the next corrector solution (; ) 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 solutions (; ) and (; ). However, the process of finding solutions can quickly reach the voltage collapse point or not depending on the step length control . The process of corrector may be ineffective; if the step length control 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 .

##### 3.2. Corrector

From Figure 2(a), in order to find the next correct solution (; ), the corrector step is performed by adding a line equation that is perpendicular to line connected between the two previous correct solutions (; ) and (; ) at the prediction solution (; ). The added line equation is expressed as follows:

The actual correct solution (; ) 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).

The 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 2: Set the initial parameters including the initial solution and the initial step length control ; note that . Step 3: Calculate the initial power flow with equation (3) using the Newton-Raphson method and update the initial solution ; note that . 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 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.

##### 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]:where is the differential change in load *j*th, is the multiplier to designate the rate of load change at bus *j*th as changes, is the power factor angle of load change at bus *j*th, and is the apparent power selected to provide appropriate scaling of . The weakest bus *h*th is determined as follows:

When the weakest bus *h*th reaches its steady-state voltage stability limit, the differential of changes closes to zero, and the ratio of will become infinite. This ratio is defined as the voltage stability index at the bus *h*th. The ratio of in equation (21) shows the magnitude comparison of the elements 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.

##### 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 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.

**(a)**

**(b)**

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. The maximum loading factor () 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 () 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.

Figure 5(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.

**(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 Figure 6. As the loading factor is increased, the voltage at load buses decreases. The most reduction in bus voltages occurs on bus 14.

**(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 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; naturally, it is an expected result since there are huge generating units.

From the simulated results in 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. 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 Figure 8 considering either installing three 20 MW DGs at buses 26, 39, and 54 or not. The numerical results are summarized in Table 3.

**(a)**

**(b)**

The 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. The maximum loading factor () 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 () 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 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.

**(a)**

**(b)**

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 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.

**(a)**

**(b)**

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.

##### 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.

**(a)**

**(b)**

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. The maximum loading factor () 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 () 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.

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.

**(a)**

**(b)**

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 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.

**(a)**

**(b)**

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. Therefore, the simulated results in Figures 13–15 could conclude that buses 33 and 44 are identified as critical ones.

#### 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.

#### Data Availability

The data used to support the study are presented in [31–33].

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The authors sincerely acknowledge the financial support provided by Industrial University of Ho Chi Minh City, Ton Duc Thang University, and Quy Nhon University, Vietnam, for carrying out this work.