A Nested-Iterative Newton-Raphson based Power Flow Formulation for Droop-based Islanded Microgrids

Electric Power Systems Research 180 (2020) 106131

Contents lists available at ScienceDirect

Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr

A Nested-Iterative Newton-Raphson based Power Flow Formulation for Droop-based Islanded Microgrids

T

Abhishek Kumara, Bablesh Kumar Jhaa,*, Dharmendra Kumar Dheerb, Rakesh Kumar Misraa, Devender Singha a b

Department of Electrical Engineering, IIT (BHU), Varanasi, India Department of Electrical Engineering, NIT Patna, India

A R T I C LE I N FO

A B S T R A C T

Keywords: Distributed Generations Droop Based Islanded Microgrid Power Flow Technique Nested-Iterative Newton-Raphson Modified Newton-Raphson Newton Trust-Region

In this paper, an iterative novel power flow technique is proposed to obtain the operating point of Droop Based Islanded Microgrid (DBIMG). The proposed technique considers system frequency as an additional variable to obtain the steady-state operating point of the system. To generalize the proposed technique, four operating modes of Distributed Generations (DGs) including droop control, isochronous, PV and PQ mode are considered. In this study, the formulated power flow problem consists of a set of non linear and linear power flow equations. To solve these set of equations, a Nested-Iterative Newton-Raphson algorithm is proposed. The proposed algorithm is implemented on several test systems including 6-bus, 22-bus, 38-bus, 69-bus, and 160-bus. To examine the robustness and effectiveness of the proposed algorithm, the power flow solutions obtained by implementing the proposed algorithm are compared with the power flow solutions obtained by implementing the existing Newton-Raphson algorithms including Modified Newton-Raphson (MNR), Newton Trust-Region (NTR) and timedomain simulator: PSCAD/EMTDC. The results show better efficiency and superior convergence of the proposed algorithm in comparison to the existing algorithms.

1. Introduction The smart grid architecture based electrical distribution system is mainly distinguished by higher penetration of the Distributed Generations (DGs) in the distribution network. In the microgrid, increasing penetration of DGs with adequate generation suffices the power (active and reactive) requirement for all or most of its loads [1,2]. The microgrid brings several advantages to the consumer as well as local distribution utilities including increment in reliability, grid support, reducing greenhouse gas emission, local employment, etc [3]. The microgrid is capable of operating in grid-connected mode, islanded mode or during the transition between these two modes of operation [4–6]. In the grid-connected mode, the voltage and frequency of the microgrid are governed by the main grid [7]. However, in the case of the islanded mode, microgrid has to update the voltage and frequency of the system according to its loading condition and the control technique applied to it. The conventional power flow algorithms (used in gridconnected mode) do not apply directly to the islanded microgrid due to the absence of the slack bus (reference bus)[8–10]. Droop features of microgrids are characterized by their load-

⁎

dependent frequency and voltage magnitude. Therefore, a DBIMG differs from conventional distribution systems in the following aspects: (i) It has no slack bus; (ii) the operating frequency is not constant and behaves as a power-flow variable; and (iii) the power output of DGs and load models are frequency-dependent. Altogether, these aspects complicate the formulation of the power-flow problem and limit the applicability of numerical techniques to find the solution for this problem, in comparison to the conventional AC systems. In [11] and [12], the unit of DG with the highest rating is treated as the infinite bus, and remaining DG units are considered as PQ/PV buses with a predefined active power and/or magnitude of bus voltage. However, the proposed model of DG for power flow analysis does not satisfy the droop control operation principle of the sources in an islanded microgrid.

• In the conventional/grid-connected power flow analysis, the slack

•

bus is assumed as an infinite bus. It is also assumed that due to infinite capacity of the slack bus, the magnitude of output voltage and the system frequency are constant. In the islanded mode of operation, DGs cannot behave/operate as a slack bus due to its smaller capacity. In the case of conventional/grid-connected power flow analysis, DG

Corresponding author.

https://doi.org/10.1016/j.epsr.2019.106131 Received 30 April 2019; Received in revised form 29 October 2019; Accepted 21 November 2019 Available online 09 December 2019 0378-7796/ © 2019 Elsevier B.V. All rights reserved.

Electric Power Systems Research 180 (2020) 106131

A. Kumar, et al.

solution of the islanded microgrid. Nevertheless, the application of these proposed BFS based algorithms is only limited to DBIMG with radial structure. Furthermore, the authors of these algorithms ignore the different operating models of DGs. In [7], [18] and [19], an algorithm based on swarm intelligence and evolutionary computation is proposed to solve the power flow problem of a droop-controlled islanded microgrid. In this paper, a generic approach to obtain the power flow solution of droop based islanded microgrids implementing a nested-iterative current-injection Newton-Rapshon approach (NICINR) is proposed. In this approach, any droop bus can be considered as an angle reference (AR) bus in the system. In every loop, the AR bus voltage magnitude, active and reactive power of other droop buses will be updated according to the droop law. In every iteration, all the bus voltage (excluding the AR bus) are updated according to the current-injection based power flow technique. Based on the output impedance seen by DGs and the R/X ratio of the network, three different droop characteristics (conventional droop, inverse droop, and mixed droop) are implemented into the proposed power flow approach. The proposed technique is implemented on several test systems including 6-bus, 22bus, 38-bus, and 69-bus system. To verify the execution of the proposed power flow approach on the power flow problem of islanded microgrids, the results are compared with simulation results obtained from a time-modal simulator PSCAD [22], MNR [16], and NTR [13]. Based on the above-mentioned limitations of the power-flow algorithms reported in the literature, the significant contribution of the paper can be summarized as follows.

units can be represented as PV/PQ buses. The voltage magnitude and/or power for these DGs units are pre-specified. In the case of an islanded mode of operation, the voltage magnitude and power are defined locally according to the droop characteristics. Hence, the DG units of DBIMG cannot be represented as PV/PQ buses. To obtain the power flow solution of the islanded microgrid, several techniques are proposed in the literature [13]-[14]. These techniques include system frequency as an additional variable to the system. In [14], a new methodology based on Nested-iterative Newton-Raphson power flow is proposed for droop based islanded microgrids. In this method, a dummy slack bus is considered, which can be connected at any bus within the system. It is assumed that, at the time of convergence, the power flow through the dummy slack bus will become zero. However, this approach does not provide convergence for the system having a high R/X ratio and thus not applicable to the network having a high R/X ratio. A Particle Swarm Optimization (PSO) based power flow solution is proposed in [15]. In [15], PSO has been utilized to calculate the optimal setting of the droop controller for optimal operation of DBIMG, and stability of the DBIMG system has also been studied. However, isochronous mode of operation of one of the DG is not considered in [15]. Modified Newton-Raphson (MNR) method is proposed in [16], which incorporates the non-linearity of the droop equation in the power flow solution. However, this algorithm is failed to provide a solution for the DBIMG system with a high R/X ratio of lines due to the ill-conditioning of that system. A Newton-trust region (NTR) based robust power flow algorithm for droop based islanded microgrid is proposed in [13], which addresses the issues of high R/X ratio and ill-conditioning of the systems. However, this algorithm requires high computational effort for the calculation and/or inversion of the Jacobian and Hessian matrices. In addition to their robustness and high computational cost, the authors of [16] and [13] have not considered the coupling between lines in their study. Direct calculation of the Jacobian and Hessian matrices of power flow equations is not possible due to the unavailability of the analytical function of the elements of the admittance matrix in the case of the coupling of lines. Furthermore, their models ignore the influence of the output impedance (shunt admittances) of the DGs on its droop characteristics and some existing models of the operation of DGs different from the PQ and PV buses. A brief review on the power flow methods of both the grid-connected and islanded MGs is presented in [1]. A method for the power flow analysis based on Swarm Intelligence (SI) is presented in [1]. Power flow method based on particle swarm optimization (PSO) with Gaussian mutation (GPSO-GM) is introduced in [3] to obtain the power flow solution of the renewable energy sources integrated DBIMGs. However, the formulation of the power flow problem is represented as an optimization problem, and the solution obtained using SI's algorithms may not be accurate for considering as a power flow solution. A forward-return forward-sweep (FRFS) power flow technique based on the system's branch structure is introduced in [4] to obtain the power flow solution of the droop controlled islanded MGs. However, the FRFS algorithm does not apply to meshed DBIMG. To calculate the power solution of small-sized DBIMGs, a Gauss-Seidel method is proposed in [17]. In [18], the drawback of the power flow technique proposed in [17] is reported. It is mentioned that the power flow technique proposed in [17] cannot work for the larger systems. In [19], an evolutionary algorithm based power flow method for the DBIMGs is proposed. Authors in [7] have proposed a generic power flow model to obtain the power flow solution of the DBIMGs. In [20], a nested backward/forward sweep algorithm is proposed and implemented to obtain the power flow solution of DBIMGs. In this power flow problem formulation, the system frequency, and voltage magnitude of the slack bus are considered as additional variables to the system. A backward-forward sweep based technique with matrix-based calculation is proposed in [21], [20], and [17] to obtain the power flow

1. To develop a generic, computationally-efficient power flow approach based on the current injection equations. 2. To develop an admittance matrix's gradient-free approach to calculate the steady-state value of the system frequency of the DBIMG. 3. To investigate the performance of droop-controlled DGs in the isochronous mode of operation of one of the DG in microgrids, and 4. To investigate the robustness and efficacy of the proposed approach for ill-conditioned DBIMG test systems. The organization of the paper is as follows: DG and load modeling is presented in section-2. Proposed power flow approach with load flow problem formulation is presented in section-3. In section-4, comparative analysis and results are presented to validate the accuracy and effectiveness of the proposed techniques in comparison to the existing techniques. In Section-5 conclusions and the scope of future work are presented. 2. Modeling of DGs and different loads The modeling of the power system and it's components is an important action that influences the power flow solutions. The DG and load models are presented in the following subsection. 2.1. Distributed generation modeling The penetration of DGs in the distribution system is increasing dayby-day. A wide category of sources (wind turbines, photovoltaic systems, energy storage systems, fuel cells, micro-turbines, etc.) is integrated with system in the form of different types of DGs. Consequently, there is a requirement of various modeling schemes to handle these DGs in the power flow analysis of DBIMG. The electrical characteristics of DGs have been determined based on the energy converter practiced in them. The different models of DG are proposed in literature based on electrical characteristics including constant voltage model (PV), voltage-dependent power model (PQ(v)), constant current model (PI), and constant power factor model (PQ), as shown in Fig. 1. The system frequency, voltage, and power-sharing among the sources in a DBIMG is controlled by the droop law applied to the DG 2

Electric Power Systems Research 180 (2020) 106131

A. Kumar, et al.

Fig. 1. Different connection of DGs and energy converters.

regardless of the connected load. Other DGs can operate according to their droop characteristics. To model isochronously controlled DG, the same droop equations can be applied by setting the droop coefficients equal to zero, i.e., |Vk| = |V0,k| and ω = ω0. In this paper, two parameter x and y are considered to transform all above-mentioned four droop equations of DGs in the form of a single composite droop equation which can be represented as,

units [23]. In DBIMGs, droop-controlled DG units are controlled to mimic the droop characteristics of synchronous generators operating in parallel. Generally, conventional droop law is applied to share proportional power among DG units in DBIMGs. Based on the output impedance seen by the DGs and R/X ratio of the line, three different droop equations (conventional droop, inverse droop, and mixed droop) are considered for the modeling of the DGs in [24,25]. For the inductive network, DGs operating in conventional droop [4,26,27] can be modeled as,

ω = ω0 − mk (xk (PG, k − P0, k ) − yk (QG, k − Q0, k )),

(7) (8)

|Vk | = |V0, k | − nk (QG, k − Q0, k ),

(1)

|Vk | = |V0, k | − nk (yk (PG, k − P0, k ) + xk (QG, k − Q0, k )).

ω = ω0 − mk (PG, k − P0, k ),

(2)

In the composite droop equations (7) and (8), the value of x and y is considered to be either 0 or 1. There can be four possible combinations of x and y, and each combination represents a different droop operation of the DGs. For an example, for the value of x = 1 and y = 0, the composite droop equations become the conventional droop equations (1) and (2). While {x, y} = {0, 1} and {x, y} = {1, 1} represents the inverse droop and mixed droop operation respectively. The last combination, {x, y} = {0, 0}, represents that the DG is acting as either an isochronously controlled or grid-connected mode.

where, PG,k and QG,k are the active and reactive power output of kth DG respectively. nk and mk are the voltage and frequency droop coefficients of kth DG respectively. Q0,k and P0,k are the power set points at kth DG respectively. For the resistive network, DGs operating in inverse droop [28,29] can be modeled as,

|Vk | = |V0, k | − nk (PG, k − P0, k ),

(3)

ω = ω0 − mk (QG, k − Q0, k ).

(4)

2.2. Load modeling In the static load model, active and reactive power absorbed by the load depends upon the bus voltage and system frequency. Voltage and frequency-based active and reactive loads can be represented using equations (9) and (10) respectively as reported in [20,30].

In a distribution system, practically the line resistance is comparable to the line reactance. The performance of the above-mentioned droop characteristics (for fully resistive and fully inductive case) is not justified in the case of comparable R and X of the line since the decoupling of the active and reactive power output of DG is impractical. Hence, a mixed droop [24,25] is considered which is expresses as,

|Vk | = |V0, k | − nk (PG, k + QG, k ),

(5)

ω = ω0 − mk (PG, k − QG, k ).

(6)

Pl, k = P0, l, k (ap + bp |Vk | + cp |Vk |2 + dp |Vk |α )(1 + ep (ω − ω0)),

(9)

Ql, k = Q0, l, k (aq + bq |Vk | + cq |Vk |2 + dq |Vk |β )(1 + eq (ω − ω0)),

(10)

where, ap + bp + cp + dp = 1 and aq + bq + cq + dq = 1. ω0 is the nominal frequency of the system, the frequency dependability coefficients ep and eq are for the active and reactive load respectively, ω is the operating system frequency, P0,l,k and Q0,l,k are the active and reactive load demand of kth-bus (when the system is operating at it's nominal

A dominant DG (typically synchronous generator based DG) can be operated in isochronous mode, thereby functioning as a non-ideal slack bus, i.e., providing constant frequency and voltage at its terminals 3

Electric Power Systems Research 180 (2020) 106131

A. Kumar, et al.

frequency and bus voltage), α and β are active and reactive power exponent respectively. Nowadays, a large number of plug-in hybrid electric vehicles (PHEVs) are connected to the distribution system. The influence of PHEV loads can be estimated by using different models based on their electrical characteristics, which mainly depends upon their different charging operation. In this work, three different voltage-dependent load models of PHEVs are considered for the power flow analysis of DBIMG. The first type of PHEV load model (PHEVI) is represented by a polynomial load model. The active and reactive power for the PHEVI can be represented as,

PEVI, k = P0,EVI, k (ap + bp |Vk | + cp |Vk

|2 ),

algorithms, unknown state variables (bus voltages) are determined by solving these node equations. Modifications in conventional algorithms are required for a power flow analysis of DBIMG because bus voltages and system frequency are regulated using droop controllers. Powers, frequency, and voltages can be calculated using droop equations (7) and (8). Therefore, droop equations should be included in the power flow formulation. The primary step of the inclusion of the droop equations in the formulation is to reclassify the buses within the system. 3.1.1. Classification of bus-types Four types of buses are proposed in this study, which are

(11)

and

QEVI, k = Q0,EVI, k (aq + bq |Vk | + cq |Vk

|2 ),

1. Droop bus: Droop based DGs are connected, where none of the state variables are pre-specified. 2. AR bus: One of the droop buses is considered as the AR bus, where the voltage angle is assumed to be zero that will serve as the reference angle for the rest of the buses in the system. 3. PV bus: Active power and voltage magnitude and pre-specified. 4. PQ bus: Active and reactive power are pre-specified.

(12)

where, the coefficients aq, bq, and cq are for the reactive power and ap, bp, and cp are for the active power. In the case of second type of PHEV load model (PHEVII), the active and reactive power are calculated according to equations (13) and (14), respectively. The modeling of a battery and a fast charging station are given in [31].

PEVII, k = P0,EVII, k (ap + bp |Vk |α ),

3.1.2. Node equations for the buses The second step is to obtain the node equations for each type of bus. For a bus, five node equations can be obtained, which can be used as mismatch equations for a bus. For k-th bus, these node equations are

(13)

where a, b, and α represent the power coefficient, exponent coefficient, and exponential index, respectively.

QEVII, k = PEVII, k tan(θ),

(14)

Nbus

∆Ir , k =

where θ is the power factor angle, and the value of the power factor is assumed to equal 0.97 [31]. The third type of PHEV load model (PHEVIII) can be represented by the constant current load model. In this case, active and reactive power are obtained using the polynomial load model. PEVIII,k of PHEVIII can be calculated as,

PEVIII, k = P0,EVIII, k |Vk |2 .

(G ki Vr , i − Bki Vm, i ) −

Vr , k Pksp + Vm, k Qksp = 0, Vr2, k + Vm2 , k

(17)

(Bki Vr , i + G ki Vm, i ) −

Vm, k Pksp − Vr , k Qksp = 0, Vr2, k + Vm2 , k

(18)

xk (ω − ω0) + yk (|Vk | − |V0, k |) = 0, mk xk + nk yk

(19)

i=1 Nbus

∆Im, k =

∑ i=1

∆Pk = (Pksp + Pl, k ) +

(15)

The reactive power is assumed to equal 0, i.e.

QEVIII, k = 0.

∑

∆Qk = (Qksp + Ql, k ) + (16)

yk (ω − ω0) + xk (|Vk | − |V0, k |) = 0, nk xk − mk yk

(20)

and

All three EV load models (PHEVI, PHEVII, and PHEVIII) can be handled using equations (9) and (10).

∆|Vk | = |Vk | −

3. Proposed Power flow algorithm for droop based islanded microgrid

Pksp /Qksp

Vr2, k + Vm2 , k = 0,

(21)

are the specified real/reactive power injected at the kth where, bus respectively. Gij/Bij are the real/imaginary part of {i, j}th element of bus-admittance matrix respectively. Vr,j/Vm,j are the real/imaginary part of voltage of kth bus respectively. Here, different values of {xk, yk} define the type of bus as discussed in section 2.1. The mismatch equations defined in equations (17)-(20) for each bus create a system of non-linear equations. This system of equations can be calculated using Newton-Raphson algorithm. The Newton-Raphson algorithm is widely implemented in [32,33,16,34,13] to obtain the power flow solution of the droop based islanded microgrid. In [13] and [16], the formation of Jacobian matrix requires derivatives of bus admittance matrix (Y-bus) with respect to system frequency because elements of Y-bus matrix change with the change of system frequency. In [13], the calculation of elements of Jacobian matrix has not been described, on the other hand, a brief description of elements of Jacobian matrix has been provided in [16]. But, the effect of mutual coupling of lines on elements of Jacobian matrix has not been addressed in [16]. In the presence of mutual coupling in lines connecting the buses, derivatives of elements of Y-bus matrix with respect to system frequency cannot be calculated analytically due to the inaccessibility of analytical function of elements of Y-bus matrix. In such situations, finite-difference based derivatives can be used in place of actual derivatives, but the performance of Newton-Raphson is highly depended upon the setting of step-size used in the finite-difference approach [35]. To choose the optimal step-size, numerous methods

In this section classification of buses, Nested-Iterative Current Injection Newton Raphson (NICINR) based power flow technique, modeling of droop bus as a PQ bus, and iterative approach to update the system frequency and voltage of Angle Reference (AR) bus are discussed. 3.1. Overview of proposed power flow algorithm Conventional power flow algorithms have been well defined for the conventional distribution system. The main objective of these algorithms is to calculate the steady-state value of bus voltages in the distribution network. The main requirement of power flow analysis using conventional algorithms are as follows. 1. Magnitude value of voltage with phase angle is known for the main feeder (slack bus). 2. Magnitude value of voltage and real power is known for the regulated buses (PV buses). 3. Real and reactive power are known for load buses (PQ buses). Using network analysis, two node equations based on power balance or current balance are estimated for each bus. In conventional 4

Electric Power Systems Research 180 (2020) 106131

A. Kumar, et al.

have been introduced in [35]. However, after employing these methods in the Newton-Raphson algorithm, the computational complexity is increased at a point where power flow analysis of large-sized DBIMG system will become impracticable using these approaches. It is worth to note that some non-derivative based load flow techniques have also been introduced in [21,17]. These algorithms employ an outer and inner loop-based approach for adjusting the value of the voltage of slack bus and system frequency outside the main power flow routine. However, the performance of these algorithms degrades in the case of a high R/X ratio of the lines and high loading conditions [20]. In the case of the DBIMG system with mutually coupled lines, adjustment of the frequency outside the routine of the Newton-Raphson power flow algorithm provides better performance in comparison with finite-difference based techniques. However, an efficient and robust technique to update the frequency has not been discovered and employed with existing derivative-based algorithms in case of a power flow analysis of the DBIMG system with mutually coupled lines. Therefore, to address this issue, a new nestediterative approach has been proposed in this paper to update the system frequency and voltage of AR bus outside the Current Injection based Newton-Raphson (CINR) [36] routine. To define the steps of the proposed algorithm NICINR, some premises need to be addressed, which are as follows.

Ql, k = QGvar ,k =

yk ω + xk |Vk | . nk xk − mk yk

(27)

Equations (24) and (25) give the pre-specified active and reactive power at droop buses, which is a necessary condition to become a PQ bus in the system. Voltage and system frequency based active and reactive load models [37] are presented in equations (26) and (27) respectively. 3.3. Current-injection based Newton-Raphson power flow In the rectangular coordinates, for the CINR power flow method, the real and imaginary part of the voltage of the AR bus is fixed at 1 and 0, respectively. All other PQ buses are initialized with a real and imaginary part of voltage with a value of 1 and 0, respectively [37]. All the initialized values will change in every iteration. The active and reactive current mismatch will be estimated in every iteration, and by means of the calculated Jacobian matrix, the real and imaginary part of voltage will be updated. The procedure to obtain ▵Ir and ▵Im for different types of buses are given in [37]. The Jacobian matrix is obtained by differentiating the equations (17) and (18) with respect to the real and imaginary parts of all the bus voltages which are represented as,

J 11 J 12 ⎤ J=⎡ , ⎣ J 21 J 22 ⎦

• The droop bus requires to be expressed as a PQ bus after considering system frequency as a constant. • A CINR algorithm needs to modify for solving the mismatch equations of each bus except the AR bus. • There is no slack bus in the system, and the voltage magnitude of the

(28)

where, J11, J12, J21, and J22 are the sub-matrices of the Jacobian matrix [37]. After calculation of the Jacobian Matrix for the (i + 1)th iteration, the real and imaginary part of the bus voltages can be calculated as,

AR bus (providing reference to the voltage angle) is not constant. In addition, system frequency is also a state variable. So, a nestediterative approach needs to define for updating the system frequency and voltage magnitude of the AR bus outside the CINR routine.

v i + 1 = v i + J −1∆I,

(29)

where,

v = [VrT VmT]T .

(30)

3.2. Modeling of droop bus as a PQ bus

Vr/Vm are the vectors of real/imaginary part of voltages of all the buses (except the AR bus). The value |V| of each bus is calculated using equation (21).

In order to express the droop buses as a PQ bus, the scheduled active and reactive power should be pre-specified. Active and reactive power output of kth DG at system frequency ω and bus voltage |Vk| can be obtained as per equations (22) and (23) respectively.

3.4. Iterative approach to update the system frequency and voltage of AR bus

PG, k =

xk (ω0 − ω) + yk (|V0, k | − |Vk |) var = PGconst , k − PG, k , mk xk + nk yk

(22)

QG, k =

yk (ω0 − ω) + xk (|V0, k | − |Vk |) var = QGconst , k − QG, k . nk xk − mk yk

(23)

The proposed nested-iterative methodology is shown in Fig. 2 and following steps are performed to obtain the power flow solution of droop based islanded microgrids using NICINR. 1. Initialize the operating condition of the system. One of the droop buses is selected as an AR bus, and for the remaining droop buses, the active and reactive power sharing are calculated according to equations (24) and (25) respectively. Other initial conditions are following the same procedure of current injection Newton-Raphson (CINR) except the initialization of system frequency, ωgrid, and voltage magnitude of AR bus, Vs. ωgrid and Vs are initialized on their nominal value ω0 and V0,s respectively. 2. Calculate Y-bus of the network, which includes system frequency (ωgrid) as a variable to the system. The frequency will be updated in every iteration until the convergence criteria are satisfied. 3. Implement the CINR to obtain the power flow solution. The main reason behind this step is to determine the total real and reactive power injected at all the droop buses, including the AR bus. The information of injected real and reactive power is needed to update the operating frequency of the system and voltage magnitude of the AR bus. 4. The operating frequency of the system is updated according to equation (31).

var Equations (22) and (23) consist two parts in active (PGconst , k and PG, k ) and const var reactive (QG, k and QG, k ) power generation of DGs. First part of both the const active (PGconst , k ) and reactive (QG, k ) power is constant. The scheduled active and reactive power (Pg,k and Qg,k) at k-th bus can be calculated as,

Pg, k = PGconst ,k =

Qg, k = QGconst ,k =

xk ω0 + yk |V0, k | , mk xk + nk yk

yk ω0 + xk |V0, k | . nk xk − mk yk

(24)

(25)

var In the similar fashion, the second part of active (PGvar , k ) and reactive (QG, k ) power of equations (22) and (23) are variable which depends upon the system frequency (ω) and the magnitude of the bus voltage (|Vk|). PGvar ,k and QGvar , k can be calculated as,

Pl, k = PGvar ,k =

xk ω + yk |Vk | , mk xk + nk yk

(26) 5

Electric Power Systems Research 180 (2020) 106131

A. Kumar, et al.

Fig. 2. Flow chart of NICINR.

i+1 wgrid = w0 −

1 ∑k = 1

1 mk (xk PG, k − yk QG, k )

state solution of DBIMG when the value of del is close to 0. In this work, condition, del ≤ 10−8, is utilized to accept the solution of NICINR as a power flow solution.

. (31)

Above equation adjusts the system frequency according to the mismatch between scheduled active power and droop characteristics of AR bus. However, the adjusted system frequency disturbs the scheduled active power of other droop buses. Therefore, equivalent adjustment of system frequency is done with minimum disturbance of scheduled active power of other droop buses using equation. (31). 5. The magnitude of voltage of AR bus will be updated according to equation (32).

V si + 1 =

(V si ) + (|V0, s | − ns (ys PG, s + xs QG, s )) 2

4. Validation of Proposed Algorithm The performance of the proposed algorithm is compared with the results obtained from time-domain simulations and other Newton-based algorithms in terms of accuracy and execution time.

4.1. Comparison of proposed Algorithm with Time-Domain Simulation

(32)

To confirm the applicability and accuracy, the proposed power flow algorithm is implemented on a six-bus droop based islanded microgrid test system [16]. Network parameters of the test system has been directly taken from [16]. Results obtained from the proposed power flow algorithm (NICINR) are compared to the results obtained from the time domain simulation of the test system in PSCAD/EMTDC. Voltages magnitudes and angles obtained from NICINR is compared with the voltage magnitudes and angles obtained in time-domain simulation in PSCAD. It is found that maximum error in voltage magnitude is 0.0081% and in the angle is 0.26%. This validates the accuracy of the NICINR to obtain the power flow solution of DBIMG. However, PSCAD takes approximately 172s to achieve the steady-state condition, while the NICINR requires 0.04s to converge. In every loop, the updated value of voltage and frequency of AR bus approaches to the steady-state value and the convergence is fast as shown in Fig. 3. From the above analysis, it can be concluded that the convergence rate of the proposed algorithm is better than the time-domain simulators in terms of the accuracy of the power flow solutions.

Above equation adjusts the voltage magnitude of AR bus acoording to the mismatch between scheduled reactive power and droop characteristics of AR bus. The average of current voltage magnitude, V si , and target voltage magnitude, |V0,s| − ns(ysPG,s + xsQG,s), of AR bus is utilized as voltage magnitude of AR bus for the next iteration. 6. Upon updating the operating frequency of the system, and voltage magnitude of the AR bus, value of the variable del is calculated (using equation (33)) to check the termination criteria. i+1 i del = max{|wgrid − wgrid |, |V si + 1 − V si |}.

(33)

7. If the obtained value of del satisfies the pre-specified tolerance criteria, the iterative loop will be terminated, and the power flow solution obtained in the last iteration will be the final power flow solution. 8. While, if the value of del is more than tolerance, the next iteration will be performed with updated values, and this procedure will continue until the convergence criteria are satisfied. The power flow solution of NICINR is only acceptable as a steady6

Electric Power Systems Research 180 (2020) 106131

A. Kumar, et al.

Table 1 Computation time required to solve power flow for different cases considering NICINR, MNR and NTR algorithm. (NC: Not Converged, NA: Not Applicable, CT: Computation Time, %: Percentage improvement in computation time.) System

case22

case38

Fig. 3. Convergence curve of system frequency and voltage of AR bus for scenario-1 with constant power load.

4.2. Comparison of Proposed Algorithm with MNR and NTR case69

In this section, four radial (case6, case22, case38, and case69) and one meshed (case160) distribution systems are considered to compare the results obtained from NICINR with the results obtained from MNR and NTR. Four different cases including: 1) conventional droop, 2) inverse droop, 3) mixed droop, and 4) isochronous mode of operation are studied to demonstrate the robustness of the proposed algorithm for the different operating modes of droop-based DGs. All the algorithms are executed on MATLAB R2017b in PC with an INTEL Core i7 @ 3.2 GHz, 8 GB of RAM. The flat start is considered as an initial solution in all algorithms (for operating frequency the flat start is 1.0 p.u.). The stopping criteria for all the algorithms are same and selected in a way that either the gradient norm or the total number of iterations does not exceed the specified value. It is to be noted that the NTR algorithm is implemented only for case 1 and MNR is not implemented for case 4 (isochronous operation). Therefore, the power flow results for these cases are not available. For case 1 (conventional droop), all the algorithms are applicable and the power flow results obtained from the algorithms are similar except MNR which fails to provide a solution in case of case69 and case160. Similarly, in cases 2, 3, and 4, NICINR and MNR provide similar power flow results for case22, and case38 but MNR fails to provide a solution for case69 and case160. For case 5, only NICINR can provide the power flow solutions because MNR and NTR are not defined for an isochronous mode of operation. To analyze the performance of MNR, NTR, and NICINR in terms of computation time, the time required to obtain the power flow solution is recorded for all the cases and reported in Table 1. It is found that expected execution time for the cases case22 and case38, NICINR is almost 10 times faster than MNR. Similarly, NICINR is approximately 3, 2, and 9 times faster than NTR for the cases case22, case38, and case69, respectively. From the above analysis, it can be concluded that the solutions obtained by NICINR converges faster in comparison to MNR and NTR without compromising the accuracy. Additionally, NICINR can also be used to perform power flow analysis for the systems having DGs with several types of droop characteristics including the isochronous mode of operation.

case160

Cases

NICINR

MNR

NTR

CT(s)

CT(s)

(%)

CT(s)

(%)

1

4.96E-03

6.35E-02

1180.16

1.77E-02

256.63

2

5.58E-03

4.84E-02

767.82

NA

NA

3

4.22E-03

4.98E-02

1080.50

NA

NA

4 1

3.81E-03 2.22E-02

NA 1.46E-01

NA 556.68

NA 3.48E-02

NA 56.83

2

1.42E-02

1.59E-01

1020.27

NA

NA

3

2.28E-02

1.64E-01

621.31

NA

NA

4 1

2.36E-02 1.33E-02

NA NC

NA NC

NA 1.43E-01

NA 970.79

2

1.72E-02

NC

NC

NA

NA

3

1.82E-02

NC

NC

NA

NA

4 1

1.81E-02 6.08E-02

NA NC

NA NC

NA 5.10E-01

NA 738.64

2

4.99E-02

NC

NC

NA

NA

3

4.47E-02

NC

NC

NA

NA

4

6.03E-02

NA

NA

NA

NA

Fig. 4. R/X ratio vs computation time (sec) for case38 test systems.

adapted from [13]. The performance of NICINR, MNR, and NTR evaluated in terms of computation time for different values of R/X ratios of line and the system loading condition are shown in Figs. (4) and (5) respectively. In the Figs. (4) and (5), it can be seen that the computation time is less for the NICINR in comparison to MNR. It is to be noted that MNR is not designed for the voltage dependent loads.

4.3. Robustness of Proposed Algorithm The performance of power flow algorithms depends on the R/X ratio of lines and the loading condition of the system. Hence, investigation of robustness of the proposed algorithm under various loading condition with different R/X ratios is performed. The radial network case38 is considered for the study which has the zone wise voltage dependent loads in addition to the lines having a high R/X ratio. The line/load data and the values of droop parameters of DGs for this system are

4.4. Application of proposed approach in NTR As discussed in the previous subsection, the robustness of NTR and NICINR is similar but the NTR algorithm cannot be applied to the systems having (i) DGs with mixed droop characteristics, (ii) DG with 7

Electric Power Systems Research 180 (2020) 106131

A. Kumar, et al.

Table 2 Computation time required to solve power flow for different cases considering NICINR, NTR-pDGm algorithm. (NC: Not Converged, CT: Computation Time, %: Percentage improvement in computation time.) System

Cases

case22

case38

case69

Fig. 5. Loading factor vs computation time (sec) for case38 test systems. case160

isochronous mode of operation and (iii) coupling between the lines. In this work, these issues of NTR are addressed and a new DG model and procedure for updating operating frequency are proposed. These procedures are in general structure and can also be employed with NTR to improve the versatility of the algorithm. In this section, to observe the versatility of the NTR, the modified algorithm is implemented on different types of DBIMGs.

NICINR

NTR-pDGm

CT(s)

CT(s)

(%)

2

5.58E-03

1.79E-02

220.84

3

4.22E-03

1.97E-02

365.86

4 2

3.81E-03 1.42E-02

1.73E-02 3.50E-02

353.31 147.17

3

2.28E-02

3.53E-02

54.77

4 2

2.36E-02 1.72E-02

4.41E-02 1.40E-01

87.21 716.45

3

1.82E-02

1.38E-01

658.41

4 2

1.81E-02 4.99E-02

1.50E-01 6.39E-01

729.12 1178.58

3

4.47E-02

5.27E-01

1078.50

4

6.03E-02

5.20E-01

763.13

Table 3 Computational effort of NTR-pr, NTR-fd, and NICINR for solving different test cases. (iter: number of iteration, CT: Computation Time) Test Cases

4.4.1. Significance of proposed DG model in NTR The proposed DG model is implemented in NTR by modifying the non-linear equations used in [13]. These non-linear equations are modified in line with the equations proposed in this work. The modified NTR, named NTR-pDMm, can perform the power flow analysis for different cases (discussed in section 4.2). The accuracy of the obtained solution of NTR-pDMm is similar to the obtained solution of NICINR for all the cases. However, computation time is still the issue in the NTRpDMm due to the high complexity of the steps of the conventional Newton-Trust algorithm. The results obtained using the proposed DG model in NTR (NTR-pDGm) and NICINR in terms of execution time are shown in Table 2. It is found that the computation time to converge on a steady-state value in NICINR is less than the NTR-pDGm due to the low complexity of the steps required in NICINR. It can be concluded from the above result and discussion that the proposed DG model can be implemented on other Newton-based algorithms to improve their versatility. Also, the convergence rate of NICINR is better than modified NTR (NTR-pDGm) due to its simple steps.

Case22 Case33 Case69 Case160 Case1458 Case3139

NTR-mod

NTR-fd

NICINR

iter

CT(s)

iter

CT(s)

iter

CT(s)

37 34 39 32 37 35

1.98E-02 2.97E-02 1.43E-01 5.48E-01 3.27E-00 7.14E-01

4 4 5 4 6 6

2.93E-01 5.68E-01 8.70E-01 2.81E-00 7.68E-00 1.67E+01

27 24 28 25 27 29

4.93E-03 6.30E-03 1.33E-02 6.08E-02 3.85E-01 8.26E-01

5. Conclusion In this paper, a nested-iterative approach, NICINR, is proposed to obtain the power flow solution of the droop-based islanded microgrid. A loop-based approach is employed to update the system frequency and voltage of the angle reference bus after every iteration. To analyze the effectiveness of NICINR, several case studies are carried out. In each case study, the load dependency and droop characteristics of DGs are considered. In each case study, outcomes obtained from NICINR show superior performance in terms of computational time and accuracy in comparison to MNR, NTR, and PSCAD. The main contributions of proposed approach are as follows.

4.4.2. Significance of proposed frequency update procedure in NTR The approach taken in the present work is further tested for its superiority by applying the modifications in NTR (NTR-mod). For the systems having coupling between the lines, the derivative of the admittance matrix cannot be calculated analytically since there is no explicit expression for the same. In this study the derivatives of the elements of admittance matrix is calculated using the finite differences method. It is assumed that the frequency is kept constant and the terms associated with the variation in frequency is eliminated. The results obtained using (i) the modifications in NTR (NTR-mod) and (ii) the elements of Jacobian are calculated using finite difference method (NTR-fd) are shown in Table (3). Significant reduction in computation time is achieved with the proposed modifications approach in comparison to the finite difference version NTR-fd as shown in Table (3). It shows that even if the NTR is adapted for the proposed framework, the NICINR shows better performance in comparison to NTR-mod.

1. The proposed algorithm, NICINR, deals with the issues and limitations related to the Newton-based algorithms, including MNR and NTR. In the Jacobian matrix, the need for gradient (with respect to system frequency) of the bus admittance matrix can also be dealt with the proposed NICINR. 2. In the proposed algorithm, the system frequency is updated in every loop according to the droop law. This update does not need the calculation of either gradient of the bus admittance matrix or the system frequency dependent variables. 3. To evaluate the system frequency and voltage magnitude of the AR bus, a closed-loop formulation is proposed, which results in fast convergence. 4. The proposed algorithm also converges for the isochronous mode of operation of DBIMG systems. 8

Electric Power Systems Research 180 (2020) 106131

A. Kumar, et al.

5. In the proposed approach, DG model and procedure for updating the system frequency are utilized in NTR to improve its versatility.

[19]

References [20] [1] M.Z. Kamh, R. Iravani, Unbalanced model and power-flow analysis of microgrids and active distribution systems, IEEE Transactions on Power Delivery 25 (4) (2010) 2851–2858. [2] D.K. Dheer, N. Soni, S. Doolla, Small signal stability in microgrids with high penetration of power electronics interfaced sources, IECON 2014-40th Annual Conference of the IEEE Industrial Electronics Society, IEEE, 2014, pp. 2272–2278. [3] D. G. Photovoltaics, E. Storage, Ieee guide for design, operation, and integration of distributed resource island systems with electric power systems. [4] N. Pogaku, M. Prodanovic, T.C. Green, Modeling, analysis and testing of autonomous operation of an inverter-based microgrid, IEEE Transactions on power electronics 22 (2) (2007) 613–625. [5] D. Dheer, S. Doolla, S. Bandyopadhyay, J.M. Guerrero, Effect of placement of droop based generators in distribution network on small signal stability margin and network loss, International Journal of Electrical Power & Energy Systems 88 (2017) 108–118. [6] A. Deihimi, B.K. Zahed, R. Iravani, An interactive operation management of a micro-grid with multiple distributed generations using multi-objective uniform water cycle algorithm, Energy 106 (2016) 482–509. [7] A. Esmaeli, M. Abedini, M.H. Moradi, A novel power flow analysis in an islanded renewable microgrid, Renewable Energy 96 (2016) 914–927. [8] D. Westermann, M. Kratz, A real-time development platform for the next generation of power system control functions, IEEE Transactions on Industrial Electronics 57 (4) (2010) 1159–1166, https://doi.org/10.1109/TIE.2009.2038403. [9] R. Majumder, G. Ledwich, A. Ghosh, S. Chakrabarti, F. Zare, Droop control of converter-interfaced microsources in rural distributed generation, IEEE Transactions on Power Delivery 25 (4) (2010) 2768–2778, https://doi.org/10. 1109/TPWRD.2010.2042974. [10] J.A. Mueller, J.W. Kimball, An efficient method of determining operating points of droop-controlled microgrids, IEEE Transactions on Energy Conversion 32 (4) (2017) 1432–1446. [11] H. Nikkhajoei, R. Iravani, Steady-state model and power flow analysis of electronically-coupled distributed resource units, IEEE Transactions on Power Delivery 22 (1) (2007) 721–728. [12] M.Z. Kamh, R. Iravani, Unbalanced model and power-flow analysis of microgrids and active distribution systems, IEEE Transactions on Power Delivery 25 (4) (2010) 2851–2858. [13] M.M.A. Abdelaziz, H.E. Farag, E.F. El-Saadany, Y.A.-R.I. Mohamed, A novel and generalized three-phase power flow algorithm for islanded microgrids using a newton trust region method, IEEE Transactions on Power Systems 28 (1) (2013) 190–201. [14] G.C. Kryonidis, E.O. Kontis, A.I. Chrysochos, K.O. Oureilidis, C.S. Demoulias, G.K. Papagiannis, Power flow of islanded ac microgrids: Revisited, IEEE Transactions on Smart Grid 9 (4) (2018) 3903–3905, https://doi.org/10.1109/TSG. 2018.2799480. [15] A. Elrayyah, Y. Sozer, M.E. Elbuluk, A novel load-flow analysis for stable and optimized microgrid operation, IEEE Transactions on Power Delivery 29 (4) (2014) 1709–1717. [16] F. Mumtaz, M. Syed, M. Al Hosani, H. Zeineldin, A novel approach to solve power flow for islanded microgrids using modified newton raphson with droop control of dg, IEEE Transactions on Sustainable Energy 7 (2) (2016) 493–503. [17] Hameed, Faisal, Al Hosani, Mohamed, Zeineldin, H. Hatem, A modified backward/ forward sweep load flow method for islanded radial microgrids, IEEE Transactions on Smart Grid 10 (1) (2017) 910–918. [18] M. Abedini, A novel algorithm for load flow analysis in island microgrids using an

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35] [36]

[37]

9

improved evolutionary algorithm, International Transactions on Electrical Energy Systems 26 (12) (2016) 2727–2743. A. Elrayyah, Y. Sozer, M.E. Elbuluk, A novel load-flow analysis for stable and optimized microgrid operation, IEEE Transactions on Power Delivery 29 (4) (2014) 1709–1717. A. Kumar, B.K. Jha, D.K. Dheer, D. Singh, R.K. Misra, Nested backward/forward sweep algorithm for power flow analysis of droop regulated islanded microgrids, IET Generation, Transmission Distribution 13 (14) (2019) 3086–3095, https://doi. org/10.1049/iet-gtd.2019.0388. G. Díaz, J. Gómez-Aleixandre, J. Coto, Direct backward/forward sweep algorithm for solving load power flows in ac droop-regulated microgrids, IEEE Transactions on Smart Grid 7 (5) (2016) 2208–2217. O. Anaya-Lara, E. Acha, Modeling and analysis of custom power systems by pscad/ emtdc, IEEE Transactions on Power Delivery 17 (1) (2002) 266–272, https://doi. org/10.1109/61.974217. M. Hossain, H.R. Pota, M.A. Mahmud, M. Aldeen, Robust control for power sharing in microgrids with low-inertia wind and pv generators, IEEE Transactions on Sustainable Energy 6 (3) (2014) 1067–1077. W. Yao, M. Chen, J. Matas, J.M. Guerrero, Z.-M. Qian, Design and analysis of the droop control method for parallel inverters considering the impact of the complex impedance on the power sharing, IEEE Transactions on Industrial Electronics 58 (2) (2011) 576–588. K. De Brabandere, B. Bolsens, J. Van den Keybus, A. Woyte, J. Driesen, R. Belmans, A voltage and frequency droop control method for parallel inverters, IEEE Transactions on power electronics 22 (4) (2007) 1107–1115. Y.A.-R.I. Mohamed, E.F. El-Saadany, Adaptive decentralized droop controller to preserve power sharing stability of paralleled inverters in distributed generation microgrids, IEEE Transactions on Power Electronics 23 (6) (2008) 2806–2816. D.K. Dheer, N. Soni, S. Doolla, Improvement of small signal stability margin and transient response in inverter-dominated microgrids, Sustainable Energy, Grids and Networks 5 (2016) 135–147. J.M. Guerrero, J. Matas, L.G. de Vicuna, M. Castilla, J. Miret, Decentralized control for parallel operation of distributed generation inverters using resistive output impedance, IEEE Transactions on industrial electronics 54 (2) (2007) 994–1004. H. Bevrani, S. Shokoohi, An intelligent droop control for simultaneous voltage and frequency regulation in islanded microgrids, IEEE transactions on smart grid 4 (3) (2013) 1505–1513. A. Kumar, B. JHA, R. Misra, D. Singh, A novel current injection based newtonraphson power flow algorithm for droop-based islanded microgrids, IET Generation, Transmission & Distribution. A.M. Haidar, K.M. Muttaqi, Behavioral characterization of electric vehicle charging loads in a distribution power grid through modeling of battery chargers, IEEE Transactions on Industry Applications 52 (1) (2015) 483–492. H. Nikkhajoei, R. Iravani, Steady-state model and power flow analysis of electronically-coupled distributed resource units, IEEE Transactions on Power Delivery 22 (1) (2007) 721–728. M.Z. Kamh, R. Iravani, Unbalanced model and power-flow analysis of microgrids and active distribution systems, IEEE Transactions on Power Delivery 25 (4) (2010) 2851–2858. G. C. Kryonidis, E.O. Kontis, A.I. Chrysochos, K.O. Oureilidis, C.S. Demoulias, G.K. Papagiannis, Power flow of islanded ac microgrids: Revisited, IEEE Transactions on Smart Grid. J.E. Dennis Jr., R.B. Schnabel, Numerical methods for unconstrained optimization and nonlinear equations, Vol. 16, Siam, 1996. P.A. Garcia, J.L.R. Pereira, S. Carneiro, V.M. da Costa, N. Martins, Three-phase power flow calculations using the current injection method, IEEE Transactions on Power Systems 15 (2) (2000) 508–514. V.M. da Costa, N. Martins, J.L.R. Pereira, Developments in the newton raphson power flow formulation based on current injections, IEEE Transactions on power systems 14 (4) (1999) 1320–1326.