A second order power flow based on current injection equations

Electrical Power and Energy Systems 27 (2005) 254–263 www.elsevier.com/locate/ijepes

A second order power flow based on current injection equations Carlos Aparecido Ferreiraa, Vander Menengoy da Costab,* b

a EletrObraS, Rio de Janeiro 22210-903, Brazil Universidade Federal de Juiz de Fora, Faculdade de Engenharia, Juiz de Fora, Minas Gerais 36036-330, Brazil

Received 19 March 2003; revised 5 August 2004; accepted 29 October 2004

Abstract This paper presents a new second order power flow methodology by using current injection equations expressed in voltage rectangular coordinates instead of using traditional rectangular power flow equations. The results presented show that the proposed method leads to a substantially faster second order power flow solution, when compared to the conventional method expressed in terms of power mismatches and written in rectangular coordinates. In order to evaluate the performance of the proposed methodology, different heavily loaded and overloaded systems, as well as ill-conditioned distribution systems have been tested. These results show the feasibility of employing this new methodology as an additional tool for power flow studies. q 2004 Elsevier Ltd. All rights reserved. Keywords: Power flow analysis; Ill-conditioned systems; Current injection equations

1. Introduction The power flow problem deals with the calculation of voltages and line flows, in a large sparse electrical network, for a given load and generation schedule. The conventional power flow solution comprises power equations expressed in terms of polar or rectangular coordinates. Many important contributions have been reported on this area along the last thirty years [1–9]. Current-versions power flow in voltage rectangular coordinates are proposed in Refs. [10,11]. In [10] the key idea is to solve an augmented system of equations in which both bus voltages and current injections appear as state variables, and both power and current mismatches are zeroed. In [11], a new dependent variable DQ is introduced for each PV bus together with an additional equation imposing the voltage constraint at these buses. Besides, except for PV buses, the Jacobian matrix has the elements of the (2!2) off-diagonal blocks equal to those of the bus

* Corresponding author. Address: Rua Domingos dos Anjos, 5, Centro, Treˆs Rios, 25815-000 Estado do Rio de Janeiro, Brazil. E-mail addresses: [email protected] (C.A. Ferreira), [email protected] (V.M. da Costa). 0142-0615/$ - see front matter q 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijepes.2004.10.005

admittance matrix expanded into real and imaginary coordinates. Similarly to [10], the elements of the (2!2) diagonal blocks are updated at each iteration according to the bus load model adopted. Basically, the main difference between these methods is that in [11] the state vector is exclusively composed of bus voltages in rectangular coordinates. Based on [11], a sparse augmented formulation for solving a set of control devices in the current injection power flow problem is described in Refs. [12,13]. This formulation has the same convergence characteristics of the conventional Newton power flow problem, expressed in terms of power mismatches written in polar coordinates. Many different flexible AC transmission system devices and controls are incorporated into the power flow problem by using this augmented methodology. Although many advances have been made in order to improve both efficiency and robustness of the power flow solution methods, some problems can still be observed in many practical situations. Firstly, the system has an operating point, but the power flow fails to provide it. This situation occurs when: (a) the system is ill-conditioned and the power flow methods cannot find the solution due to numerical difficulties. (b) The system becomes

C.A. Ferreira, V.M. da Costa / Electrical Power and Energy Systems 27 (2005) 254–263

255

Nomenclature h iteration counter n number of buses Ek voltage phasor at bus k DPkCjDQk complex power mismatch at bus k PGk C j QGk generated complex power at bus k PkCjQk net complex injected power at bus k DIrk C jDImk complex current mismatch at bus k Irk C jImk injected current at bus k Vrk C j Vmk complex voltage at bus k Vk voltage magnitude at bus k qk voltage angle at bus k Vksp scheduled voltage magnitude at bus k ill-conditioned due to its constant load increase without the corresponding increase in generation and transmission facilities. Secondly, there is no operating point. This situation may arise after severe disturbances. Finally, it is also possible that the system admits multiple solutions. In order to overcome some difficulties for convergence and to provide additional information on the iterative process, in Ref. [14] the conventional Newton Raphson power flow calculation, in voltage rectangular coordinates, is treated as a non-linear programming problem, which determines the direction and the magnitude of the solution so that a certain cost function is minimized. The cost function value becomes zero if there is a solution from the initial estimate, and stays at a positive value if no solution exists. In this method the fully expansion in Taylor’s series up to the third term of the load flow equations in rectangular coordinates is explored. The methodology proposed in Ref. [14] is named in this paper as Iwamoto’s method. This paper proposes a new second order power flow, following the main ideas presented in Ref. [14], but considering the current injection formulation [11] expressed only in terms of voltage rectangular coordinates. This new methodology not only presents a quite satisfactory convergence characteristics, but also reduces the computational effort for the second order power flow solution. In this paper, firstly the current injection method [11–13] and the conventional second order power flow [14] are briefly described, in order to make the paper self-contained. Secondly, the proposed methodology, results and conclusions are presented.

2. Current injection formulation 2.1. PQ buses representation This formulation uses a set of 2n current injection equations written in rectangular coordinates [11]. The real and imaginary parts of these equations, denoted by Irk and

calculated voltage magnitude at bus k Vkcalc GkjCjBkj (k,j)th element of bus admittance matrix DVk voltage magnitude correction at bus k DVrk ; DVmk real and imaginary voltage component correction at a bus k fk set of buses directly connected to bus k, including the bus k Uk set of buses directly connected to bus k J Jacobian matrix H Hessian matrix Matrices are shown in bold and vectors are shown in bold underlined Imk respectively, are given by 0Z

Vrk Pk C Vmk Qk X K ðGki Vri K Bki Vmi Þ Vr2k C Vm2 k i2f

(1)

Vmk Pk K Vrk Qk X K ðGki Vmi C Bki Vri Þ Vr2k C Vm2 k i2f

(2)

k

0Z

k

The Newton Raphson solution algorithm applied to Eqs. (1) and (2), considering all system buses as being of the PQ type, is given by 3 2 0 0 3 2 B11 G11 B12 G12 / B1n G1n 2 DVr 3 DIm1 1 7 6 6 DI 7 6 G00 B00 G12 KB12 / G1n KB1n 76 DV 7 11 76 m1 7 6 r1 7 6 11 76 7 6 7 6 0 0 76 7 6 DI 7 6 B 6 m2 7 6 21 G21 B22 G22 / B2n G2n 76 DVr2 7 76 7 6 7 6 6 DIr 7 Z 6 G KB G00 B00 / G KB 76 DVm 7 21 2n 2n 76 27 22 22 6 2 7 6 21 76 7 6 7 6 6 « 7 6 « 6 « 7 « « « 1 « « 7 76 7 6 7 6 76 7 6 7 6 0 0 DV 7 5 4 DImn 5 6 Bn1 Gn1 Bn2 Gn2 / Bnn 4 r G n nn 5 4 DIrn DVmn G KB G KB / G00 B00 n1

n1

n2

n2

nn

nn

(3) where [11,12] 0 Bkk Z

vImk ZBkk Kak vVrk

(4)

0 Z Gkk

vImk ZGkk Kbk vVmk

(5)

00 Z Gkk

vIrk ZGkk Kck vVrk

(6)

00 Z Bkk

vIrk ZKBkk Kdk vVmk

(7)

Vmk DPk KVrk DQk Vr2k CVm2 k

(8)

Vrk DPk CVmk DQk Vr2k CVm2 k

(9)

DImk Z DIrk Z

256

C.A. Ferreira, V.M. da Costa / Electrical Power and Energy Systems 27 (2005) 254–263

hand, the diagonal (2!2) blocks must be updated at each iteration according to Eqs. (4)–(7).

The parameters ak, bk, ck, and dk, are presented in Appendix A of [12].

2.3. Bus voltage updating

2.2. PV buses representation

After determining the bus voltage corrections by using (13), the updated voltage values in rectangular coordinates, at a generic iteration (hC1), are given by

The computation of real and imaginary current mismatches is straightforward for PQ buses, because of real and reactive power mismatches are known. For PV buses, the reactive power mismatches are unknown and treated as a state variable. With the insertion of this new variable, it must be introduced one additional equation for each PV bus. The equation adopted for imposing the voltage constraint at a PV bus k in rectangular coordinates is given by Vk2 Z Vr2k C Vm2 k

(10)

vVk2 vV 2 DVrk C k DVmk vVrk vVmk

Z 2Vrk DVrk C 2Vmk DVmk

(14)

VmðhC1Þ Z Vmh C DVmh

(15)

At this point, it is very important to emphasize the voltage updating process presented in this paper. In Refs. [11–13], the current equations are written in voltage rectangular coordinates, but the voltage corrections are performed by using polar coordinates, through a rectangular-polar transformation. On the other hand, in this paper there is no need to proceed any type of transformation, since the updating process uses the rectangular coordinates.

Eq. (10) becomes DVk2 Z

VrðhC1Þ Z Vrh C DVrh

(11)

where DVk2 Z ðVksp Þ2 K ðVkcalc Þ2

3. Description of Iwamoto’s method

(12)

For this situation, the current injection Newton-Raphson algorithm is shown in Eq. (13). Note that bus k is of the PV type. 3 2 Vm 1 Vr1 DP K DQ 1 1 7 2 0 6 V12 V12 0 7 6 B11 G11 B12 G12 / 7 6 Vr 1 Vm1 7 6 6 DP C DQ 1 17 00 6 V2 6 00 V12 7 6 G11 B11 G12 KB12 / 6 1 7 6 Vm 2 6 V r2 0 0 7 6 6 6 V 2 DP2 K V 2 DQ2 7 6 B21 G21 B22 G22 / 7 6 2 6 2 00 00 7 6 6 Vr 6 2 DP C Vm2 DQ 7 6 G21 KB21 G22 B22 / 2 27 6 V2 6 2 V2 7 6 « 6 2 « « « « 7 6 6 « 7 6 6 7 6 6 7 Z 6 Bk1 Gk1 Bk2 Gk2 / 6 Vm k 7 6 6 DP k 7 6 6 2 Vk 7 6 6 7 6 Gk1 KBk1 Gk2 KBk2 / 6 Vrk 7 6 6 DPk 7 6 6 2 Vk 7 6 « 6 « « « « 7 6 6 « 7 6 6 7 6 Bn1 Gn1 Bn2 Gn2 / 6 7 6 6 Vm 6 n DP K Vrn DQ 7 6 n n7 6 V2 6 2 Vn 7 4 Gn1 KBn1 Gn2 KBn2 / 6 n 7 6 Vr V m 6 n DP C n DQ 7 0 0 0 0 / n n5 4 V2 Vn2 n DVk2

This section contains a brief description of the method proposed in Ref. [14]. The conventional power flow equations in rectangular coordinates can be written in the

B1k G1k B2k G2k

G1k

/ B1n

G1n

KB1k / G1n KB1n G2k

/ B2n

G2n

KB2k / G2n KB2n

«

«

«

«

«

0 Bkk

0 Gkk

/ Bkn

Gkn

00 Gkk

00 Bkk

/ Gkn

KBkn

«

«

Bnk

Gnk

0 0 0 0 «

«

Vr k Vk2 KVmk Vk2 «

0 / Bnn

0 Gnn

0

00 Gnn

00 Bnn

0

0

0

0

«

Gnk

KBnk /

2Vrk

2Vmk

/

«

3

3 2 7 DVr 7 1 76 7 76 DVm1 7 76 7 76 7 76 DVr2 7 76 7 76 7 76 DVm2 7 76 7 76 « 7 76 7 76 7 76 DVr 7 k 7 76 76 7 76 DVmk 7 76 7 76 7 76 « 7 76 7 76 DV 7 76 rn 7 76 7 76 DV 7 74 mn 5 7 5 DQk

(13)

following form It can be easily observed that the off-diagonal (2!2) blocks are equal to the corresponding blocks of the bus admittance matrix, remaining constant during all the iterative process. This is one of the most important feature regarding the current injection Jacobian matrix. On the other

2

3 PðxÞ   5 S Z4 QðxÞ   

(16)

C.A. Ferreira, V.M. da Costa / Electrical Power and Energy Systems 27 (2005) 254–263

where g1 Z

t

x Z ½Vr1 Vm1 Vr2 Vm2 /Vrn Vmn  

2n X

g2 Z 3

(18)

where J is the conventional rectangular Jacobian matrix and SðDxÞ is the third term of the series which represents the   power equations calculated for the state vector increment at a given iteration. The solution of Eq. (18) can be obtained through a nonlinear programming method. Let us denote by m, the step size optimization factor. Thus, applying this factor to the correction vector in Eq. (18) results Ssp Z SðxÞ C JðmDxÞ C SðmDxÞ      

(29)

(19)

or

2n X

bi c i

(30)

c2i

(31)

iZ1

(17)

Eq. (16) is a quadratic function of the real and imaginary parts of the voltages, and its expansion in terms of Taylor’s series can still be given by Ssp Z SðxÞ C JDx C SðDxÞ      

ðb2i C 2ai ci Þ

iZ1

By expanding Eq. (16) in terms of Taylor’s series results 1 Ssp Z SðxÞ C JDx C HDx2   2   

257

g3 Z 2

2n X iZ1

For well-conditioned systems, m assumes values close to one. For systems with no feasible operating point, m tends to assume very low values, indicating that the actual voltage vector cannot be changed further in order to minimize function F. The state vector for the next iteration is given by xhC1 Z xh C mh Dxh   

(32)

4. Proposed second order power flow methodology 4.1. Basic principles: PQ buses

sp

2

S Z SðxÞ C mJðDxÞ C m SðDxÞ      

(20)

Eq. (20) can be re-written as a C mb C m2 c Z 0    

(21)

where a Z Ssp K SðxÞ    

(22)

b Z KJDx Z Ka   

(23)

c Z KSðDxÞ   

(24)

Imspk Z Imk C X C

vImk vImk DVrk C DVmk vVrk vVmk

ðBki DVri C Gki DVmi Þ C

i2Uk

The factor m is computed in order to minimize the following quadratic mismatch based function 2n 1X FZ ða C mbi C m2 ci Þ2 2 iZ1 i

vF Z0 vm

!

1 2

X  v2 Im v2 Imk k DVrk DVri C DVrk DVmi vVrk vVri vVrk vVmi i2F k

(25)

The optimal solution is given by

v2 Imk v2 I m k C DVmk DVri C DVmk DVmi vVmk vVri vVmk vVmi

(26)

By applying the factor m to the correction vector in Eq. (33), results in Eq. (21), where aImk Z DImk

2

3

g0 C g1 m C g2 m C g3 m Z 0

(33)

which leads to (27)

vImk vImk DVrk K DVmk vVrk vVmk X K ðBki DVri CGki DVmi Þ Z KDImk

(34)

bImk Z K

where g0 Z

The proposed methodology uses the power flow formulation based on current injection equations expressed in terms of voltage rectangular coordinates. To find the factor m is used the same idea presented in Ref. [14]. Let us consider the imaginary current injection Eq. (2) related to bus k. Its expansion in function of Taylor’s series, neglecting all terms of order higher than two, yields

2n X iZ1

ai bi

(28)

i2Uk

(35)

258

C.A. Ferreira, V.M. da Costa / Electrical Power and Energy Systems 27 (2005) 254–263

cImk Z K

 2 v Imk v2 I m k 1 X DVrk DVri C DVrk DVmi 2 i2F vVrk vVri vVrk vVmi k

v2 Imk v2 I m k C DVmk DVri C DVmk DVmi vVmk vVri vVmk vVmi



The reactive power mismatch can be written as X  Gki V 2 B V2 k DVri K ki k DVmi DQk Z Vmk Vmk i2U k

C

(36) Now, taking into account the real current injection Eq. (1) and proceeding in a similar way, the expansion in Taylor’s series and the application of factor m results aIrk Z DIrk

(37)

vIr vIrk bIrk Z K k DVrk K DVmk vVrk vVmk X K ðGki DVri K Bki DVmi Þ Z KDIrk 2

(38)

1 2

Substituting this value in another equation related to PV bus and effecting proper mathematical manipulations one has X DPk Z ðBki Vmk C Gki Vrk ÞDVri C ðGki Vmk i2Uk

0 00 C ðGkk Vmk C Bkk Vrk ÞDVmk

(42)

0 00 Bkk Vmk C Gkk Vrk Z Bkk Vmk C Gkk Vrk C Irk

(43)

0 00 Gkk Vmk C Bkk Vrk Z Gkk Vmk K Bkk Vrk C Imk

(44)

But

2

v Irk v Irk DVrk DVri C DVrk DVmi vVrk vVri vVrk vVmi i2Fk v2 Irk v2 Irk C DVmk DVri C DVmk DVmi ð39Þ vVmk vVri vVmk vVmi

cIrk Z K

(41)

0 00 K Bki Vrk ÞDVmi C ðBkk Vmk C Gkk Vrk ÞDVrk

i2Uk

X

00 2 Vr Gkk Vk B 00 V 2 DVrk C kk k DVmk K k DPk Vmk Vmk Vm k

Thus, in the proposed methodology, the components of vectors a, b and c related to a generic PQ bus k, are given by    ak Z ½DImk DIrk t t

bk Z K½DImk DIrk 

Finally, substituting Eqs. (43) and (44) in Eq. (42) one has DPk Z

vPk vPk DVrk C DVmk vVrk vVmk X  vPk vP C DVri C k DVmi vVri vVmi i2U

(45)

k

ck Z ½cImk cIrk t

The expansion of Eq. (45) in terms of Taylor’s series results in the following components of the vectors

4.2. Inclusion of PV buses

aPk Z DPk

(46)

For a generic PV bus k, it can be seen from Eq. (13) that the reactive power mismatch DQk is now a dependent variable. Besides, there are three variables, namely, DVr, DVm, DQ, associated with each PV bus. As previously stated, the main objective is to develop a new method for solving the second order power flow, using the current injection formulation only expressed in terms of voltage rectangular coordinates Vr and Vm. Thus, it becomes necessary to eliminate the variable DQ from the set of equations related to each PV bus, in order to make possible the comparison of results presented by this new method with those ones produced by Ref. [14], which solely uses voltage rectangular coordinates. From Eq. (13) for a PV bus one has

bPk Z KaPk

(47)

X Vr k 00 DPk Z ðGki DVri K Bki DVmi Þ C Gkk DVrk 2 Vk i2U k

00 C Bkk DVmk K

Vmk DQk Vk2

(40)

cPk Z K

X

DVrk ðDVri Gki C DVmi Bki Þ

i2fk

C DVmk ðDVmi Gki K DVri Bki Þ Z KPk ðDVr ; DVm Þ

(48)

In the same way, the expansion of Eq. (11) leads to aDVk Z ðVksp Þ2 K ðVkcalc Þ2 Z DVk2

(49)

bDVk Z KaDVk

(50)

cDVk Z KðDVr2k C DVm2 k Þ Z KVk2 ðDVr ; DVm Þ

(51)

Hence, in the proposed methodology the components of vectors a, b and c related to a generic PV bus k, are given by    2 t ak Z ½DPk DVk 

C.A. Ferreira, V.M. da Costa / Electrical Power and Energy Systems 27 (2005) 254–263 Table 1 Values of m—IEEE 300 buses/base case h

Iwamoto

Proposed

0 1 2 3

0.8815 1.1025 0.9980 1.0000

0.94247 1.0438 1.0039 1.0000

259

† Step 7: Solve Eq. (27) to determine the step size factor mh; † Step 8: Update the voltages 2 hC1 3 2 h 3 2 3 DVh V V 4  r 5 Z 4  r 5 C m h ,4  r 5 (52) VhC1 Vh DV h m m m † Step 9: If the solution process converged then stop, else go back to step 2. At this point, it is important to clarify that the insertion of the step optimization m does not alter the conventional procedures to consider limits and controls in the power flow problem.

5. Results

Fig. 1. Cost function F evolution—IEEE 300 buses/base case.

bk Z K½DPk DVk2 t Z Kak ck Z K½Pk ðDVr ; DVm Þ Vk2 ðDVr ; DVm Þt 4.3. Solution methodology

The conventional rectangular-based Newton Raphson, the Iwamoto’s method and the proposed second order methodology have been coded in Matlab using the same programming techniques. The results displayed have been obtained through tests with the IEEE-300 and the Brazilian Southeastern practical system. These systems have 300 buses and 401 circuits, and 730 buses and 1146 circuits, respectively. In addition, 11-bus and 43-bus radial distribution systems have also been tested. The convergence of iterative process is assumed when the maximum power mismatches at all buses are less than 1!10K5 p.u. The constant power model is adopted for the active and reactive loads. The threshold for m is 0.1; that is, if m is smaller than the threshold, the process is interrupted and it is assumed that the best possible solution has been reached.

The proposed algorithm can be summarized in the following steps:

5.1. IEEE-300 bus test system

† Step 1: Assemble the bus admittance matrix Y; † Step 2: Calculate the nodal current injections IZ Y,E;   † Step 3: Determine the power mismatches at all buses of the electrical network; † Step 4: If the maximum power mismatch is greater than a given tolerance, then determine the voltage corrections using Eq. (13) and increment the iteration counter, hZ hC1; † Step 5: Determine the vectors a, b and c considering the    different types of buses; † Step 6: Calculate the scalars gh0 , gh1 , gh2 and gh3 using Eqs. (28)–(31), respectively;

5.1.1. Base case For the base-case condition, the power flow solution is achieved in four iterations, by using the conventional Newton-Raphson method. The same number of iterations is required if the step size optimization factor m is used as shown in Table 1, which also presents a comparison for values of m obtained through Iwamoto’s method and the proposed one. Fig. 1 presents the cost function F for the three methods under analysis. It can be observed that the second order methods present a better performance than the conventional method. From Table 2, it can be noted

Table 2 Power flow results—IEEE 300 buses/base case Method

V9033 (p.u.)

q9033 (8)

V9031 (p.u.)

q9031 (8)

V9038 (p.u.)

q9038 (8)

PG7049 (p.u.)

QG7049 (p.u.)

Conventional Iwamoto Proposed

0.9288 0.9288 0.9288

K25.331 K25.331 K25.331

0.9317 0.9317 0.9317

K25.016 K25.016 K25.016

0.93916 0.93916 0.93916

K24.411 K24.411 K24.411

4.5595 4.5595 4.5595

0.3883 0.3883 0.3883

260

C.A. Ferreira, V.M. da Costa / Electrical Power and Energy Systems 27 (2005) 254–263

Table 3 Values of m—IEEE 300 buses/overloaded h

Iwamoto

Proposed

0 1 2 3 4 5 6

0.2231 0.2061 0.1942 0.1800 0.1594 0.1317 0.0993

0.2269 0.2045 0.1906 0.1761 0.1566 0.1304 0.0997

Fig. 2. Cost function F evolution—IEEE 300 buses/overloaded.

that the voltage at a certain buses and the power generated by the slack bus are the same, regardless of the methodology used. 5.1.2. Overloaded case The point of maximum loading occurs with an additional load increase of 3.64% at all load buses. In order to characterize the overloaded condition, a load increase of 5% has been taken into account. For this situation, the conventional Newton method fails to converge, whereas in the second order methods the best possible solution is achieved. Table 3 presents a comparison for values of m. Both methodologies find the best solution in seven iterations. Fig. 2 presents the cost function F for the three methods. It can be observed that the proposed second order methodology has the same convergence characteristics of the Iwamoto’s method. From Table 4, it can be noted that the best possible solution obtained by each method are slightly

different. Even though the operating point is not feasible, valuable information can be provided by this result in order to restore the power flow solution. 5.2. Brazilian southeastern test system 5.2.1. Heavily loaded case The point of maximum loading for this practical Brazilian system occurs with an additional load increase of 5.33% at all load buses. In order to characterize the heavily loaded condition, a load increase of 5.30% has been taken into account. For this situation, the conventional Newton method converges in six iterations. The step size optimization factor m can be utilized in order to improve the convergence of the iterative process. Table 5 shows the values of m obtained through both second order methodologies. The convergence was obtained with five iterations. From Table 6 it can be noted that the voltage at a certain buses and the power generated by the slack bus are the same, regardless of the methodology used. Fig. 3 presents the cost function F for the three methods. Again, it can be observed a very similar behavior related to the second order methods. 5.2.2. Overloaded case By increasing the loads at 6%, the conventional Newton method fails to converge. By using the step size optimization m the best possible solution is achieved. Table 7 presents the values of m obtained through both second order methods. From this table it is very interesting to note that the Iwamoto’s method converges in six iterations, whereas just five iterations are necessary using the proposed approach. From Table 8, it can be noted that the best possible solution obtained by each method are not equal. As for IEEE-300 operating at overloading condition, this slight variation stems from the fact that the values of m at the end of each iteration are different when calculated by both second order methodologies. For this operating condition, the convergence through power mismatches has not been achieved, being m (m!0.1) the parameter used to interrupt the iterative process. The evolution of function F is shown in Fig. 4. 5.2.3. Contingency case The modal analysis indicates that the transmission line connecting the buses 1135 and 1137 is critical with regard to voltage stability. Let us suppose a contingency involving the loss of this line. For this situation, by using the conventional Newton method the convergence is not achieved. By using both Iwamoto and proposed methods, the best possible

Table 4 Power flow results—IEEE 300 buses/overloaded Method

V9033 (p.u.)

q9033 (8)

V9031 (p.u.)

q9031 (8)

V9038 (p.u.)

q9038 (8)

PG7049 (p.u.)

QG7049 (p.u.)

Iwamoto Proposed

0.8085 0.8102

K69.054 K68.622

0.8124 0.8141

K68.602 K68.171

0.8219 0.8236

K67.718 K67.293

15.451 16.162

5.1911 6.081

C.A. Ferreira, V.M. da Costa / Electrical Power and Energy Systems 27 (2005) 254–263

261

Table 5 Values of m—Brazilian heavily loaded southeastern system

Table 8 Power flow results—Brazilian overloaded southeastern system

h

Iwamoto

Proposed

Method

0 1 2 3 4

0.9015 1.0434 1.0058 1.0029 1.0830

1.0323 1.0320 1.0019 1.0076 1.0670

V1137 (p.u.)

Iwamoto Proposed

0.5829 0.7394 0.6891 7.1293 9.7262 9.7262 44.452 7.6345 0.5762 0.4691 0.6845 6.9689 9.5921 9.5921 44.448 7.6295

q1137 (8)

q1139 (8)

V1139 (p.u.)

V1135 (p.u.)

q1135 (8)

PG501 (p.u.)

QG501 (p.u.)

Table 6 Power flow results—Brazilian heavily loaded southeastern system q1137 (8)

q1139 (8)

q1135 (8)

QG501 (p.u.)

Method

V1137 (p.u.)

Conventional Iwamoto Proposed

0.6238 6.1251 0.7266 11.682 0.7572 14.033 41.669 5.3288

V1139 (p.u.)

V1135 (p.u.)

PG501 (p.u.)

0.6238 6.1252 0.7266 11.682 0.7572 14.033 41.669 5.3288 0.6238 6.1263 0.7266 11.682 0.7572 14.033 41.669 5.3288

solution is achieved in four and three iterations, respectively, as shown in Table 9. Table 10 shows the best possible solution found through the two methods. Fig. 5 presents the cost function F for the three methods. 5.3. Computational effort Table 11 shows the ratio between the computational times to run both proposed and Iwamoto second order power

Fig. 4. Cost function F evolution—Brazilian overloaded southeastern system.

flow formulations, denoted by Tp and TI, respectively, under different operating conditions. The CPU times have been obtained using a personal computer ATHLON 1.3 GHz. It can be noted that the proposed second order formulation in terms of current injection equations is faster than the Iwamoto’s method. 5.4. 11-Bus and 43-bus test systems Table 12 presents a comparison between the values of m for ill-conditioned test systems. It can be noted that for Table 9 Values of m—Brazilian southeastern system/contingency

Fig. 3. Cost function F evolution—Brazilian heavily loaded southeastern system.

h

Iwamoto

Proposed

0 1 2 3 4

0.9997 0.9967 0.4114 0.1747 0.0261

1.0001 0.9985 0.2634 0.0020 –

Table 7 Values of m—Brazilian overloaded southeastern system h

Iwamoto

Proposed

0 1 2 3 4 5

0.8691 1.0477 1.0121 0.8588 0.6283 0.0070

1.0230 1.0429 1.0429 0.6951 0.0561 –

Table 10 Power flow results—Brazilian southeastern system/contingency q1137 (8)

q1139 (8)

q1135 (8)

QG501 (p.u.)

Method

V1137 (p.u.)

Iwamoto Proposed

0.5905 22.927 0.8018 33.097 0.8627 35.937 22.728 K5.019 0.5824 22.644 0.7976 33.008 0.8598 35.871 22.720 K5.024

V1139 (p.u.)

V1135 (p.u.)

PG501 (p.u.)

262

C.A. Ferreira, V.M. da Costa / Electrical Power and Energy Systems 27 (2005) 254–263

11-bus system m tends to zero in both methodologies. This fact means that the solution does not exist from the initial estimate. On the other hand, for 43-bus system Iwamoto’s and proposed methods converge easily for the solution at the end of the sixth iteration, being the cost function F equal to 2.799!10K9 and 4.180!10K9, respectively. Table 13 shows the voltages produced by second order and conventional rectangular methodologies at some buses of the ill-conditioned test systems. The conventional rectangular method does not converge for 11-bus system, and converges in seven iterations for 43-bus system. As can be seen, the final results are essentially the same, validating the proposed approach.

6. Conclusions Fig. 5. Cost function F evolution—Brazilian southeastern system/contingency case.

Table 11 Computational time System IEEE-300

Brazilian Southeastern

Operating condition

T1/TP

Base case Heavily loaded Overloaded Base case Heavily loaded Overloaded Contingency

1.238 1.355 1.368 1.315 1.396 1.694 1.652

Table 12 Values of m—Ill-conditioned test systems h

Step size optimization factor 11-Bus system

0 1 2 3 4 5 6 7 8

43-Bus system

Iwamoto

Proposed

Iwamoto

Proposed

1.0322 0.9202 0.6631 0.6910 0.9386 1.2400 0.7664 0.1752 0.0004

0.9649 1.3167 0.7808 0.8140 1.0047 1.1340 0.8191 0.1737 0.0021

0.6522 1.1043 1.0740 1.1333 1.0006 1.0011 0.9876 1.0011 1.0001

0.7456 1.0653 1.0526 1.2453 1.0428 1.0001 – – –

This paper proposes a new algorithm for the second order power flow calculation, using the current injection formulation in rectangular coordinates. The main advantage of this formulation lies on the calculation of Jacobian matrix, because its off-diagonal elements are exactly the terms of bus admittance matrix and the diagonal elements are easily calculated. This new methodology uses the current injection Jacobian matrix and the current mismatches for all the iterative process. As a consequence, the numerical values of the step size optimization factor, at each iteration, are different from those yielded by the Iwamoto’s method. On the other hand, both methodologies produce the same final results. The extra steps inserted for calculating the second order power flow solution represent a small additional computational effort to the original problem. These steps are easily incorporated into the conventional Newton Raphson power flow program. Besides, the proposed methodology not only maintains a good performance in terms of number of iterations and convergence characteristic, but also reduces the computational time for the solution process. By using step size optimization factor it is possible to find a solution for the load flow equations even when the conventional method fails. When there is no real solution, it provides the best possible solution and becomes possible to bring the system back to feasibility. In fact, this approach can be a useful tool for voltage instability and voltage collapse studies.

Table 13 Power flow results—Ill-conditioned test systems Method

Voltage magnitudes and angles 11-Bus system

Conventional Iwamoto Proposed

43-Bus system o

V5 (p.u.)

q5 ( )

V11 (p.u.)

q11 (8)

V20 (p.u.)

q20 (8)

V43 (p.u.)

q43 (8)

– 1.0337 1.0337

– K4.876 K4.876

– 1.0311 1.0314

– K25.217 K25.218

1.0938 1.0938 1.0938

K13.665 K13.665 K13.665

1.0694 1.0694 1.0694

K11.250 K11.250 K11.250

C.A. Ferreira, V.M. da Costa / Electrical Power and Energy Systems 27 (2005) 254–263

7. Future work The authors are actually working on the application of the step size optimization factor, obtained through current injection formulation, with the objective of defining control strategies in order to restore the power flow solution. This issue is a subject of another paper, which is under preparation by the authors.

Acknowledgements The authors would like to thank the Brazilian agencies CAPES and FAPEMIG for financial support.

References [1] Stott B. Review of load flow calculation methods. Proc IEEE 1974;62: 916–29. [2] Stott B, Alsac O. Fast decoupled load flow. IEEE Trans Power Syst 1974;93(3):859–69. [3] Bacher R, Tinney WF. Faster local power flow solutions: the zero mismatch approach. IEEE Trans Power Syst 1989;4(4):1345–54.

263

[4] Semlyen A. Fundamental concepts of a Krilov subspace power flow methodology. IEEE Trans Power Syst 1996;11(3):1528–37. [5] Fuerte-Esquivel CR, Acha E. A Newton-type algorithm for the control of power flow in electrical power networks. IEEE Trans Power Syst 1997;12(4):1474–80. [6] Gotham DJ, Heydt GT. Power flow control and power flow studies for systems with FACTS devices. IEEE Trans Power Syst 1998;13(1):60–5. [7] Fuerte-Esquivel CR, Acha E, Ambriz-Pe´rez H. A comprehensive Newton Raphson UPFC model for the quadratic power flow solution of practical power networks. IEEE Trans Power Syst 2000;15(1):102–9. [8] Wang L, Li XR. Robust fast decoupled power flow. IEEE Trans Power Syst 2000;15(1):208–15. [9] Semlyen A, De Leo´n F. Quasi-Newton power flow using partial Jacobian updates. IEEE Trans Power Syst 2001;16(3):332–9. [10] Expo´sito AG, Ramos ER. Augmented rectangular load flow model. IEEE Trans Power Syst 2002;17(2):271–6. [11] Da Costa VM, Martins N, Pereira JLR. Developments in the NewtonRaphson power flow formulation based on current injections. IEEE Trans Power Syst 1999;14(4):1320–6. [12] Da Costa VM, Pereira JLR, Martins N. An augmented NewtonRaphson power flow formulation based on current injections. Int J Electrical Power Energy Syst 2001;23(4):305–12. [13] Variz AM, Da Costa VM, Pereira JLR, Martins N. Improved representation of control adjustments into the Newton-Raphson power flow. Int J Electrical Power Energy Syst 2003;25(7):501–13. [14] Iwamoto S, Tamura Y. A load flow calculation method for illconditioned power systems. IEEE Trans Power Syst 1981;100(4): 1736–43.