Advanced
Quick and Accurate Computation of Voltage Stability Margin
Quick and Accurate Computation of Voltage Stability Margin
Journal of Electrical Engineering and Technology. 2016. Jan, 11(1): 1-8
Copyright © 2016, The Korean Institute of Electrical Engineers
This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.
  • Received : March 01, 2015
  • Accepted : August 06, 2015
  • Published : January 01, 2016
Download
PDF
e-PUB
PubReader
PPT
Export by style
Article
Author
Metrics
Cited by
TagCloud
About the Authors
Farid Karbalaei
Corresponding Author: Faculty of Electrical Engineering, Shahid Rajaee Teacher Training University (SRTTU), Tehran, Iran. (f_karbalaei@Srttu.edu)
Shahriar Abasi
Electrical Engineering Department, Faculty of Engineering, Razi University, Kermanshah, Iran. (shahriarabasi@gmail.com)

Abstract
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 maximum-loading 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 maximum-loading point. This method is tested on four IEEE test systems.
Keywords
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 measurement-based methods [1 , 2] and model-based ones [3 - 16] . The formers are based on the fact that maximum-loading point is detectable using the local measurement of voltage and current phasors. Since the system model is not used, the measurement-based methods are suitable for the analysis of the current power system only, and are unusable for the examination of contingency states.
The model-based 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 maximum-loading point or present indices that can provide a measure of the MW-distance from the current operating point to the maximum-loading 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 L-index [7] is one other well-known 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 single-one, single-branch 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 maximum-loading 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 maximum-loading 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 second-order 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 second-order 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 118-bus 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 maximum-loading 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:
PPT Slide
Lager Image
Typical PV curve
PPT Slide
Lager Image
The coefficients ai , bi and ci are determined by only one power flow solution in each loading factor. To do this, the first-order and second-order derivatives dVi / and d2Vi / 2 are calculated from the power flow equations and substituted in Eqs. (2) and (3):
PPT Slide
Lager Image
PPT Slide
Lager Image
The coefficients ai , bi and ci can be found by solving Eqs. (1), (2) and (3) as:
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
After the determination of the coefficients ai , bi and ci , the estimated collapse point voltage V c,i and the estimated critical loading factor λ c,i are obtained as follows:
PPT Slide
Lager Image
PPT Slide
Lager Image
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 ai , bi and ci 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.
PPT Slide
Lager Image
Collapse point voltages at different contingency states in 14-bus test system
PPT Slide
Lager Image
Collapse point voltages at different contingency states in 30-bus 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 Nc + 1. Nc is the number of contingency states. The proposed quadratic function is as:
PPT Slide
Lager Image
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:
PPT Slide
Lager Image
PPT Slide
Lager Image
where Pi0 ( Qi0 ) is the net-injected 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 second-order derivative d2Vi / 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 118-bus 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.
PPT Slide
Lager Image
Collapse point voltages at different contingency states in 57-bus test system
PPT Slide
Lager Image
Collapse point voltages at different contingency states in 118-bus 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 118-bus 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 30-bus test systems, from 0.47 to 0.52 pu in the 57-bus test system, and from 0.77 to 0.82 pu in the 118-bus 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 ai and bi , the first-order derivative dVi / is obtained using (9) as follow:
PPT Slide
Lager Image
In each iteration, with a given Vi , λ and dVi / is calculated by solving the power flow equations. Having Vi , λ and dVi / , the coefficients ai and bi can be found from Eqs. (9) and (11) as:
PPT Slide
Lager Image
PPT Slide
Lager Image
After determining the coefficients ai and bi , the estimated collapse point voltage V c,i and the estimated critical loading factor λc are obtained as follows:
PPT Slide
Lager Image
PPT Slide
Lager Image
To explain how to calculate dVi / , it is assumed that the power flow equations are as the compact form:
PPT Slide
Lager Image
where X is the system-state vector (voltage magnitudes and phase angles) and λ is the loading factor. A Taylor series expansion of Eq. (16) provides the relation:
PPT Slide
Lager Image
where gX 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:
PPT Slide
Lager Image
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 Nc + 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 dVi / using (18).
6) Calculate ai and bi 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 maximum-loading 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 dVi / using (18).
11) Calculate ai and bi 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 .
PPT Slide
Lager Image
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 118-bus 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 118-bus 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.
PPT Slide
Lager Image
The number of required iterations at different contingency states in 14-bus test system
PPT Slide
Lager Image
The number of required iterations at different contingency states in 30-bus test system
PPT Slide
Lager Image
The number of required iterations at different contingency states in 57-bus test system
PPT Slide
Lager Image
The number of required iterations at different contingency states in 118-bus 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 118-bus 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
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
The initial and critical loading factors for different contingency states in 14-bus test system
PPT Slide
Lager Image
The initial and critical loading factors for different contingency states in 30-bus test system
PPT Slide
Lager Image
The initial and critical loading factors for different contingency states in 57-bus test system
PPT Slide
Lager Image
The initial and critical loading factors for different contingency states in 118-bus 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
PPT Slide
Lager Image
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
PPT Slide
Lager Image
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 118-bus test systems. For instance, in 118-bus 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.
References
Verbic G. , Gubina F. 2003 “Fast voltage-collapse line-protection algorithm based on local phasors” IEE Proc. Gener.Transm.Distrib. 150 (4) 482 - 486    DOI : 10.1049/ip-gtd: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 PWRD-1 (3) 346 - 354
Lee CY. , Tsai SH. , Wu YK. 2010 “A new approach to the assessment of steady-state voltage stability margins using the P-Q-V 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 real-time 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/ip-gtd:20045156
Pama A. , Radman G. 2009 “A new approach for estimating voltage collapse point based on quadratic approximation of PV-curves” 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