It is well known that the behavior of PV curves is similar to a quadratic function. This is used in some papers to approximate PV curves and calculate the maximumloading point by minimum number of power flow runs. This paper also based on quadratic approximation of the PV curves is aimed at completing previous works so that the computational efforts are reduced and the accuracy is maintained. To do this, an iterative method based on a quadratic function with two constant coefficients, instead of the three ones, is used. This simplifies the calculation of the quadratic function. In each iteration, to prevent the calculations from diverging, the equations are solved on the assumption that voltage magnitude at a selected load bus is known and the loading factor is unknown instead. The voltage magnitude except in the first iteration is selected equal to the one at the nose point of the latest approximated PV curve. A method is presented to put the mentioned voltage in the first iteration as close as possible to the collapse point voltage. This reduces the number of iterations needed to determine the maximumloading point. This method is tested on four IEEE test systems.
1. Introduction
The quick and accurate determination of voltage collapse proximity for normal and contingency states is needed for secure operation of power systems. Methods have been proposed in the literature for voltage stability analysis can be classified into measurementbased methods
[1
,
2]
and modelbased ones
[3

16]
. The formers are based on the fact that maximumloading point is detectable using the local measurement of voltage and current phasors. Since the system model is not used, the measurementbased methods are suitable for the analysis of the current power system only, and are unusable for the examination of contingency states.
The modelbased methods use dynamic and static techniques. The formers are time consuming and unsuitable when the investigation of a wide range of system conditions and a large number of contingencies is required
[3]
. This paper focuses on the static techniques that use equilibrium equations especially the power flow equations instead of differential equations. The main objective of too many papers that are based on the static methods is to accurately calculate the maximumloading point or present indices that can provide a measure of the MWdistance from the current operating point to the maximumloading point. Some information embedded in the Jacobian and the Hessian matrix of the load flow equations, such as the minimum singular value and the minimum magnitude of eigenvalue, are the commonly used voltage stability indices
[4

6]
. The Lindex
[7]
is one other wellknown index that uses the admittance matrix and voltage phasors. In
[8]
, a PQ curve that shows the boundary of the maximum active / reactive power is used to define a new voltage stability index. This index provides information about identifying weak buses and areas. Another method to identify weak buses has been proposed in
[9]
based on creating a set of decoupled singleone, singlebranch equivalent circuits. The decoupled circuits can reveal important characteristics of a power system. None of voltage stability indices can accurately estimate the proximity to voltage collapse because of their nonlinear behavior.
Besides the voltage stability indices, many approaches have been proposed to accurately compute the distance to the maximumloading point as the voltage stability margin. A large number of them are based on iterative methods. The authors in
[10
,
11]
propose an iterative method to calculate the shortest distance from the current operating point to loadability boundary. In each iteration, the direction of load increase is determined using the normal vector on the loadability boundary obtained in the previous iteration.
In other methods, unlike the methods proposed in
[10
,
11]
, the direction of load increase in each load bus is usually assumed to be constant and proportional to the initial load
[12

16]
. The focus is on the determination of the step size of load increase in each iteration. In some proposed methods, the step size is so determined that no variables such as reactive power generations and load bus voltages reach their limits. This may be calculated using the sensitivity of these variables with respect to loading factor
[12
,
13]
. In
[14
,
15]
, the authors assume a quadratic relation between loading factor and some variables that are zero in the maximumloading point. These variables are the inverse of the sensitivities of the voltage magnitude at PQ buses
[14]
and the reactive power generated at PV buses
[15]
with respect to loading factor. The step size of load increase is determined based on this quadratic relation. The coefficients of quadratic functions are successively calculated using the two last power flow solutions. In each iteration, the loading factor is selected equal to the critical loading factor predicted from the latest calculated quadratic function.
It is well known that PV curves can be approximated by a quadratic function. This is used in
[16]
for quick calculation of the maximum loading point. The three coefficients of the quadratic function are determined by only one power flow solution using first and secondorder derivatives of voltage magnitude with respect to loading factor. The loading factor in the nose point of the approximated PV curve is selected as the loading factor in the next power flow run. The iteration continues until the difference between the maximum loading points calculated in two successive iterations becomes smaller than a given tolerance. Among the presented methods, the one proposed in
[16]
needs the least number of power flow solutions to calculate the maximum loading point.
It is important to note that voltage stability margin must be calculated for a large number of contingency states. This work needs to be repeated every some minutes because of change in system’s operating point. So a very limited time may be available for each contingency. This is why the voltage stability margin needs to be determined as fast as possible. The main objective of this paper is to improve the efficiency and reduce the computational effort and the number of iterations of the method proposed in
[16]
. To do this, a quadratic function with two constant coefficients, instead of the three ones, is used. This simplifies the calculation of the quadratic function because the determination of the secondorder sensitivities is no longer needed. In each iteration, to prevent the calculations from diverging, the equations are solved on the assumption that voltage magnitude at a selected load bus is known and the loading factor is unknown instead. The voltage magnitude except in the first iteration is selected equal to the one at the nose point of the latest approximated PV curve. A method is presented to put the mentioned voltage in the first iteration as close as possible to the collapse point voltage. This reduces the number of iterations needed to determine the maximum loading point. The simulations have been performed on the IEEE 14, 30, 57 and 118bus test systems. The results show that the voltage stability margin can be accurately calculated after two to five iterations of power flow run. The iterations stop when the difference between the maximum loading points calculated in two successive iterations becomes less than 0.0005.
2. A Brief Description of the Method Proposed in[16]
The method is described using
Fig. 1
that is a typical PV curve for a generic bus (bus
i
). In this Fig.,
λ_{c}
is the critical loading factor (the loading factor in the maximumloading point),
V
_{c,i}
is the collapse point voltage at bus
i
(the voltage magnitude at bus
i
in the critical loading factor),
λ_{co}
and
V
_{co,i}
are the loading factor and the voltage magnitude in the current operating point, respectively. In
[16]
, a quadratic approximation for the PV curve has been defined as:
Typical PV curve
The coefficients
a_{i}
,
b_{i}
and
c_{i}
are determined by only one power flow solution in each loading factor. To do this, the firstorder and secondorder derivatives
dV_{i}
/
dλ
and
d^{2}V_{i}
/
dλ^{2}
are calculated from the power flow equations and substituted in Eqs. (2) and (3):
The coefficients
a_{i}
,
b_{i}
and
c_{i}
can be found by solving Eqs. (1), (2) and (3) as:
After the determination of the coefficients
a_{i}
,
b_{i}
and
c_{i}
, the estimated collapse point voltage
V
_{c,i}
and the estimated critical loading factor
λ
_{c,i}
are obtained as follows:
The critical bus is the bus with the least value of
λ
_{c,i}
. This
λ
_{c,i}
is an estimate of the critical loading factor. In the next step, the loading factor is replaced by
λ
_{c,i}
and the entire process is carried out again (i.e. first the coefficients
a_{i}
,
b_{i}
and
c_{i}
are calculated again, then
V
_{c,i}
and
λ
_{c,i}
are determined). The iteration continues until the difference between the critical loading factors calculated in two successive iterations becomes smaller than a given tolerance. The method needs 5 iterations for the normal loading and far fewer iterations for the stressed loading (based on
Figs. 2
and
3
in
[16]
). The closer the initial loading factor is to the critical loading factor, the fewer iterations are needed. The paper does not propose any method to bring the initial loading factor close to the critical one. Also, it is not clear what decision is made when the power flow calculations fail to converge. This occurs if the estimated critical loading factor becomes larger than the actual one.
Collapse point voltages at different contingency states in 14bus test system
Collapse point voltages at different contingency states in 30bus test system
3. The Proposed Method
This paper, also, uses a quadratic approximation of the PV curve. The aim is to calculate the critical loading factor for a set of contingency states with minimum number of power flow runs. The normal state is numbered as contingency state 1. Other contingencies are numbered arbitrarily from 1 to
N_{c}
+ 1.
N_{c}
is the number of contingency states. The proposed quadratic function is as:
i
refers to the selected bus number (the bus where the PV curve is drawn). The loading factor
λ
is a scalar factor used to move system operating point to the loadability boundary. This is included in the power flow equations as follows:
where
P_{i0}
(
Q_{i0}
) is the netinjected active (reactive) power at bus
i
when
λ
= 1. An iterative method is used for each contingency state. In each iteration, the two constant coefficients of the quadratic function are calculated. Since only two coefficients exist, the determination of the secondorder derivative
d^{2}V_{i}
/
dλ^{2}
is no longer needed. Several simulations confirm that the proposed quadratic function can nicely approximate PV curves. In all iterations, to prevent the power flow calculations from diverging, it is assumed that the voltage magnitude at the considered bus is known and loading factor is unknown instead. The voltage magnitude is selected equal to the latest estimate of the collapse point voltage obtained for that contingency state. The initial voltage magnitude (the voltage magnitude at the considered bus in the first iteration) is selected based on the variation range of the actual collapse point voltage. The magnitude of this voltage depends on the level of reactive power compensation
[17]
. The results of a large number of simulations performed on the IEEE 14, 30, 57 and 118bus test systems for contingency states show that the collapse point voltage is within a limited range for most of the states. The results for some load buses are shown in
Figs. 2
to
5
. The similar results are obtained for other load buses.
Collapse point voltages at different contingency states in 57bus test system
Collapse point voltages at different contingency states in 118bus test system
The contingencies are single line or generator outages. Considering the normal state, the number of selected contingencies for the IEEE 14, 30, 57 and 118bus test systems are 21, 41, 76, and 228, respectively. As can be seen from the Figs., most of the calculated collapse point voltages in the considered buses range from 0.57 to 0.62 pu in the 14 and 30bus test systems, from 0.47 to 0.52 pu in the 57bus test system, and from 0.77 to 0.82 pu in the 118bus test system.
Therefore, after determining the actual collapse point voltage for few contingency states, the range of variations in this voltage for other contingencies can be guessed. The selection of the initial voltage magnitude in the within this range put the initial operating point close to the nose point of the PV curve, may be on the upper part or lower part of the PV curve. This reduces the number of iterations needed to determine the maximum loading point. It is proposed that for each contingency state, the initial voltage magnitude is selected equal to the average value of the collapse point voltage calculated for the previous contingency states. For the first contingency state, i.e. the normal state, the initial voltage magnitude is selected equal to 0.8 pu.
To determine the coefficients
a_{i}
and
b_{i}
, the firstorder derivative
dV_{i}
/
dλ
is obtained using (9) as follow:
In each iteration, with a given
V_{i}
,
λ
and
dV_{i}
/
dλ
is calculated by solving the power flow equations. Having
V_{i}
,
λ
and
dV_{i}
/
dλ
, the coefficients
a_{i}
and
b_{i}
can be found from Eqs. (9) and (11) as:
After determining the coefficients
a_{i}
and
b_{i}
, the estimated collapse point voltage
V
_{c,i}
and the estimated critical loading factor
λ_{c}
are obtained as follows:
To explain how to calculate
dV_{i}
/
dλ
, it is assumed that the power flow equations are as the compact form:
where
X
is the systemstate vector (voltage magnitudes and phase angles) and
λ
is the loading factor. A Taylor series expansion of Eq. (16) provides the relation:
where
g_{X}
is the Jacobian matrix and
g_{λ}
is the vector of derivatives of load flow equations with respect to
λ
. Therefore, the sensitivity of the state variables
X
with respect to the loading factor
λ
is given by:
The determination of the critical loading factor can be done at each of weak load buses. In each iteration, the main computational effort is related to the modified power flow solution with voltage magnitude at one of load buses being known and the loading factor being unknown. Other variables, such as the sensitivity vector in Eq. (18), are easily available from the power flow solution. So, the number of iterations means the number of power flow solutions. The steps of the proposed method can be described as bellow:
1) Select one of weak load buses (based on voltage stability indices or the voltage magnitude in the base load). The determination of the weakest load bus is not required.
2) Number the selected contingency states from 1 to
N_{c}
+ 1. The contingency state 1 refers to the normal state. Assume that
k
represents the contingency number.
3) For the normal state (
k
= 1 ), select the initial voltage magnitude at the considered load bus equal to 0.8 pu.
4) Perform a power flow solution to calculate the loading factor. The calculated loading factor is considered as an estimate of critical loading factor.
5) Calculate
dV_{i}
/
dλ
using (18).
6) Calculate
a_{i}
and
b_{i}
using (12) and (13).
7) Calculate the new estimated collapse point voltage
V
_{c,i}
and the new estimated critical loading factor
λ_{c}
using (14) and (15). If the difference between the new estimated critical loading factor and the last one is less than 0.0005, the critical loading factor for the normal state has been determined. In this case, go to step 8 for the determination of the critical loading factor for other contingency states. If the difference is more than 0.0005, replace the initial voltage magnitude in the considered load bus by the latest estimated collapse point voltage and repeat the steps 4 to 7 for next iteration. The iterations continue until the difference between the estimated critical loading factors obtained in two successive iterations becomes less than 0.0005 (Note that, the critical loading factor obtained by the proposed method is the same as determined by the continuation power flow method
[18]
that is known as an accurate method to determine voltage stability margin. In this method load is changed along a predetermined direction to derive the system’s operating point from base case to maximumloading point. The continuation power flow method needs several power flow runs and is not suitable when the examination of a large number of contingencies is needed. This is shown in the next section).
8) Set
k
=
k
+1 , For contingency state
k
, select the initial voltage magnitude equal to the average value of the collapse point voltage calculated for contingency states 1 to
k
1.
9) Perform a power flow solution to calculate the loading factor. The calculated loading factor is considered as an estimate of critical loading factor for contingency state
k
.
10) Calculate
dV_{i}
/
dλ
using (18).
11) Calculate
a_{i}
and
b_{i}
using (12) and (13).
12) Calculate the new estimated collapse point voltage
V
_{c,i}
and the new estimated critical loading factor
λ_{c}
using (14) and (15). If the difference between the new estimated critical loading factor and the last one is less than 0.0005, the critical loading factor for the contingency state has been determined. In this case, go back to step 8 for the determination of the critical loading factor for other contingency states. If the difference is more than 0.0005, replace the initial voltage magnitude in the considered load bus by the latest estimated collapse point voltage and repeat the steps 9 to 12 for next iteration. The iterations continue until the difference between the estimated critical loading factors obtained in two successive iterations becomes less than 0.0005.
For better realization of the proposed method, the flowchart of the above mentioned steps is presented in
Fig. 6
.
The flowchart of the proposed method
4. Simulation Results
The simulations have been performed for normal and contingency states on the IEEE 14, 30, 57 and 118bus test systems. The contingencies are single line or generator outages. Considering the normal state, the number of selected contingencies for the IEEE 14, 30, 57 and 118bus test systems are 21, 41, 76, and 228, respectively. The generator reactive power limits have been enforced.
Figs. 7
to
10
show the performance of the proposed method for the four test systems. The figures show the number of iterations needed for determination of the critical loading factor in each contingency state.
The number of required iterations at different contingency states in 14bus test system
The number of required iterations at different contingency states in 30bus test system
The number of required iterations at different contingency states in 57bus test system
The number of required iterations at different contingency states in 118bus test system
Table 1
shows the results in percentage terms. The results confirm that the proposed method considerably reduces the number of required iterations. For instance, in 118bus test system, the critical loading factor can be determined by only two iterations for more than 90 percent of contingency states.
The results in percentage terms
The results in percentage terms
The reduction in required iterations is because of putting initial loading factor close to the critical one.
Figs. 11
to
14
show the initial and critical loading factors for different contingency states. In these figures, ‘*’ denotes the initial loading factor and ‘○’ denotes the critical one. As can be seen, most of the initial loading factors are very close to the critical ones.
The initial and critical loading factors for different contingency states in 14bus test system
The initial and critical loading factors for different contingency states in 30bus test system
The initial and critical loading factors for different contingency states in 57bus test system
The initial and critical loading factors for different contingency states in 118bus test system
As mentioned in the last section, all the critical loading factors calculated using the proposed method are equal to the ones obtained by the continuation power flow method. Due to the limited space in the paper it is not possible to show all results, but for some contingency states, the results of both methods are tabulated in
Table 2
. The continuation power flow method needs several iterations and is not suitable when the examination of a large number of contingencies is needed. The proposed method considerably reduces the needed number of iterations. This is shown in
Table 3
.
The critical loading factors at some contingency states, calculated by the proposed method and the continuation power flow method
The critical loading factors at some contingency states, calculated by the proposed method and the continuation power flow method
The comparison of the needed number of iterations in the proposed method and the continuation power flow method
The comparison of the needed number of iterations in the proposed method and the continuation power flow method
5. Conclusion
The aim of this paper is to calculate the critical loading factor for a set of contingency states with minimum number of power flow runs. The method is based on quadratic approximation of the PV curves. An iterative method is used for each contingency state. In each iteration, the two constant coefficients of the quadratic function are calculated. In all iterations, to prevent the power flow calculations from diverging, it is assumed that the voltage magnitude at the considered bus (the bus where the PV curve is drawn) is known and loading factor is unknown instead. The voltage magnitude is selected equal to the latest estimate of the collapse point voltage obtained for that contingency state. For each contingency state, the initial voltage magnitude is selected equal to the average value of the collapse point voltage calculated for the previous contingency states. This selection put the initial operating point close to the nose point of the PV curve, and consequently reduces the number of iterations needed to determine the critical loading factor. This is confirmed by the simulations performed on the IEEE 14, 30, 57 and 118bus test systems. For instance, in 118bus test system, the critical loading factor can be determined by only two power flow runs for more than 90 percent of contingency states.
BIO
Farid Karbalaei He received B.Sc. degree in power engineering from K. N. Toosi University of Technology, Tehran, Iran, in 1997, and M.Sc. and Ph.D. degrees in power engineering from Iran University of Science and Technology, Tehran in 2000 and 2009, respectively. Currently he is an assistant professor in Shahid Rajaee Teacher Training University. His research interests are power system dynamics and control, voltage stability and collapse and reactive power control.
Shahriar Abasi He received B.Sc. and M.Sc. degree in power engineering from Shahid Rajaee Teacher Training University, Tehran, Iran, in 2007 and 2011, respectively. Currently he is a Ph.D. candidate in power engineering in Razi University, Kermanshah, Iran. His research interests are voltage stability Studies, reactive power planning and transmission network expansion planning.
Verbic G.
,
Gubina F.
2003
“Fast voltagecollapse lineprotection algorithm based on local phasors”
IEE Proc. Gener.Transm.Distrib.
150
(4)
482 
486
DOI : 10.1049/ipgtd:20030461
Parniani M.
,
Vanouni M.
2010
“A fast local index for online estimation of closeness to loadability limit”
IEEE Trans. Power Systems
25
(1)
584 
585
DOI : 10.1109/TPWRS.2009.2036460
Morison G. K.
,
Gao B.
,
Kundur P.
1993
“Voltage stability analysis using static and dynamic approaches”
IEEE Trans. Power Systems
8
(3)
1159 
1171
DOI : 10.1109/59.260881
Gao B.
,
Morison G. K.
,
Kundur P.
1992
“Voltage stability evaluation using modal analysis”
IEEE Trans. Power Systems
7
(4)
1529 
1541
DOI : 10.1109/59.207377
Berizzi A.
,
Finazzi P.
,
Dosi D.
,
Marannino P.
,
Corsi S.
1998
“First and second order methods for voltage collapse assessment and security enhancement”
IEEE Trans. Power Systems
13
(2)
543 
551
DOI : 10.1109/59.667380
Lof P. A.
,
Smed T.
,
Andersson G.
,
Hill D. J.
1992
“Fast calculation of a voltage stability index”
IEEE Trans. Power Systems
7
(1)
54 
64
Kessel P.
,
Glavitsch H.
1986
“Estimating the voltage stability of a power system”
IEEE Trans. Power Delivery
PWRD1
(3)
346 
354
Lee CY.
,
Tsai SH.
,
Wu YK.
2010
“A new approach to the assessment of steadystate voltage stability margins using the PQV curve”
Int. J. Electrical Power and Energy Systems
32
1091 
1098
DOI : 10.1016/j.ijepes.2010.06.005
Xu W.
,
Pordanjani I. R.
,
Wang Y.
,
Vaahedi E.
2012
“A network decoupling transform for phasor data based voltage stability analysis and monitoring”
IEEE Trans. Smart Grid
3
(1)
261 
270
DOI : 10.1109/TSG.2011.2163175
Dobson I.
,
Lu L.
1993
“New method for computing a closest saddle node bifurcation and worst case load power margin for voltage collapse”
IEEE Trans. Power Systems
8
(3)
905 
913
Alvarado F.
,
Dobson I.
,
Hu Y.
1994
“Computation of closest bifurcations in power systems”
IEEE Trans. Power Systems
9
(2)
918 
928
DOI : 10.1109/59.317655
Flatabo N.
,
Fosso O. B.
,
Ognedal R.
,
Carlson T.
,
Heggland K.R.
1993
“A method for calculation of margins to voltage instability applied on the Norwegian system for maintaining required security level”
IEEE Trans. Power Systems
8
(3)
920 
928
DOI : 10.1109/59.260910
Lemaltre C.
,
Paul J. P.
,
Tesseron J. M.
,
Harmand Y.
,
Zhao Y. S.
1990
“An indicator of the risk of voltage profile instability for realtime control applications”
IEEE Trans. Power Systems
5
(1)
154 
161
DOI : 10.1109/59.49100
de Souza A.C.Z.
,
Canizares C.A.
,
Quintana V.H.
1997
“New techniques to speed up voltage collapse computations using tangent vectors”
IEEE Trans. Power Systems
12
(3)
1380 
1387
DOI : 10.1109/59.630485
Zarate L. A. Li.
,
Castro C. A.
2006
“Fast computation of security margins to voltage collapse based on sensitivity analysis”
IEE Proc.Gener. Transm. Distrib.
153
(1)
35 
43
DOI : 10.1049/ipgtd:20045156
Pama A.
,
Radman G.
2009
“A new approach for estimating voltage collapse point based on quadratic approximation of PVcurves”
Int. J. Electric Power Systems Research
79
653 
659
DOI : 10.1016/j.epsr.2008.09.018
Cutsem T. V.
,
Vournas C.
1998
“Voltage stability of electric power systems”
Kluwer Academic
Ajjarapu V.
,
Christy C.
1992
“The continuation power flow: a tool for stead state voltage stability analysis”
IEEE Trans. Power Systems
7
(1)
416 
423