An improved hybrid load flow calculation algorithm for weakly-meshed power distribution system

Electrical Power and Energy Systems 74 (2016) 437–445

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

An improved hybrid load flow calculation algorithm for weakly-meshed power distribution system Hongwei Li a,⇑, Yong Jin b, Anan Zhang a, Xia Shen a, Chao Li a, Bing Kong a a b

School of Electrical Engineering and Information, Southwest Petroleum University, Chengdu 610500, PR China School of Electrical Engineering and Information, Sichuan University, Chengdu 610065, PR China

a r t i c l e

i n f o

Article history: Received 20 March 2014 Received in revised form 5 August 2015 Accepted 7 August 2015

Keywords: Load flow Weakly meshed network Link branch power PV node Slack bus

a b s t r a c t An improved hybrid load flow calculation algorithm and its simplified and approximate version have been presented in this paper for weakly meshed distribution systems with PV nodes. The link branches, slack buses and PV nodes can be dealt with together with a combined matrix. The calculation accuracy of the proposed methods equals to that of Newton–Raphson method and loop theory based method. The method has relatively faster calculation speed because of the real-number matrix operations. Furthermore, because the coefficient matrices are constants for meshed networks, the proposed simplified and approximate version has less calculation time than other methods. The numerical tests proved that the new method provides a simpler, more robust, and potentially faster solution to handle complex meshed distribution networks with a high amount of dispersed generators. Ó 2015 Elsevier Ltd. All rights reserved.

Introduction As one of the most common computational procedures used in distribution management system, load flow calculation is a fundamental tool in planning, operation and control of distribution systems to analyze the steady-state performance of the systems. The increasing penetration of dispersed generators (DG), which are most based on the renewable energy sources, has great influence on the voltage quality, load flow and power loss etc [1]. So the analysis of distribution systems plays even more critical role than before and the complexity of analysis has significantly increased as well. Currently, Newton–Raphson method [2,3] is the most common method in the commercial and educational power system software for solving load flow problems. In comparison with other load flow methods such as Gauss–Seidel and Fast Decoupled Power Flow [4], Newton–Raphson method stands out with faster convergence speed (compared to Gauss–Seidel) and/or better applicability and accuracy (compared to Fast Decoupled Power Flow) [5]. So up to today, based on its high execution efficiency and good applicability to varied network structures and different component models, the Newton–Raphson method still holds significant advantages over other known approaches as the most mature load flow algorithm. However, the deficiencies of Newton–Raphson method lie in its

⇑ Corresponding author. E-mail address: [email protected] (H. Li). http://dx.doi.org/10.1016/j.ijepes.2015.08.004 0142-0615/Ó 2015 Elsevier Ltd. All rights reserved.

poor performance for ill-conditioned networks and heavy reliance on a good guess of initial condition [5]. The forward/backward sweep based method [6–11] is another most commonly used load flow solutions in distribution system. The method utilizes the radial or weakly-meshed topological feature of distribution networks to reduce dimension of calculations, and the programming is relatively simple for no large matrix calculation. First proposed in [6], this method was adapted and modified from current-iteration to power-iteration in [7], revised with better loop handling in [8], and added with more active component models in [9]. As the extension from [6,7], Refs. [10,11] present two improved forward backward/sweep based radial load flow method and were more effective and fast, but they did not consider the weakly meshed network. For meshed network handling, the major drawback with backward/forward method comes from the adoption of loop breakpoints—this operation makes it only applicable to simple topologies. Aiming at overcoming this drawback, some new solutions have been proposed. Ref. [12] presents a novel forward/backward sweep algorithm for meshed network based on loop analysis theorem. The method has excellent convergence characteristics even when the network becomes more meshed. Two load flow methods have been presented in [13,14] for weakly meshed distribution respectively. The concise formula to formulate the relationship between branch currents and node injections is developed and can be used to solve load flow problem directly. The two methods all have higher calculation efficiency and speed, but the method in [14] has a clear theory foundation and a more universal form.

438

H. Li et al. / Electrical Power and Energy Systems 74 (2016) 437–445

Nomenclature N l b aij A At Al Tt T tk Bt dU da dP dQ P Q P It

number of the independent nodes in the network number of link branches (independent loops) in the network number of the total branches in the network, b ¼ N þ l e element for row i and column j in A reduced node-branch incidence matrix in the network reduced node-branch incidence matrix corresponding to the tree branches in A reduced node-branch incidence matrix corresponding to the link branches in A path matrix corresponding to tree branches in the network node k corresponding row vectors(path vector) in T t loop matrix corresponding to tree branches in the network branch voltage drop branch voltage angle difference branch active power loss branch reactive power loss node consumption active power vector node consumption reactive power vector tree branch input active power vector

Most forward/backward sweep based load flow calculations belong to the category of the current iteration method. In [7,15], two power iteration methods were presented. The methods reduce the related computational effort by using the branch power flows instead of the branch complex currents as flow variables. But the method in [7] needs to set the loop breakpoints and this make it have more iteration and become less efficient when node voltage is lower or when network become more meshed. Compared with the method in [7], the method in [15] realized the decoupled calculation of the active and reactive power flows, so it has better convergence and faster calculation speed even when network becomes more meshed and has more PV nodes. However, when handling meshed network, the load flow results have errors for the link branch power calculation is not direct. In this paper, as an improved version of the method in [15], the calculation of link branch powers, which includes the loop-links, multiple slacks and PV nodes, were redeveloped. The improved method retains all the advantages of the original method, and it solves the error problem fundamentally. The features of the improved method are summarized as follows: (1) The model of tap changing transformer was introduced and can be dealt with like the distribution line. (2) Based on Kirchhoff’s voltage law and the accurate formulas to estimate the voltage drop and angle difference, the calculation of link branches powers has been re-deduced. (3) By constructing the artificial loops, the multiple slacks and PV nodes can be handled as loop links. The improved method realized the seamless integration of the loop links, multiple slacks and PV nodes into load flow calculation by combining them in one matrix form. (4) The method uses real-number matrix operations and has good convergence and accuracy.

Q It P Il Q Il P Isl Q Isl P IPV Q IPV Un

tree branch input reactive power vector link branch input active power vector link branch input reactive power vector slack bus injection active power vector slack bus injection reactive power vector PV node injection active power vector PV node injection reactive power vector node voltage magnitude vector node voltage angle vector branch voltage drop vector branch voltage angle difference vector loop-link voltage drop vector loop-link voltage angle difference vector artificial branch voltage drop vector for slack bus artificial branch voltage angle difference vector for slack bus artificial branch voltage drop vector for PV node the tree-branch resistance diagonal matrix the loop-link resistance diagonal matrix the tree-branch reactance diagonal matrix the loop-link reactance diagonal matrix

an

dU da dU l dal dU sl dasl dU PV Rt Rl Xt Xl

distribution line can be modeled as a pure p element which shown in Fig. 1 with corresponding electrical parameters [15]. In the networks, two nodes may be connected by a tap changing transformer. For a tap changing transformer with off-nominal turns ratio, the series impedance Z T (p.u.) is different from the viewpoint of transformer terminals. Therefore it must consider the effect of off-nominal transformation ratio and Z T must be modified. General model of such a transformer can be modeled as shown in Fig. 2. This figure represents a transformer with one series impedance and an ideal transformer (the core loss of transformer could be evaluated by the shunt core loss functions of the terminal voltage [16] and will not be discussed here). And impedance Z T is derived based on nominal turn ratio and b is off-nominal tap position in terms of p.u. From Fig. 2, the following equations in matrix form are true,

"

I_p
I_s

#

 ¼


1=ðbZ T Þ

The modeling of the serial components in distribution network is used as a basis of most load flow calculation algorithms. A typical

ð1Þ

Then from (1), the p model for tap changing transformer, in which series branch is modeled with impedance and equivalent shunt admittances ðY T ¼ 1=Z T Þ, can be illustrated as Fig. 3. For the distribution line model, the shunt admittance G
jB can be regarded as the constant admittance load. Their expended power (both active and reactive) can be obtained with terminal voltage as PGi ¼ GU 2i ; Q Gi ¼
BU 2i and PGo ¼ GU 2o ; Q Go ¼
BU 2o . Similarly, aiming at the transformers model, the equivalent shunt admittances in both sides are also considered as constant admittance loads, so a tap changing transformer can be dealt with like the distribution line. During calculating the load flow, those

Pi + jQi U i ,αi

Networks modeling

" _ # Up 1=ðb2 Z T Þ U_ s


1=ðbZ T Þ

1=Z T

G-jB

Fig. 1. Equivalent

I R+jX

Po + jQo

G-jB

U o ,αo

p circuit of distribution lines.

439

H. Li et al. / Electrical Power and Energy Systems 74 (2016) 437–445

1: β

ZT

Up

Designating one tree in the weakly meshed network and numbering the tree branches in the front and the link branches in the back, then A can be expressed as,

Us

Ip

A ¼ ½At ; Al 

Is

Fig. 2. Tap changing transformer model with a series impedance and an ideal transformer.

ð7Þ

The order of invertible At is N  N and the order of Al is N  m. The p-models as Figs. 1 and 3 are used during calculating load flow. Let superscript ‘T’ indicates matrix transpose. Define     Q It P It QI ¼ ¼ ½Q i1 ; Q i2 ; . . . ; Q iðNþmÞ T and P I ¼ ¼ ½P i1 ; Pi2 ; . . . ; Q Il PI l PiðNþmÞ T be the input reactive and active powers of branches,

U p ,α p

β ZT

β −1 YT β

Fig. 3. Equivalent

1− β

β2

P ¼ ½P 1 ; P 2 ; . . . ; PN T and Q ¼ ½Q 1 ; Q 2 ; . . . ; Q N T be the node consumption active and reactive powers. The node consumption powers include the loads, the constant admittance loads and the serial branch power loss in the p-models as Figs. 1 and 3. Then there exist [15],

U s ,αs

YT

p model for a tap changing transformer.

powers that consumed by constant admittance loads can be merged into the loads at the corresponding nodes. Based on Fig. 1, define voltage drop and angle difference as,



dU ¼ U i
U o

ð2Þ

da ¼ ai
ao

It is possible to deduce the accurate formulas to calculate differences in voltage magnitude across a branch through knowledge of terminal power (both active and reactive). Based on Fig. 1 there exist [15],

dU ¼

Rð2Pi
dPÞ þ Xð2Q i
dQ Þ Ui þ Uo

sinðdaÞ ¼

Xð2Pi
dPÞ
Rð2Q i
dQ Þ 2U o U i  
1 Xð2P i
dPÞ
Rð2Q i
dQÞ ) da ¼ sin 2U o U i Xð2P i
dPÞ
Rð2Q i
dQÞ  2U o U i

ð3Þ

dP ¼ P i
Po ¼ RI ¼

RðP2i þ Q 2i Þ

dQ ¼ Q i
Q o ¼ XI2 ¼

U 2i XðP2i þ Q 2i Þ U 2i

ð8Þ


1 T T Q It ¼ A
1 t Q
At Al Q Il ¼ Tt Q þ Bt Q Il

ð9Þ

where, Tt ðTTt ¼ A
1 t Þ is the path matrix corresponding to tree T branches and Bt (BTt ¼
A
1 t Al ¼
Tt Al ) is the loop matrix corresponding to tree branches [14]. From (8) and (9), the tree branch input powers P It ðQ It Þ consist of

two parts: P I1 ¼ TTt P (Q I1 ¼ TTt Q ) and P I2 ¼ BTt P Il ðQ I2 ¼ BTt Q Il Þ.P I1 ðQ I1 Þ can be viewed as the contribution of the node consumption powers P (Q), which is equivalent to do a backward sweep calculation for a radial network by breaking the link; And P I2 ðQ I2 Þ can be viewed as the contribution of the link branch input powers P Il ðQ Il Þ. Let U 0 be the power supply (reference node) voltage magnitude and U n be the node voltage magnitude vector of N  1. Define T

dU t ¼ ½dU 1 ; dU 2 ; . . . ; dU N  . So the voltage magnitude difference between power supply and any other node k equals to the summation of the tree branch voltage drop along the path of node k, that is,

ð4Þ

where, dP is the branch active power loss and dQ is the branch reactive power loss. Eq. (3) is an accurate formula to estimate the voltage drop. The dP and dQ in (3) and (4) can be calculated by 2


1 T T P It ¼ A
1 t P
At Al P Il ¼ Tt P þ Bt P Il

ð5Þ

ð6Þ

U nk ¼ U 0
Ttk dU t ðk ¼ 1; 2; . . . ; NÞ

ð10Þ

where, Ttk is the row vector corresponding to node k in Tt . Let the angle of power supply (reference node) be zero and an be the node voltage angle vector of N  1. Define dat ¼ ½da1 ; da2 ; . . . ; daN  . Similarly there exists, T

ank ¼
Ttk dat ðk ¼ 1; 2; . . . ; NÞ

ð11Þ

Supposing P Il and Q Il be known, the calculation steps of the proposed load flow algorithm can be found in [15]. For weakly meshed network, one key to calculate the load flow is to calculate P Il and Q Il . In [15], in order to calculate P Il and Q Il , the loop currents were calculated first and then P Il and Q Il were obtained indirectly. But this brings errors to the load flow results.

Realization of load flow algorithm

Calculation of link branch power

For a weakly meshed distribution network with N þ 1 nodes and m link branches (loops), there are N independent nodes when the power supply is referred as the reference node (the first node). The number of the total branches b equals N þ m. The reduced node-branch incidence matrix A can be used to describe the network and the elements of A are defined as:

It need to reconsider the calculation of P Il and Q Il . And here, the loop-links, multiple slacks and PV nodes are considered. A simple network is shown in Fig. 4 and discussed below.

8
1 Node i is the starting point of branch j: > > > < 1 Node i is the ending point of branch j: aij ¼ > 0 Node i is neither the starting > > : point nor the ending point of branch j:

Handling of loops Define the loop-links voltage drop vector and angle difference vector be dU l ¼ ½dU l1 ;dU l2 ;...;dU lm  and dal ¼ ½dal1 ;dal2 ;...;dalm  . Based on Kirchhoff’s voltage law for loops, it gives [14] T

Bt dU t þ dU l ¼ 0

T

ð12Þ

440

H. Li et al. / Electrical Power and Energy Systems 74 (2016) 437–445

ðkPl þ Bt kPt BTt ÞP Il
ðkQl þ Bt kQ t BTt ÞQ Il ¼
Bt ðkPt TTt P
kQt TTt Q
kCt Þ þ kCl ¼
Bt S 00t þ kCl

00

write as

! X 00 P Il
R00 Q Il ¼ b

ð19Þ

where S 0t ¼ cPt TTt P þ cQ t TTt Q
cCt and S 00t ¼ kPt TTt P
kQ t TTt Q
kCt . Rearranging (18) and (19) in matrix form has,



Fig. 4. A simple network with loops, PV node and slack node.

Bt dat þ dal ¼ 0

From (3) and (4), the voltage drop and angle difference for the kth branch can be derived as,

Rk ð2Pik
dPik Þ þ Xk ð2Q ik
dQ ik Þ U ik þ U ok 2Rk 2Xk Rk dPik þ Xk dQ ik ¼ Pik þ Q
U ik þ U ok U ik þ U ok ik U ik þ U ok

dU k ¼

write as

ðk ¼ 1; 2; . . . ; N þ mÞ ð14Þ

Xk ð2Pik
dPik Þ
Rk ð2Q ik
dQ ik Þ dak  2U ik U ok Xk Rk Xk dPik
Rk dQ ik ¼ Pik
Q
U ik U ok U ik þ U ok ik 2U ik U ok ! dak ¼ kPk Pik
kQk Q ik
kCk write as

cPk ¼ Uik2RþUk ok ,

cQk ¼ Uik2XþUk ok ,

ik



¼

dU slk ¼ U slk
U 0

0

b

#

ð20Þ

00

b

ðk ¼ 1; 2; . . . ; Nsl Þ

ð21Þ

dU sl ¼
Tsl dU t ¼
Tsl ½cPt ðTTt P þ TTsl P Isl Þ þ cQt ðTTt Q þ TTsl Q Isl Þ
cCt  ) ðTsl cPt TTsl ÞP Isl þ ðTsl cQt TTsl ÞQ Isl

ok

dU l ¼ cPl P Il þ cQl Q Il
cCl ¼
Bt dU t

¼
Tsl ðcPt TTt P þ cQ t TTt Q
cCt Þ
dU sl ¼
Tsl S 0t
dU sl 0

write as

! R0sl P Isl þ X 0sl Q Isl ¼ bsl

ð22Þ

dasl ¼
Tsl dat ¼
Tsl ½kPt ðTTt P þ TTsl P Isl Þ
kQt ðTTt Q þ TTsl Q Isl Þ
kCt  ) ðTsl kPt TTsl ÞP Isl
ðTsl kQt TTsl ÞQ Isl ¼
Tsl ðkPt TTt P
kQ t TTt Q
kCt Þ
dasl ¼
Tsl S 00t
dasl

00

write as

00 00 ! X sl P Isl
Rsl Q Isl ¼ bsl

"

ð23Þ

Rearranging (18) and (19) in matrix form has,

R0sl

X 0sl

#

P Isl

"



0

bsl

#


cCt 

ð16Þ

¼
Bt ½kPt ðTTt P þ BTt P Il Þ
kQ t ðTTt Q þ BTt Q Il Þ
kCt 

ð17Þ

Handling of PV nodes

ð18Þ

Power system generators are known to be operated under two basic modes: PQ and PV. A PV generator is characterized by an active power set point and a voltage magnitude set point, thus it can be seen as a mixture of slack and negative load. Due to the absence of reactive power set point, a PV generator cannot be modeled as a normal PQ node element. However, the missing voltage angle set point also makes it possible to transform a PV

T Pt ðTt P

þ

BTt P Il Þ

þc

"



For load flow calculations of transmission networks, it is quite often that multiple slack nodes need to be modeled. For each slack, a set point for voltage magnitude and a set point for voltage angle are respectively given, thus all the slacks can be seen as mutually correlated to each other by the differences in their set points. By arbitrarily choosing one slack node as total injection point, the remaining N sl slacks can be dealt with artificial loops by adding artificial lines that stems from the total injection point to the remaining slacks that have fixed values for voltage magnitude ðU sl Þ and voltage angle ðasl Þ. As shown in Fig. 4, one additional slack is transformed into one artificial loop (dashed line). Define (let the voltage angle of total injection point be zero),

kPk ¼ U XUk ,

kQk ¼ U RUk , kCk ¼ Xk dP2Uik
RU k dQ ik . ik ok  ik ok  c Define cP ¼ Pt c ¼ diagð½cP1 ; cP2 ; . . . ; cPðNþmÞ Þ, cQ ¼ Pl   cQ t c ¼ diagð½cQ 1 ; cQ 2 ; . . . ; cQ ðNþmÞ Þ,  Ql  c cC ¼ cCt ¼ ½cC1 ; cC2 ; . . . ; cCðNþmÞ T ,  Cl  kPt ¼ diagð½kP1 ; kP2 ; . . . ; kPðNþmÞ Þ, kP ¼ kPl   k ¼ diagð½kQ 1 ; kQ 2 ; . . . ; kQ ðNþmÞ Þ, kQ ¼ Q t kQl   kCt ¼ ½kC1 ; kC2 ; . . . ; kCðNþmÞ T . Then based (12) and (13), the kC ¼ kCl following equations can be obtained in vector and matrix forms,

¼
Bt ½c


R00

P Il Q Il

Except the slack bus of total injection point, extracting all other N sl slacks buses corresponding row vectors from path matrices Tt makes up a new matrix Tsl ðN sl  NÞ, then,

ðk ¼ 1; 2; . . . ; N þ mÞ

k dQ ik cCk ¼ Rk dPUikik þX , þU ok

X 00



daslk ¼ aslk

ð15Þ where

X0

Handling of multiple slacks

ð13Þ

! dU k ¼ cPk Pik þ cQk Q ik
cCk

R0

T Q t ðTt Q

þ

BTt Q Il Þ

X 00sl


R00sl

Q Isl

¼

00

bsl

ð24Þ

dal ¼ kPl P Il
kQl Q Il
kCl ¼
Bt dat

Rearranging (16) and (17) gives,

ðc

c

c c

c c

T T Pl þ Bt Pt Bt ÞP Il þ ð Ql þ Bt Q t Bt ÞQ Il T T ¼
Bt ð Pt Tt P þ Qt Tt Q
Ct Þ þ Cl ¼
Bt S 0t þ Cl

c

c

c

write as

0

0

! R P Il þ X 0 Q Il ¼ b

441

H. Li et al. / Electrical Power and Energy Systems 74 (2016) 437–445

generator into an artificial loop by adding one artificial line as shown in Fig. 4 (dashed line). But the artificial link contains only reactive power flow and is limited to the physical constraint from the voltage magnitude set point. That is,

dU PVk ¼ U PVk
U 0

ðk ¼ 1; 2; . . . ; NPV Þ

ð25Þ

Extracting all N PV PV nodes corresponding row vectors from path matrices Tt makes up a new matrix TPV ðN PV  NÞ. For the known active powers PPV from PV node is invariable, that is P IPV ¼ const, so

dU PV ¼
TPV dU t ¼
TPV ½cPt ðTTt P þ TTPV P IPV Þ þ cQ t ðTTt Q þ TTPV Q IPV Þ
cCt  ) ðTPV cQ t TTPV ÞQ IPV ¼
TPV ðcPt TTt P þ cQ t TTt Q
cCt Þ
TPV cPt TTPV P IPV
dU PV ¼
TPV S 0t
TPV cPt TTPV P IPV
dU PV

Fig. 6. A 69-bus test system with five loops and six PV nodes.

write as

! X PV Q IPV ¼ bPV

ð26Þ

During calculating the load flow, it need note that the output reactive power is limited at each PV node. If Q PVk at kth PV node violates the upper or lower limit, it will be set to upper or lower limits and then the PV node will be converted to PQ node during next iteration. Combination

kPt ðTTt P þ TTPV P IPV Þ
kQ t TTt Q
kCt , then the following equilibrium equations can be obtained,

R0R

X 00R

X 0R


R00R

#

Case no.

Closed loops

Case no.

Closed loops

No0 No1 No2

No loop Loop 1 Loop 1–2

No3 No4 No5

Loop 1–3 Loop 1–4 Loop 1–5

Table 2 Results for closed five loops (No5) with four methods (p.u.).

The loop-links, multiple slacks and PV nodes can be combined 2 3   Q Il P Il and dealt with together. Define P IlR ¼ , Q IlR ¼ 4 Q Isl 5, P Isl Q IPV 2 3 2 3 2 3T      T Bt Bt Bt c 0 B c 0 0 0 t Ql RR ¼ Pl þ 4 Tsl 5cPt þ 4 Tsl 5cQ t 4 Tsl 5 , , XR ¼ 0 0 0 0 Tsl TPV TPV TPV 2 3T      T     Bt Bt Bt c 0 B kPl 0 00 00 t Ql þ þ , RR ¼ c 4 Tsl 5 , k XR ¼ Tsl Pt Tsl Tsl Q t 0 0 0 0 TPV 2 3 2 3     cCl Bt B kCl 0 00 , S 0tn ¼ S 0t þ bR ¼
4 Tsl 5S 0tn þ 4
dU sl 5, bR ¼
t S 00tn þ Tsl
dasl TPV
dU PV  cPt TTPV P IPV S 0t ¼ cPt TTt P þTTPV P IPV þ cQ t TTt Q
cCt , S 00tn ¼ S 00t þ kPt TTPV P IPV ¼

"

Table 1 System configuration with different loops.

P IlR Q IlR

"

 ¼

0

bR 00 bR

#

 )

P IlR Q IlR

"

 ¼

R0R

X 00R

X 0R

#
1 "


R00R

0

bR 00 bR

#

ð27Þ So P IlR and Q IlR can be calculated with (27), and then,

Bus

Newton–Raphson method

Loop theory based method

New method ONE

New method TWO

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

0.997092 0.986238 0.982551 0.979100 0.971050 0.970077 0.968956 0.965663 0.965234 0.965234 0.965365 0.961966 0.960760 0.960406 0.958598 0.955067 0.953959 0.995332 0.980742 0.976657 0.972928 0.980736 0.970003 0.962650 0.970051 0.968786 0.963629 0.960136 0.956945 0.953827 0.953280 0.953498

0.997092 0.986238 0.982551 0.979100 0.971050 0.970077 0.968956 0.965663 0.965234 0.965234 0.965365 0.961966 0.960760 0.960406 0.958598 0.955067 0.953959 0.995332 0.980742 0.976657 0.972928 0.980736 0.970003 0.962650 0.970051 0.968786 0.963629 0.960136 0.956945 0.953827 0.953280 0.953498

0.997092 0.986238 0.982551 0.979100 0.971050 0.970077 0.968956 0.965663 0.965234 0.965234 0.965365 0.961966 0.960760 0.960406 0.958598 0.955067 0.953959 0.995332 0.980742 0.976657 0.972928 0.980736 0.970003 0.962650 0.970051 0.968786 0.963629 0.960136 0.956945 0.953827 0.953280 0.953498

0.997092 0.986236 0.982548 0.979097 0.971046 0.970075 0.968953 0.965660 0.965231 0.96523 0.965362 0.961963 0.960757 0.960403 0.958596 0.955061 0.953954 0.995331 0.980742 0.976655 0.972924 0.980733 0.969999 0.962645 0.970047 0.968781 0.963624 0.960131 0.956938 0.953820 0.953273 0.953493

2

P It ¼

Fig. 5. A 33-bus test system with five loops.

TTt P

3T Bt   6 7 P IlR ; þ 4 Tsl 5 P IPV TPV

2

Bt

3T

6 7
1 T Q It ¼ A
1 t Q
At Al Q Il ¼ Tt Q þ 4 Tsl 5 Q IlR TPV

ð28Þ

442

H. Li et al. / Electrical Power and Energy Systems 74 (2016) 437–445

Table 3 The convergence maximum voltage magnitude error under different loops with the proposed method. (10
6 , p.u.). System configuration Maximum voltage error

33-bus system 69-bus system

No0

No1

No2

No3

No4

No5

0

7

7

7

12

7

0

7

7

8

18

13

0

00

0

00

0

It need note that during calculating b ; b ; bsl ; bsl ; bPV ; bR 00

bR ,

and S 0t

¼

Rt ðTTt P

X t TTt Q

there þ

exist

TTPV P IPV Þ

cCk ¼

Rk dP ik þXk dQ ik 2

þ X t TTt Q
cCt and

k dQ ik and kCk ¼ Xk dPik
R , 2 00 T T S t ¼ Rt ðTt P þ TPV P IPV Þ



kCt .

Numerical tests Accuracy and efficiency analysis of load flow algorithm under good voltage conditions

Table 4 The convergence results for the 33-bus system. Load flow method

Item

System configuration with different loops No0

No1

No2

No3

No4

No5

The loop theory based method in [12]

Iteration number Calculation time (ms)

6

5

5

5

5

5

0.32

0.56

0.63

0.73

0.85

0.98

The method in [15]

Iteration number Calculation time (ms)

5

4

4

4

4

4

0.16

0.36

0.44

0.54

0.64

0.72

New method ONE

Iteration number Calculation time (ms)

5

4

4

4

5

5

0.22

0.93

1.08

1.18

1.46

1.61

New method TWO

Iteration number Calculation time (ms)

5

4

4

4

5

5

0.18

0.32

0.30

0.29

0.35

0.36

Simplifying of coefficient matrices If distribution network works under a good voltage condition, the voltage amplitude at each node is close to the rated value. So the terminal voltage amplitudes for the nodes across a branch can be approximatively set as rated value, that is U ik  U ok  1. Then the coefficients in (14) and (15) can be simplified as cPk  kQk  Rk ; cQk  kPk  Xk . Let Rt ¼ diagð½R1 ; R2 ; . . . ; RN Þ and Rl ¼ diagð½RNþ1 ; RNþ2 ; . . . ; RNþm Þ be the tree-branch and loop-link resistance diagonal matrix respectively, X t ¼ diagð½X 1 ; X 2 ; . . . ; X N Þ and X l ¼ diagð½X Nþ1 ; X Nþ2 ; . . . ; X Nþm Þ be the tree-branch and loop-link reactance diagonal matrix respectively. And then the coefficient matrices in (20), (24) and (26) are R0 ¼ R00 ¼ Rl þ Bt Rt BTt , X 0 ¼ X 00 ¼ X l þ Bt X t BTt , R0sl ¼ R00sl ¼ Tsl Rt TTsl , X 0sl ¼ X 00sl ¼ Tsl X t TTsl , X PV ¼ TPV X t TTPV which all are transformed into theconstant matrices. Moreover, the coefficient matrices in (27) can be also transformed into the constant matrices as following,

R0R ¼

00



Rl 0



RR ¼

Rl 0



Xl X 0R ¼ 0 X 00R ¼



Xl 0

3 Bt  T Bt 6 7 ¼ R0Rcon ; þ 4 Tsl 5Rt 0 Tsl TPV 2 3T Bt    0 Bt T 6 7 Rt 4 Tsl 5 ¼ R00Rcon ¼ ðR0Rcon Þ ; þ 0 Tsl TPV 2 3 2 3T Bt B  t 0 6 7 6 7 þ 4 Tsl 5X t 4 Tsl 5 ¼ X 0Rcon ; 0 TPV TPV     T 0 Bt Bt Xt ¼ X 00Rcon þ 0 Tsl Tsl 0



2

First, a 33-bus, five loops distribution system [17] and a 69-bus, five loops distribution system [18] are used to test accuracy of the proposed method (see Figs. 5 and 6). As listed in Table 1, six cases with different loops are used to test the ability of the different methods to deal with the meshed networks. For the propose method, the initial values of the voltage magnitudes of all nodes are set to be U 0 and the initial values of the branch power losses are set to be zero. The results from the proposed method with non-constant coefficient matrices (new method ONE, as described in Section ‘Co mbination’) and constant coefficient matrices (new method TWO, as described in Section ‘Simplifying of coefficient matrices’) are compared with those from Newton–Raphson method and loop theory based method [12]. The convergence accuracy is 10
6 p.u. (the same below) and the results are rounded to 6 digits after the decimal point. The results for closed five loops (No5) are listed in Table 2. In Table 2, there are consistent results for Newton–Raphson method, loop theory based method and new method ONE (proposed method with non-constant coefficient matrices). It can draw the same conclusion for six cases in Table 1 for the 33-bus system and the 69-bus system. Since Newton–Raphson method and loop theory based method have been generally acknowledged to be the best performing algorithms for networks in good condition, the accuracy of new method ONE has been validated. When the results from the loop theory based method are used as the benchmarks below, however, the results from new method TWO (the proposed method with constant coefficient matrices) have errors because it assuming all bus voltages being close to 1.0 p.u. Compared with the benchmarks, the maximum voltage magnitude error (p.u.) for new method TWO are listed in Table 3. The error terms are defined as absolute values of the differences between calculation results from two algorithms (the same below).

Table 5 The convergence results for the 69-bus system. Load flow method

item

system configuration with different loops No0

No1

No2

No3

No4

No5

The loop theory based method in [12]

Iteration number Calculation time (ms)

6

6

6

6

5

5

0.67

1.19

1.31

1.51

1.46

1.63

The method in [15]

Iteration number Calculation time (ms)

5

5

5

5

4

4

0.35

0.74

0.86

1.01

0.95

1.02

Iteration number Calculation time (ms)

5

5

5

5

4

5

0.53

1.77

2.19

2.44

2.01

2.36

Iteration number Calculation time (ms)

5

5

5

5

4

5

0.39

0.71

0.63

0.65

0.59

0.65

New method ONE

New method TWO

443

H. Li et al. / Electrical Power and Energy Systems 74 (2016) 437–445 Table 6 The convergence results for the 210 bus rural MV distribution network. Load flow method

Item

System configuration with different loops No loop

5 closed loops

10 closed loops

15 closed loops

The loop theory based method in [12]

Iteration number Calculation time (ms)

9 2.33

9 5.29

8 6.81

7 7.58

The method in [15]

Iteration number Calculation time (ms) Maximum voltage error (p.u.)

8 1.45 0

8 3.28 0.000696

6 3.34 0.004870

6 4.15 0.003189

New method ONE

Iteration number Calculation time (ms) Maximum voltage error (p.u.)

8 1.84 0

8 18.25 0

6 11.59 0

6 13.01 0

New method TWO

Iteration number Calculation time (ms) Maximum voltage error (p.u.)

8 1.53 0

8 2.46 0.000002

6 2.48 0.000003

6 2.08 0.000003

Table 7 The algorithm’s convergence properties under heavy load for six cases (see Table 1) in 69-bus system. Case

No0 No1 No2 No3 No4 No5

Load ratio

The loop theory based method

New method ONE

New method TWO

Iteration number

Calculation time (ms)

Iteration number

Calculation time (ms)

Iteration number

Calculation time (ms)

3.2 3.3 3.3 3.4 5.7 7.5

73 32 32 93 40 49

8.77 6.65 7.39 25.55 12.12 16.91

73 32 32 94 42 50

5.42 10.58 13.39 44.05 19.74 38.97

73 32 32 97 48 70

4.76 4.18 3.83 11.42 5.58 8.10

From Table 3, it shows that the maximum voltage magnitude error (p.u.) is 0.000 018 for the original loads. The error is very small and the new method TWO is capable of offering an enough accuracy level for networks under good voltage conditions. Compared with the maximum voltage magnitude error (p.u.), 0.0014 with the method in [15], the accuracy of new method TWO has been improved nearly 100 times. In order to examine the algorithm’s robustness and efficiency, loop theory based method and the method in [15] are used to compare with the proposed methods. The 33-bus system and the 69-bus system were solved by the four algorithms under different loops respectively. The total iteration numbers and calculation times are listed in Tables 4 and 5. It can be seen from Tables 4 and 5 that the four methods have the good and similar convergence property. But, the calculation times for four methods are very different. Except the radial network, new method TWO has the minimal calculation time, the second is method in [15], the third is loop theory based method in [12], and new method ONE has the maximum time. New method TWO and method in [15] use real reactive and active power rather than complex currents as flow variables and real operation is in place of complex operation, so they have faster calculation speed compared with loop theory based method. However, although new method ONE also uses powers as flow variables, it takes longer calculation time for meshed network because it need update the coefficient matrices and perform matrix inversion operation in each iteration. But for redial network, new method ONE has faster calculation speed than loop theory based method because of the real-number matrix operations. For new method TWO, because the coefficient matrices are constant, so the inverse matrices can be solved before the start of the iteration during the load flow calculation. Therefore the calculation time of new method TWO decreases greatly. In order to further examine the algorithm’s robustness and accuracy in large-scale distribution systems, the above four

Maximum voltage error(p.u.)

minimum voltage/ maximum voltage

0.000001 0.000196 0.000193 0.000858 0.008535 0.014931

0.501891 0.560235 0.560790 0.496553 0.515878 0.5112281

Fig. 7. A real 210 buses MV distribution network.

methods was applied to a MV utility distribution system shown as Fig. 7. The MV network is a real rural network with 2 distinct areas with different voltage levels: 15 kV and 30 kV and a

444

H. Li et al. / Electrical Power and Energy Systems 74 (2016) 437–445

Table 8 System configuration in 69-bus system with PV node. Case

Closed loops

PV nodes

Case1 Case2 Case3 Case4 Case5 Case6

No loop No loop Loop 1&4 Loop 1&4 Loop 1–5 Loop 1–5

PV1-3 PV1-6 PV1-3 PV1-6 PV1-3 PV1-6

30 kV/15 kV transformer at bus 2063. The 30 kV feeder has 73 buses and the 15 kV feeder has 137 buses including 15 link branches. The latter three methods have similar and better convergence characteristic than loop theory based method. The analysis of calculation time is accordant with the above conclusion. The results from the loop theory based method are used as the benchmark table too, and the maximum voltage magnitude error (p.u.) of new method TWO is no larger than 3e-6 for four cases listed in Table 6. However, the maximum voltage magnitude error (p.u.) is 0.004870 for the method in [15]. From the above analysis, it can draw the conclusion that new method TWO is the best choice for weakly meshed network under good voltage conditions because of its faster calculation speed and enough accuracy. Robustness and accuracy analysis of load flow algorithm under bad voltage conditions To further validate the algorithm’s ability to deal with ill-conditioned network, 69-bus system under heavy load for six network configurations was tested and the results are listed in Table 7.

In Table 7, the load values are the maximum ratio between the new load values and their original values. With even a very small increase (by 0.1 times) in the load value, the load flow calculation will diverge. From Table 7, the three methods all have a good ability to deal with ill-conditioned network. The new method ONE has similar convergence and accuracy characteristic compared with the loop theory based method. And although the convergence is less good, new method TWO has much faster calculation speed and need less calculation time. The maximum voltage magnitude error (p.u.) for new method TWO is less than 1.5% and is completely acceptable because of the ill-conditioned network. So it can come to the same conclusion that the new method TWO is the best choice. Efficiency analysis considering PV nodes In 69-bus test system, six PV nodes are added at node 88, 46, 34, 52, 14 and 23 respectively (see Fig. 6). Their rated active powers are 200 kW, 300 kW, 250 kW, 300 kW, 200 kW and 250 kW in turn. The lower and upper limits of output reactive power at each PV node are set to the value of the rated active power. The convergence accuracy is 10
6 p.u. The method in [15], new method ONE and new method TWO proposed in this paper are selected to calculate load flow. Six cases in Table 8 will be analyzed and discussed. The total iteration numbers and calculation time with the three methods are listed in Table 9. It can also be seen that the proposed two methods have better convergence, but new method TWO has minimum calculation time among three methods. The convergence output reactive powers at PV nodes solved with the three methods are shown in Table 10 for Case1 and Case4. Because the reactive power has upper/lower limits at each PV node, some PV nodes that violated the limits have been converted to PQ node during the iteration calculation. From Table 10, it can be seen that the numerical

Table 9 The convergence results in the 69-bus system with PV nodes. Load flow method

Item

System configuration with different loops Case1

Case2

Case3

Case4

Case5

Case6

The method in [15]

Iteration number Calculation time (ms)

6 0.60

6 0.55

5 0.93

8 1.45

5 1.30

5 1.35

New method one

Iteration number Calculation time (ms)

5 1.50

4 1.56

5 2.71

4 3.80

5 3.92

5 3.80

New method two

Iteration number Calculation time (ms)

5 0.45

4 0.44

4 0.85

5 0.90

5 1.13

5 1.07

Fig. 8. The voltage profiles for the 69-bus test system under Case4.

H. Li et al. / Electrical Power and Energy Systems 74 (2016) 437–445 Table 10 Part of output reactive power of PV nodes from convergence results with three methods/kVar. Case no.

PV nodes

The method in [15]

New method ONE

New method TWO

Case1

PV1 PV2 PV3

4.2 300 250


2.6 300 250


2.6 300 250

Case4

PV1 PV2 PV3 PV4 PV5 PV6

200 300 250 300
200 150.7

200 300 196.4 300
200
72.9

200 300 196.1 300
200
72.9

results with three methods are close except for individual bus. But for new method TWO is the simplified and approximate version of new method ONE, so the results of the two methods keeps consistent (see Table 10). The voltage amplitude profile from above three methods is shown in Fig. 8 for Case4. The voltage at each node has been improved with PV nodes and loops. It can be seen from Fig. 8 that the voltage level calculated with three methods is basically the consistent. Although the voltage level calculated with loop theory based method is slightly better, it needs a larger reactive power (total of leading and lagging reactive powers is 1400.4 kVar vs 1269.3 kVar). There are similar conclusions for other cases. Conclusions As an extension of the method in [15], this paper proposed an efficient and robust load method and its simplified and approximate version that overcome the relatively low precision of method in [15]. The link branches, slack bus and PV node can be processed together with a combined matrix. The proposed method with nonconstant coefficient matrices has the calculation accuracy which equals Newton–Raphson method and loop theory based method. And its simplified and approximate version with the constant coefficient matrices is accurate to one in 100000 in terms of p.u. under good voltage conditions. The methods have relatively faster calculation speed because of the real-number matrix operations, and the simplified and approximate version has minimum calculation time. So taking all those factors, the recommended better method is the proposed method with constant coefficient matrices. Some numerical tests verified the above conclusions. In a word, compared with existing main-stream methods (Newton–Raphson, forward/backward sweep based method etc.), the new load flow algorithms provide a simpler, more robust, and faster solution to handle complex distribution networks with a high amount of dispersed generators. The proposed load flow

445

algorithm is essentially still belongs to the loop-analysis based method and has a strong ability to deal with meshed network. It reduces the iteration number and has a faster calculation speed even when network becomes more meshed and has more PV nodes. The methods are also easy to extend to optimal power flow and probabilistic power flow problems. Acknowledgements This work is supported financially by the Chinese National Key Basic Research and Development Program(973) with item #2013CB228203, the Chinese National Nature Science Foundation with item #51107107 and the Chinese Postdoctoral Science Foundation with item #2014M562335. References [1] Pecas Lopes JA, Hatziargyriou N, Mutale J, et al. Integrating distributed generation into electric power systems: a review of drivers, challenges and opportunities [J]. Electr Power Syst Res 2007;77(9):1189–203. [2] Zhang F, Cheng CS. A modified Newton method for radial distribution system power flow analysis[J]. IEEE Trans Power Syst 1997;12(1):389–97. [3] Thukaram D, Banda HMW, Jerome J. A robust three-phase power flow algorithm for radial distribution systems [J]. Electr Power Syst Res 1999;50 (3):227–36. [4] Zimmerman RD, Chiang HD. Fast decoupled power flow for unbalanced radial distribution systems [J]. IEEE Trans Power Syst 1995;10(4):2045–51. [5] Keyhani A, Abur A, Hao S. Evaluation of power flow techniques for personal computers [J]. IEEE Trans Power Syst 1989;4(2):817–826,. [6] Shirmohammadi D, Hong HW, Semlyen A, et al. A compensation-based power flow method for weakly meshed distribution and transmission networks [J]. IEEE Trans Power Syst 1988;3(2):753–61. [7] Luo GX, Semlyen A. Efficient load flow for large weakly meshed networks [J]. IEEE Trans Power Syst 1990;5(4):1309–16. [8] Dragoslav R, Risto A, Rubin T. Voltage correction power flow [J]. IEEE Trans Power Syst 1994;9(2):1056–62. [9] Cheng CS, Shirmohammadi D. A three-phase power flow method for realtime distribution system analysis [J]. IEEE Trans Power Syst 1995;10(2): 671–9. [10] Chang GW, Chu SY, Wang HL. An improved backward/forward sweep load flow algorithm for radial distribution systems [J]. IEEE Trans Power Syst 2007;22 (2):882–4. [11] Singh Sachin, Ghose T. Improved radial load flow method [J]. Electr Power Energy Syst 2013;44:721–7. [12] Wu WC, Zhang BM. A three-phase power flow algorithm for distribution system power flow based on loop-analysis method [J]. Int J Elect Power Energy Syst 2008;30(1):8–15. [13] Teng Jen-Hao. A direct approach for distribution system load flow solutions [J]. IEEE Trans Power Deliv 2003;18(3):882–7. [14] Hongwei Li, An’an Zhang. Three-phase power flow solution for weakly meshed distribution system including PV type distributed generation [J]. Proc CSEE 2012;32(4):128–35. [15] Hongwei Li, Anan Zhang, Xia Shen, Jin Xu. A load flow method for weakly meshed distribution networks using powers as flow variables [J]. Int J Elect Power Energy Syst 2014;58:291–9. [16] Chen TH, Chen MS, Inoue T, Kotas P, Chebli EA. Three phase cogenerator and transformer models for distribution system analysis [J]. IEEE Trans Power Deliv 1991;6(4):1671–81. [17] Wang Shouxiang, Wang Chengshan. Modern distribution system analysis. Beijing: Higher Education Press; 2007. [18] Civanlar S, Grainger JJ, Yin H, et al. Distribution feeder reconfiguration for loss reduction [J]. IEEE Trans Power Deliv 1988;3:1217–23.