A decoupled extended power flow analysis based on Newton-Raphson method for islanded microgrids

Electrical Power and Energy Systems 117 (2020) 105705

Contents lists available at ScienceDirect

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

A decoupled extended power flow analysis based on Newton-Raphson method for islanded microgrids

T

⁎

Amir Ali Nazari, Reza Keypour , M.H. Beiranvand, Nima Amjady Faculty of Electrical and Computer Engineering, Semnan University, Semnan, Iran

A R T I C LE I N FO

A B S T R A C T

Keywords: Power Flow (PF) Decoupled Extended Newton-Raphson (DENR) method Droop Control Microgrid (MG)

Recently, microgrid (MG) power flow (PF) studies have gained a lot of attention due to the emergence of autonomous MGs which feature distributed generations (DGs). There are some inherent limitations in the islanded operation mode of MGs which cannot be addressed by conventional PF methods. In this regard, recent studies have proposed some updates to the conventional approaches by inclusion of network frequency as a variable in their modeling and omission of the slack bus from the grid. These considerations specifically in the NewtonRaphson (NR) method change the Jacobean matrix (JM) formulations. This paper concentrates on improving the previous NR methods for MGs which include droop controlled DGs. For this purpose, the partial derivatives of both calculated and scheduled powers are considered in the Taylor series expansion of bus power injections. Thus, the convergence and accuracy of the PF method are enhanced by adding these derivatives in modeling of generations, loads and losses. Moreover, these extended PF equations are decoupled by reformulating the JM based on the consideration that the lines in MGs are mostly resistive, which results in simplified JM calculations and improved convergence speed. The effectiveness of the proposed decoupled extended NR (DENR) method for MG PF analysis is illustrated in several case studies including 6-bus and 38-bus networks. Moreover, two convergence enhancement methods are also incorporated into the proposed approaches and their merits are investigated.

1. Introduction PF analysis is one of the most common and basic studies of power system which determines the steady state operating point of any network. Moreover, it is also required within the generation scheduling problems and generation expansion programs of any power system [1,2]. Two well-known and reliable PF analysis methods are Gauss–Seidel (GS) and NR methods [3]. With the emergence of small MGs which enable local generation and demand supply, the power system structure is being reshaped. The expansion of electrical MGs is favorable for decision makers as it leads to less operational costs and environmental pollution [4,5]. These MGs also lower the power losses and increase the penetration of renewable energy sources. Not only PF is a basic analysis tool for power systems and MGs, but also it is the base of some other more sophisticated planning and operation tools that require a large number of PF simulations. In the case of MGs, as they grow in size or form coalitions in multi-MG systems, the complexity and computation burden of the PF analysis also increase [6]. Moreover, in some operation tools, like optimal PF (OPF), where successive solution candidates are evaluated continuously, the capability of the

⁎

incorporated PF methods becomes crucial [7,8]. Another example is power distribution planning which consist of huge PF runs in its decision-making evaluations [9]. Thus, even in MGs, reducing the PF computation burden is an essential issue for the planning and operation. Typical MGs inherit the capability to operate in both islanded and grid-connected operation modes. In grid-connected mode, the voltage and frequency of the MG are determined by the main grid. Thus, the MG components act together as a whole in the point of common coupling to inject/absorb a specified quantity of power to/from the grid. In this case, the MG operates as a virtual power plant [10–12]. On the contrary, in islanded mode, it is the task of individual DG units to regulate the voltage and frequency of the MG by their controlling components [13,14]. One effective solution to this problem is to include a droop control strategy in DGs for voltage and frequency stabilization which is an imitation of the synchronous generators’ performance in this regard [15,16]. In conventional power systems, one bus which usually holds the largest power plant is considered as the slack bus. The generator in this bus is able to maintain the frequency of the system fixed by compensating the excess/deficit active power. In an MG, all DGs must operate together to avoid power mismatch in the MG. Thus,

Corresponding author. E-mail address: [email protected] (R. Keypour).

https://doi.org/10.1016/j.ijepes.2019.105705 Received 24 June 2019; Received in revised form 30 September 2019; Accepted 9 November 2019 0142-0615/ © 2019 Elsevier Ltd. All rights reserved.

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al.

speed of convergence in NR-based PF algorithms.

there is no slack bus in the grid and frequency is also a variable to be controlled. In [17,18], some models for PF analysis in islanded MGs have been proposed. In their models, the DG with the highest generation is considered as the slack bus generator and other DGs act as P-V or P-Q buses. These models are not realistic and verifiable in practical MGs, since the usual small DGs within MGs are not able to act as a regulating unit. PF analysis in islanded MGs has been carried out by meta-heuristic optimization algorithms in [19,20]. Frequency of the grid has been considered as a variable and different cases for DGs as P-V buses, P-Q buses, or droop controlled buses have been taken into account. A generalized PF algorithm intended for islanded MGs with balanced three phase loads has been introduced in [21]. A system of non-linear equations is formed and solved by Newton trust region method. However, the method can be sensitive to the initial parameters’ settings and the computational burden is increased enormously by the addition of more DGs. In [22], a backward/forward sweep method is proposed for solving PF in islanded MGs. In their approach, the voltage of one bus associated with a DG is selected as reference voltage so that the reactive powers of other DGs can be obtained accordingly. Its practicality can hardly be proved since the reactive power injections of DGs are to be determined according to the local voltage quantities. In [23], a revised GS method is presented for PF analysis in islanded MGs. Their framework includes the operational features of DGs with droop control strategy. However, the GS method shows poor convergence performance compared to NR alternatives. A NR-based PF method for islanded MGs is modeled in [24]. In that model, the frequency variability is added to the formulations and the slack bus is omitted. Two cases of P-V or P-Q bus types are considered for DG buses. However, the dependency of line impedances and loads on the frequency variations is not considered in their study. As a complement to the former research, another modification of the conventional NR method has been presented by the authors of [25]. In their modified NR approach, all DGs are considered as droop controlled units, and the variability of frequency and its effect on other system variables are taken into account. The main contributions of this research work can be summarized as follows:

The rest of this paper is organized as follows. In Section 2, the load model, Ybus matrix, and the DGs’ droop models are presented. In Section 3, the proposed PF analysis method with its components are introduced. The suggested convergence enhancement methods to be incorporated into the proposed PF approaches are presented in Section 4. The simulation results of the paper and the comparative studies are depicted in Section 5. Finally, Section 6 concludes the paper. 2. Problem modeling Different system components including loads, line impedances and DGs have to be modeled for PF analysis. In this section, the formulations required for PF implementation are briefly presented. 2.1. Load model Every load bus of the system has predefined or forecasted power consumption at the rated bus voltage and grid frequency. However, loads naturally respond to variations in the voltage or frequency. Since in PF modeling of MG, the system frequency is no longer kept fixed by a generator in a slack bus, and unexpected voltage variations are also viable, it becomes essential to consider loads’ responses to these variations. Therefore, the following formulations are defined for the active and reactive powers of loads. α

|V | PLK = PLK 0 ⎛ k ⎞ (1 + Kpf (ω − ωr )) ⎝ |Vr | ⎠ ⎜

⎟

(1)

β

|V | QLK = QLK 0 ⎛ k ⎞ (1 + K qf (ω − ωr )) ⎝ |Vr | ⎠ ⎜

⎟

(2)

where |Vk| is the voltage magnitude of k th bus and ω is the system frequency; |Vr| and ωr are respectively the rated bus voltage and the rated system frequency of the load; PLK and QLK are the active and reactive powers of k th bus at voltage magnitude |Vk| and system frequency ω ; PLK 0 and QLK 0 are the active and reactive powers of k th bus at voltage magnitude |Vr| and system frequency ωr ; Kpf and K qf are frequency sensitivity coefficients for respectively the active and reactive parts of the load. The Kpf and K qf values can be in the range of [0,3] and [−2,0], respectively [26]. Moreover, α and β depend on the load type and their values have been presented in [27] for specific loads (resistive, inductive and capacitive).

1. In this paper, the PF analysis method, presented in the aforementioned study, is improved. First by taking a deeper look into the model and the JM, the effect of slack bus omission and droop scheme on the behavior of islanded MGs is addressed by taking the variability of scheduled power injections into account. In other words, in this study, the partial differential equations for active and reactive powers, not only for the calculated powers but also for the scheduled ones, are appropriately modeled. This modeling approach is applied to the scheduled powers of both droop controlled buses and P-Q load buses, each of which responds to the deviations of frequency and local voltages from the rated values. To the best of the authors’ knowledge, the proposed PF analysis method, incorporating the additional effective partial derivatives, has not been addressed in the literature. Thus, this expansion in the NR method, which boosts iterative changes to system variables, improves the convergence of the algorithm by finding more optimal steps towards global optimum resulting in a lower computation burden. 2. The increased size of the JM leads to higher complexity and computation time. To remedy this problem, in this paper, by considering the resistive nature of the distribution system impedances, the decoupling of the JM is carried out to alleviate the sizable calculations in each iteration. Combining these two features leads to the proposed decoupled extended NR (DENR) method for MG PF analysis. 3. Two convergence enhancement methods called acceleration factor method (AFM) and Levenberg-Marquardt method (LMM) are incorporated into the proposed PF approaches to increase the convergence speed of PFs. These methods modify the obtaining of the system variables’ mismatch vector in specific ways to boost the

2.2. Admittance matrix model The system admittance matrix is a function of ω as:

⎡ Y11 (ω) ⋯ Y1N (ω) ⎤ ⋱ ⋮ ⎥ Ybus (ω) = ⎢ ⋮ ⎢Y (ω) ⋯ Y (ω) ⎥ NN ⎣ N1 ⎦

Ykn (ω) =

(3)

−1 − Zkn (ω) ∀ k ≠ n ⎧ ⎪ N −1 ⎨∑ Zkn (ω) ∀ k = n ⎪ k = 1, k ≠ n ⎩

where Zkn is the per-unit impedance of the branch connecting buses to each other.

(4)

k th

and nth

2.3. Steady-state model of DGs In a droop controlled system, every DG takes a portion of the overall power mismatch to be compensated. Every DG monitors the local signals of the system. Therefore extra communicational devices are not required which lowers the costs of the system and improves its reliability [28]. The droop controlled generation of each DG is dependent on 2

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al.

Fig. 1. Droop characteristics for resistive MGs.

3.1. Conventional NR formulations

its local voltage magnitude and the overall system frequency. It can be shown as the following equations:

ω = ω0 − mp . PG

(5)

First the conventional NR method is discussed and formulated in this section. Modifications will be made on it in the next sections. The problem variables are defined as the following vector:

|Vk| = |V0| − nq. QG

(6)

X = [δ T |V |T ]T

where δ is the voltage angle vector of all buses except the slack bus; |V | is the voltage magnitude vector of P-Q buses. The power injection values of network buses can be calculated according to the latest state of the system variables as shown in the following equations:

where PG and QG are DG’s active and reactive generations; mp and nq are respectively the droop coefficients for P -ω control and Q -V control schemes; ω0 and |V0| are no load frequency and no load voltage, respectively. These equations are based on the IEEE 1547.7 standard which has been defined for islanded MGs with inductive inverters [29]. However, their realization is based on an assumption that the output impedance seen from the inverter side is inductive [30,31] or a virtual impedance has been added to the inverter [32]. Considering the fact that the impedances of the lines within MGs are mostly resistive, these conventional droop controlling equations are not able to properly compensate MG power mismatches. For resistive MGs, droop control strategy may be defined according to the diagrams in Fig. 1. According to the droop control scheme portrayed in Fig. 1, the active and reactive power adjustments of the system should be reformulated to address the resistive nature of the MG structure [33]:

ω = ω0 + mq . QG

(7)

|Vk| = |V0| − np. PG

(8)

(9)

F (x ) = [PcT (X )QcT (X )]T

(10)

N

Pck (X ) = |Vk |

∑ |Ykn ||Vn| cos(δk − δn − θkn) n=1

(11)

N

Qck (X ) = |Vk |

∑ |Ykn ||Vn | sin (δk − δn − θkn) n=1

(12)

where F (X ) is the vector of calculated active and reactive powers of network; Pck and Qck are the active and reactive powers injected from k th bus to grid; |Ykn| and θkn are the magnitude and angle of the k th row- nth column element of Ybus matrix. Eqs. (11) and (12) yield the calculated power values at each iteration. These values can be compared with the predetermined or scheduled power values to obtain the PF mismatches or errors. The aim of the PF solution method is to minimize the mismatches. The G vector

Eqs. (7) and (8) indicate that the output active powers control the voltages of the inverters and the output reactive powers control the frequency of the MG.

G = [PsT

QsT ]T

(13)

consists of scheduled active and reactive powers of the system. Using the first-order Taylor series expansion of (10), the power mismatches can be calculated by defining a JM. Thus, the vector of the system variables X can also be updated iteratively.

3. Proposed PF analysis model The first step to model a proper PF is to determine bus types in the system. Load buses are modeled as P-Q buses in which the active and reactive powers are predetermined. Generator buses can be modeled as either droop-controlled buses or simple P-V buses. In Table 1, the known and unknown variables of each bus type are portrayed. As it is seen in this table, in a droop controlled bus, all δ , |V |, P , Q parameters are unknown. Table 1 Different bus types in MGs. Bus type PQ

PV

– P, Q, V, δ

P, Q V, δ

P, V Q, δ

(14)

X i + 1 = X i + J −1. ΔC (X )

(15)

∂F (X ) ⎤ J=⎡ ⎣ ∂X ⎦

(16)

ΔC (X ) = [G − F (X )]

(17)

where ΔC and J represent the mismatch vector and the JM, respectively. The JM can further be divided into some sub-matrices, each of which forms a matrix of partial derivatives with respect to |V | or δ :

Variable state

Droop

∂F (X ) ⎤ F (X ) + ⎡ . ΔX = G ⎣ ∂X ⎦

J J J = ⎡ 11 12 ⎤ ⎣ J21 J22 ⎦

Known Unknown

(18)

The formulations of the components in (18) have been given in 3

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al.

loads are constant, here these values are updated in each iteration as it is shown in (1) and (2). The total power losses of the system are also calculated as follows.

Appendix A. Eqs. (9)–(18) depict the conventional NR method in which one bus is considered as the slack bus whose voltage magnitude and angle remain unchanged. This bus conventionally includes a large-capacity generator with adequate rated power which can hold the system frequency constant by compensating the excess/deficit power of the whole system. In the conventional power systems this assumption usually holds. However, islanded MGs have a different situation as: (1) Due to the lack of a slack bus with a large unit, the excess/deficit power of an MG cannot be compensated by a single DG. (2) The frequency deviations due to the mismatch between supply and demand should be taken into account. (3) A new bus type called droop-controlled bus should be defined and modelled. In the next section, the proposed extended form of the conventional NR method for PF analysis in islanded MGs is introduced. Thereafter, its decoupled version is presented as the proposed DENR approach.

N

Ploss =

k=1 n=1

N

Qloss = − ∑

X =

T

|V |T ω |V1 |]

G (X ') = [PsT (X ')QsT (X ') Ptot (X ') Qtot (X ')]T

1 (|V0| − |Vk|) npk

QGk =

1 (ω − ω0) mqk

∂ (F (X ') − G (X ')) ⎤ (F (X ') − G (X ')) + ⎡ . ΔX ' = 0 ⎢ ⎥ ∂X ' ⎣ ⎦

Psys =

d

(19)

∂ (F (X ') − G (X ')) ⎤ J' = ⎡ ⎢ ⎥ ∂X ' ⎣ ⎦

∑ PGk = ∑ k=1

k=1

ΔC (X ) = [PsT − PcT QsT − QcT Ptot − Psys Qtot − Qsys ]T

d

Qsys =

⎡ J11 ⎢ ' ⎢ J21 J' = ⎢ ' ⎢ J31 ⎢ J' ⎣ 41

(21)

1 (ω − ω0) mqk

(23)

d

∑ QGk = ∑ k=1

k=1

(32)

' J12 ' J22 ' J32 ' J42

' ' J13 J14 ⎤ ' ' ⎥ J23 J24 ⎥ ' ' ⎥ J33 J34 ⎥ ' ' ⎥ J43 J44 ⎦

(33)

The detailed formulations of the JM sub-matrices, shown in (33), are presented in Appendix B. These additions and modifications build the proposed extended NR (ENR) model. In order to schematically highlight the specifications of the proposed ENR compared to the conventional NR, the structures of these two methods are compared in Fig. 2.

The two variables in (22) and (23), which indicate the overall power supply of the system, are added to the vector F :

F (X ') = [PcT (X ')QcT (X ') Psys (X ') Qsys (X ')]T

(31)

The third and fourth terms in (32) indicate the supply adequacy constraints of the whole system. The purpose is that, in any stable frequency and voltage values in the system, the power mismatches between supply and demand are kept close to zero. In other words, these two added terms simulate the function of slack bus, with a difference that here the frequency of the network is variable. The submatrices of the new JM matrix, shown in (31), have to be defined. Moreover, the formulations of the already existing sub-matrices must be modified to address the new partial derivatives added to the framework: '

(22)

(30)

The power mismatch vector ΔC is modified as: '

(20)

1 (|V0| − |Vk|) npk

(29)

Moreover, JM is also modified to address the added partial derivatives of vector G :

Using these droop equations, the total generated powers of the whole system can be obtained. d

(28)

Eq. (29) indicates that the components of vector G are now functions of system variables X ' . Hence, in the proposed extended NR method, the partial derivatives of power mismatches in the Taylor series expansion apply for both F and G components, rather than only for F components which was the case in the conventional NR method. The required modification is as:

Thus, |V1| and ω are added to the unknown system variables. In the conventional NR method, |V1| represents the voltage magnitude of slack bus. Therefore, due to being fixed at its predefined value, it is not considered in the vector X . On the contrary, in the extended NR model, |V1| is considered in the vector of system variables, since an islanded MG lacks slack bus feature. In order to have a criterion for determining the voltage angles, one arbitrary bus is chosen to have a δ value equal to zero. This bus can be a P-Q, P-V or droop-controlled bus. Here, bus 1 is selected as the reference bus for voltage angle. The voltage angles of other buses will be calculated with respect to this bus. In order to attain the injected power values in droop-controlled buses, by considering the resistive nature of MG lines and droop formulations provided in (7) and (8), the following equations are obtained for DGs:

PGk =

∑ |Ykn ||Vn ||Vk|{sin (θkn) cos (δn − δk )}

The vector G , which consists of the scheduled power values, has to be reformulated as well. In this regard, G(X ') is defined as:

The authors of [25] have suggested some changes in their modified NR (MNR) method to make the PF approach compatible for islanded MGs. Our proposed extended NR method is built upon their MNR formulation. To introduce the proposed extended NR method, the vector of variables X ' is redefined as:

[δT

(27)

N

k=1 n=1

3.2. The proposed extended form of the conventional NR method

'

N

∑ ∑ |Ykn ||Vn ||Vk|{cos (θkn) cos (δn − δk )}

3.3. Decoupling the extended NR

(24) Most PF approaches which are based on NR method yield sufficiently accurate results [34]. However, their main drawback is their slowness and high computation burden. This occurs due to the complexity of inverse JM calculations at each iteration. The ENR method, presented in the previous section, has added 12 more sub-matrices to the already large JM. Thus, the computation burden has increased considerably. Therefore, in this sub-section, a decoupled ENR (DENR) method is proposed to alleviate this problem and lower the JM size by decoupling it. In order to do so, the low voltage and resistive nature of

After obtaining Psys and Qsys , the total power demand of the system should be determined:

Ptot = Pload + Ploss

(25)

Qtot = Qload + Qloss

(26)

where Ptot and Qtot are total active and reactive power consumptions of the system, respectively, including the power losses of all lines. Unlike the conventional NR method in which the scheduled power values for 4

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al.

Fig. 2. Comparative structures of the proposed ENR and conventional NR.

P=

Vi Vj R

Q=−

cosδ −

Vi Vj R

V j2 (36)

R

sinδ

(37)

Moreover, another assumption can be made on the voltage angle differences. If δ is considered to be small, one can assume that sinδ = δ and cosδ = 1. Thus, the transferring powers can be further simplified as:

Fig. 3. A typical line connecting two buses in MG.

islanded MGs are taken into account. Fig. 3 shows a connecting line between two buses of a MG [35]. According to this figure the active and reactive powers transferring from bus i to bus j can be written as follows:

P≈

Vj R

Q≈−

V j2 ⎞

Vi Vj ⎛ Vi Vj cosδ − P=⎜ ⎟ cosθ + Z sinδ sinθ Z Z ⎝ ⎠

(34)

V j2 ⎞ Vi Vj ⎛ Vi Vj cosδ − Q=⎜ ⎟ sinθ − Z sinδ cosθ Z Z ⎝ ⎠

(35)

(Vi − Vj )

Vi Vj R

(38)

δ

(39)

According to the above equations, it can be deduced that in MGs which consist of resistive feeders, the active (reactive) power flows in lines are only affected by the voltage magnitudes (angles) of buses and are independent of voltage angles (magnitudes) of buses. Thus, by taking this simplification into account, the JM of the extended NR method can be decoupled. To do so, all ∂P and ∂Q terms are neglected. ∂δ

where δ indicates the angle difference between busi and bus j . The other variables in (34) and (35) are as shown in Fig. 3. In low voltage networks, line resistance is considerably higher than line reactance [36]. Thus, in MGs, it is reasonable to neglect the imaginary part of the impedances and consider them as pure resistive impedances (i.e., Z = R, θ = 0 ). By taking this assumption into account, the power transferring through lines can be reformulated as:

∂ | V|

' ' ' ' ' ' ' ' , J13 , J22 , J24 , J31 , J33 , J42 , J44 Thus, J11 sub-matrices can be neglected. Accordingly, the final decoupled form of the proposed DENR method can be portrayed as the following system of equations.

5

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al.

⎡⎡ ∂ (Pc2 − Ps2) ⎢⎢ ∂ | V2 | ⋮ ⎢⎢ ⎢⎢ ∂ (PcN − PsN ) ⎢⎣ ∂ | V2 | ⎢ ⎢ ⎡ ∂ (Psys − Ptot ) ⎢ ⎣ ⎣ ∂ | V2 |

The choice of λ value depends on the initial guess quality. If the initial point is close to the global optimum, higher values of λ are more appropriate [38]. In (46), I is the identity matrix.

∂ (Pc 2 − Ps2) ∂ | VN |

∂ (P − P ) ⎡ ∂c2| V |s2 ⎤⎤ 1 ⎢ ⎥⎥ ⎡ Δ |V2 | ⎤ ⎡ Ps2 − Pc 2 ⎤ ⋮ ⎢ ⎥⎥ ⎢ ⋮ ⎥ ⎢ ⋮ ⎥ ⎢ ∂ (PcN − PsN ) ⎥⎥ ⎢ ⎥ = ⎢ PsN − PcN ⎥ V Δ | | ⎥ | | N ∂ V 1 ⎣ ⎦ ⎢ ⎥ ⎥ ⎢ ⎥ Ptot − Psys ∂ (Psys − Ptot ) ⎥ ⎣ Δ |V1| ⎦ ∂ (Psys − Ptot ) ⎣ ⎦ ⎤ ⎡ ⎤ ⋯ ∂ | VN | ⎦ ⎣ ∂ | V1 | ⎦⎥ ⎦

⋯

⎤ ⎥ ⋱ ⋮ ⎥ ∂ (PcN − PsN ) ⎥ ⋯ ∂ | VN | ⎦

5. Simulation results In order to analyze the capability and effectiveness of the proposed DENR method, three different case studies for PF simulation are considered in this section. The results of the proposed DENR method are compared with the results of its non-decoupled variant (ENR), the results of the MNR method presented in [25], and the results of the Newton trust region (NTR) method [21], in terms of accuracy and runtime. Regarding the NTR method as a general minimization approach, it is noted that it relies on approximating the target function with a second order Taylor series expansion. Therefore, in its modelling, both the Jacobean and the Hessian matrices are required [39–41]. Its implementation as PF analysis [21] is carried out to minimize the mismatch power injections of each bus separately in a sequential order. Regarding the merit of convergence enhancement techniques, discussed in Section 4, their effects on the performance of the DENR, ENR and MNR methods are investigated. The source codes for the simulations have been implemented in the MATLAB software package using an Intel Core i7 processor with 12 GB of RAM. The stopping criterion for all cases is considered as 10−5 or lower variation in all system variables in successive iterations.

(40)

⎡⎡ ∂ (Qc2 −Qs2) ∂δ2 ⎢⎢ ⋮ ⎢⎢ ⎢⎢ ∂ (QcN − QsN ) ⎢⎣ ∂δ2 ⎢ ⎢⎡ ∂ (Qsys − Qtot ) ∂δ2 ⎢ ⎣⎣

∂ (Qc 2 −Qs2) ∂δN

∂ (Q − Q ) ⎡ c2∂ω s2 ⎤⎤ ⎢ ⎥⎥ ⎡ Δδ2 ⎤ ⎡ Qs2 ⋮ ⎢ ⎥⎥ ⎢ ⋮ ⎥ ⎢ ⎢ ∂ (QcN − QsN ) ⎥⎥ ⎢ ⎥=⎢ ∂ω ⎣ ⎦⎥ ⎢ ΔδN ⎥ ⎢ QsN ⎥ Q ⎢ ∂ (Qsys − Qtot ) ⎥ ⎣ Δω ⎦ ∂ (Qsys − Qtot ) ⎣ tot ⎤ ⎡ ⎤ ⋯ ∂ ω ∂δN ⎥ ⎦⎦ ⎦ ⎣

⋯

⎤ ⎥ ⋱ ⋮ ⎥ ⋯ ∂ (QcN − QsN ) ⎥ ∂δN ⎦

− Qc 2 ⎤ ⋮ ⎥ − QcN ⎥ ⎥ − Qsys ⎥ ⎦ (41)

In the proposed DENR method, the voltage angle deviations Δδ and frequency deviation Δω are controlled only by reactive power injections. Similarly, the voltage magnitude deviations Δ |V | are only controlled by active power Injections. Since the system variables’ vector X ' is not changed, the convergence capability of the DENR method is not significantly changed with respect to the non-decoupled variant, while the PF computation burden is significantly decreased in the DENR method. 4. Convergence enhancement

5.1. Case 1

There are several reported methods for boosting the speed of NR approach convergence. These methods are used to achieve a pre-specified deviation in the system variables within fewer iterations and shorter runtimes. Two of them are chosen in this paper and their performances are compared with each other.

In this case, a 6-bus islanded MG is simulated. Its nominal voltage magnitude, frequency and apparent power are respectively 127 V, 60 Hz and 10 KVA. The topology of this MG is illustrated in Fig. 4. This islanded MG has three identical DGs which operate as droop controlled units. The parameters of the MG, regarding the DGs, lines and loads, have been extracted from [21]. The specifications of DGs are listed in Table 2. In this case, α and β coefficients, presented in (1) and (2), are considered equal to zero which models the demands as constant loads. The Kpf and K qf parameters are also chosen equal to 1 and −1, respectively. The comparative PF results for case 1 are shown in Table 3. It is seen in the table that the proposed DENR, the proposed ENR, and the benchmark MNR and NTR methods have obtained the same or very close results for different PF variables. However, while all of these four methods have been implemented within the same software package and have been run on the same hardware set (i.e., have had the same conditions), the runtimes of the proposed DENR and ENR are lower than the runtime of the benchmark MNR and NTR. Moreover, the

1. AFM In this method, an acceleration factor is used to scale the displacement of system variables proportional to ΔX vector obtained at each iteration [37]. Eq. (15) is then modified as follows:

X i + 1 = X i + ρi . J (X i )−1 . ΔC (X i ) ρi = ρini + (ρfin − ρini )(

i ) imax

(42)

(43)

ρ is the acceleration factor. Values of ρ higher than 1 result in a faster convergence in the initial iterations and values lower than 1 boost the convergence at the final iterations. Therefore, the value of the acceleration factor changes adaptively during the algorithm progression. In (43), ρini and ρfin indicate the initial and final values of ρ , respectively, and imax denotes the maximum number of iterations. 2. LMM Considering ΔC (X ) as the error in power injections, ε is defined as the error vector at each iteration.

ε i + 1 = ΔC (X i + ΔX ) = ε i + J(X i ). ΔX

(44)

The sum of squared errors is defined as S = εT ε . The aim of LMM is to minimize such term plus deviation from current iteration as shown in the following equation:

SM = εT ε + λ (X i + 1 − X i )

(45)

where SM is the function to be minimized and λ is called the damping factor. Minimization of SM in (45) yields ΔX according to the following equation:

ΔX = [J(X i )T . J(X i ) + λ. I ]−1 . J(X i )T . ΔC (X i )

(46)

Fig. 4. The 6-bus islanded MG. 6

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al.

Table 2 Specification of DGs in the six-bus MG.

Table 4 Power injections of the MG obtained by DENR in Case 1.

DG location

np

mq

|V0|(p.u)

Pmax (p.u)

Qmax (p.u)

DG number

DG1

DG2

DG3

4 5 6

0.0059 0.0059 0.0059

0.0025 0.0025 0.0025

1.01 1.01 1.01

2 2 2

0.9 0.9 0.9

P (p.u) Q (p.u)

0.1409 0.8619

0 0.8619

1.7050 0.8619

Overall variables (p.u)

computation time of the DENR is lower than the computation time of the ENR. These comparisons illustrate the computational efficiency of the proposed DENR method. As it is seen from the table, the AFM and LMM have not been applied to the NTR method. It is due to the fact that the speed of convergence in the NTR, is controlled by updating the trust region radius at each iteration. The power injections of the MG, obtained from DENR method for case 1, are presented in Table 4. It is seen from Table 4 that all DGs generate active and reactive powers within their operational limits. Moreover, since the network frequency is constant in the MG and the droop parameters are the same for all DGs, the resulting reactive power injections are equal (see Eq. (7)).

Psys

Qsys

Pload

Qload

Ploss

Qloss

1.8459

2.5856

1.6648

2.4123

0.1749

0.1729

Their corresponding parameters are given in Table 7. Different load types including Residential (R), Commercial (C) and Industrial (L) loads have been considered for the load buses. The required data for these loads along with line impedances are taken from [42]. The results of PF analysis in this case are presented in Table 8. In this larger test case, again, all methods have obtained the same or close results, while the proposed DENR has the lowest computation time. The resulted power injections in the MG obtained by DENR in this case are shown in Table 9. The resulted power injections by DGs into the MG depend on both the droop coefficients (npk , mqk ) and the steady-state system variables. The reactive power shares of DGs are only affected by the droop coefficients (mqk ) considering that the frequency (ω ) is constant throughout the network. Thus, according to Table 9, the unit DG2 takes the largest share of reactive power injection. However, overall active power requirement of network is shared between DGs depending on both the local voltage magnitudes and the droop coefficients (npk ). Thus, by considering the results of voltage magnitudes for Busses 34–38 in Table 5 and the droop coefficients of DGs in Table 9, it can be deduced that the unit DG1 injects more active power compared to other DGs. The very close results of the three methods in cases 1–3 show that the proposed DENR method obtains reliable PF analysis results in different system sizes and with different types of loads. Moreover, the system frequency in the three cases is stabilized at certain values, usually a bit above 1 p.u, i.e., higher than ω0 (Fig. 1). The superiority of the proposed DENR method relies on its faster convergence speed compared to the MNR method presented in [25], i.e., it is 27–29% faster in 6-bus system and 15–33% faster in 38-bus system, depending on the convergence enhancement method used. The NTR method, due

5.2. Case 2 In this case, the dependency of loads in P-Q buses on bus voltage variations are taken into account as constant impedance loads by considering α and β coefficients equal to 2 [25]. All other assumptions are the same as in case 1. The PF results for this case are shown in Table 5. Again it is seen that, while the results of the four methods are the same or very close, the computation time of the proposed DENR method is slightly lower than the computation times of the ENR and MNR methods and considerably lower than that of the NTR under the same software package and hardware set. The power injection results obtained from DENR method in case 2 are presented in Table 6. 5.3. Case 3 In this case, a 38-bus microgrid, consisting of five separate DGs, is chosen for PF analysis. The configuration of this MG is depicted in Fig. 5. Again it is supposed that All DGs operate in droop controlled mode. Table 3 Comparative results for case 1 (α = 0, β = 0) .

Bus variables Bus number

1 2 3 4 5 6

System frequency (p.u)

Voltage angle (deg)

Voltage magnitude (p.u)

DENR

ENR

MNR

NTR

DENR

ENR

MNR

NTR

0 −0.1939 6.4237 −2.8788 −2.0608 6.2942

0 −0.1961 6.4207 −2.8832 −2.0657 6.2915

0 −0.1913 6.4087 −2.8729 −2.0521 6.279

0 −0.1941 6.43 −2.8804 −2.0611 6.2965

0.9785 1.0391 0.979 1.0092 1.0543 0.9999

0.9787 1.0391 0.9789 1.0092 1.0542 0.9999

0.9782 1.0391 0.979 1.0092 1.0545 0.9999

0.978 1.0385 0.9785 1.0093 1.0544 0.9999

DENR

ENR

MNR

Overall Variables NTR

1.0022

1.0022

1.0022

1.0022 Runtime (s)

Default AFM Applied LMM Applied

DENR

ENR

MNR

NTR

0.185 0.116 0.110

0.204 0.138 0.134

0.259 0.162 0.150

0.391 – –

7

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al.

Table 5 Comparative results for case 2 (α = 2, β = 2) . Bus variables Bus number

Voltage angle (deg)

1 2 3 4 5 6

Voltage magnitude (p.u)

DENR

ENR

MNR

NTR

DENR

ENR

MNR

NTR

0 −0.1908 6.1396 −2.757 −1.9836 6.0157

0 −0.193 6.1481 −2.7608 −1.9906 6.0249

0 −0.1792 6.1424 −2.7511 −1.9647 6.0188

0 −0.1907 6.1336 −2.7595 −1.9828 6.0192

0.9799 1.0381 0.9803 1.0092 1.0527 1.0004

0.9801 1.0381 0.9802 1.0092 1.0525 1.0003

0.9797 1.0382 0.9803 1.0092 1.053 1.0004

0.978 1.0385 0.9785 1.0093 1.0544 1.0008

Overall variables

System frequency (p.u)

DENR

ENR

MNR

NTR

1.0021

1.0021

1.0021

1.0021 Runtime (s)

Default AFM Applied LMM Applied

DENR

ENR

MNR

NTR

0.194 0.118 0.111

0.210 0.143 0.129

0.268 0.171 0.158

0.427 – –

Table 6 Power injections of the MG obtained by DENR in Case 2. DG number

DG1

DG2

DG3

P (p.u) Q (p.u)

0.1342 0.8252

0 0.8252

1.6350 0.8252

Overall variables (p.u)

Psys

Qsys

Pload

Qload

Ploss

Qloss

1.7692

2.4755

1.5988

2.3171

0.1610

0.1585

to its inherent complex formulations and sequential solution approach, is not competitive with NR based PFs in terms of runtime. For instance, its runtime is between three to five times higher than the DENR runtimes as the fastest method. In terms of convergence enhancement, both AFM and LMM significantly boost the convergence speed of the implemented NR methods. Meanwhile, the reduction of runtimes is slightly superior in the LMM compared to the AFM. 6. Conclusions In this paper, a novel PF analysis method based on NR approach for islanded MGs has been presented. The proposed method enhances the already existing modified NR approaches presented for islanded MGs. In the proposed PF analysis method, a droop control scheme has been implemented for DGs to control their injected active and reactive powers by monitoring the MG variables consisting of local voltage magnitudes and overall system frequency. Considering the fact that a droop controlled scheme leads to variable scheduled power generations for DGs, new partial derivatives have been added to the Taylor series expansion of bus power injections. These additional equations and partial derivatives improve the accuracy of the JM to estimate the required changes of variables in each iteration. Although the computation time of each iteration is increased in this approach, the overall runtime of the PF is decreased due to considerably reducing the number of iterations. For this reason, the computation time of the proposed ENR method is lower than the computation time of the previous MNR method presented for islanded MGs. To further decrease the runtime of the proposed PF analysis method,

Fig. 5. . The configuration of the 38-bus MG in Case 3. Table 7 DG Parameters in the 38-bus MG.

8

DG location

np

mq

|V0|(p.u)

Pmax (p.u)

Qmax (p.u)

34 35 36 37 38

0.02 0.0333 0.02 0.05 0.05

0.0051 0.0015 0.0045 0.0022 0.0022

1.01 1.01 1.01 1.01 1.01

2 2 2 2 2

0.9 0.6 0.9 0.3 0.3

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al.

Table 8 Comparative results for case 3. Bus variables Bus number

Voltage angle (deg)

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 33 34 35 36 37 38

Voltage magnitude (p.u)

DENR

ENR

MNR

NTR

DENR

ENR

MNR

NTR

0 0.0004 0.0121 0.0169 0.0135 0.0393 0.1116 0.0393 0.0221 −0.0727 −0.0984 −0.0992 −0.2069 −0.2754 −0.3131 −0.3312 −0.3913 −0.4101 0.0025 0.0406 0.0729 0.1668 0.0118 0.0137 0.0617 0.0503 0.0667 0.1194 0.1703 0.2453 0.1802 0.1611 0.1565 0.3394 0.2893 0.0539 0.2048 0.134

0 0.0001 0.0103 0.0165 0.0134 0.0326 0.1106 0.0363 −0.0256 −0.0786 −0.0909 −0.1163 −0.206 −0.2779 −0.313 −0.3363 −0.4072 −0.4175 0.0021 0.0394 0.071 0.1629 0.0094 0.0094 0.0556 0.0451 0.0633 0.1275 0.1863 0.2655 0.1831 0.1612 0.1542 0.3387 0.2833 0.0513 0.2093 0.1364

0 0.0002 0.0109 0.0172 0.0143 0.0339 0.112 0.0374 −0.0243 −0.0772 −0.0895 −0.1149 −0.2044 −0.2761 −0.3112 −0.3345 −0.4052 −0.4155 0.0023 0.0402 0.072 0.1642 0.0101 0.0105 0.057 0.0466 0.0651 0.1301 0.1895 0.2688 0.1866 0.1647 0.1578 0.3388 0.2876 0.0526 0.2108 0.138

0 0.0004 0.0110 0.0189 0.0165 0.0393 0.1125 0.0393 0.0121 −0.027 −0.0884 −0.0122 −0.2068 −0.2855 −0.3132 −0.331 −0.4011 −0.4001 0.0025 0.0406 0.0729 0.1667 0.0118 0.0121 0.0547 0.0503 0.0668 0.1195 0.1803 0.2451 0.1801 0.1610 0.1566 0.3592 0.2795 0.0839 0.2151 0.1350

0.9692 0.9692 0.9688 0.9693 0.9698 0.9718 0.974 0.9779 0.9793 0.9812 0.9815 0.9822 0.9773 0.9755 0.9743 0.9732 0.9716 0.971 0.9696 0.9739 0.9756 0.9793 0.9681 0.9675 0.9699 0.9713 0.9708 0.9691 0.9681 0.9648 0.9609 0.9601 0.9598 0.9935 0.9849 0.9968 0.9818 0.9731

0.9697 0.9697 0.9692 0.9697 0.9702 0.972 0.9741 0.9778 0.9789 0.9804 0.9807 0.9813 0.9758 0.9738 0.9725 0.9713 0.9695 0.969 0.9701 0.9744 0.9761 0.9797 0.9686 0.9679 0.9703 0.9716 0.9712 0.9696 0.9688 0.9657 0.962 0.9612 0.9609 0.9934 0.9856 0.996 0.9824 0.9735

0.9697 0.9697 0.9692 0.9696 0.9702 0.972 0.9741 0.9778 0.9789 0.9804 0.9807 0.9813 0.9758 0.9738 0.9726 0.9713 0.9695 0.969 0.9701 0.9744 0.9761 0.9797 0.9686 0.9678 0.9702 0.9716 0.9712 0.9696 0.9688 0.9656 0.962 0.9612 0.9609 0.9934 0.9856 0.996 0.9824 0.9735

0.9698 0.9700 0.9681 0.9701 0.9701 0.971 0.9736 0.978 0.9802 0.9821 0.9808 0.9831 0.9782 0.9755 0.9749 0.9725 0.9714 0.9718 0.9702 0.9748 0.9759 0.9784 0.9688 0.9683 0.9702 0.9718 0.9713 0.9689 0.9684 0.9642 0.9613 0.9592 0.9594 0.9926 0.9843 0.9972 0.9824 0.9728

Overall variables

System Frequency (p.u)

DENR

ENR

MNR

NTR

1.0021

1.0021

1.0021

1.0021

Runtime (s)

Default AFM Applied LMM Applied

DENR

ENR

MNR

NTR

0.204 0.168 0.162

0.242 0.183 0.175

0.297 0.205 0.191

0.814 – –

inverse. This leads to the superior performance of the proposed DENR method in resistive MGs with high R ratios. For this reason, the comX putation times of the proposed DENR method for all MG test cases of the paper are considerably lower than the computation times of both the proposed ENR and the benchmark MNR methods. Even in non-resistive networks, the non-decoupled variant of the proposed method, i.e., ENR, can still outperform the previous MNR method. The comparative results indicate that the proposed PF analysis methods save the accuracy of similar methods, while improving the convergence speed.

Table 9 Power injections of the MG obtained by DENR in Case 3. DG number

DG1

DG2

DG3

DG4

DG5

P (p.u) Q (p.u)

0.8248 0.4172

0.7556 0.6

0.6588 0.4724

0.5647 0.3

0.7373 0.3

Overall variables (p.u)

Psys

Qsys

Pload

Qload

Ploss

Qloss

3.5412

2.0896

3.4888

2.0413

0.0525

0.0489

Declaration of Competing Interest the resistive nature of MG has been taken into account to decouple the JM and to reduce the computation effort required for computing the JM

The authors declared that there is no conflict of interest. 9

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al.

Appendix A In the conventional NR method, the PF mismatch vector is as:

Δc = [PsT − PcT QsT − QcT ]T

(47)

⎡ Qc 2 ⎤ ⎡ Qs2 ⎤ ⎡ Pc 2 ⎤ ⎡ Ps2 ⎤ Ps = ⎢ ⋮ ⎥, Pc = ⎢ ⋮ ⎥, Qs = ⎢ ⋮ ⎥, Qc = ⎢ ⋮ ⎥ ⎢QcN ⎥ ⎢QsN ⎥ ⎢ PcN ⎥ ⎢ PsN ⎥ ⎣ ⎦ ⎣ ⎦ ⎦ ⎣ ⎣ ⎦

(48)

where Psk and Qsk are the scheduled active and reactive power injections of injections of k th bus. The sub-matrices of the JM are: ∂P ⎡ ∂δc2 ⋯ ⎢ 2 J11 = ⎢ ⋮ ⋱ ⎢ ∂PcN ⋯ ⎣ ∂δ2 ∂Q ⎡ ∂δc2 ⋯ ⎢ 2 J21 = ⎢ ⋮ ⋱ ⎢ ∂QcN ⋯ ⎣ ∂δ2

∂Pc 2 ∂δN

∂P ⎤ ⎡ ∂ | cV22 | ⋯ ⎥ ⎢ ⋮ ⎥, J12 = ⎢ ⋮ ⋱ ∂PcN ⎥ ⎢ ∂PcN ⋯ ∂δN ⎦ ⎣ ∂ | V2 |

k th

bus. Similarly, Pck and Qck are the calculated active and reactive power

∂Pc 2 ⎤ ∂ | VN |

⎥ ⋮ ⎥

∂PcN ⎥ ∂ | VN | ⎦

∂Qc 2 ∂δN

∂Q ⎤ ⎡ ∂ | Vc22 | ⋯ ⎥ ⎢ ⋮ ⎥, J22 = ⎢ ⋮ ⋱ ∂QcN ⎥ ⎢ ∂QcN ⋯ ∂δN ⎦ ⎣ ∂ | V2 |

(49)

∂Qc 2 ⎤ ∂ | VN |

⎥ ⋮ ⎥

∂QcN ⎥ ∂ | VN | ⎦

(50)

The partial derivatives existing in the sub-matrices of the JM can be written for diagonal and off-diagonal elements as follows [43]:

∂Pck = ∂δk

∑

|Vk ||Vn ||Ykn| sin(θkn − δk + δn ) (51)

n≠k

∂Pck = −|Vk ||Vn ||Ykn| sin(θkn − δk + δn ) ∂δn ∂Pck = 2 |Vk ||Ykk| cosθkk + ∂ |Vk |

∑

∑

(52)

|Vn ||Ykn| cos(θkn − δk + δn ) (53)

n≠k

∂Pck = |Vk ||Ykn| cos(θkn − δk + δn ) ∂ |Vn | ∂Qck = ∂δk

n≠k

n≠k

(54)

|Vk ||Vn ||Ykn| cos(θkn − δk + δn ) (55)

n≠k

∂Qck = −|Vk ||Vn ||Ykn| cos(θkn − δk + δn ) ∂δn ∂Qck = −2 |Vk ||Ykk| sinθkk − ∂ |Vk |

∑

n≠k

(56)

|Vn ||Ykn| sin(θkn − δk + δn ) (57)

n≠k

∂Qck = −|Vk ||Ykn| sin(θkn − δk + δn ) ∂ |Vn |

n≠k

(58)

Appendix B Compared to the JM of the conventional NR method in (18), which has only four sub-matrices, the JM of the proposed ENR method includes 16 sub-matrices as shown in (33). Each of these 16 sub-matrices is formulated in detail in the following: ∂ (P − P ) ⎡ c2∂δ s2 ⋯ 2 ⎢ ' ⋮ ⋱ J11 =⎢ ⎢ ∂ (PcN − PsN ) ⋯ ∂δ2 ⎣ ∂ (Q − Q ) ⎡ c2∂δ s2 ⋯ 2 ⎢ ' J21 =⎢ ⋮ ⋱ ⎢ ∂ (QcN − QsN ) ⋯ ∂δ2 ⎣ ' J31 =⎡ ⎣

∂ (Psys −Ptot )

' J41 =⎡ ⎣

∂ (Qsys − Qtot )

∂δ2

∂δ2

⋯

⋯

∂ (Pc 2 − Ps2) ∂δN

∂ (P − P ) ⎡ ∂c2| V |s2 ⋯ ⎤ 2 ⎢ ⎥ ' ⋮ ⋱ ⋮ ⎥, J12 = ⎢ ∂ (PcN − PsN ) ⎥ ⎢ ∂ (PcN − PsN ) ⋯ ∂δN ⎣ ∂ | V2 | ⎦

∂ (Pc 2 − Ps2) ∂ | VN |

⎤ ⎥ ⎥ ∂ (PcN − PsN ) ⎥ ∂ | VN | ⎦ ⋮

∂ (Qc 2 − Qs2) ∂δN

∂ (Q − Q ) ⎤ ⎡ ∂c2| V | s2 ⋯ 2 ⎥ ' ⎢ ⋮ ⋮ ⋱ ⎥, J22 = ⎢ ∂ (QcN − QsN ) ⎥ ⎢ ∂ (QcN − PsN ) ⋯ ∂δN ⎦ ⎣ ∂ | V2 | ∂ (Psys −Ptot ) ∂δN

' ⎤, J32 =⎡ ⎦ ⎣

∂ (Qsys − Qtot ) ∂δN

∂ (Psys −Ptot ) ∂ | V2 |

' ⎤, J42 =⎡ ⎦ ⎣

⋯

∂ (Qsys − Qtot ) ∂ | V2 |

∂ (Qc 2 − Qs2) ∂ | VN |

⎤ ⎥ ⋮ ⎥ ∂ (QcN − PsN ) ⎥ ∂ | VN | ⎦

∂ (Psys −Ptot ) ∂ | VN |

⋯

(59)

(60)

⎤ ⎦

∂ (Qsys − Qtot ) ∂ | VN |

(61)

⎤ ⎦

(62)

10

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al. ∂ (PcN − PsN ) T ' ⎤ , J14 ∂ω ⎦

' J13 = ⎡ ∂ (Pc2∂−ω Ps2) ⋯ ⎣

∂ (QcN −QsN ) T ⎤ , ∂ω ⎦

' J23 = ⎡ ∂ (Qc2∂ω−Qs2) ⋯ ⎣

=⎡ ⎣

∂ (Pc 2 − Ps2) ∂ | V1 |

∂ (PcN − PsN ) T ⎤ ∂ | V1 | ⎦

⋯

' = ⎡ ∂ (Q∂c2| −VQ| s2) ⋯ J24 1 ⎣

∂ (QcN − QsN ) ∂ | V1 |

(63) T

⎤ ⎦

(64)

∂ (Psys −Ptot ) ⎤ ∂ (Psys −Ptot ) ⎤ ' ' , J34 = ⎡ J33 =⎡ ⎥ ⎢ ⎢ ∂ ω ⎦ ⎣ ∂ |V1 | ⎥ ⎦ ⎣

(65)

∂ (Qsys − Qtot ) ⎤ ' ∂ (Qsys − Qtot ) ⎤ ' , J44 = ⎡ =⎡ J43 ⎥ ⎢ ⎥ ⎢ ∂ω ∂ |V1 | ⎦ ⎣ ⎦ ⎣

(66) ∂Pc , ∂Pc , ∂Qc ∂δ ∂ | Vk | ∂δ

∂Qc ∂ | Vk |

and are similar to those of the conventional NR method. However, the other The equations required for obtaining the terms partial derivatives in the aforementioned sub-matrices are new, which are mathematically formulated as:

∂ (Pc − Ps ) ∂X '

∂Pc

=

∂Pck = |Vk | ∂ω

−

∂X ' N

∂Ps ∂X '

,

∂ (Qc − Qs )

∂Qc

=

∂X '

−

∂X '

∂Qs ∂X '

(67)

∂ |Ykn | ∂θ |Vn| cos (δk − δn − θkn ) + kn |Ykn ||Vn| sin (δk − δn − θkn )⎤ ∂ω ⎣ ∂ω ⎦

(68)

∑⎡ n=1 N

∂Qck = |Vk | ∂ω

∂ |Ykn | ∂θ |Vn| sin (δk − δn − θkn ) − kn |Ykn ||Vn| cos (δk − δn − θkn )⎤ ∂ω ⎣ ∂ω ⎦

∑⎡ n=1

(69)

2 Xkn /ω ∂θ Xkn /(ωRkn ) ∂ |Ykn | =− , kn = − 2 2 3 ∂ω 1 + (Xkn / Rkn )2 (Rkn + Xkn ) 2 ∂ω

(70)

∂Psk ∂Qsk = =0 ∂δ ∂δ

(71) −1

∂P

= n ⎧ ∂ | Gk Vk | ⎪ pk ∂Psk = ⎨ ∂PLK = P ∂ |Vk | LK 0 ∗ ⎪ ∂ | Vk | ⎩

if k is a droop bus α |Vr |

|Vk | α − 1 (1 |Vr |

( )

+ Kpf (ω − ωr )) if k is a PQ bus

(72)

∂P

Gk = 0 if k is a droop bus ⎧ ∂ω ∂Psk = ∂P |V | α LK ⎨ ∂ω = PLK 0 |Vk| Kpf if k is a PQ bus ∂ω r ⎩

( )

⎧ ⎪ ∂Qsk = ⎨ ∂QLK = Q ∂ |Vk | LK 0 ∗ ⎪ ∂ | Vk | ⎩ ∂Q

∂QGk ∂ | Vk | β |Vr |

=0

(73)

if k is a droop bus

|Vk | β − 1 (1 |Vr |

( )

+ K qf (ω − ωr ))

if k is a PQ bus

1

Gk if k is a droop bus = m ⎧ ∂ω ⎪ qk ∂Qsk = |Vk | β ⎨ ∂QLK = Q ∂ω K qf if k is a PQ bus LK 0 |V | ⎪ ∂ω r ⎩

( )

∂ (Psys − Ptot ) ∂x

'

=

∂ (Qsys − Qtot )

∂PLoss = ∂ω

N

∂Psys ∂x

=

∂x '

'

−

∂Qsys ∂x '

∂ (Ploss + Pload ) ∂x

−

'

∂x '

N

k=1 n=1 N

=

∂ (Qloss + Qload )

∑ ∑ |Vn ||Vk| cos(δn − δk )

∂QLoss =−∑ ∂ω k=1

(74)

{

N

∂PLoss = 2 |Yii ||Vi | cos(θii ) + 2 ∂ |Vi |

∂x '

=

−

∂Qsys ∂x '

∂Ploss ∂x '

−

−

∂Qloss ∂x '

∂Pload ∂x '

−

(76)

∂Qload ∂x '

∂ |Ykn | ∂θ cos(θkn) − |Ykn | kn sin(θkn) ∂ω ∂ω

∑ |Vn ||Vk| cos(δn − δk ) n=1

∂Psys

(75)

{

(77)

}

(78)

}

∂ |Ykn | ∂θ sin(θkn ) + |Ykn | kn cos(θkn ) ∂ω ∂ω

(79)

N

∑

|Yki ||Vk| cos(θki )cos(δi − δk )

k=1 k≠i

∂QLoss = −2 |Yii ||Vi | cos(θii ) − 2 ∂ |Vi |

(80) N

∑

|Yki ||Vk| sin(θki )cos(δi − δk )

k=1 k≠i

(81)

11

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al. N

∂PLoss =2 ∂δi

∑

|Yki ||Vi ||Vk| cos(θki )sin(δk − δi )

k=1 k≠i

∂QLoss = −2 ∂δi

(82) N

∑

|Yki ||Vi ||Vk| sin(θki )sin(δk − δi )

k=1 k≠i

(83)

∂PLoad α ⎛ |Vk | ⎞ = PLK 0 ∗ |Vr | ⎝ |Vr | ⎠ ∂ |Vk | ⎜

α−1

(1 + Kpf (ω − ωr ))

⎟

β ⎛ |Vk | ⎞ ∂QLoad = QLK 0 ∗ |Vr | ⎝ |Vr | ⎠ ∂ |Vk | ⎜

(84)

β−1

(1 + K qf (ω − ωr ))

⎟

(85)

∂PLoad ∂QLoad = =0 ∂δi ∂δi N

(86) α

∂PLoad = ∂ω

∑ PLK 0 ⎛ |Vk | ⎞

|V |

∂QLoad = ∂ω

∑ QLK 0 ⎛ |Vk | ⎞

⎜

⎟

⎝

k=1

⎠

r

N

Kpf β

|V |

⎜

⎝

k=1

⎟

r

(87)

⎠

K qf

(88)

The partial derivatives of Psk and Qsk , presented in (71)–(75), are based on the droop control scheme applied to DGs which has been illustrated in (7) and (8) for resistive MGs. The same applies to Psys and Qsys derivatives which can be extracted according to the following equations.

∂Qsys ∂ |Vk | ∂Psys ∂δk

= 0,

= 0,

∂Psys ∂ |Vk |

=

−1 ⎧ npk

⎨ ⎩

∂Qsys ∂ |V1 | ∂Psys ∂ω

= 0,

∂Qsys ∂ω

d

=

∑ k=1

1 , mqk

∂Qsys ∂δk

=0 (89)

=0

(90)

if k is a droop bus 0

otherwise

(91)

Finally, it must be noted that the assumption of modelling feeders as resistive impedances in the decoupled variant of the proposed PF method was considered to yield simplified equations of bus power injections and consequently neglect some sub-matrices in the modified JM. However, according to Eqs. (11), (12) and Eqs. (46)–(49), the complex values of line admittances |Ykn | ∠θkn remain intact in both ENR and DENR methods.

[14] Rocabert J, Luna A, Blaabjerg F, Rodriguez P. Control of power converters in AC microgrids. IEEE Trans Power Electron 2012;27(11):4734–49. [15] De Brabandere K, Bolsens B, Van den Keybus J, Woyte A, Driesen J, Belmans R. A voltage and frequency droop control method for parallel inverters. IEEE Trans Power Electron 2007;22(4):1107–15. [16] Vovos PN, Kiprakis AE, Wallace AR, Harrison GP. Centralized and distributed voltage control: Impact on distributed generation penetration. IEEE Trans Power Syst 2007;22(1):476–83. [17] Kamh MZ, Iravani R. A sequence frame-based distributed slack bus model for energy management of active distribution networks. IEEE Trans Smart Grid 2012;3(2):828–36. [18] Kamh MZ, Iravani R. Unbalanced model and power-flow analysis of microgrids and active distribution systems. IEEE Trans Power Delivery 2010;25(4):2851–8. [19] Moradi MH, Foroutan VB, Abedini M. Power flow analysis in islanded Micro-Grids via modeling different operational modes of DGs: A review and a new approach. Renew Sustain Energy Rev 2017;69:248–62. [20] Abedini M. A novel algorithm for load flow analysis in island microgrids using an improved evolutionary algorithm. Int Trans Electr Energy Syst 2016;26(12):2727–43. [21] Abdelaziz MMA, Farag HE, El-Saadany EF, Mohamed YARI. A novel and generalized three-phase power flow algorithm for islanded microgrids using a newton trust region method. IEEE Trans Power Syst 2013;28(1):190–201. [22] Díaz G, Gómez-Aleixandre J, Coto J. Direct backward/forward sweep algorithm for solving load power flows in AC droop-regulated microgrids. IEEE Trans Smart Grid 2016;7(5):2208–17. [23] Mumtaz F, Syed MH, Al Hosani M, Zeineldin HH. A simple and accurate approach to solve the power flow for balanced islanded microgrids. In: Environment and electrical engineering (EEEIC), 2015 IEEE 15th international conference on 2015 Jun 10. IEEE. p. 1852–6. [24] Rese L, Costa AS, e Silva AS. A modified load flow algorithm for microgrids operating in islanded mode. InInnovative smart grid technologies Latin America (ISGT LA), 2013 IEEE PES conference on 2013 Apr 15. IEEE. p. 1–7. [25] Mumtaz F, Syed MH, Al Hosani M, Zeineldin HH. A novel approach to solve power flow for islanded microgrids using modified newton raphson with droop control of dg. IEEE Trans Sustain Energy 2016;7(2):493–503.

References [1] Bhowmik PS, Bose SP, Rajan DV, Deb S. Power flow analysis of power system using power perturbation method. Power engineering and automation conference (PEAM). IEEE; 2011. p. 380–4. [2] Brown HE, Carter GK, Happ HH, Person CE. Z-matrix algorithms in load-flow programs. IEEE Trans Power Apparat Syst 1968;3:807–14. [3] Grainger John J, Stevenson WD. Power system analysis. New York: McGraw-Hill; 1994. [4] Venkata SS, Pahwa A, Brown RE, Christie RD. What future distribution engineers need to learn. IEEE Trans Power Syst 2004;19(1):17–23. [5] Basak P, Saha AK, Chowdhury S, Chowdhury SP. Microgrid: control techniques and modeling. Universities power engineering conference (UPEC), 2009 proceedings of the 44th international. IEEE; 2009. p. 1–5. [6] Kargarian A, Falahati B, Fu Y, Baradar M. Multiobjective optimal power flow algorithm to enhance multi-microgrids performance incorporating IPFC. IEEE power and energy society general meeting; 2012. [7] Critical review of recent advances and further developments needed in AC optimal power flow. Electric Power Syst Res 136:57–68. [8] Abdi H, Beigvand SD, La Scala M. A review of optimal power flow studies applied to smart grids and microgrids. Renew Sustain Energy Rev 2017;71:742–66. [9] Khator SK, Leung LC. Power distribution planning: a review of models and issues. IEEE Trans Power Syst 1997;12(3):1151–9. [10] Pudjianto D, Ramsay C, Strbac G. Virtual power plant and system integration of distributed energy resources. IET Renew Power Gener 2007;1(1):10–6. [11] Braun M, Strauss P. A review on aggregation approaches of controllable distributed energy units in electrical power systems. Int J Distributed Energy Resources 2008;4(4):297–319. [12] Werner TG, Remberg R. Technical, economical and regulatory aspects of Virtual Power Plants. In: Electric utility deregulation and restructuring and power technologies, 2008. DRPT 2008. Third international conference on. IEEE; 2008. p. 2427–33. [13] Zeng Z, Yang H, Zhao R. Study on small signal stability of microgrids: a review and a new approach. Renew Sustain Energy Rev 2011;15(9):4818–28.

12

Electrical Power and Energy Systems 117 (2020) 105705

A.A. Nazari, et al.

[26] Kundur P, Balu NJ, Lauby MG. Power system stability and control. New York: McGraw-hill; 1994. [27] Payasi RP, Singh AK, Singh D. Effect of voltage step constraint and load models in optimal location and size of distributed generation. Power, energy and control (ICPEC), 2013 international conference on. IEEE; 2013. p. 710–6. [28] Han H, Hou X, Yang J, Wu J, Su M, Guerrero JM. Review of power sharing control strategies for islanding operation of AC microgrids. IEEE Trans Smart Grid 2016;7(1):200–15. [29] IEEE Standards Coordinating Committee. IEEE standard for interconnecting distributed resources with electric power systems. IEEE Std1547-2003; 2009. [30] Mohamed YARI, El-Saadany EF. Adaptive decentralized droop controller to preserve power sharing stability of paralleled inverters in distributed generation microgrids. IEEE Trans Power Electron 2008;23(6):2806–16. [31] Pogaku N, Prodanovic M, Green TC. Modeling, analysis and testing of autonomous operation of an inverter-based microgrid. IEEE Trans Power Electron 2007;22(2):613–25. [32] He J, Li YW. Analysis, design, and implementation of virtual impedance for power electronics interfaced distributed generation. IEEE Trans Ind Appl 2011;47(6):2525–38. [33] Guerrero JM, Matas J, de Vicuna LG, Castilla M, Miret J. Decentralized control for parallel operation of distributed generation inverters using resistive output impedance. IEEE Trans Ind Electron 2007;54(2):994–1004. [34] Lagacé PJ, Vuong MH, Kamwa I. Improving power flow convergence by Newton

[35] [36] [37] [38]

[39] [40]

[41] [42] [43]

13

Raphson with a Levenberg-Marquardt method. IEEE power and energy society general meeting-conversion and delivery of electrical energy in the 21st century; 2008. Bergen AR. Power system. Englewood Cliffs, NJ: AnalysisPrentice-Hall Inc; 1986. p. 529. Klaus Heuk, Klaus-Dieter Dettmann, Elektrische Energieversorgung, Vieweg, 3rd ed. Li G, Zhang XP, Wang XF. Accelerated newton–raphson power flow. Eur Trans Electr Power 2012;22(4):504–17. Lagacé PJ, Vuong MH, Kamwa I. Improving power flow convergence by Newton Raphson with a Levenberg-Marquardt method. In 2008 IEEE power and energy society conference; 2008. Yuan YX. A review of trust region algorithms for optimization. Iciam 2000;99(1):271–82. Conn AR, Gould N, Sartenaer A, Toint PL. Global convergence of a class of trust region algorithms for optimization using inexact projections on convex constraints. SIAM J Optim 1993;3(1):164–221. Byrd RH, Gilbert JC, Nocedal J. A trust region method based on interior point techniques for nonlinear programming. Math Program 2000;89(1):149–85. Singh D, Misra RK, Singh D. Effect of load models in distributed generation planning. IEEE Trans Power Syst 2007;22(4):2204–12. Saadat H. Power system analysis. 2nd ed. McGraw-Hill Higher Education; 2009.