Engineering, Technology & Applied Science Research Vol. 8, No. 2, 2018, 2758-2763 2758 www.etasr.com Ngo et al.: Assessing Power System Stability Following Load Changes and Considering Uncertainty Assessing Power System Stability Following Load Changes and Considering Uncertainty Van Duong Ngo University of Danang Danang, Vietnam nvduong@ac.udn.vn Van Kien Pham Faculty of Electrical Engineering The University of Danang—University of Science and Technology, Vietnam pvkien@dut.udn.vn Dinh Duong Le Faculty of Electrical Engineering The University of Danang—University of Science and Technology, Vietnam ldduong@dut.udn.vn Kim Hung Le Faculty of Electrical Engineering The University of Danang—University of Science and Technology, Vietnam lekimhung@dut.udn.vn Van Ky Huynh University of Danang Danang, Vietnam hvky@ac.udn.vn Abstract—An increase in load capacity during the operation of a power system usually causes voltage drop and leads to system instability, so it is necessary to monitor the effect of load changes. This article presents a method of assessing the power system stability according to the load node capacity considering uncertainty factors in the system. The proposed approach can be applied to large-scale power systems for voltage stability assessment in real-time. Keywords-stability; power system; power plane; uncertainty; Gaussian elimination I. INTRODUCTION Power system stability could be defined as “the ability of an electric power system, for a given initially operating condition, to regain a state of operating equilibrium after being subjected to a physical disturbance with most system variables bounded, so that practically the entire system remains intact” [1]. In [1], power system stability is broadly classified into three groups: frequency stability [2-4], rotor angle stability [5-7], and voltage stability [8-14]. Voltage stability could be defined as the ability of the system to maintain steady voltages at all buses in the system after being subjected to a disturbance from a given initial operating condition [1, 8, 11]. Voltage stability is divided into two subcategories: large and small disturbance stability. System faults, tripping transmission lines or dropping in generators are considered large disturbances whereas load changes are considered small disturbances. Static voltage stability is considered in the present paper in particular. The system is assumed to be operated in an equilibrium state and static voltage stability analysis assesses the feasibility of the operating point to provide the system operators with a permissible region in which the system can operate normally. Many techniques have been used in static voltage stability analysis in the literature. The well-known P-V and Q-V curves are widely used [8, 15-20] to determine the maximum permissible loading of the system. In particular, a P-V curve provides the relationship between the real power load and bus voltage and it is depicted with a constant power factor [17]. Contrariwise, a Q-V curve, plotted for a constant power [17], gives the change of bus voltages with respect to reactive power injection or absorption. Procedures for constructing both P-V and Q-V curves are time-consuming because a large number of power flows is needed to be executed using conventional methods and models [21]. Due to this drawback, they are not suitable to be used for online analysis. Moreover, they could be used only for certain increasing modes and not for providing the whole view of the analysis. Generally, they examine an individual bus by stressing the considered bus independently so they are not able to fully reflect the real stability condition of the system. In fact, voltage collapse occurs when the system load, i.e., real power and/or reactive power load, grows over a certain limit. Hence, voltage stability boundary needs to be plotted on a power plane [17]. It will form a P-Q curve [8, 22-26] that is very useful to determine boundary and operating regions for the system. Once the boundary is determined, the distance from the operating point to the voltage collape of the system can be directly assessed. For this work, stability reserve ratios are widely used in practice. Figure 1 shows, as an example, a P-Q curve (stability boundary) that divides power plane into two regions: normal region and impossible operating region. From Figure 1, stability reserve ratios Pk and Qk can be calculated as follows: lim 0 P 0 lim 0 Q 0 P - P k = 100% P Q - Q k = 100% Q where, P0 and Q0 are real and reactive power of the considered load at the operating point; Plim and Qlim are the limitations of real and reactive power, recpectively. A P-Q curve can be determined by various approaches. In [8, 22-24], the curve is characterized by a parabolic equation, while it is assumed to follow a circle in [25]. Nevertheless, such assumptions on the shape of the curve make it unrealistic. In addition, those techniques usually require high computational time so they are not suitable for online voltage stability ass by cur acc loa out sta nu env sim ana cur pro po Se sys sta cur loa gro nu rem adm (1) wh is t Engineerin www.etasr sessment. In [2 -point without rve. However, counted for so ads, uncertain tages of branc Fig. 1. In this pape ability conside mber of attrac vironment, w mplify a comp alytical techn rves. Thanks ovide results q wer systems condly, in or stem, a large ability is asse rves built are ad is used to co II. SIMPL Considering ounding bus, mbered from maining buse mittance equa ). . . . 11 121 . . . 21 221 . . . F1 F21 Y V Y Y V Y ... .. Y V Y ... ..       . . . i1 i 21 . . . N1 N 21 Y V Y ... .. Y V Y                   here: . ii . 1 ij 1N i i j Y Z   the impedance ng, Technology r.com 26], P-Q curv t using any as , all of the abo ources of unce nties from re ches and gener P-Q curve and er, we propose ering uncertai ctive features. we adopt Gau plex diagram t nique develope to these techn quickly, hence for voltage s rder to take i number of P- essed probabi realistic, since onstruct the cu LIFIED POWER S a power netw denoted 0) 1 to F, loads a es are interm ations to presen . . . 2 1F2 F . . . 2 112 F . . . 2 FF2 F V ... Y V V ... Y V . ... ... V ... Y V . ... ..    . . . 2 iF2 F . . . 2 NF2 F . V ... Y V . ... ... V ... Y V   :Y = ∑ e of branch bet y & Applied Sci Ngo e e can be fastly ssumption abo ove mentioned ertainty in pow enewable sou rating units. operating region e an approach inties in pow Firstly, for wo ussian elimina to a simple on ed in [26] fo niques, the pro e it can be ap stability assess into account -Q curves is p ilistically. In e no assumptio urves. SYSTEM DIAGR work with N+1 in which ge are represented mediate buse nt steady-state . . 1i i . . 2i i . . Fi i ... Y V ... ... Y V ... ... ... ... ... Y V ...       . . ii i . . Ni i ... ... ... ... Y V ... ... ... ... ... Y V ...    self admittan tween i and j; ience Research et al.: Assessing y constructed out the shape d approaches a wer systems, su urces and ra s on power plane h to assess v wer systems w orking in a rea ation techniq ne, then we u or constructing oposed metho pplied to large sment in real randomness i plotted and v addition, the on on the chan RAM ALGORITH 1 buses (includ neration buse d by admittanc es, we can e of the system . . . 1N N 1 . . . 2 N N 2 . . . FN N F Y V J Y V J ... ... ... Y V J       . . iN N . . NN N ... ... ... Y V 0 ... ... ... Y V 0     nce ( . ij ijZ = R + i = 1÷N; j = 1÷ h g Power System point- of the are not uch as andom . oltage with a al-time que to use the g P-Q od can e-scale l-time. in the oltage e P-Q nge of HM ding a es are ces so form m as in (1) ij+ jX ÷N), j. mak foll . N . NN Y Y and wil i.e., adm last equ mak gen 2. gro ord load Vol. 8, No. 2, 2 Stability Follo . ij . ij 1 Z Y : Y From (1), ad ke a simplifi lows. First, devide b . .N1 1. NN Y Y - V ...- Y Y Next, multiply . . .N1 NF iN 1 . NN NN Y Y V ...- Y Eventually, su d the Nth comp l be eliminated After the abov , the last equa mittances for t bus. We rep uations from th ke an equivale neration buses Fig. 2. Equiv considered. In Figure 2, ounding bus 0 der to assess th d, electromoti 2018, 2758-276 owing Load Cha = mutual a dopting Gaussi ied equivalen both sides of t . . .NF Ni F. . NN NN Y Y V ... Y Y y both sides o . . .F Ni iN F . N NN Y Y V ... Y ubtract (3) from *   . . iFiFY Y ponent in each d. ve steps, we c ation in (1). E the equivalen peat the abov he initial N eq ent network fo and a load bu . 1E . 2E . iE . FE . 1( F 1)Z  . 10Y . 2(F 1)Z  . 20Y . i ( F 1)Z  . i0Y . F( F 1)Z  . F0Y valent network w admittance Y included admi he ability accor ive force Εi (i 63 anges and Cons admittance bet ian eliminatio nt diagram fo the last equatio . . i N N V ...+V = 0 f each equatio . . . iN iNiY V ...+Y V m the ith equat . . NF iN . NN Y Y Y h equation in can eliminate o Equation (4) i nt network aft ve process un quations in (1). or the consider us considered a . ( F 1)V  . ( F 1Y  with F generation Y(F+1)0 between ittance of the rding to chang i=1÷F) is con 2759 sidering Uncer tween bus i an on method, we or the networ on in (1) by . Y 0 (2) on in (2) with Y . NV = 0 (3) tion in (2) yiel (4) (2) (i.e., iNY . one equation in s used to calc fter eliminating ntill obtaining . Therefore, w red network w as shown in F 1)0 n buses and a loa n bus F+1 and load at bus F+ ge at the consid nsiderd as con rtainty nd bus e can rk as NN . iNY lds NV . . ) n (1), culate g the g F+1 we can with F Figure ad bus d the +1. In dered nstant Engineering, Technology & Applied Science Research Vol. 8, No. 2, 2018, 2758-2763 2760 www.etasr.com Ngo et al.: Assessing Power System Stability Following Load Changes and Considering Uncertainty (so Yi0 can be eliminated) and the load is separated as in Figure 3. . 1E . ( F 1)V  . 2E . iE . FE SL = PL+jQL . 1( F 1)Z  . 2( F 1)Z  . i ( F 1)Z  . F( F 1)Z  . 0Y Fig. 3. Simplified equivalent diagram. In Figure 3, Y = G + jB is parallel with the admittance of the load Y = G + jB , where 0 ( F 1)0 LG G G  0 ( F 1)0 LB B B  2 L L F 1,ratedG P / V  2 L L F 1,ratedB Q / V  and F 1,ratedV  is the rated voltage at bus F+1. Equation (1) can be represented in a matrix form as in (5), in which the first matrix is denoted as Y (called admittance matrix): . . . . . 11 12 1F 1i 1N . . . . . 21 22 2F 2i 2N . . . . . F1 F2 FF Fi FN . . . . . i1 i2 iF ii iN +Y - Y ... - Y ... -Y ... - Y - Y +Y ... +Y ... - Y ... -Y ... ... ... ... ... ... ... ... - Y - Y ... +Y ... - Y ... - Y ... ... ... ... ... ... ... ... - Y - Y ... -Y ... +Y ... - Y ... ... ... ... ... . . . 1 1 . . 2 2 . . F F . . i i . . . . . . . N1 N2 NF Ni NN N N V J V J ... ... V J= ... ... V J .. ... ... ... ... - Y - Y ... -Y ... - Y ... +Y V J (5) It is worth noting that when we need to assess voltage stability at a certain load bus, that bus is numbered F+1 before forming either (1) or (5). Alternatively, at first, the equations are formed, then in order to consider any load bus i, rows i and F+1, and columns i and F+1 are exchanged. Based on the above analysis, we developped an algorithm for making a simplified equivalent diagram for the considered network, as in Figure 4. III. REPRESENTING UNCERTAIN FACTORS IN ELECTRICAL POWER SYSTEMS In a power system, random factors related to the load, failures of elements in the system, renewable energy sources can be represented by probability distribution functions [27- 29]. These functions describe the inherent nature of uncertainty factors and need to be integrated into the analysis of the system. In practice, such functions can be estimated based on historical data measured at loads, data on failures of lines, generating units, etc. In the literature load is usually expressed by a normal distribution function, while renewable power generation is usually represented by a generic function such as Beta, Gamma, Weibull functions, etc [29]. Random outage of an element such as transmission line, transformer, generating unit can be described by a 0-1 distribution function (0: outage state, 1: working state; an operating element can be failed with a certain probability) [29]. If all generating units in a power plant are the same, the combination of 0-1 distribution function for each unit could be represented by a binomial distribution function [29]. Compute and form admittance matrix Y of the system Start Power system components libraries Calculate admittances of all branches of the system Check information of power system (System parameters and operating parameters) S a ve m a tr ix Y Select the load bus i to calculate Calculate admittances for simplifying the network diagram using (4) and k = N-1 Exchange row i and row (F+1), column i and column (F+1) in matrix Y k < (N - F-1) Set value k = N End True False True False Calculate admittances of the simplified equivalent diagram with (F+1) buses in Fig. 3 Save result Input data of the system (System parameters and operating parameters) Fig. 4. Algorithm for making a simplified equivalent diagram based on Gaussian elimination method. IV. STABILITY ASSESSING ALGORITHM DEVELOPMENT For an electrical power system with N buses (excluding a grounding bus, with F generation buses, bus 1 is considered as the slack bus), in order to examine the stability of the system according to the change at a certain load bus, first, we use the algorithm discussed in Section II to obtain the simplified equivalent diagram for the considered network. When the system is operating at a certain state with known network configuration, generating power, load, ect., we increase load at the considered bus and use the pragmatic criterion dQ/dV presented in [26, 30] to construct the permissible operating region on a power plane under the conditions of a static stability limit as in Figure 5. In Figure 5, M1 is a certain operating point of the load considered belonging to stable region, while M2 belongs to unstable region of the system (see Figure 1). Suppose that the system is operating at M1, its stability and the dangerous level of increasing load can be evaluated according to the distance from M1 to the stability boundary. When random factors, as discussed in Section III, are considered, stability boundary is not represented by a single curve but by a set of curves as shown in Figure 6. Uncertainties in the system, related to load, renewable sources, random outages of lines, generating units, etc., can be represented by probability functions and set of random samples representing the sta sam are (se sys lin and sys loa and reg po wo If wit gre car avo ass con (Fi nor (de be bra Engineerin www.etasr ese functions c atistics [31, 32 mple. Assume e plotted) and ee Figure 6). stem in terms ne connecting d the probab stem can be ca Therefore, b ad on power p d the system gions 1, 2, 3 o int M1, the orking at M2, it is working th probability eatter than a ce rry out approp oid instability sessing stabilit nsidering unce F Fig. 6. Set uncertainty in th We consider igure 8). In th rmal distribu eterministic va equal to 10% anch is assume ng, Technology r.com can be genera 2]. In Figure 6 e that N1 samp d the considere In order to of static volt the origin O a bility of unsta alculated as in 2 1 N 1 N p = ased on the o lane, the stabi falls into one on power plan system is cer the system is at M3, the sy y of instabilit ertain value, th priate remedi y for the syste ty of power sy ertainty is show 0 P [MW] Q1 P1 M1 Q P2 M Fig. 5. Operat 0 P [MW] Q1 P1 M1 Q P2 M P3 of P-Q curves he system. V. MODEL r the IEEE 39 he system, loa ution charact alue [33]) and of its expecte ed to have pro y & Applied Sci Ngo e ated using sam 6, each curve ples are genera ed load is ope evaluate dang tage stability, and M3 that i ability (dange (6):  100 % operating poin lity of the syst of three case ne): (1) if the rtainly stable certainly unst ystem may be ty p (0%