# A Parameterization Technique for the Continuation Power Flow Developed from the Analysis of Power Flow Curves

Authors: Elisabete de Mello Magalhães; Alfredo Bonini Neto; Dilson Amancio Alves
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This paper presents an efficient geometric parameterization technique for the continuation power flow. It was developed from the observation of the geometrical behavior of load flow solutions. The parameterization technique eliminates the singularity of load flow Jacobian matrix and therefore all the consequent problems of ill-conditioning. This is obtained by adding equations lines passing through the points in the plane determined by the loading factor and the total real power losses that is rewritten as a function of the real power generated by the slack bus. An automatic step size control is also provided, which is used when it is necessary. Thus, the resulting method enables the complete tracing ofP-V curves and the computation of maximum loading point of any electric power systems. Intending to reduce the CPU time, the effectiveness caused by updating the Jacobian matrix is investigated only when the system undergoes a significant change. Moreover, the tangent and trivial predictors are compared with each other. The robustness and simplicity as well as the simple interpretation of the proposed technique are the highlights of this method. The results obtained for the IEEE 300-bus system and for real large systems show the effectiveness of the proposed method.


---

## Body

## 1. Introduction

The power flow problem (PF) consists of an algebraic analysis of power system under steady-state operating conditions. In this analysis, the electric power system is represented by a set of nonlinear algebraic equations that are used for computing the operating points of the electrical power system for various loading conditions. Its solution provides the magnitudes and angles of voltages, tap changer setting values for the on-load tap changers (OLTCs) transformers, and the real and reactive power flows and losses on each branch of power network (line transmission and transformer). The PF is used in steady-state stability analysis for assessing voltage stability margins and the areas prone to voltage collapse. It is important to know whether the system has a feasible and secure operating point when either a sudden change in system loading or branch outages occur. When the PF equations have no solution for a given loading condition, it is concluded that the generation and network are not physically able to meet the demand required. In this situation, modifications are necessary in the generation dispatch and/or in the electrical network topology.Among the three types of load representation (constant power, constant current, and constant impedance) for steady-state stability analysis, constant power typically results in the most pessimistic MLP and in the smallest voltage collapse margin [1, 2]. Unless more accurate load representations are known, it is recommended to use the constant power representation because it approximates the action of distribution system voltage regulating devices [2]. This model will result in a more secure operational system condition and should be used if the systems security is evaluated through the maintenance of a minimum voltage stability margin [1, 2]. However, for systems where the constant PQ loads models are used, the gradual load increment will lead to a saddle-node bifurcation, which corresponds to the maximum loading point (MLP) [3–5]. The conventional power flow presents convergence problems to obtain the MLP of electric power systems, because the Jacobian matrix is singular at this point. Therefore, the use of conventional PF to obtain the power flow curves (P-V curve, which is the curve of bus voltage profiles as a function of their loading) is restricted to its upper part. The MLP indirectly defines the boundary between the stable and unstable operating regions and it is important in the application of modal analysis, since the most important information regarding effective remedial measures to enhance system voltage stability is obtained at the MLP or in its adjacent. Besides the MLP determination, these curves have also been used to determine the maximum power transfer among system areas, to adjust margins, and to compare planning strategies [6].By reformulation of the power flow equations, the continuation methods eliminate the singularity of Jacobian matrix and the related numerical problems. Usually this is done by adding parameterized equations [5–13]. Due to robustness provided by these methods in solving nonlinear algebraic equations [9], they have been widely used in the analysis of electric power systems for obtaining multiple solutions, contingency analysis, real power losses reduction, the tracing of loading curves (P-V curves), and the determination of MLP [5–8, 10, 11, 13, 14]. Such publications, including the latest books on the subject [15–19], show that there is a growing interest in the power industry even in small improvements of CPF methods, which provide better performance for tracing the whole P-V curves. The most common parameterization techniques used by CPF to remove the singularity of the PF Jacobian matrix are the geometric [5, 8, 10, 12, 20] and local ones [7, 9].The continuation power flow traces the completeP-V curves by automatically changing the value of a parameter. In the local parameterization technique [7], a parameter change always occurs close to MLP. Generally, the loading factor (λ) is an initially chosen parameter. Close to the MLP, it changes to the voltage magnitude that presents the largest variations and after a few points, it changes back to λ. The voltage magnitudes and angles may also be chosen as parameters, but, in these cases, the new Jacobian matrix can become singular at the MLP, or in the lower or upper part of the P-V curve [8].The addition of the equation of total real power losses to the PF equations has proposed in [10]. In this case, instead of specifying the loading factor and geting the converged state, it specifies the desired amount of total real power losses, and the solution provides the operating point, including the loading factor, for those that losses occur. Adopting a fixed step size for the value of new parameter and through successive solutions of the new system of equations, one can determine all the other points on the P-V curve. The advantage given to the use of this technique was that, in most cases examined, the local parameterization was only needed to points located just after the MLP. Later it was found that for many large systems, the singularities of both Jacobian matrices are practically coincident, that is, the noses are also coincident. Therefore, in many cases still remains the difficulty to know the real cause of the divergence. It is consequent of a poor initial guess, of physical limitation of electric power system, or of numerical problems related to power flow algorithms. To overcome these limitations, in [13], a line equation was added to the power flow problem, which passes through a point in the plane determined by the loading factor (λ) and total real power losses (Pa). The angular coefficient of line is the only used parameter, but in order to avoid the singularity of the Jacobian matrix of the new parameter and thus determination of the MLP, it was needed to define an automatic procedure to switch from a set of line to other. Nevertheless, the method fails in determining the MLP of some systems, such as the 904-bus, even with the selection of very small step sizes. The method proposed in [13] fails to obtain the MLP because it is a real large, heavily loaded, electrical power system with voltage instability problems that have the strong local characteristic. For systems like this, the P-V curve of most buses presents a sharp nose, that is, the loading factor and the voltage magnitude present a simultaneous reversion in its variation tendency, and they reach their maximum value at the same point. In other words, the curve noses are coincident and both Jacobian matrices are singular at the MLP.Aiming to improve the technique proposed by Garbelini et al. in [13], this paper proposes some modifications consisting in rewriting the equation of total real power losses as a function of the real power generated at the slack bus; using the coordinates of the set of lines of a point located between two points close to the MLP, in case of divergence; using the evolution of total power mismatch rather than a fixed step to change the coordinates of the set of lines; considering the total real power losses normalized by its base case value. With these modifications, the proposed method allows the complete tracing of P-V curves, obtaining of MLP, and afterward the assessment of voltage stability margin.The proposed geometric parameterization technique shows the robustness and also is simple and easy to implement and interpret. It is applied to obtain the wholeP-V curves of the IEEE 300-bus system and from three real large, heavy loaded systems such as 638-bus and 787-bus systems that correspond to parts of South-Southeast Brazilian system, and of a 904-bus Southwestern American system. The results show that the method presents good convergence characteristics and the MLP can be computed with any specified precision, without the numerical problems related to the singularity.


## 2. Formulation of the Proposed Continuation Power Flow

TheP-V curve is obtained by tracing a bus voltage profile as a function of its loading. To automate this procedure, the load flow equations are reformulated to include the loading factor (λ), which is used to gradually increase load and the generation level. The new set of PF equations is written as
(2.1)ΔP(θ,V,λ)=Psp(λ)-Pcal(θ,V)=λ(Pgsp-Pcsp)-Pcal(θ,V)=0ΔQ(θ,V,λ)=Qsp(λ)-Qcal(θ,V)=λQcsp-Qcal(θ,V)=0,
where V and θ are, respectively, the vectors of voltage magnitudes and phase angles, Psp(λ) is the difference between the vectors of real power generation (Pgsp) and consumption (Pcsp) specified at the load (PQ) and generator (PV) buses, and Qsp(λ) is the vector of reactive power consumption (Qcsp) specified at the PQ buses. The superscripts sp and cal mean specified and calculated values, respectively. The real and reactive powers at bus k are given by
(2.2)Pkcal(θ,V)=Vk∑m=1NBVm[Gkmcos(θk-θm)+Bkmsin(θk-θm)]Qkcal(θ,V)=Vk∑m=1NBVm[Gkmsin(θk-θm)-Bkmcos(θk-θm)],
where NB is the total number of system buses, Vk and θk are the voltage magnitude and angle at the bus k, and Gkm and Bkm terms represent the conductance and susceptance of (k, m) element in the nodal admittance matrix Y=[G]+j[B].The system of (2.1) considers a network loading proportional to the base case and a constant power factor. The specified real Psp(λ) and reactive Qsp(λ) power vectors can be defined as being equal to λ(kPgPgsp+kPcPcsp) and λkQcQcsp, respectively. The vectors kPg, kPc, and kQc are parameters used to characterize a specific load scenario. Using the aforementioned parameters, it is possible to simulate different variations of real and reactive power for each bus.In general, the continuation power flow (CPF) consists of a parameterization procedure, a predictor step with a step size control, and a corrector step.


### 2.1. Parameterization Techniques and the Problems Related to the Choice of the Continuation Parameter

The parameterization provides a way to identify each solution along with the trajectory to be obtained. In the local parameterization technique [7], a parameter change always occurs close to MLP. Generally, the loading factor (λ) is the parameter initially chosen. Close to the MLP, the parameter changes to the voltage magnitude that presents the largest variations in tangent vector, and after a few points, it changes back to λ. However, as will be shown in Figures 1(a), 1(b), 2(a), 2(c), 2(d), and 2(f), the use of this technique for the automatic choice of the parameter can present difficulty, because the set of buses whose voltage magnitude can be used as the continuation parameter can be considerably constrained, particularly in systems with large number of generation buses (PV), and those which have problems of local voltage instability. In these systems, the P-V curve of most buses presents a sharp nose. The loading factor and the voltage magnitude present a simultaneous reversion in its variation tendency, that is, the noses are coincident and then both Jacobian matrices are singular at the MLP [20], as the normalized determinants (|Jλ|, |JV9|) presented in Figures 2(a) and 2(c). As stated in [21], even the arc length parameterization technique [8, 14] fails to obtain the MLP for P-V curves with sharp nose.Figure 1
IEEE 118-bus system: (a)P-V curves for the base case, (b) total real power losses (Pa) as function of loading factor (λ) for the base case, (c) P-V curves for a contingency case (outage of transmission line between the buses 11 and 13), (d) λ-Pa curve for a contingency case (outage of transmission line between the buses 11 and 13).Figure 2
IEEE 118-bus system: (a) effect of limits on theP-V curves, (b) total real power losses (Pa) as function of loading factor (λ), that is, λ-Pa curve, (c) normalized determinants, (d) voltage magnitudes at critical bus (9) and at PV buses (4, 8, and 15), (e) reactive powers generated by the PV buses, (f) voltage magnitudes at tap controlled buses, (g) transformers taps variations.The goal of the following figures is to show in detail the difficulties that are present during the choice of continuation parameter. The explanation may be helpful to better understanding of the most relevant difficulties to overcome and also to develop an efficient and automatic procedure to traceP-V curves of electric power systems. In these figures, one can observe the curvature of the P-V curve of some variables that are candidates to be used as the continuation parameter. Figure 1(a) shows that close to the MLP, at least four voltage magnitudes (V1, V8, V9, and V10) could be used as continuation parameter in the local parameterization technique, while, in Figure 1(c), only the voltage magnitude of bus 13 (V13) could be used as parameter.Figure2 shows the MLP (“A,” “B,” and “C”) of three operating conditions: without limit control of reactive power and tap, with tap limits, and with control of both limits of reactive power and tap. Figure 3 shows the pre- and postcontingency P-V curves for the outage of a transmission line. It is important to note that in general, these curves are not previously known and their curvatures may be very different from each other due to changes in operating conditions such as the limits of reactive power (Qg) of generators and synchronous condensers, the limits of transformers tap, and transmission line outages. From Figures 2(d) and 2(e), one can verify that the voltage magnitudes of buses 4 and 8 (V4 and V8) are maintained constant as long as their respective generated reactive powers (Qg4 and Qg8) are within their respective limits. When their upper limits are hit, their types are switched to PQ. Then, as the system is progressively loaded, the voltage magnitudes begin to fall and the generated reactive powers will be kept constant. Note from Figures 2(f) and 2(g) that, while the transformer tap (t8) remains within its upper and lower limits (1.05 and 0.95, resp.), its voltage magnitude (V5) is maintained constant. When it hits its tap range limit, the voltage magnitude will begin to fall. The same does not happen with the transformer tap (t51) that is kept constant at its maximum value throughout the procedure and so it practically does not regulate the voltage magnitude at the bus 37 (V37). Hence, such variables (V4, V5, V8, V37, t51, and t8) are not appropriated to trace the complete P-V curves since, in many cases, or their values remain constant over a large portion of the curve, or both Jacobian matrices are singular at MLP, or as the case of voltage magnitudes of the buses 9, 75, and 118, which are appropriated as continuation parameter at a given operating condition (Figure 3), but are not in other conditions, as for the contingency of transmission line presented at Figure 3, because there is often a coincidence of noses. The same can be said about the variables V5 and V44 shown in Figure 2(a). The angle variables can have singularities not only at the maximum loading point, but also at the top of the curve, see point U in the λ-θ9 curve shown in Figure 3(c).Figure 3
Performance of IEEE 118-bus system for the base case and for the contingency of the branch 116 (transmission line between buses 69 and 75): (a)P-V curves, (b) λ-Pa curve, (c) voltage angles.Considering all the problems aforementioned, it can become very difficult to choose among all possible parameters that would allow the complete tracing ofP-V curve. In general, an approach to define the parameter changes during the computation process will be needed. Furthermore, it may be also necessary to switch the parameter a few times during the P-V curve tracing process. Despite all such care, quite often the singularity is not removed. Several global parameterization techniques have been proposed to overcome these difficulties and provide algorithms which have good convergence characteristics and computational efficiency [8, 10, 12, 13, 21]. The techniques that use the arc length [8], reactive power of a PV bus [12], and total real power losses [10] as a parameter are some examples of global parameterization techniques. The use of these techniques is interesting when the curvature of the solution trajectory in all the analyzed systems is similar, because this characteristic will simplify the steps needed for success of the method. However, not all resulting trajectories will be similar for all systems. For example, once the reactive power equation of a PV bus is already included in a conventional power flow program, it is advantageous to use it as parameter because it only requires a simple modification, which is the substitution of a column corresponding to the new variable (λ). To obtain a new operating point, it is only necessary to specify the reactive power value (Q) of any PV bus that is within its generated reactive power bounds.  Nevertheless, the use of this parameter for obtaining the MLP was only possible for few buses of some systems. In general, the PV buses reach their limits before or very close to the MLP, as shown in Figure 2(e). Initially, both reactive powers are within their bounds and so any one of them can be used as the parameter. However, as the system approaches the MLP, but before this, both the generators hit their generated reactive power limits. So, an automatic procedure for choosing that one that is within its bounds will be needed. For some systems, this never happens; all the buses hit their generated reactive power limits before or at the MLP. Also note for these two PV  buses that the resulting trajectories are not similar. Therefore, the reactive power generated at PV  bus may not be the most appropriate parameter for obtaining the MLP.As it can be seen from Figures1(b), 1(d), 2(b), and 3(b), the total real power losses curve presents the similar and so desired curvature for all operating conditions and therefore it is a strong candidate to be used as alternative parameter. Consequently, to avoid the exchange of parameter along the P-V curve tracing, in the parameterization technique proposed in [10], the following equation is added to (2.1):
(2.3)W(θ,V,λ,μ)=μF1-F(θ,V)=0,
where F(θ,V) and F1 correspond to the total real power losses (Pa) equation and its respective value calculated at base case. The total real power losses equation is given by
(2.4)Pa(θ,V)=∑kandm∈Ωgkm(Vk2+Vm2-2VkVmcosθkm),
where Ω is the set of all network buses. As one new equation is added to the problem, λ  can now be treated as a dependent variable, while the new variable (μ) is considered as a parameter. As a consequence, to obtain a new operating point, including the  λ, it is sufficient to specify the desired amount of Pa by presetting a value for μ. For example, to μ=1, the converged solution should result in λ=1. The other points on the P-V curve can be determined through successive solutions of the system formed by (2.1) and (2.3) and adopting a fixed step size for the value of parameter μ. Remember that the modified zero-order polynomial (or trivial predictor) predictor technique uses a fixed step size for the parameter and the current solution as an estimate for the next solution [9]. However, this method has only succeeded in obtaining the MLP of small systems because often, for real larger systems, the Jacobian matrices have singularities close to the MLP. Therefore, it is not possible to obtain MLP by increasing any of these parameters (λ or μ). Once the numerical problems still prevail in the region of the MLP, in [13] has been proposed the addition of a line equation (2.5) that passes through a chosen point (λ0, Pa0) in the plane determined by λ and Pa:
(2.5)W(θ,V,λ,α)=α(λ-λ0)-(Pa(θ,V,λ)-Pa0)=0,
where α is the new parameter added to the problem. So, after computing the initial angular coefficient value (α) from the coordinates of solved PF base case and a chosen point, the P-V curve is traced by successive increments in  α. When the method fails to converge, reduce the step size. When it fails again, the coordinates of the set of lines are switched to new coordinates located at abscissa axis and that corresponds to a midpoint between the base case loading and the last loading point before the divergence. This step is essential to overcome the singularity of the modified Jacobian matrix at MLP. During this procedure, it is not necessary to switch the parameter, but only change the coordinates of the set of lines. Despite the addition of the equation has allowed the determination of the MLP of various test systems, the method has failed to determine the MLP of some large realistic systems, particularly when they have local voltage instability problems.In this work, some modifications are proposed to the method presented in [13], which enables the P-V curves to trace any electric power system. The first modification is to rewrite the total real power losses (Pa) as a function of real power generated by the slack bus (Pgs). Here the goal is to avoid the degradation of the Jacobian matrix sparsity. The derivatives of (2.5) introduce a row of nonzero elements into the augmented Jacobian matrix. On the other hand, the use of the real power at a slack bus does not affect the Jacobian matrix sparsity. Table 1 shows the size of the augmented Jacobian matrix and the total number of elements added by the proposed method and the one proposed in [13]. Note that the number of elements added by the proposed method is always much less than which is added by the proposed in [13]. Also note that these elements must be updated at each iteration need to compute each point of P-V curve.Table 1
Comparison between total number of elements added by the methods.
System
Size of matrix
Number of elements added to the matrix
Proposed method
Method proposed in [13]
300
551× 551
3
519
638
1190× 1190
3
1036
787
1477× 1477
3
1168
904
1701× 1701
19
1055The total real power losses(Pa) can be computed by the summation in the right hand side of (2.4), or it can be rewritten as a function of Pgs. The slack bus voltage angle is used as the reference for all systems buses. The slack bus is also used to balance the real and reactive power (i.e., the summation of generation and demand powers over all buses should be equal the summation of the summation of power loss over all the transmission lines) in the electric power systems [22]. Therefore, one may write the total real power generated by the slack bus s as
(2.6)Pgs(θ,V,λ)=λ∑j=1NCkPcjPcjsp-λ∑i=1i≠sNGkPgiPgjsp+∑kandm∈Ωgkm(Vk2+Vm2-2VkVmcosθkm)orPgs(θ,V,λ)=λ∑j=1NCkPcjPcjsp-λ∑i=1i≠sNGkPgiPgjsp+Pa(θ,V),
where NC and NG are the numbers of load and generation buses, respectively. Solving the last equation for Pa yields
(2.7)Pa(θ,V,λ)=λ∑i=1i≠sNGkPgiPgjsp-λ∑j=1NCkPcjPcjsp+Pgs(θ,V,λ),
which now is also function of λ.  Using (2.1), Pgs may also be written as
(2.8)Pgs(θ,V,λ)=λkPcsPcssp+Pgscal(θ,V),
where Pcssp is the real power consumed at the slack bus and Pgscal(θ,V) is calculated by
(2.9)Pgscal(θ,V)=Vs∑m=1NBVm(Gsmcos(θs-θm)+Bsmsin(θs-θm)),
which may be rewritten as
(2.10)Pa(θ,V,λ)=λ∑j=1j≠sNB(kPgjPgjsp-kPcjPcjsp)-λkPcsPcssp+Pgs(θ,V,λ).
Substitution of  (2.8) into (2.10) yields
(2.11)Pa(θ,V,λ)=λ∑j=1j≠sNB(kPgjPgjsp-kPcjPcjsp)-λkPcsPcssp+λkPcsPcssp+Pgscal(θ,V).
Finally, the equation of the total real power losses becomes
(2.12)Pa(θ,V,λ)=λ∑j=1j≠sNB(kPgjPgjsp-kPcjPcjsp)+Pgscal(θ,V).In [23], it is recommended that before the notions of continuation methods are applied, the two axes must have the same scale. This recommendation is important not only to avoid curves with poor scaling which could lead to the development of inefficient techniques in consequence of a misinterpretation of their curvatures, but also to facilitate and simplify the definition of an efficient procedure of control of step size. Despite being in per-unit, the numerical values of Pa are very different from those of loading factor. Besides, the numerical values of Pa can also present a large variation for different electric power systems. For this reason, it is proposed to normalize the total real power losses values by their base case values. Using the normalization, the values of variables λ and Pa remain within the same range of numerical values and the axes have the same scales, see Figure 4. Thus, one can use the same step size for tracing the P-V curves for all operating conditions of any electric power systems. Consequently, we divide (2.12) by its value of the base case, before substituting it in (2.5). The resulting system of equations is given by
(2.13)ΔP(θ,V,λ)=Psp(λ)-Pcal(θ,V)=λ(Pgsp-Pcsp)-Pcal(θ,V)=0ΔQ(θ,V,λ)=Qsp(λ)-Qcal(θ,V)=λQcsp-Qcal(θ,V)=0ΔR(θ,V,λ,α)=α(λ-λ0)-[Pa(θ,V,λ)-Pa0Pa1]=0,
where the parameter α is the angular coefficient of the line and Pa1 is the value of the total real power losses calculated in the base case. With the addition of this new equation, λ can be treated as a dependent variable and α is considered as an independent variable, that is, it is chosen as the continuation parameter (its value is prefixed). The necessary condition to solve the above equation will remain satisfied, while the number of unknown variables and equations remains the same, that is, while the augmented Jacobian matrix has full rank and is not singular.Figure 4
P
a as function of λ (λ-Pa curve), obtained using α as parameter.


### 2.2. The Predictor Step and the Step Size Control

From the solution of the base case (pointP(θ1, V1, Pa1, λ1) in Figure 4) obtained by using a conventional PF, the value for α is computed by
(2.14)α1=(Pa1-Pa0)/Pa1(λ1-λ0),
where point O(λ0,Pa0) is an initial chosen point. Next, the PCPF is used to compute the other solutions through successive increments (Δα) in the value of α. For α=α1+Δα, the solution of (2.13) will provide the new operating point (θ2,V2,Pa2,λ2) corresponding to the intersection of the solution trajectory (curve λ-Pa) with the line whose new angular coefficient value (α1+Δα) was specified. For α=α1, the converged solution should provide λ=1. After calculating an initial value for α, a predictor step is carried out in order to find an estimate for the next solution point of P-V curve. The tangent and the secant are the most used predictors. The tangent predictor finds the estimate by taking an appropriate step size in the direction of the tangent to the P-V curve at the current solution [7–9]. In the PCPF, the tangent vector is computed by taking the derivative of (2.13):
(2.15)[∂ΔP(θ,V,λ)∂θ∂ΔP(θ,V,λ)∂V∂ΔP(θ,V,λ)∂λ∂ΔP(θ,V,λ)∂α∂ΔQ(θ,V,λ)∂θ∂ΔQ(θ,V,λ)∂V∂ΔQ(θ,V,λ)∂λ∂ΔQ(θ,V,λ)∂α∂ΔR(θ,V,λ,α)∂θ∂ΔR(θ,V,λ,α)∂V∂ΔR(θ,V,λ,α)∂λ∂ΔR(θ,V,λ,α)∂α0001][dθdVdλdα]=[000±1],
which may also be written as
(2.16)[-∂P(θ,V)∂θ-∂P(θ,V)∂VJ-∂Q(θ,V)∂θ-∂Q(θ,V)∂VPspQsp00∂ΔR(θ,V,λ,α)∂θ∂ΔR(θ,V,λ,α)∂Vα-[∑j=1j≠sNB(Pgjsp-Pcjsp)]Pa1(λ-λ0)0001][dθdVdλdα]=[000±1],
where J is the Jacobian matrix of the conventional power flow and t=[dθTdVTdλdα]T is the tangent vector. An estimate for θe, Ve, λe, and αe for the next solution is obtained by:
(2.17)[θeVeλeαe]=[θjVjλjαj]+σ[dθdVdλdα],
where σ is a scalar that defines the predictor step size and the subscript “j” means current solution.The two secant methods most widely used have been introduced in [8, 9]: the first-order polynomial predictor that uses the current and previous solutions to estimate the next one, and the modified zero-order polynomial predictor or trivial predictor, which uses the current solution and a fixed increment in the parameter (θk, Vk, λ or α in the case of proposed method) as an estimate for the next solution. In this work, the performance of PCPF will be compared considering the tangent and the trivial predictors.Another issue considered as the most critical for the success of a continuation power flow is related to the choice and control of the step size. As discussed in [23], the selection of a single step size does not guarantee that it will always work, regardless of the characteristics of the electric system under study, as well as its operating condition. Moreover, a single step size may not be adequate even for the analysis of a single operating condition. In general, a large step-size can be used on light load condition, but on heavy load condition, it should be smaller. The step size control should be flexible to adapt to the behavior of the power system in order to obtain good global convergence. In this sense in the PCPF is proposed another important modification, which is the change of the coordinates in the set of lines to the coordinates of a midpoint (MP) located between the two last points obtained. This change is only used in case of divergence. This procedure has been shown to be sufficient for the success of the method. The changes always take place near the MLP. Although the proposed method uses a single step size to whole λ-Pa curve, the change of the coordinates of initial point to the MP introduces an automatic step size control around the MLP. This occurs because of the proximity of its coordinates with those of the MLP.


### 2.3. The Corrector Step

Since the estimate is only an approximate solution, after the predictor step, it is necessary to perform the corrector step to avoid error accumulation. In most cases, the estimate is close to the correct solution, and therefore a few iterations are needed in the corrector step to obtain a solution within a required precision. Usually the Newton method is used in the corrector step. The following equation is added to (2.13):
(2.18)y-ye=0,
where y and ye correspond to the variable selected as the continuation parameter and its estimated value, obtained by the predictor step. In case of zero-order predictor (or trivial predictor), the value of the parameter can be simply fixed at ye. The expansion of the system (2.13) in Taylor series, including only terms of first order, considering the parameter value α calculated for the base case, resulting in
(2.19)-[-∂P(θ,V)∂θ-∂P(θ,V)∂VPsp-∂Q(θ,V)∂θ-∂Q(θ,V)∂VQsp∂ΔR(θ,V,λ,α)∂θ∂ΔR(θ,V,λ,α)∂Vα-[∑j=1j≠sNB(Pgjsp-Pcjsp)]Pa1][ΔθΔVΔλ]=-[JPspQsp∂ΔR(θ,V,λ,α)∂θ∂ΔR(θ,V,λ,α)∂Vα-[∑j=1j≠sNB(Pgjsp-Pcjsp)]Pa1][ΔθΔVΔλ]=-Jm[ΔθΔVΔλ]=[ΔPΔQΔR],
where J and Jm are Jacobian matrices of the PF and the PCPF, respectively. The symbols ΔP, ΔQ, and ΔR represent the correction factors (mismatches) of respective functions in the system of (2.13). For a solved PF base case, the mismatches should be equal to zero (or at least close to zero, that is, less than the adopted tolerance). In this case, only ΔR will be different from zero due to the change in α, that is, due to the increment Δα.


### 2.4. General Procedure for Tracing theP-V Curve of Bus k

The tracing of any desiredP-V curve is carried out with the corresponding desired values of voltage magnitude and loading factor, which are stored during the procedure used for tracing the λ-Pa curve. The following steps are needed to trace the λ-Pa curve.Step 1.
After solving the base case operating point “P”(λ1, Pa1) using a conventional PF, compute, by using (2.14), the angular coefficient value (α1) of the first line that passes through the initial chosen point “O”(λ0, Pa0) and point “P” (see Figure 4).Step 2.
The other points of theλ-Pa curve are obtained by using the PCPF and applying a gradual increment (Δα) to the continuation parameter α (angular coefficient of the line), αi+1=αi+Δα.Step 3.
When the PCPF fails to find a solution, it changes the coordinates of the set of lines to the midpoint MP((λa+λb)/2, (Paa+Pab)/2) located between the last two points obtained, points “a” and “b” (see detail in Figure 4(b)). Then, it computes the remaining points of the curve using the line equation that passes through the MP and the last point “b” obtained.Step 4.
When the value ofPa (point “d” in Figure 4(b)) is less than the value of the previous point (point “c”), the coordinates of the set of lines are changed to the point “D.” Then, it completes the tracing of the upper part of the λ-Pa curve considering Δα=-Δα and the line that passes through the initial point and the last point obtained in the previous step (Figure 4(a)).Aiming to increase the efficiency of the method, only a few points of the curve are computed with the set of lines that pass through the MP point, whereas all others are computed by using the set of line through the pointD. The set of lines passing through MP point is essential for the robustness of the method, once it is needed to overcome the singularity of augmented Jacobian matrix. Additionally, its use provides an automatic step size control around the MLP, see Figure 4(b).An advantage of the PCPF is that all the systems present a similar curvature for the solution trajectory, that is, theλ-Pa curvature and the next parameter are known a priori. Since the same parameter is used for all the P-V curve, the changing of one set of lines to another requires alterations only on the value of parameter α but not in matrix structure (see (2.16)), that is, the changes do not introduce new nonzero elements in the matrix.


## 2.1. Parameterization Techniques and the Problems Related to the Choice of the Continuation Parameter

The parameterization provides a way to identify each solution along with the trajectory to be obtained. In the local parameterization technique [7], a parameter change always occurs close to MLP. Generally, the loading factor (λ) is the parameter initially chosen. Close to the MLP, the parameter changes to the voltage magnitude that presents the largest variations in tangent vector, and after a few points, it changes back to λ. However, as will be shown in Figures 1(a), 1(b), 2(a), 2(c), 2(d), and 2(f), the use of this technique for the automatic choice of the parameter can present difficulty, because the set of buses whose voltage magnitude can be used as the continuation parameter can be considerably constrained, particularly in systems with large number of generation buses (PV), and those which have problems of local voltage instability. In these systems, the P-V curve of most buses presents a sharp nose. The loading factor and the voltage magnitude present a simultaneous reversion in its variation tendency, that is, the noses are coincident and then both Jacobian matrices are singular at the MLP [20], as the normalized determinants (|Jλ|, |JV9|) presented in Figures 2(a) and 2(c). As stated in [21], even the arc length parameterization technique [8, 14] fails to obtain the MLP for P-V curves with sharp nose.Figure 1
IEEE 118-bus system: (a)P-V curves for the base case, (b) total real power losses (Pa) as function of loading factor (λ) for the base case, (c) P-V curves for a contingency case (outage of transmission line between the buses 11 and 13), (d) λ-Pa curve for a contingency case (outage of transmission line between the buses 11 and 13).Figure 2
IEEE 118-bus system: (a) effect of limits on theP-V curves, (b) total real power losses (Pa) as function of loading factor (λ), that is, λ-Pa curve, (c) normalized determinants, (d) voltage magnitudes at critical bus (9) and at PV buses (4, 8, and 15), (e) reactive powers generated by the PV buses, (f) voltage magnitudes at tap controlled buses, (g) transformers taps variations.The goal of the following figures is to show in detail the difficulties that are present during the choice of continuation parameter. The explanation may be helpful to better understanding of the most relevant difficulties to overcome and also to develop an efficient and automatic procedure to traceP-V curves of electric power systems. In these figures, one can observe the curvature of the P-V curve of some variables that are candidates to be used as the continuation parameter. Figure 1(a) shows that close to the MLP, at least four voltage magnitudes (V1, V8, V9, and V10) could be used as continuation parameter in the local parameterization technique, while, in Figure 1(c), only the voltage magnitude of bus 13 (V13) could be used as parameter.Figure2 shows the MLP (“A,” “B,” and “C”) of three operating conditions: without limit control of reactive power and tap, with tap limits, and with control of both limits of reactive power and tap. Figure 3 shows the pre- and postcontingency P-V curves for the outage of a transmission line. It is important to note that in general, these curves are not previously known and their curvatures may be very different from each other due to changes in operating conditions such as the limits of reactive power (Qg) of generators and synchronous condensers, the limits of transformers tap, and transmission line outages. From Figures 2(d) and 2(e), one can verify that the voltage magnitudes of buses 4 and 8 (V4 and V8) are maintained constant as long as their respective generated reactive powers (Qg4 and Qg8) are within their respective limits. When their upper limits are hit, their types are switched to PQ. Then, as the system is progressively loaded, the voltage magnitudes begin to fall and the generated reactive powers will be kept constant. Note from Figures 2(f) and 2(g) that, while the transformer tap (t8) remains within its upper and lower limits (1.05 and 0.95, resp.), its voltage magnitude (V5) is maintained constant. When it hits its tap range limit, the voltage magnitude will begin to fall. The same does not happen with the transformer tap (t51) that is kept constant at its maximum value throughout the procedure and so it practically does not regulate the voltage magnitude at the bus 37 (V37). Hence, such variables (V4, V5, V8, V37, t51, and t8) are not appropriated to trace the complete P-V curves since, in many cases, or their values remain constant over a large portion of the curve, or both Jacobian matrices are singular at MLP, or as the case of voltage magnitudes of the buses 9, 75, and 118, which are appropriated as continuation parameter at a given operating condition (Figure 3), but are not in other conditions, as for the contingency of transmission line presented at Figure 3, because there is often a coincidence of noses. The same can be said about the variables V5 and V44 shown in Figure 2(a). The angle variables can have singularities not only at the maximum loading point, but also at the top of the curve, see point U in the λ-θ9 curve shown in Figure 3(c).Figure 3
Performance of IEEE 118-bus system for the base case and for the contingency of the branch 116 (transmission line between buses 69 and 75): (a)P-V curves, (b) λ-Pa curve, (c) voltage angles.Considering all the problems aforementioned, it can become very difficult to choose among all possible parameters that would allow the complete tracing ofP-V curve. In general, an approach to define the parameter changes during the computation process will be needed. Furthermore, it may be also necessary to switch the parameter a few times during the P-V curve tracing process. Despite all such care, quite often the singularity is not removed. Several global parameterization techniques have been proposed to overcome these difficulties and provide algorithms which have good convergence characteristics and computational efficiency [8, 10, 12, 13, 21]. The techniques that use the arc length [8], reactive power of a PV bus [12], and total real power losses [10] as a parameter are some examples of global parameterization techniques. The use of these techniques is interesting when the curvature of the solution trajectory in all the analyzed systems is similar, because this characteristic will simplify the steps needed for success of the method. However, not all resulting trajectories will be similar for all systems. For example, once the reactive power equation of a PV bus is already included in a conventional power flow program, it is advantageous to use it as parameter because it only requires a simple modification, which is the substitution of a column corresponding to the new variable (λ). To obtain a new operating point, it is only necessary to specify the reactive power value (Q) of any PV bus that is within its generated reactive power bounds.  Nevertheless, the use of this parameter for obtaining the MLP was only possible for few buses of some systems. In general, the PV buses reach their limits before or very close to the MLP, as shown in Figure 2(e). Initially, both reactive powers are within their bounds and so any one of them can be used as the parameter. However, as the system approaches the MLP, but before this, both the generators hit their generated reactive power limits. So, an automatic procedure for choosing that one that is within its bounds will be needed. For some systems, this never happens; all the buses hit their generated reactive power limits before or at the MLP. Also note for these two PV  buses that the resulting trajectories are not similar. Therefore, the reactive power generated at PV  bus may not be the most appropriate parameter for obtaining the MLP.As it can be seen from Figures1(b), 1(d), 2(b), and 3(b), the total real power losses curve presents the similar and so desired curvature for all operating conditions and therefore it is a strong candidate to be used as alternative parameter. Consequently, to avoid the exchange of parameter along the P-V curve tracing, in the parameterization technique proposed in [10], the following equation is added to (2.1):
(2.3)W(θ,V,λ,μ)=μF1-F(θ,V)=0,
where F(θ,V) and F1 correspond to the total real power losses (Pa) equation and its respective value calculated at base case. The total real power losses equation is given by
(2.4)Pa(θ,V)=∑kandm∈Ωgkm(Vk2+Vm2-2VkVmcosθkm),
where Ω is the set of all network buses. As one new equation is added to the problem, λ  can now be treated as a dependent variable, while the new variable (μ) is considered as a parameter. As a consequence, to obtain a new operating point, including the  λ, it is sufficient to specify the desired amount of Pa by presetting a value for μ. For example, to μ=1, the converged solution should result in λ=1. The other points on the P-V curve can be determined through successive solutions of the system formed by (2.1) and (2.3) and adopting a fixed step size for the value of parameter μ. Remember that the modified zero-order polynomial (or trivial predictor) predictor technique uses a fixed step size for the parameter and the current solution as an estimate for the next solution [9]. However, this method has only succeeded in obtaining the MLP of small systems because often, for real larger systems, the Jacobian matrices have singularities close to the MLP. Therefore, it is not possible to obtain MLP by increasing any of these parameters (λ or μ). Once the numerical problems still prevail in the region of the MLP, in [13] has been proposed the addition of a line equation (2.5) that passes through a chosen point (λ0, Pa0) in the plane determined by λ and Pa:
(2.5)W(θ,V,λ,α)=α(λ-λ0)-(Pa(θ,V,λ)-Pa0)=0,
where α is the new parameter added to the problem. So, after computing the initial angular coefficient value (α) from the coordinates of solved PF base case and a chosen point, the P-V curve is traced by successive increments in  α. When the method fails to converge, reduce the step size. When it fails again, the coordinates of the set of lines are switched to new coordinates located at abscissa axis and that corresponds to a midpoint between the base case loading and the last loading point before the divergence. This step is essential to overcome the singularity of the modified Jacobian matrix at MLP. During this procedure, it is not necessary to switch the parameter, but only change the coordinates of the set of lines. Despite the addition of the equation has allowed the determination of the MLP of various test systems, the method has failed to determine the MLP of some large realistic systems, particularly when they have local voltage instability problems.In this work, some modifications are proposed to the method presented in [13], which enables the P-V curves to trace any electric power system. The first modification is to rewrite the total real power losses (Pa) as a function of real power generated by the slack bus (Pgs). Here the goal is to avoid the degradation of the Jacobian matrix sparsity. The derivatives of (2.5) introduce a row of nonzero elements into the augmented Jacobian matrix. On the other hand, the use of the real power at a slack bus does not affect the Jacobian matrix sparsity. Table 1 shows the size of the augmented Jacobian matrix and the total number of elements added by the proposed method and the one proposed in [13]. Note that the number of elements added by the proposed method is always much less than which is added by the proposed in [13]. Also note that these elements must be updated at each iteration need to compute each point of P-V curve.Table 1
Comparison between total number of elements added by the methods.
System
Size of matrix
Number of elements added to the matrix
Proposed method
Method proposed in [13]
300
551× 551
3
519
638
1190× 1190
3
1036
787
1477× 1477
3
1168
904
1701× 1701
19
1055The total real power losses(Pa) can be computed by the summation in the right hand side of (2.4), or it can be rewritten as a function of Pgs. The slack bus voltage angle is used as the reference for all systems buses. The slack bus is also used to balance the real and reactive power (i.e., the summation of generation and demand powers over all buses should be equal the summation of the summation of power loss over all the transmission lines) in the electric power systems [22]. Therefore, one may write the total real power generated by the slack bus s as
(2.6)Pgs(θ,V,λ)=λ∑j=1NCkPcjPcjsp-λ∑i=1i≠sNGkPgiPgjsp+∑kandm∈Ωgkm(Vk2+Vm2-2VkVmcosθkm)orPgs(θ,V,λ)=λ∑j=1NCkPcjPcjsp-λ∑i=1i≠sNGkPgiPgjsp+Pa(θ,V),
where NC and NG are the numbers of load and generation buses, respectively. Solving the last equation for Pa yields
(2.7)Pa(θ,V,λ)=λ∑i=1i≠sNGkPgiPgjsp-λ∑j=1NCkPcjPcjsp+Pgs(θ,V,λ),
which now is also function of λ.  Using (2.1), Pgs may also be written as
(2.8)Pgs(θ,V,λ)=λkPcsPcssp+Pgscal(θ,V),
where Pcssp is the real power consumed at the slack bus and Pgscal(θ,V) is calculated by
(2.9)Pgscal(θ,V)=Vs∑m=1NBVm(Gsmcos(θs-θm)+Bsmsin(θs-θm)),
which may be rewritten as
(2.10)Pa(θ,V,λ)=λ∑j=1j≠sNB(kPgjPgjsp-kPcjPcjsp)-λkPcsPcssp+Pgs(θ,V,λ).
Substitution of  (2.8) into (2.10) yields
(2.11)Pa(θ,V,λ)=λ∑j=1j≠sNB(kPgjPgjsp-kPcjPcjsp)-λkPcsPcssp+λkPcsPcssp+Pgscal(θ,V).
Finally, the equation of the total real power losses becomes
(2.12)Pa(θ,V,λ)=λ∑j=1j≠sNB(kPgjPgjsp-kPcjPcjsp)+Pgscal(θ,V).In [23], it is recommended that before the notions of continuation methods are applied, the two axes must have the same scale. This recommendation is important not only to avoid curves with poor scaling which could lead to the development of inefficient techniques in consequence of a misinterpretation of their curvatures, but also to facilitate and simplify the definition of an efficient procedure of control of step size. Despite being in per-unit, the numerical values of Pa are very different from those of loading factor. Besides, the numerical values of Pa can also present a large variation for different electric power systems. For this reason, it is proposed to normalize the total real power losses values by their base case values. Using the normalization, the values of variables λ and Pa remain within the same range of numerical values and the axes have the same scales, see Figure 4. Thus, one can use the same step size for tracing the P-V curves for all operating conditions of any electric power systems. Consequently, we divide (2.12) by its value of the base case, before substituting it in (2.5). The resulting system of equations is given by
(2.13)ΔP(θ,V,λ)=Psp(λ)-Pcal(θ,V)=λ(Pgsp-Pcsp)-Pcal(θ,V)=0ΔQ(θ,V,λ)=Qsp(λ)-Qcal(θ,V)=λQcsp-Qcal(θ,V)=0ΔR(θ,V,λ,α)=α(λ-λ0)-[Pa(θ,V,λ)-Pa0Pa1]=0,
where the parameter α is the angular coefficient of the line and Pa1 is the value of the total real power losses calculated in the base case. With the addition of this new equation, λ can be treated as a dependent variable and α is considered as an independent variable, that is, it is chosen as the continuation parameter (its value is prefixed). The necessary condition to solve the above equation will remain satisfied, while the number of unknown variables and equations remains the same, that is, while the augmented Jacobian matrix has full rank and is not singular.Figure 4
P
a as function of λ (λ-Pa curve), obtained using α as parameter.


## 2.2. The Predictor Step and the Step Size Control

From the solution of the base case (pointP(θ1, V1, Pa1, λ1) in Figure 4) obtained by using a conventional PF, the value for α is computed by
(2.14)α1=(Pa1-Pa0)/Pa1(λ1-λ0),
where point O(λ0,Pa0) is an initial chosen point. Next, the PCPF is used to compute the other solutions through successive increments (Δα) in the value of α. For α=α1+Δα, the solution of (2.13) will provide the new operating point (θ2,V2,Pa2,λ2) corresponding to the intersection of the solution trajectory (curve λ-Pa) with the line whose new angular coefficient value (α1+Δα) was specified. For α=α1, the converged solution should provide λ=1. After calculating an initial value for α, a predictor step is carried out in order to find an estimate for the next solution point of P-V curve. The tangent and the secant are the most used predictors. The tangent predictor finds the estimate by taking an appropriate step size in the direction of the tangent to the P-V curve at the current solution [7–9]. In the PCPF, the tangent vector is computed by taking the derivative of (2.13):
(2.15)[∂ΔP(θ,V,λ)∂θ∂ΔP(θ,V,λ)∂V∂ΔP(θ,V,λ)∂λ∂ΔP(θ,V,λ)∂α∂ΔQ(θ,V,λ)∂θ∂ΔQ(θ,V,λ)∂V∂ΔQ(θ,V,λ)∂λ∂ΔQ(θ,V,λ)∂α∂ΔR(θ,V,λ,α)∂θ∂ΔR(θ,V,λ,α)∂V∂ΔR(θ,V,λ,α)∂λ∂ΔR(θ,V,λ,α)∂α0001][dθdVdλdα]=[000±1],
which may also be written as
(2.16)[-∂P(θ,V)∂θ-∂P(θ,V)∂VJ-∂Q(θ,V)∂θ-∂Q(θ,V)∂VPspQsp00∂ΔR(θ,V,λ,α)∂θ∂ΔR(θ,V,λ,α)∂Vα-[∑j=1j≠sNB(Pgjsp-Pcjsp)]Pa1(λ-λ0)0001][dθdVdλdα]=[000±1],
where J is the Jacobian matrix of the conventional power flow and t=[dθTdVTdλdα]T is the tangent vector. An estimate for θe, Ve, λe, and αe for the next solution is obtained by:
(2.17)[θeVeλeαe]=[θjVjλjαj]+σ[dθdVdλdα],
where σ is a scalar that defines the predictor step size and the subscript “j” means current solution.The two secant methods most widely used have been introduced in [8, 9]: the first-order polynomial predictor that uses the current and previous solutions to estimate the next one, and the modified zero-order polynomial predictor or trivial predictor, which uses the current solution and a fixed increment in the parameter (θk, Vk, λ or α in the case of proposed method) as an estimate for the next solution. In this work, the performance of PCPF will be compared considering the tangent and the trivial predictors.Another issue considered as the most critical for the success of a continuation power flow is related to the choice and control of the step size. As discussed in [23], the selection of a single step size does not guarantee that it will always work, regardless of the characteristics of the electric system under study, as well as its operating condition. Moreover, a single step size may not be adequate even for the analysis of a single operating condition. In general, a large step-size can be used on light load condition, but on heavy load condition, it should be smaller. The step size control should be flexible to adapt to the behavior of the power system in order to obtain good global convergence. In this sense in the PCPF is proposed another important modification, which is the change of the coordinates in the set of lines to the coordinates of a midpoint (MP) located between the two last points obtained. This change is only used in case of divergence. This procedure has been shown to be sufficient for the success of the method. The changes always take place near the MLP. Although the proposed method uses a single step size to whole λ-Pa curve, the change of the coordinates of initial point to the MP introduces an automatic step size control around the MLP. This occurs because of the proximity of its coordinates with those of the MLP.


## 2.3. The Corrector Step

Since the estimate is only an approximate solution, after the predictor step, it is necessary to perform the corrector step to avoid error accumulation. In most cases, the estimate is close to the correct solution, and therefore a few iterations are needed in the corrector step to obtain a solution within a required precision. Usually the Newton method is used in the corrector step. The following equation is added to (2.13):
(2.18)y-ye=0,
where y and ye correspond to the variable selected as the continuation parameter and its estimated value, obtained by the predictor step. In case of zero-order predictor (or trivial predictor), the value of the parameter can be simply fixed at ye. The expansion of the system (2.13) in Taylor series, including only terms of first order, considering the parameter value α calculated for the base case, resulting in
(2.19)-[-∂P(θ,V)∂θ-∂P(θ,V)∂VPsp-∂Q(θ,V)∂θ-∂Q(θ,V)∂VQsp∂ΔR(θ,V,λ,α)∂θ∂ΔR(θ,V,λ,α)∂Vα-[∑j=1j≠sNB(Pgjsp-Pcjsp)]Pa1][ΔθΔVΔλ]=-[JPspQsp∂ΔR(θ,V,λ,α)∂θ∂ΔR(θ,V,λ,α)∂Vα-[∑j=1j≠sNB(Pgjsp-Pcjsp)]Pa1][ΔθΔVΔλ]=-Jm[ΔθΔVΔλ]=[ΔPΔQΔR],
where J and Jm are Jacobian matrices of the PF and the PCPF, respectively. The symbols ΔP, ΔQ, and ΔR represent the correction factors (mismatches) of respective functions in the system of (2.13). For a solved PF base case, the mismatches should be equal to zero (or at least close to zero, that is, less than the adopted tolerance). In this case, only ΔR will be different from zero due to the change in α, that is, due to the increment Δα.


## 2.4. General Procedure for Tracing theP-V Curve of Bus k

The tracing of any desiredP-V curve is carried out with the corresponding desired values of voltage magnitude and loading factor, which are stored during the procedure used for tracing the λ-Pa curve. The following steps are needed to trace the λ-Pa curve.Step 1.
After solving the base case operating point “P”(λ1, Pa1) using a conventional PF, compute, by using (2.14), the angular coefficient value (α1) of the first line that passes through the initial chosen point “O”(λ0, Pa0) and point “P” (see Figure 4).Step 2.
The other points of theλ-Pa curve are obtained by using the PCPF and applying a gradual increment (Δα) to the continuation parameter α (angular coefficient of the line), αi+1=αi+Δα.Step 3.
When the PCPF fails to find a solution, it changes the coordinates of the set of lines to the midpoint MP((λa+λb)/2, (Paa+Pab)/2) located between the last two points obtained, points “a” and “b” (see detail in Figure 4(b)). Then, it computes the remaining points of the curve using the line equation that passes through the MP and the last point “b” obtained.Step 4.
When the value ofPa (point “d” in Figure 4(b)) is less than the value of the previous point (point “c”), the coordinates of the set of lines are changed to the point “D.” Then, it completes the tracing of the upper part of the λ-Pa curve considering Δα=-Δα and the line that passes through the initial point and the last point obtained in the previous step (Figure 4(a)).Aiming to increase the efficiency of the method, only a few points of the curve are computed with the set of lines that pass through the MP point, whereas all others are computed by using the set of line through the pointD. The set of lines passing through MP point is essential for the robustness of the method, once it is needed to overcome the singularity of augmented Jacobian matrix. Additionally, its use provides an automatic step size control around the MLP, see Figure 4(b).An advantage of the PCPF is that all the systems present a similar curvature for the solution trajectory, that is, theλ-Pa curvature and the next parameter are known a priori. Since the same parameter is used for all the P-V curve, the changing of one set of lines to another requires alterations only on the value of parameter α but not in matrix structure (see (2.16)), that is, the changes do not introduce new nonzero elements in the matrix.


## 3. Test Results

For all tests, the tolerance adopted for mismatches was 10−4 p.u. The initial coordinates of the set of lines, point “O,” were λ0=0.0 and Pa0=0.0 p.u., that is, the origin. For the trivial predictor, a fixed step size (Δα) of 0.05 is adopted to obtain all points on the λ-Pa curve. A fixed value of 0.05 was also adopted for σ in the tangent predictor. So a step size reduction or control is not used in case the singularities associated to all parameters coincide, but only a change to the coordinates of MP. The first point of each curve is computed by a conventional PF. The reactive power limits (Q) in PV 's buses are the same used in the conventional PF. In each iteration, the reactive generations for all PV buses are compared to their respective limits. The PV bus is switched to type PQ when its generated reactive power hits its upper or lower limit. It can also be switched back to PV in future iterations. During the iterations, if its voltage magnitude value is equal or greater than the specified value, its type is switched back to PV and its voltage magnitude is fixed in the specified value. The generated reactive power is changed, and while the generated reactive power remains within its upper and lower limits, its voltage magnitude is maintained constant. The loads are modeled as constant power and the parameter λ is used to simulate real and reactive increments of load, considering constant power factor. Each load increase is followed by an equivalent increase of generation by using  λ. The purpose of the tests is to highlight the efficiency and robustness of the proposed methods to trace the P-V curve of real large and heavy load electrical power systems.Figure5 shows the results from the PCPF for the IEEE-300 bus system considering the trivial and tangent predictors. Figure 5(a) shows the λ-Pa curve and lines used during the tracing process. Figure 5(b) shows the details of the tracing process of the MLP region. Note that the loading factor (λ) and the total real power losses (Pa) show a simultaneous reversal in its variation tendency, where the noses are coincident. This characteristic is commonly found in curves of the real large electric power systems. Consequently, these variables cannot be used as parameters to obtain the MLP because the PF method will present numerical difficulties in their vicinity. It can be seen from this figure that, despite the sharp nose presented by the λ-Pa curve, the change for the coordinates of midpoint allows overcoming the singularity of the Jacobian matrix at the MLP.Performance of the PCPF for the IEEE-300: (a)λ-Pa curve, (b) detail of the convergence process, (c) P-V curve of the critical bus, (d) number of iterations for the trivial and tangent predictors, (e) convergence pattern for point “b,” and (f) convergence pattern for the next point (dashed line “S”) for which the iterative process does not converge.
(a)
(b)
(c)
(d)
(e)
(f)At the MLP, the values ofλ and Pa found by the proposed method were 1.0553 p.u. and 1.1624 p.u, respectively. Figure 5(c) shows the P-V curve of critical bus (V526), whose corresponding voltage magnitudes values were stored while obtaining the λ-Pa curve. The values of the λ and V526 at the MLP were 1.0553 and 0.7302, respectively. Figure 5(d) presents a comparison of the number of iterations for convergence of each point of λ-Pa curve for the trivial and tangent predictors. As it can be seen from Figure 5(d), the proposed method succeeds in finding each point of the curve, including the MLP, with a small number of iterations. Also note that, even using a fixed step size (Δα), an automatic step size control occurs around the MLP as a consequence of the proximity between the midpoint coordinates and the curvature of the solution trajectory (λ-Pa curve).


### 3.1. Analysis of the Total Mismatch

Among several possible candidates for convergence criterion of an iterative procedure, the simplest is the one that uses a predefined number of iterations. Despite the simplicity of programming, it is unknown how close the solution is. Moreover, the analysis of the total mismatch behavior is a good indicator of the possibility of ill-conditioning of the Jacobian matrix. The total mismatch is defined as the sum of the absolute values of the real and reactive power mismatches. Therefore, the analysis of the total mismatch behavior associated with a predefined number of iterations is the criterion used to change the coordinates of the set of lines to the MP. This criterion prevents the process spending too much time on cases that do not converge or diverge. Figures5(e) and 5(f) present the evolution of the mismatches for two points (see Figures 5(a) and 5(b)), the point “b” and the next point (dashed line “S”) corresponding with the last point before the divergence happens. Figure 5(e) depicts the evolution of the mismatch for the last point “b” obtained using the first set of line. Note from it that the convergence occurs in only four iterations using a tolerance of 10−4 p.u. However, for dashed line “S,” the process shows an oscillatory behavior and slow divergence, as shown in Figure 5(f). This occurs because, as it can be seen from Figure 5(a), there is no intersection of the dashed line (S) and the λ-Pa curve, and therefore the problem actually has no solution. It can be seen that, after the third iteration, it is already possible to verify this behavior. Thus, it is proposed to compare the magnitudes of the last two mismatches after the third iteration. If the last value is higher than the previous one, then the iterative process is interrupted and the λ-Pa curve tracing is resumed from the previous point. In this case, it changes the coordinates of the set of lines to the midpoint (MP), see Figure 5(b). Otherwise, if the last value is smaller than the previous one, then the iterative process continues. Table 2 shows the number of iterations needed to carry out the change of coordinates of the set of lines to the MP. It can be seen that the number of iterations is always less than ten.Table 2
Number of iterations needed to change the coordinates in the set of lines to MP through total mismatch criterion.
Predictor
System
P1
P2
Trivial
300
4
4
638
5
4
787
5
5
904
7
6
300(1)
4
7
Tangent
300
5
4
638
4
4
787
5
5
904
6
7
Without reactive limits
(
1
).Figure6 shows the performance of the PCPF for the real large systems: 638 and 787-bus systems corresponding to part of South-Southeast Brazilian system. Figure 6(a) shows the λ-Pa curve for both systems, and Figure 6(b) presents the P-V curves of critical bus (V150 and V576) of each system obtained by storing, during the tracing of the curve λ-Pa, the corresponding values of voltage magnitudes. In Figure 6(c), a comparison of the number of iterations needed for each point of λ-Pa curve by using the trivial and tangent predictors. In general, the number of iterations for tangent predictor is less than the trivial predictor is shown. The proposed method showed a good overall performance, as it can be confirmed from the low number of iterations needed for most points.Figure 6
Performance of the PCPF for the South-Southeast Brazilian systems: (a)λ-Pa curves, (b) voltage magnitudes at the critical buses (150 and 576) as function of  λ, (c) number of iterations for tangent and trivial predictors.Figure7 shows the results of the PCPF for the 904-bus Southwestern American system. It is a large, heavy load, with 904 buses and 1283 branches. From Figure 7(a), which illustrates the P-V curves of all buses of the system, one can see that this is a system with voltage instability problems that have the strong local characteristic. Note that, except the P-V curve of critical bus (V138), all the others present a sharp nose, or a straight line, that is, their voltage magnitudes remain fixed at the specified values. Figure 7(b) shows the λ-Pa curve, while Figure 7(c) presents the P-V curve of critical bus. In Figure 7(c), the numbers of iterations needed to obtain the curve λ-Pa by trivial predictor are shown.Performance of the PCPF for the 904-bus Southwestern American system: (a)P-V curves, (b) λ-Pa curve, (c) voltage magnitude at the critical bus (138) as function of λ, (d) number of iterations.
(a)
(b)
(c)
(d)


## 3.1. Analysis of the Total Mismatch

Among several possible candidates for convergence criterion of an iterative procedure, the simplest is the one that uses a predefined number of iterations. Despite the simplicity of programming, it is unknown how close the solution is. Moreover, the analysis of the total mismatch behavior is a good indicator of the possibility of ill-conditioning of the Jacobian matrix. The total mismatch is defined as the sum of the absolute values of the real and reactive power mismatches. Therefore, the analysis of the total mismatch behavior associated with a predefined number of iterations is the criterion used to change the coordinates of the set of lines to the MP. This criterion prevents the process spending too much time on cases that do not converge or diverge. Figures5(e) and 5(f) present the evolution of the mismatches for two points (see Figures 5(a) and 5(b)), the point “b” and the next point (dashed line “S”) corresponding with the last point before the divergence happens. Figure 5(e) depicts the evolution of the mismatch for the last point “b” obtained using the first set of line. Note from it that the convergence occurs in only four iterations using a tolerance of 10−4 p.u. However, for dashed line “S,” the process shows an oscillatory behavior and slow divergence, as shown in Figure 5(f). This occurs because, as it can be seen from Figure 5(a), there is no intersection of the dashed line (S) and the λ-Pa curve, and therefore the problem actually has no solution. It can be seen that, after the third iteration, it is already possible to verify this behavior. Thus, it is proposed to compare the magnitudes of the last two mismatches after the third iteration. If the last value is higher than the previous one, then the iterative process is interrupted and the λ-Pa curve tracing is resumed from the previous point. In this case, it changes the coordinates of the set of lines to the midpoint (MP), see Figure 5(b). Otherwise, if the last value is smaller than the previous one, then the iterative process continues. Table 2 shows the number of iterations needed to carry out the change of coordinates of the set of lines to the MP. It can be seen that the number of iterations is always less than ten.Table 2
Number of iterations needed to change the coordinates in the set of lines to MP through total mismatch criterion.
Predictor
System
P1
P2
Trivial
300
4
4
638
5
4
787
5
5
904
7
6
300(1)
4
7
Tangent
300
5
4
638
4
4
787
5
5
904
6
7
Without reactive limits
(
1
).Figure6 shows the performance of the PCPF for the real large systems: 638 and 787-bus systems corresponding to part of South-Southeast Brazilian system. Figure 6(a) shows the λ-Pa curve for both systems, and Figure 6(b) presents the P-V curves of critical bus (V150 and V576) of each system obtained by storing, during the tracing of the curve λ-Pa, the corresponding values of voltage magnitudes. In Figure 6(c), a comparison of the number of iterations needed for each point of λ-Pa curve by using the trivial and tangent predictors. In general, the number of iterations for tangent predictor is less than the trivial predictor is shown. The proposed method showed a good overall performance, as it can be confirmed from the low number of iterations needed for most points.Figure 6
Performance of the PCPF for the South-Southeast Brazilian systems: (a)λ-Pa curves, (b) voltage magnitudes at the critical buses (150 and 576) as function of  λ, (c) number of iterations for tangent and trivial predictors.Figure7 shows the results of the PCPF for the 904-bus Southwestern American system. It is a large, heavy load, with 904 buses and 1283 branches. From Figure 7(a), which illustrates the P-V curves of all buses of the system, one can see that this is a system with voltage instability problems that have the strong local characteristic. Note that, except the P-V curve of critical bus (V138), all the others present a sharp nose, or a straight line, that is, their voltage magnitudes remain fixed at the specified values. Figure 7(b) shows the λ-Pa curve, while Figure 7(c) presents the P-V curve of critical bus. In Figure 7(c), the numbers of iterations needed to obtain the curve λ-Pa by trivial predictor are shown.Performance of the PCPF for the 904-bus Southwestern American system: (a)P-V curves, (b) λ-Pa curve, (c) voltage magnitude at the critical bus (138) as function of λ, (d) number of iterations.
(a)
(b)
(c)
(d)


## 4. Performance of the PCPF by Using Constant Jacobian

The robustness and effectiveness are some of the main features required for a CPF that is used in the steady-state voltage stability analysis. In these cases, the Newton-Raphson algorithm is the most appropriate one. From the analysis presented at the previous section, it is verified that the PCPF method is robust and effective for determining the MLP and the completeP-V curve. On the other hand, in these analyses, it is often necessary to evaluate the loading margin of the system for a large number of configurations and of operating conditions, and consequently it is also necessary to have computational efficiency.Usually, the CPF needs a few iterations to obtain each point on theP-V because the iterative process begins with a good initial estimate of the solution. For the tangent and first-order polynomial secant predictors, the next operation point is obtained from an estimate and from a previous point for the modified zero-order polynomial predictor. For each point on the curve, successive updates and inversions of the Jacobian matrix are required. Besides, the convergence of the power flow methods is also affected by adjustment of solutions due to, for example, reactive limit violations. Therefore, the CPU time required to get all points of the curve may be very high. For planning applications, this computing demand can be acceptable, but not for system control applications where, for example, the effects of several outages on the system loading margin must be available in real time [24].In order to reduce the computational burden, several updating procedures have been proposed [22, 24–27]. A procedure known as the “dishonest Newton method” is commonly used. In this procedure, the Jacobian matrix is not updated at each iteration, but only occasionally. Sometimes, it is proposed to make its update only once, that is, at the first iteration. From several studies conducted, it is concluded that the Jacobian matrix is important for the convergence of the process but does not influence the final solution, because at each iteration the function value is accurately calculated, despite the corrections are approximated. Generally, a relatively small increase in the number of iterations is the only effect observed. Therefore, it is possible to use approximated values for Jacobian matrix without losing the global convergence.It has been proposed to update the Jacobian matrix only when the system undergoes a significant change, for example, when a voltage controlled bus (PV) is converted to a load bus (PQ) as a consequence of a violation of one of their reactive limits, or when the number of iterations exceeds a predefined threshold value. So, in the first procedure (P1), updating the whole Jacobian matrix is performed at every iteration and, in the second one (P2), only when the system undergoes a significant change. When the limits of reactive powers are not taken into account, the Jacobian matrix updating is performed at every iteration in procedure P1, and P2 only when the number of iterations exceeds the seventh one.Table3 shows for the two procedures the loading factor at the MLP (λmax, at 3th and 5th columns) and the corresponding voltage magnitude of the critical bus (at 4th and 6th columns), computed by the trivial and tangent predictors. Note that the values obtained by each of the procedures are practically the same. For the IEEE-300, the table also shows that the reactive power limits have significant influence on the value of  λmax.Table 3
Maximum loading point (λmax) and critical voltage of the analyzed systems.
Predictor
System
P1
P2
λ
max
Critical voltage (p.u.)
λ
max
Critical voltage (p.u.)
Trivial
300
1.0553
0.7302
1.0553
0.7305
638
1.0080
0.8086
1.0080
0.8086
787
1.1273
0.7465
1.1273
0.7465
904
1.1993
0.6508
1.1993
0.6508
300(1)
1.4124
0.4643
1.4125
0.4663
Tangent
300
1.0553
0.7302
1.0553
0.7304
638
1.0080
0.8086
1.0080
0.8086
787
1.1273
0.7464
1.1273
0.7466
904
1.1993
0.6507
1.1993
0.6507
without reactive limits
(
1
).The results presented in Tables4 and 5 allow us to compare the performance of PCPF for both procedures, considering the trivial and the tangent predictors. Table 4 shows for both procedures the total number of iterations (IC) to trace the complete P-V curve, and, for P2, the total number of iterations (ACo) for which there is a matrix updating. The computational times (CPU time) required for procedures P1 and P2 are shown at the fourth and seventh columns, respectively. Their values were normalized by the respective CPU times of the procedure P1. Although the procedure P2 requires a larger total number of iterations than the P1, it requires less computational time, as shown at the seventh column. As it can be confirmed at the eighth column, the procedure P2 presents a better performance than the P1 for both predictors. Therefore, it is possible to obtain an overall CPU time reduction without losing robustness. The efficiency improvement is achieved by a simple change in the procedure, which is to update the Jacobian matrix only when the system undergoes a significant change.Table 4
Performance of the PCPF for procedures P1 and P2.
Predictor
System
P1
P2
CPU ratio [%]
IC
CPU time (p.u.)
IC
ACo
CPU time (p.u.)
Trivial
300
50
1.000
59
23
0.435
56.45
638
86
1.000
147
18
0.509
49.13
787
72
1.000
106
22
0.494
50.61
904
115
1.000
126
46
0.413
58.73
300(1)
151
1.000
378
23
0.326
67.38
Tangent
300
40
1.000
47
18
0.421
57.87
638
46
1.000
103
18
0.572
42.78
787
53
1.000
93
16
0.485
51.46
904
93
1.000
114
53
0.416
58.40
without reactive limits
(
1
), IC: iteration count, ACo: actualization count.Table 5
Comparison between PCPF considering the trivial predictor and tangent predictor.
Procedure
System
Tangent predictor
Trivial predictor
CPU ratio [%]
CPU time (p.u.)
CPU time (p.u.)
P1
300
1.000
0.857
14.22
638
1.000
0.862
13.84
787
1.000
0.864
13.54
904
1.000
0.866
13.35
P2
300
1.000
0.886
11.32
638
1.000
0.766
23.39
787
1.000
0.879
12.02
904
1.000
0.859
14.03Table5 compares the overall CPU times needed for the tangent and trivial predictors. As it can be seen in the last column, the overall CPU time of tangent predictor is higher than the trivial predictor.


## 5. Conclusion

This paper shows the most relevant difficulties that are present during the choice of continuation parameter to trace theP-V curve in steady-state voltage stability analysis. It also presents significant modifications for the continuation method proposed by Garbelini et al. in [13], which was developed from the geometrical behavior of the solutions trajectories of the power flow equations. The proposed modifications not only allow obtaining the maximum loading point (MLP) and, subsequently, assessment of voltage stability margin, but also make possible the complete tracing of P-V curve of any power systems with a low number of iterations, including those which have problems with local voltage instability. The PCPF combines robustness with simplicity and ease of interpretation. It also shows a very attractive option and easy computational implementation, since it only requires few changes on power flow program currently in use.Theλ-Pa curve has been chosen because it has a behavior (curvature) similar to all the systems, which greatly facilitates defining an overall procedure able to trace the P-V curves of any systems. Consequently, the parameterization technique proposes the addition of a line equation which passes through a point in the plane determined by the variables loading factor (λ) and the total real power losses (Pa). The angular coefficient of the line (α) is used as parameter for tracing the whole P-V curve. So it is not necessary to switch the parameter, but only change the coordinates of the set of lines. The changing of one set of lines to another requires only alterations on the value of parameter α but not in matrix structure. The degradation of the Jacobian matrix sparsity is avoided by rewriting the total real power losses (Pa) equation as a function of real power generated by the slack bus. This modification provides a smaller number of elements than the method proposed in [13]. Another proposed modification is the normalization of the total real power losses values by its base case value, which makes possible to use the same step size for tracing the P-V curves for all operating conditions of all electric power systems. The proposal of using the coordinates of the midpoint located close to the MLP is essential for the robustness of the method, once it is needed to overcome the singularity of augmented Jacobian matrix and all the consequent problems of ill-conditioning. Furthermore, it provides the advantage of an automatic step size control around the MLP, despite of using a single step size for tracing the whole P-V curve. The changing of coordinates is based on the analysis of the total mismatch behavior. This criterion leads to a low overall number of iterations.To reduce the computational burden, it is also investigated to update the Jacobian matrix only when the system undergoes a significant change (changes in the system’s topology). This simple change of procedure increases the efficiency of the proposed technique and proves that it is possible to obtain a reduction in computational time without losing robustness. The results confirm the efficiency of the proposed method, including its application feasibility in studies related with the assessment of static voltage stability.



---

*Source: 762371-2012-10-01.xml*