Eigenanalysis-based small signal stability of the system of coupled sustainable microgrids

Electrical Power and Energy Systems 91 (2017) 42–60

Contents lists available at ScienceDirect

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

Eigenanalysis-based small signal stability of the system of coupled sustainable microgrids Farhad Shahnia ⇑, Ali Arefi School of Engineering and Information Technology, Murdoch University, Australia

a r t i c l e

i n f o

Article history: Received 3 July 2016 Received in revised form 14 February 2017 Accepted 6 March 2017

Keywords: Coupled microgrids Distributed energy resources Battery energy storage Eigenanalysis-based stability

a b s t r a c t Until now, microgrids were assumed to operate either in grid-connected or islanded modes. As a third alternative, neighboring microgrids may interconnect to support each other during faults in a section of one MG, or in the course of overloading. They may also interconnect to reduce the cost of electricity generation in their systems. Thereby, a microgrid will experience a significant transformation in its structure, when coupled with one or more neighboring microgrids. Before the formation of the system of coupled microgrids, the stability of the new system is vital to be cautiously examined to intercept the transformation, if instability is to occur. An eigenanalysis-based small signal stability evaluation technique is developed in this paper which can be used to allow/deny the interconnection of the microgrids. The analysis also defines the suitable range of control parameters for the energy resources of the CMG to guarantee the stability of the new system, if it was determined unstable. Numerical analyses, realized in MATLAB, are performed for the interconnection of a few sample microgrids. The accuracy of the developed technique is validated by comparing its results with the time-domain performance of the similar system, realized in PSCAD/EMTDC. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction A microgrid (MG) is a cluster of distributed energy resources (DERs), battery energy storages (BESes) and loads and can operate in grid-connected as well as islanded modes [1]. An MG can run in islanded mode temporarily due to any maintenance or emergency conditions in the utility while the MGs of remote and rural areas operate in islanded mode permanently (if the utility feeder is not available due to cost-benefit concerns) [2]. When an MG operates in islanded mode, the voltage, and frequency of the MG needs to be coordinated by one of the dispatchable DERs or by a group of them. The MG control is achieved under centralized, decentralized, or distributed schemes [3–5]. Fuel cost and its transportation difficulties increase the levelized cost of electricity for MG operators when they use diesel or liquid petroleum gas-based DERs [2]. If a renewable energybased source is found to be ample in an area, a larger portion of the electricity demand can be supplied from such resources. Thereby, the need for diesel/gas-based generators is expected to be minimized or removed; thus, all of the electricity demand can be supplied by the renewable energy resources. In such a case,

⇑ Corresponding author. E-mail address: [email protected] (F. Shahnia). http://dx.doi.org/10.1016/j.ijepes.2017.03.003 0142-0615/Ó 2017 Elsevier Ltd. All rights reserved.

energy storage systems may be required to store and return the energy in accordance with the intermittencies and availabilities of the resources [6]. Thus, the ultimate vision is to have a sustainable MG, which is composed of renewable energy-based, converter-interfaced DERs [2]. Very recently, the concept of coupled MGs (CMGs) has been proposed in which two or more neighboring MGs can interconnect temporarily [7]. Neighboring MGs can interconnect to support each other during emergencies (e.g. when a fault occurs in a section of an MG, leading to the outage of one or some of its DERs [8], or when an MG is overloaded [9]) or during normal conditions to minimize the levelized cost of electricity [10]. An MG may also interconnect to another MG to import power if the neighboring MG temporarily observes a high power generation by one of its DERs and thus offers electricity for trade with a lower price. Let us consider the network of Fig. 1 which shows a group of neighboring sustainable MGs (that are stretched over a few kilometers). Under the concept of coupling the MGs, an MG in Fig. 1 may be connected to one or more of its neighboring MGs. A transformative architecture is proposed for coupling the nearby MGs in [11] to improve the system resiliency during faults. A decision-making-based approach is proposed in [12] to determine the most suitable MG(s) to be coupled with an overloaded MG, which considers different criteria such as available surplus

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

43

Nomenclature Abbreviations BES Battery energy storage CMG Coupled microgrids DER Distributed energy resource IGBT Insulated gate bipolar transistor LPF Low pass filter LQR linear quadratic regulator MG Microgrid PI proportional-integral RoD Rate of discharge SoC State of charge SSS Small signal stability VSC Voltage source converter

q, Q S Tk u uc |Vc|

Parameters a turns ratio of transformer b bandwidth of hysteresis control Cf capacitance of LC filter f frequency m, m’ P-f droop coefficient n Q-V droop coefficient N number of DERs, BESes, loads, lines J objective function k state-feedback gain Lf inductance of LC filter terminal inductance LT NoC number of cycles R control cost variable Rf power losses of transformer and H-bridge

f

Variables E i1 ic i2 p, P

vc vt

d k

r t x xn

Matrices and vectors A, B, C, F, G, H system matrices identity matrix of order 2 I2
2 k matrix of state-feedback gains O2
2 zero matrix of order 2 PF participation factor matrix Q state weighting matrix u uc matrix in direct-quadrature frame v voltage matrix in direct-quadrature frame z state vector in time domain x, Z state vectors in direct-quadrature frame C right eigenvectors matrix Indices, subscripts, and superscripts a, b, c Phase-a, b, and c d, q components of direct-quadrature frame h discrete-time index i, j, k, l, kk, ll, l numbering Indices t continuous-time index operating point o

stored energy in BES current flowing though Lf current flowing through Cf current flowing through LT instantaneous and average active power

power, electricity cost, reliability and the distance of the neighboring MGs as well as the voltage/frequency deviation in the CMG. Ref. [9] presents the conditions based on which the overloading of an MG and the availability of excess power in the neighboring MG can be detected. The trade of power among MGs of a CMG is discussed in [10] while [13] presents an interactive control of CMGs to guarantee adequate load sharing and system-wide stability. The dynamic operation of DERs within CMGs is investigated in [14] whereas [15] examines the dynamic security of the CMGs. The interaction among the DERs of the MGs in a CMG is

MG-2 MG-3

MG-1 S BE ad Lo

Load

Lo ad

DER

DE R

BES

ER

BE S

D

MG-N

Load

DER

Large Remote Area/Town

BES

Load

MG-N _1

DER

instantaneous and average reactive power apparent power oscillation period of an eigenvalue switching function of IGBTs continuous-time approximation of u magnitude of vc voltage across Cf terminal voltage angle of voltage eigenvalue real part of an eigenvalue imaginary part of an eigenvalue angular frequency undamped natural frequency damping ratio

BES

Fig. 1. Illustration of a large remote town supplied by isolated MGs.

investigated in [16]. Ref. [17] analyses the reliability aspects of a CMG, while the voltage and current controllability in CMGs is analyzed in [18]. Coupling of MGs can be realized by back-to-back converters [19] or by static switches [9] between the adjacent MGs. The ultimate vision is that an MG can be interconnected to any MG (and not necessarily an adjacent MG) if a general link is available to act as a power exchange highway. Thus, [20] has developed an optimal planning of the interconnecting lines between the MGs. The above literature has presented and evaluated the planning and operational aspects of transforming MGs into a CMG as a beneficial strategy for their operators. Without a further discussion on the benefits and operational aspects of CMGs, this paper has developed a small signal stability (SSS) evaluation tool for CMGs, as determining the stability status of the new system is very important because at the CMG formation time, one of the MGs may be unstable or very close to its instability boundary [9]. The developed SSS tool can determine whether the new system will become stable and can define its stability margins. Thus, its outcome can be used to accept or deny the CMG formation. Alternatively, if instability is observed, a range of droop control coefficients should be determined for the DERs to enable the interconnection. Fig. 2 illustrates schematically the flowchart of forming a CMG in which the primary research focus of this paper (i.e., development of the necessary stability evaluation tool) is highlighted. The application of the developed tool to stabilize an unstable CMG and the relevant mechanism will be addressed in a future publication.

44

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

Fetch System Data Dynamically No

CMG Criteria Satisfied?

Research Scope

Yes Developed Eigenanalysis-based SSS Evaluation Technique

Modify the Droop No Coefficients for the DERs and/or the BES Systems

Will the CMG be Stable? Yes Proceed to Interconnect the MGs

Fig. 2. Operation flowchart for coupling MGs and the research scope of this paper.

The SSS of a multi-synchronous machine-based power system is extended to MGs, and a variety of approaches such as eigenanalysis [21], frequency-response (impedance-based) methods of Middlebrook criterion [22], gain-phase margin criterion [23], and a micro-based approach [24] are used. A review of different techniques such as transformation into a rotating direct-quadrature reference frame, dynamic phasors, reduced-order modeling, and harmonic linearization is presented in [25]. Ref. [26] has also carried out a review of different SSS analyses for an MG in gridconnected and islanded modes. In eigenanalysis, the system is represented by a set of statespace equations and its eigenvalues are calculated from the state matrix [21]. An eigenanalysis-based SSS model of an MG composed of converter-interfaced DERs is presented in [21,25–33] while [34– 36] present the model for an MG consisting of both inertial and converter-interfaced DERs. The internal dynamics of the DERs are also considered in the SSS analysis in [37,38]. A reduced order SSS analysis models for an MG is developed in [39] in which the internal dynamics of the VSCs of the DERs are neglected. Another reduced order technique is developed in [40] using the singular perturbation technique, in which the filters and power filters are not considered. Moreover, the impact of communication delays on the stability of an MG [41], in addition to the effect of constant power loads (such as rectifiers, dc-dc supplied loads, and inverterdriven loads) [42] and induction motor loads [43] as well as the ramp rates of the DERs [44] have been also evaluated for an MG. The above studies have focused on an MG with a known structure and have not considered any changes in the structure of the MG, which can happen due to the interconnection of the MG to one or more MGs. This is the research gap that is focused in this paper and is required to be addressed when analyzing the stability of CMGs. Thereby, the main novelty of this research is the development of a generalized eigenanalysis-based SSS evaluation tool for a group of MGs, which some of them may be interconnected temporarily. In summary, the main contributions of the paper are:  to develop an algorithm for eigenanalysis of a CMG, that is composed of multiple interconnected sustainable MGs, and  to validate the accuracy of the developed tool by comparing its results with the system’s time-domain performance.

be easily converted to a functional code for a microprocessor, using the existing features in MATLAB, with little knowledge on microprocessor architectures. 2. Network under consideration Consider the group of neighboring MGs of Fig. 1. Each MG operates independently from the other MGs. Different NDER converterinterfaced DERs, NBES BESes, Nload loads, Nline lines and Nbus buses are considered in each MG. In this research, a BES is not categorized under a DER; although both of them are converter-interfaced. The operation of DERs and BESes are discussed below: 2.1. Operation of DERs Different control methods for the DERs of an MG are reviewed in [45]. In this research, the voltage magnitude and frequency at the converter output of DER-i, i 2 {1, . . . , NDER} is thought to be regulated by frequency droop control as [14]

f DERi ¼ f

max

 mDERi PDERi

jV c jDERi ¼ V max  nDERi Q DERi

ð1Þ

where m is the P-f droop coefficient in Hz/W, n is the Q-V droop coefficient in V/Var, Vmax and fmax are respectively the highest allowed voltage magnitude and frequency, and P and Q are respectively the average active and reactive powers injected by the DER to its terminal. Alternatively, the generalized droop control can be used [46]. If the cost of electricity generation by the DERs is different, the MG owners may adopt a cost-prioritized droop technique [47]. When a DER increases its output active and reactive power from zero to maximum, the DER frequency and voltage magnitude reduces from fmax to fmin and from Vmax to Vmin where max and min respectively show the maximum and minimum allowable limits in the MG. Thus, m and n of DER-i, are derived from [14]

mDERi ¼ ðf

max

f

min

Þ=Pcap DERi

nDERi ¼ ðV max  V min Þ=2Q cap DERi

cap

ð2Þ max

where illustrates the capacities of the DERs. In this research, f , fmin Vmax, and Vmin of all MGs are considered the same. Thus, the DERs of all MGs, albeit their capacities, have the same Df = f max  f min and the same DV = Vmax  Vmin. The ratio of the active and reactive powers generated by any two DERs in an MG (e.g. DER-i and DER-j) is desired to be equal to the proportion of their nominal powers (capacities), and is controlled by the reciprocal of the ratio of their droop coefficients [14]; i.e., cap PDERj =PDERi  Pcap DERj =P DERi ¼ mDERi =mDERj cap Q DERj =Q DERi  Q cap DERj =Q DERi ¼ nDERi =nDERj

ð3Þ

In this research, a DER refers to a combination of energy resource with its smoothening energy storage; thus, the DERs are dispatchable. Also, the DERs are thought to be capable of generating the powers defined by (3). If due to environmental reasons, the output power of any of the DERs is saturated below the desired value of (3), the other DERs will share the load based on the ratios given by (3). 2.2. Operation of BESes

It is to be highlighted that although some existing commercial power system simulation software tools have some sort of stability analysis tools, the developed algorithm can be custom-developed in MATLAB. Thus, it has more flexibility for the code developers (in comparison with existing commercial software tools) when expanding it to include adaptive control or optimization techniques which are the later stages of interconnecting the MGs and operating CMGs. Also, the developed algorithm in MATLAB can

A different operation mechanism is used for the BESes, which guarantees that different BESes in an MG will operate such that the state of charge (SoC) of all BESes reduce to the minimum acceptable SoC (SoCmin) at the same time, regardless of their capacity and initial charge [48]. For this, a BES with lower SoC is controlled to inject less active power compared to a BES with a higher SoC, regardless of the nominal capacities or nominal rate

45

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

of discharges (RoD). To this end, the RoD of BES-i, (RoDBES-i), i 2 {1, . . . , NBES} is defined as

RoDBESi ¼ SoC dBESi Pcap BESi where

SoC dBESi

ð4Þ

8 1 SoC ¼ 1 > > > > < 0:9 0:9 6 SoC < 1 ¼ . .. > > > > : 0:2 0:2 6 SoC < 0:3

From (4), the P-f droop coefficient of BES-i (mBES) is defined as

mBESi ¼ Df =RoDBESi

ð5Þ

From (4) and (5), it is seen that the output active power of BES-i reduces as SoCd becomes smaller. Thereby, a lower portion of the capacity of the converter of the BES ðScap BES Þ is utilized for active power flow. The unused part of the converter capacity can then be used for reactive power exchange. Therefore, in this research, in addition to active power exchange between the BES and the MG, a reactive power exchange also occurs. The available capacity in the converter of BES-i ðQ av BES Þ is defined dynamically as [49]

Q av BESi ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðScap BESi Þ  ðRoDBESi Þ

ð6Þ

and is used to defined the V-Q droop coefficient of BES-i (nBES) as

nBESi ¼ DV=2Q av BESi

ð7Þ

The reactive power exchange between the converter of the BES and the MG does not affect the stored energy and the SoC of the BES. 3. Considered VSC structure and control for DERs/BESes The DERs in this research are connected to the MG through a voltage source converter (VSC), with the structure of Fig. 3. The considered VSC is composed of three parallel H-bridges, where each H-bridge is comprised of 4 switches. Each switch is thought to be an insulated gate bipolar transistor (IGBT) with an antiparallel diode. A single-phase boosting transformer, with a turns ratio of 1:a is used at the output of each H-bridge. These transformers are wye-connected at their secondary and provide voltage boosting capability as well as galvanic isolation. A Lf-Cf-LT filter is used at the secondary side of the transformers. A voltage control technique is utilized in abc frame to track the desired voltage at the output of the VSC (across capacitors Cf). This structure provides a path for the zero sequence component of the current, in the case of unbalanced loads. Alternatively, a three-phase, neutral-clamped

i2

v

+ c

ic i1

_

vtn

LT

_

_ zðtÞ ¼ Az zðtÞ þ Bz uc ðtÞ þ Cz v t ðtÞ

z½h þ 1 ¼ Fz z½h þ Gz uc ½h þ Hz v t ½h R Ts

Lf Rf

S2

S2

S1

1:a

S1

Fig. 3. Considered VSC-filter structure for each DER and BES.

ð9Þ R Ts

where Fz ¼ eA:T s , Gz ¼ 0 eAt B dt and Hz ¼ 0 eAt C dt, assuming h as the discretized sample index and Ts as the sampling time. Using a state feedback control, uc[h] can be computed as

uc ½h ¼  k ðz½h  zref ½hÞ

T

ð10Þ

where k = [k1 k2 k3]. Assuming a full control over uc[h], an infinite-time linear quadratic regulator (LQR) can be designed to define k. In general, this controller has a gain margin of infinity and a phase margin of 60° and is categorized under optimal controls [51,53]. Thus, it is more stable than proportionalintegral (PI)-based controls and prevents the instability among parallel DERs when two MGs interconnect. In addition, PIbased controllers are not very effective when utilized in abc frame. In discrete LQR, an objective function (J) is chosen as [53]

JðhÞ ¼

1 X T T ½ðz½h  zref ½hÞ Q ðz½h  zref ½hÞ þ uc ½h Ruc ½h

ð11Þ

h¼0

where 0 < R < 1 is the control cost variable (e.g. R = 0.1), Q is the state weighting matrix that reflects the importance of each controlling parameter in z an J(1) represents the objective function at infinite-time (steady-state). Eq. (11) is then minimized to obtain the optimal values for k through a solution of Riccati equations while satisfying constraints of (9). k is defined once only as it solely depends on the fixed parameters of Vdc, a, Rf, Lf, Cf, and LT. Eq. (10) shows the input of the control system of the VSC, which is composed of the output tracking errors. The tracking error can be minimized by limiting this error within a small bandwidth of b > 0. To this end, switching function u is generated using a hysteresis function, denoted by hys(.), from [51]

u ¼ hys ðuc ½h; bÞ ¼

Cf

ð8Þ

where uc is the continuous-time approximation of the switching function u and vt is the voltage of the DER terminal. Eq. (8) is represented in the discrete-time domain as [52]



+

Vdc

vtc

vtb

vta

VSC or a three-phase, four-leg VSC can be used [50]. Lf and Cf suppress the harmonic components of the current and voltage while LT controls the power injected to the terminal. Resistance Rf represents the power losses of the transformer and the H-bridge. In Fig. 3, the output voltage of each H-bridge is a square waveform with a peak of uaVdc where u is the switching function. Assuming a bipolar switching, u = +1 if switches S1 are on while u = 1 if switches S2 are on. More details on converter structure are given in [51]. Consider the state vector z(t) = [i2(t) ic(t) vc(t)]T for each phase of Fig. 3, where i2 is the current of LT while vc and ic are respectively the voltage and current of Cf and T is the transpose operator. Then, the single-phase circuit of Fig. 3 is represented in continuous-time domain with state-space description of

þ1 uc ½h > þb 1 uc ½h < b

ð12Þ

The switching frequency and therefore the VSC switching losses depend on b. Small b values result in better converter output current and voltage tracking while increasing the switching frequency and losses. Thus, b should be defined to achieve an acceptable tracking in the output waveforms while keeping the losses small. In this study, it is assumed that b = 104. Alternatively, a pulse-width modulation function can be used.

46

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

A similar closed-loop control is applied for the other two Hbridges. The considered closed-loop control block diagram for each phase of Fig. 3 is shown in Fig. 4. This figure illustrates the inner-loop (switching control) and outer-loop (droop control) of a DER. In this figure, the instantaneous active and reactive powers, denoted respectively by p and q are calculated from [51]

pðtÞ ¼ i2 ðtÞv ac ðtÞ þ i2 ðtÞv bc ðtÞ þ i2 ðtÞv cc ðtÞ 1 a b qðtÞ ¼ pffiffiffi ½i2 ðtÞðv bc ðtÞ  v cc ðtÞÞ þ i2 ðtÞðv cc ðtÞ  v ac ðtÞÞ 3 c þ i2 ðtÞðv ac ðtÞ  v bc ðtÞÞ a

b

d

as expanded in Appendix-A. In (15), xDER ¼ ½i1

q

i1

d

i2

q

i2

T

v dc v qc 

T

is different from ZDER ¼ ½i2 i2 ic ic v dc v qc  , i.e., the equivalent of z in the direct-quadrature frame for which the LQR gains of [k1 k2 k3] were calculated in (11). However, the differences between xDER and ZDER (i.e., the states of i1 and ic) are correlated as ic = i1  i2 according to Kirchhoff’s current law in Fig. 3. Thus, the state-feedback control of (10) can be rewritten in directquadrature domain as d

c

q

d

q

ref

DuiDER ¼ Ki ½DZDER  ðDZiDER Þ  ¼ Gi DxiDER þ Hi ðDZiDER Þ ð13Þ

ref

ð16Þ

in which G applies appropriately the LQR gains, that were initially calculated for zDER, to xDER. It is to be noted that in (16), it is also assumed that by a perfect tracking, uc can be represented by u. From (15) and (16), the state-space equation of the VSC-filter system of DER-i with its inner-loop control becomes

where a, b and c respectively represent phase-a, b, and c. The average active and reactive powers can then be calculated by passing their instantaneous values through a low pass filter (LPF) with a Laplace domain transfer function of xc/(s + xc) where xc is its cut-off frequency and s is the Laplace operator. In Fig. 4, the references of vc(t), ic(t), and i2(t) of phase-a for DER-i are calculated from [31]

Dx_ iDER ¼ AiDER DxiDER þ BiDER ðDZiDER Þ

ref

þ CiDER Dv iDER

ð17Þ

where AiDER ¼ Ai þ Bi Gi and BiDER ¼ Bi Hi . Since the BESes have the same VSC-filter structure and closedloop control, the state-space equations for BES-i system is represented as

pffiffiffi 2jV c jDERi cosð2p f DERi tÞ pffiffiffi refa ic ðtÞDERi ¼ 2 2p f DERi C f jV c jDERi cosð2p f DERi t þ 90 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi    2 2 2 PDERi þ Q DERi 1 Q DERi i2 ðtÞrefa ¼ cos 2 p f t  tan DERi DERi jV c jDERi PDERi 3

v c ðtÞrefa DERi ¼

Dx_ iBES ¼ AiBES DxiBES þ BiBES ðDZiBES Þ

þ CiBES Dv iBES

ref

ð18Þ

If a different structure and control are used for the VSCs of the DERs and/or BESes, the corresponding state-space description will be utilized instead of (17) and (18). The loads in this study are assumed as constant impedance RL loads. Thus, the state-space equation of load-j connected to bus-jj of the MG is represented as

ð14Þ while the references for phase-b and c have respectively 120° and +120° phase shifts to those of phase-a. More details on converter control are given in [51]. A similar structure and control is used for the converters of the BESes. The only dissimilarity is that droop coefficients given by (2) are utilized in (1) for DERs while the droop coefficients given by (5) and (7) are used in (1) for BESes.

j j j j Dx_ load ¼ Aload Dxload þ Cload Dv jjload

ð19Þ

as expanded in Appendix-A. If the loads are constant power or induction motor type, the technique of [42,43] can be utilized. Likewise, the lines in this study are thought to be constant impedance RL branches; thereby, the state-space equation of line-k between adjacent buses kk and ll is represented as

4. Small signal model and eigenanalysis This Section provides the detailed small signal modeling for a CMG. The equations are developed in a comprehensive fashion such that the program can quickly evaluate the stability of a CMG, composed of any two or more MGs. The state-space description of DER-i, with the VSC-filter structure of Fig. 3, is represented in direct-quadrature frame as

as expanded in Appendix-A. From (17)–(20), the equivalent state-space description of an MG becomes

Dx_ iDER ¼ Ai DxiDER þ Bi DuiDER þ Ci Dv iDER

Dx_ 1MG ¼ A1MG Dx1MG þ B1MG ðDZ1MG Þ

þ C1MG Dv MG

vcc vcb vca

a

i2 i2 i2

Droop Control (1)

f

LQR Gains

|Vc |

Reference Generation (14)

uc

k3

vref-a c

VSC and LCL Filter System (9)

a

k1

k2

a

ic

a

i2

c (10)

Outer-loop Control

u

i1

icref-a i2ref-a

Hysteresis Function (12)

+

P

+

p

Power Calculation q LPF Q (13) b

ref

a

vc vc vc

c

ð20Þ

+

b

ð15Þ

+

c

Dx_ kline ¼ Akline Dxkline þ Ckline Dv MG

Inner-loop Control

Fig. 4. Considered closed-loop state-feedback-based control for the VSCs.

phase - a

b

ð21Þ

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

where

x1MG

2 3T Loads Lines DERs BESes zfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{ 6 1 Nload Nline 7 NBES NDER 1 1 1 ¼ 4xDER . . . xDER xBES . . . xBES xload . . . xload xline . . . xline 5

Z1MG

2 3T DERs BESes zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{ 6 1 NDER NBES 7 1 ¼ 4ZDER . . . ZDER ZBES . . . ZBES 5 0

DERs

In addition, the angle of the voltage at the output of the VSC of each DER can be calculated from the angular frequency of the converter (x = 2pf) versus a common angular frequency (xcom) as [31–33]

Z

d¼

Dd_ ¼ Dx  Dxcom

1 Loads Lines zfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ Nload Nline C A1load ; . . . Aload ; A1line ; . . . Aline A

x ¼ xmax  m0 P

3 BESes DERs zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ 6diagðB1DER ; . . . BNDER ; B1BES ; . . . BNBES Þ7 DER BES 7 ¼6 O2ðN 4 5 load þNline Þ
6ðNDER þN BES Þ 3 1 6 zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ 7 6 7 6diag@C1 ; . . . CNDER ; C1 ; . . . CNBES ; C1 ; . . . CNload A7 6 7 DER BES load DER BES load 6 7 6 7 ¼6 7 1 6 7 C line 6 7 6 7 .. 6 7 . 4 5 Nline Cline

C1MG

BESes

" d 3 i2 ¼ 2 iq q 2

q

i2 d

i2

#"

v dc v qc

ð29Þ

_ ¼ m0 DP Dx

Loads

ð30Þ

Replacing (30) in (28) yields

Dd_ iDER ¼ AidDER DpiDER þ DidDER DPcom

ð31Þ

as expanded in

where O represents a matrix of zeros, while diag(x, y, . . .) creates a diagonal matrix with diagonal elements of respectively x, y, . . .. The instantaneous output active/reactive power of each DER is expressed in direct-quadrature frame as [31–33]

  p

ð28Þ

where m’ = 2pm in rad/s.W. Differentiating (29) and linearizing it yields

2

DERs

ð27Þ

On the other hand, from (1), it can be said that

2

0

ðx  xcom Þ dt

xcom = mcomPcom is usually selected as the angular frequency of one of the converters, e.g. xcom = x1 = m1P1. Differentiating (27) and linearizing it around the operating point yields

BESes

zfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ B DER BES A1MG ¼ diag@A1DER ; . . . ANDER ; A1BES ; . . . ANBES ;

B1MG

47

#

ð22Þ

Dd_ ¼ m0 DP þ m0com DPcom

ð32Þ

Again, (31) is also valid for BES-i as

Dd_ iBES ¼ AidBES DpiBES þ DidBES DPcom

ð33Þ

BES is characterized by its SoC, which is calculated as [49]

Z

Assuming an LPF with a cutoff frequency of xc, the linearized average active/reactive powers around the operating (equilibrium) point of DER-i is expressed as

SoC ¼ SoC int 

pdt=Ecap

ð34Þ

Dp_ iDER ¼ AiPDER DpiDER þ BiPDER DxiDER

ð23Þ

where p, Ecap, and SoCint are respectively the instantaneous output power, storage capacity in Whr and the initial SoC of the BES. Linearizing (34) around the operating point of BES-i yields

Dq_ iDER ¼ AiQDER DqiDER þ BiQDER DxiDER

ð24Þ

_ iBES ¼ BiSoCBES DxiBES DSoC

as expanded in Appendix-A. Equations and (24) are also valid for BES-i as

Dp_ iBES ¼ AiPBES DpiBES þ BiPBES DxiBES

ð25Þ

DqiBES ¼ AiQBES DqiBES þ BiQ BES DxiBES

ð26Þ

x2MG

Combining (23)–(26), (31), (33), and (35) yields 2 2 2 1 Dx_ 2MG ¼ A02 MG DxMG þ BMG DxMG þ DMG DP com

2 3T BESes BESes BESes BESes DERs DERs DERs zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 6 DER BES DER BES DER BES BES 7 ¼ 4p1DER . . . pNDER p1BES . . . pNBES q1DER . . . qNDER q1BES . . . qNBES d1DER . . . dNDER d1BES . . . dNBES SoC1BES . . . SoCNBES 5 1 DERs BESes DERs BESes DERs BESes zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflBESes ffl}|fflfflfflfflfflffl{C B 1 N DER 1 N BES 1 N DER 1 N BES 1 N DER 1 N BES ¼ diag@APDER ; . . . ; APDER ; APBES ; . . . ; APBES ; AQ DER ; . . . ; AQ DER ; AQBES ; . . . ; AQ BES ; AdDER ; . . . ; AdDER ; AdBES ; . . . ; AdBES ; ONBES
NBES A 0

B2MG

1 DERs BESes DERs BESes BESes DERs & BESes zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl }|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl { B C N BES N BES 1 1 DER BES ¼ diag@B1PDER ; . . . ; BNPDER ; B1PBES ; . . . ; BNPBES ; B1Q DER ; . . . ; BNQDER DER ; BQBES ; . . . ; BQ BES ; OðN DER þN BES Þ
6ðN DER þN BES Þ ; BSoCBES ; . . . ; BSoCBES A 2

D2MG

ð36Þ

where

0

A02 MG

ð35Þ

3T DERs BESes zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 6 7 N DER NBES 1 1 ¼ 4O2ðNDER þNBES Þ
1 DdDER . . . ; DdDER DdBES . . . ; DdBES ONBES
1 5

48

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

Since DP1 is selected as DPcom in the above discussions and it is the first element in

x2MG ,

thereby, vector

D2MG

will be added to the rele-

2 vant vector (column-1) of A02 MG to form AMG . Thus, (36) becomes

Dx_ 2MG ¼ A2MG Dx2MG þ B2MG Dx1MG

ð37Þ

Assuming a large virtual resistor of Rv (e.g. Rv = 10+4) between each bus in the MG and the ground, the voltage of bus-jj is expressed as

"

Dv

jj bus

¼D

Expressing the time-domain references of (14) in directquadrature frame and then converting them to the common ref

direct-quadrature frame and linearizing them, ðDZiDER Þ [33]

ðDZiDER Þ

ref

¼ MiDER D½ piDER

diDER 

qiDER

as expanded in Appendix-A. Likewise, ðDZiBES Þ system of BES-i, becomes

ðDziBES Þ

ref

¼ MiBES D½ piBES

qiBES

ref

for the VSC-filter

T diBES 

Mi2DER

Mi3DER 

MiBES

Mi2BES

Mi3BES 

¼

½ Mi1BES

¼ Rv
bus

ð40Þ

Combining (38) and (39), ðDZMG Þ Dx2MG as

l

D

" #j d i2 q

i2

þ

N BES X

BESj

l

D

" #j d i2

DER

q

j¼1

as expanded in Appendix-A. Calculating (44) for each bus of the MG, DvMG is represented as a function of DxMG as

Dv MG ¼ NMG DxMG

ð45Þ

N

T

bus NMG ¼ ½N1bus . . . Nbus  2 DERNDER DER1 zfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 6 jj DER1 Nbus ¼ 4O2
2 l O2
2 . . . O2
2 lDERNDER O2
2

BESNBES

zfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{
O2
2 lBES1 O2
2 . . . O2
2 lBESNBES O2
2

is expressed as a function of

ðDZMG Þref ¼ M Dx2MG

: j¼1

DERj

i2 BES " #j " #k 9 N N = d d load line X X i i þ lloadj D q þ llinek D q ð44Þ i load k¼1 i line ; j¼1

BES1

ref

8
where

ð39Þ

Let us reshape MiDER and MiBES of (38) and (39) as

MiDER ¼ ½ Mi1DER

#jj

becomes

ð38Þ

T

vd vq

3 loads lines zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 7
lload1 . . . lloadNload lline1 . . . llineNline O2
ð3NDER þ4NBES Þ 5

ð41Þ

where

2

Replacing (45) in (43) yields the homogeneous form of the MG as

BES Systems

DERs

zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 6 DER BES M ¼ 4diag ðM11DER ; :::; MN1DER Þ diag ðM11BES ; :::; MN1BES Þ

Dx_ MG ¼ AHMG DxMG

where AHMG ¼ AMG þ CMG NMG . From (46), the eigenvalues of the MG ðkMG Þ are defined as the roots of [53]

BES Systems

DERs

zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ DER BES
diag ðM12DER ; :::; MN2DER Þ diag ðM12BES ; :::; MN2BES Þ 3 BES Systems DERs zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 7 DER BES
diag ðM13DER ; :::; MN3DER Þ diag ðM13BES ; :::; MN3BES Þ5

det ðkMG I  AHMG Þ ¼ 0

Replacing (41) in (21) yields

Dx_ 1MG ¼ A1MG Dx1MG þ B1MG M Dx2MG þ C1MG Dv MG

ð42Þ

Combining (37) and (42) the linearized state-space description of an MG, consisting of all states, is

Dx_ MG ¼ AMG DxMG þ CMG Dv MG where

xMG ¼

" ½x1MG "

CMG ¼

T x2MG  ;

AMG ¼

C1MG Oð3NDER þ4NBES Þ
2Nbus

#

ð43Þ

A1MG

B1MG M

B2MG

A2MG

# ;

ACMG

0 1 MGs Tielines zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ ¼ diag @AMG1 ; . . . AMGNMG ; ATieline1 ; . . . ATielineNTieline A

ðN MG k Þ is given by

NMG ¼ 9NDER þ 10NBES þ 2Nload þ 2Nline k

where

3T MGs Tielines zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ ¼ 4xMG1 . . . xMGNMG xTieline1 . . . xTielineNTieline 5

v CMG

where I is the identity matrix with the same size of AHMG and det(.) is the determinant function. As can be seen from the expanded equations in Appendix-A, each DER has 9 states (i.e., 4 states of current, two states of voltage, two states of power and one angle state), each BES has 10 states (i.e., the same states mentioned for the DERs as well as the SoC state), each load and line have two current states, the number of the eigenvalues for an MG

Dx_ CMG ¼ ACMG DxCMG þ CCMG Dv CMG

:

3T MGs zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ ¼ 4v MG1 . . . v MGNMG 5

ð47Þ

ð48Þ

For a CMG consisting of NMG  2 MGs that are interconnected through NTie-line tie-lines and static switches, the linearized statespace equation is calculated as

2

xCMG

ð46Þ

3 0 1 MGs 6 zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 7 7 6 6diag @CMG1 ; CMG2 ; . . . CMGNMG A7 7 6 7 6 7 ¼6 7 6 C Tieline1 7 6 7 6 . 7 6 . . 5 4 CTielineNTieline 2

2

CCMG

ð49Þ

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

while ATie-line and CTie-line are calculated similar to (20). It is to be noted that in (49), the state variables of all coupled MGs need to be expressed in a single reference frame. To this end, one of the converters of the CMG can be selected as CMG’s reference for the angular frequency (similar to the procedure discussed earlier for one MG), and the time-domain references of the states of all converters in the CMG will be first converted to the common direct-quadrature frame of the CMG, and then linearized. This will be reflected in the relevant column of each AMG in ACMG. Similar to (45), DvCMG is expressed as a function of DxCMG as

Dv CMG ¼ NCMG DxCMG

ð50Þ

where NCMG ¼ diagðNMG1 ; NMG2 ; . . . ; NMGNMG Þ in which the corresponding arrays of the buses of CMG to which the tie-lines are connected are updated with lTie-line-k, defined similar to lline-k. Replacing (50) in (49) yields the homogeneous form of CMG as

Dx_ CMG ¼ AHCMG DxCMG

ð51Þ

where AHCMG ¼ ACMG þ CCMG NCMG . From (51), the eigenvalues of the CMG (kCMG) are defined as the roots of

det ðkCMG I  AHCMG Þ ¼ 0 where I has the same size of

ð52Þ AHCMG .

The number of the eigenvalues

Þ is given by for a CMG ðN CMG k MG NCMG ¼ NMG1 þ NMG2 þ ::: þ NMGN þ 2NTieline k k k k

ð53Þ

For a stable CMG, all eigenvalues must be on the left side of the imaginary axis of s-plane, i.e., Re(k) = r < 0 where k = r ± jt and Re (.) represents the real function [53,54]. The eigenvalues located closer to the imaginary axis, are the dominant eigenvalues of the system. The CMG becomes unstable if any of the dominant eigenvalues are relocated to the right side of the imaginary axis; i.e., r > 0 [31,32]. Thus, the stability of the CMG will be evaluated by observing the trajectory of its dominant modes on the s-plane. Each eigenvalue of the system is affected by a group of its states. The contribution of each of the states of xCMG on each of the eigenvalues is calculated from the participation factor matrix (PF) defined as [54,55] 1 T

PF ¼ ðC Þ  C

ð54Þ

49

where C is the right eigenvectors matrix of AHCMG and s indicates the Hadamard production. PF(x, y) represents the contribution of state x on eigenvalue y. The states that have the highest impact (i.e., the largest contribution) on an eigenvalue are referred to as the effective states of that eigenvalue in the rest of this paper. The normalized weightings of each state on each eigenvalue are calculated from

kPFðx; yÞk ¼ jPFðx; yÞj =

X

jPFðl; yÞj

ð55Þ

l

where l shows all states of eigenvalue y. The damping ratio (f), undamped natural frequency (xn), and period (Tk) of each mode as well as the number of cycles for a mode during which its amplitude is damped by 50% (NoC50%) are defined from [55,56]

qffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f ¼ 1= 1 þ ðt=rÞ2 xn ¼ r=f T k ¼ 2p=xn 1  f2 qffiffiffiffiffiffiffiffiffiffiffiffiffi NoC 50% ¼ lnð2Þ 1  f2 =2pf

ð56Þ

where ln(.) represents the natural logarithm function. To define the range of the droop coefficients in which the CMG becomes unstable, using a sensitivity analysis, droop coefficient m is typically varied between 0.1mo and 10mo around the operating point of mo. Likewise, droop coefficient n is ranged from 0.1no and 10no around the operating point of no. Let us denote the very first m and n in which a stable system becomes critically (marginally) stable as mcritical and ncritical respectively. If all DERs have the same nominal powers, mcritical and ncritical are valid for all DERs. However, if they have different nominal powers, their mcritical and ncritical have a ratio as given by (3). Thereby, two new stability margins of Dfcritical and DVcritical are used here, which are derived from (2) and are expressed as

Df critical % ¼ Pcap DERi mcriticali =f nominal
100 DV critical % ¼ 2Q cap DERi ncriticali =V nominal
100

ð57Þ

which can be calculated from any DERs and is fixed for CMG. Fig. 5 illustrates schematically the flowchart of the algorithm that assesses the stability of a CMG prior to its formation.

Read Fixed Input Data and Fetch Dynamic Input Data of the DERs, BES Systems, Loads, Lines of the Coupling MGs Calculate the Linearized State-space Description for each DER, BES system, Load and Line of the Coupling MGs (18), (19), (20) and (22) Calculate the System Matrix for each Coupling MG (24) Run Dynamic Analysis and Calculate Operating Points for all States of each MG Calculate the Generalized State-space Description for each MG (48) Calculate the Generalized State-space Description for the CMG System (54)

Define the CMG System in the Homogenous Form (56) Define Eigenvalues of the CMG System (57) and fcritical , Vcritical (62) Fig. 5. Flowchart of the eigenanalysis-based SSS evaluation technique for CMGs.

50

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

Fig. 7), are to be evaluated by the developed tool, realized in MATLAB. Later, a comparison is carried out between the results of the study with system time-domain performance, realized in PSCAD/EMTDC. Although the internal parameters of multiple DERs, BESes, lines and loads in the MGs may be different in reality, for simplicity, they are considered the same in this numerical analysis. These analyses are to illustrate the performance of the evaluation technique and are used to verify its accuracy. Based on the developed

5. Numerical analysis results The developed eigenanalysis-based SSS evaluation tool is utilized to evaluate the stability margins of some sample CMGs. Let us consider MG-1, MG-2, and MG-3 of Fig. 1, with the detailed structures shown in Fig. 6 and the technical data provided in Table B1 in Appendix-B. Let us assume it is desired to couple MG-1 with either of MG-2 or MG-3, or both of them. Thus, the stability of these CMGs, denoted by CMG-A, CMG-B, and CMG-C (see

DER-1 Bus-2

Bus-1

Bus-5 Bus-4

Load-1 Bus-6

MG-1

DER-2

Bus-4

Bus-3

BES-1 DER-3

MG-2

Load-2 Bus-1

Bus-7 Bus-3 Load

Bus-2

DER-1

Bus-9

Bus-7

Bus-8

DER-2

Bus-6 Bus-5

Load-3

Load-1

Load-2

Bus-10

BES-2

MG-3

BES-1

Bus-2

Bus-1 DER-1

Bus-3

DER-2

Bus-4

DER-3

DER-4

Fig. 6. Considered structures for MG-1, MG-2 and MG-3.

CMG-A

MG2

MG1 MGN

MG2

MG3

MG_ N 1

Large Remote Area/ Town

MG1

CMG-B MGN

MG_ N 1

MG3 Large Remote Area/ Town

CMG-C

MG2

MG1 MGN

MG_ N 1

MG3 Large Remote Area/ Town

Fig. 7. Considered CMG-A, CMG-B and CMG-C systems.

1

DER-2

48%

(i2)MG-1 DER-1

Fig. 8. Eigenvalue trajectory of MG-1 and the participation factors of its eigenvalues.

48%

51

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

technique, another future publication will be focused on investigating the impact of the differences between the internal parameters of DERs, BESes, loads and lines among the MGs on the stability of the CMGs. 5.1. Stability margins of the considered MGs and CMGs First, the stability of each MG of Fig. 6 is evaluated. Let us assume that the operating point droop coefficients for all systems are mo = 0.0333 Hz/kW and no = 2.3 V/kVAr. To evaluate the system stability, droop coefficient m is varied in the range of 0.0033 m 0.333. Fig. 8 illustrates the eigenvalue trajectory of MG-1 in this scenario. This figure only shows the dominant eigenvalues of the system. Asterisks and respectively show the eigenvalues of the system at the operating point and the marginal

stability point while shows the direction of the increase of m. This trajectory is composed of one complex-conjugate pair of critical eigenvalues (k1,10 ) where the trajectory of k1 is shown in the zoomed sections of Fig. 8. The characteristics of this eigenvalue are provided in Table 1. A similar sensitivity analysis has been carried out for variation of droop coefficient n, which its eigenvalue trajectory is not presented here. The carried out sensitivity analyses yield that the stability margins of MG-1 are Dfcritical = 19.4468% and DVcritical = 7.7102%. Using these values in (2) gives mcritical and ncritical for every DER and BES in this MG. As similar DERs are considered in this study, the stable range of the droop coefficients of both DERs is determined as 0.1mo m < mcritical (i.e., 0.0033 m < 3.2411). Likewise, the stability region can be defined for the range of n as0.1no n < ncritical (i.e., 0.023 n < 3.5581). The donut chart of Fig. 8 also shows the

Table 1 Characteristics of the critical eigenvalues of MG-1 for m and n variation.

k1,10

r [1/s]

t [rad/s]

f [%]

xn [rad/s]

Tk [ms]

NoC50%

Critical for

8.408

±309.839

2.71

309.953

20

4.06

m, n

Table 2 Results of the SSS evaluation of the considered MGs and CMGs. System

Dfcritical [%]

DVcritical [%]

Nk

N kcritical [complex-conjugate pair(s)]

MG-1 MG-2 MG-3 CMG-A CMG-B CMG-C

19.4468 17.3083 12.3984 17.2865 12.3548 12.3548

7.7102 6.5085 6.2142 6.5085 6.2142 6.2142

30 57 84 89 116 175

1 4 7 5 8 15

2

1

BES-1

MG-2

(i2)

(i2)

48% DER-1

59% DER-3

31%

DER-2

BES-1

MG-2

DER-2

DER-1 19%

10%

8%

19%

16%

4

3 10%

20%

line-1

line-5

line-4

MG-2

47%

(i2) BES-1

DER-3 40%

(i)

MG-2

20% line-7

line-2 20% line-6

20% Fig. 9. Eigenvalue trajectory of MG-2 and the participation factors of its eigenvalues.

52

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

2

1

4

3

10% BES-1

13%

DER-2

38% DER-1

(i2)MG-3

DER-2 15%

12%

DER-3

DER-2

(i2)MG-3 BES-2

DER-4

28%

DER-4

DER-1

DER-1

22%

20%

36%

(i2)MG-3

DER-2

18%

(i)MG-3

10%line-10 line-9

10%

23%

line-3

line-1

33%

17%

4%

20%

line-2

10%

10%

line-6

line-8 line-7

10%

line-5

10%

10%

17%

(i2)MG-3

42%

BES-2

7

10%

10%

DER-4

53%

6

BES-2

BES-1

(i2)MG-3

17%

5

DER-3

BES-2 17%

37%

BES-1

39%

MG-3 DER-4

MG-3 DER-1 MG-3 BES-1

5% 6% MG-3 8% BES-1

MG-3 10% DER-4

MG-3 DER-2 MG-3 BES-2 MG-3 DER-2

MG-3 BES-2

10%

15%

Fig. 10. Eigenvalue trajectory of MG-3 and participation factors of its eigenvalues.

participation factor of different states of the system on its critical eigenvalues (k1,10 ). From this figure, it is seen that k1 is significantly affected by the output current states of DERs, i.e., (i2)DER-1 and (i2)DER-2, while the other states of the system have a total impact of 4%. The effective states of k1 and their weightings are the same with those for its complex conjugate eigenvalue. The above analysis is repeated for MG-2, MG-3 as well as the considered CMGs. Table 2 summarizes the outcomes of these studies while the eigenvalue trajectory and participation factor of the dominant eigenvalues of each system are shown in Figs. 9–13. The characteristics of the dominant eigenvalues of each of these systems are also listed in Tables 3–7. These analyses result in several findings as summarized below:  The eigenvalue trajectories of Figs. 8–13 illustrate that each mode of a system is composed of either one eigenvalue (such as the one in Figs. 8) or multiple eigenvalues (such as those shown by different colors in Figs. 9–13). Interestingly, it can be seen that the closer operating point eigenvalue to the imaginary axis of the s-plane is not necessarily the first critical eigenvalue that makes the system marginally stable. An example can be seen in the subset of the eigenvalue trajectory of Fig. 9 in

which the first critical eigenvalue (k4) is outside the zoomed section and farther than some other eigenvalues (k1, k2, and k3) at the operating point of the system.  From Table 2, it can be seen that the number of the complexconjugate critical eigenvalues of each CMG is equal or more than the sum of the critical complex-conjugate eigenvalues of its participating MGs. As an example, CMG-A has 5 complex-conjugate critical eigenvalues (i.e., equal to the sum of the critical complex-conjugate eigenvalues of MG-1 and MG-2) while CMG-C has 15 critical complex-conjugate eigenvalues (larger than the sum of the critical complex-conjugate eigenvalues of MG-1, MG-2, and MG-3). Also, comparing the characteristics of the dominant eigenvalues of the CMGs with the ones of participating MGs from Tables 1 and 3–7, it can be seen that the dominant eigenvalues of each CMG is a combination of the dominant eigenvalues of every participating MG in that CMG with very minor differences in the real and imaginary values (less than 1% displacement in the studies) which results in observing more or less the same f, xn, T and NoC50% in the dominant eigenvalues of the CMG when compared with those of the dominant eigenvalues of the participating MGs.

53

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

2

1

8%

7%

9%

(i2)MG-2 11% BES-1 (i2)MG-2 MG-1 DER-1 (i2)

DER-2

29%

BES-1 16%

(i2) DER-1

(i2)MG-2 DER-2

DER-3

26%

DER-2

(i2)MG-1

16%

(i2)MG-2

DER-1

(i2)MG-1 (i )MG-1 DER-1 2 DER-2 18%

DER-1 25%

27%

28%

5 8%

20%

(i)MG-1,2 tie-line

DER-2

(i2)MG-2 BES-1

BES-1

(i2)MG-2

4

48%

(i2)MG-2

29%

DER-2

MG-2

49%

3

DER-3 40%

20%

(i)MG-2 line-6 (i)MG-2 line-8

MG-2 (i) line-10

20%

20%

(i)MG-2 line-9 20%

Fig. 11. Eigenvalue trajectory of CMG-A and participation factors of its eigenvalues.

 Through the eigenvalue trajectories of Figs. 8–13, it can be seen that the dominant eigenvalues of each CMG at the operating point is approximately in the same location within the s-plane location of the dominant eigenvalues of the participating MGs in that CMG at their operating point. However, extra analyses yield that if the output power of the DERs and power flow through the lines in an MG change significantly after its coupling with other MGs, the location of the dominant eigenvalues of the CMG at the operating point might be more different compared to the location of the dominant eigenvalues of the participating MGs, prior to the coupling.  The donut charts of Figs. 8–10 show that each dominant eigenvalue of an MG is affected mostly by either the output currents of the DERs and BESes, or the line currents, or the output instantaneous active power and voltage angles. Comparing the donut charts of Figs. 11–13 with those of Figs. 8–10, it can be seen that each dominant mode of a CMG can be affected by the current, voltage, power, or angle states of one or more MGs. Through this analysis, no further systematic method was observed of their correlations. The studies also yield that the internal current and voltage states of DERs and BESes, as well as their out-

put reactive powers and the energy states of the BESes, do not have a strong effect on the dominant eigenvalues of the MGs or CMGs. These states mostly affect the non-dominant eigenvalues of the system (which do not fall within the demonstrated zone of the s-planes of Figs. 8–13 and damp very quickly).  By comparing the stability margins of Dfcritical and DVcritical between each CMG and its participating MGs from Table 2, it can be seen that these two quantities are almost equal to the respected minimum value of its participating MGs. Extended analyses yield that in general, if one MG is much stronger than the other one, the stability margins of the resulting CMG can settle somewhere between the stability margins of the participating MGs but closer to the stability margins of the less stable MG. As the primary focus of this paper was to explain the developed stability evaluation tool, the results of the studied cases are not discussed further. A future publication will concentrate on the impact of the parameters of the MGs on the stability of a CMG composed of those MGs based on which a mechanism is developed to make an unstable CMG stable.

54

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

2

1

3

4

8% BES-2

33% DER-2

12%

13%

BES-1

BES-1

BES-1

DER-1

38% MG-3

MG-3

(i2)

MG-3

(i2)

52%

43%

DER-4 21%

DER-3

16%

BES-2

(i2)MG-1

51%

BES-2

45%

DER-2

DER-1

DER-2

25%

30%

6

5

35%

DER-2

20% DER-2

BES-2

17%

32%

line-4

(i2)MG-3

DER-1

10%

10% 24% line-1

(i) MG-1

DER-4 DER-4

8

25%

BES-1

35%

(i2)MG-3

7

8%

10% BES-2

DER-3

(i2)

19%

line-2

25%

10%

line-1 line-5 line-2

(i) MG-3

10% line-9

line-3

20%

line-7 line-3

24%

10%

line-8

line-6

10%

21%

10%

line-4

10%

Fig. 12. Eigenvalue trajectory of CMG-B and participation factors of its eigenvalues.

5.2. Application of the developed analysis for decision-making during interconnection As discussed earlier, the main application of the developed SSS technique is to accept/deny an interconnection. To this end, two examples are illustrated here. Let us consider MG-1 and MG-3 introduced above in which they are intended to be coupled and form CMG-B. According to Table 2, Dfcritical of MG-1 and 3 are respectively 19.4468 and 12.3984% while this quantity for CMG-B is 12.3548%. Assuming Df as 10%, Fig. 14a illustrates the output active power of one of the DERs (i.e., DER-1) in MG-1 and MG-3, from which it is seen that these MGs are at their stable operating points. It is also seen that CMG-B is stable, if formed. Similar results can be observed for other DERs and BESes in each system. Thus in this example, the developed SSS analysis allows the interconnection as Df is found to be smaller than Dfcritical of CMG-B. As the second example, let us assume Df as 12.35%. Fig. 14b illustrates that although MG-1 and MG-3 are respectively at their stable and marginally stable operating points, CMG-B becomes unstable, if formed. To prevent an unstable CMG, the developed SSS analysis

will not allow the interconnection of these MGs as the operating conditions as Df is found to be larger than Dfcritical of CMG-B. This is because the droop control parameters of the DERs and BESes are beyond the limit corresponding to Dfcritical. A future publication will be dedicated to developing an adaptive control technique based on the outcomes of the developed SSS analysis, that can change the operating conditions of the DERs (by modifying their droop control parameters) to guarantee the CMG stability. 5.3. Verification of the developed SSS analysis To verify the accuracy of the developed eigenanalysis technique for the CMG, the findings of marginal stability of the CMG is compared with the results of the same system in the time domain when modeled in PSCAD/EMTDC. The CMGs are analyzed for different values of m in which they are stable (0.1mo m < mcritical), critically stable (m = mcritical), and unstable (mcritical < m 10mo). Since the CMGs are linearized around the operating points of mo, a small accuracy error is expected in case mcritical is closer to the upper range of m variations. For this example, a maximum error

55

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

2

1 7%

12%

DER-2

13%

BES-1

BES-1

DER-1

4

8%

BES-2

34%

3

BES-1 17%

(i2)MG-3

(i2)MG-2

48%

DER-4 22%

38%

BES-1

(i2)MG-3

52%

43%

(i2)MG-3BES-2 16%

DER-3

DER-1 BES-2

DER-3

25%

DER-2

DER-2

30%

24%

6

5

7

7%

(i2)MG-2

DER-2

42%

(i2)MG-2

DER-1

DER-2

(i2)MG-1

line-1

line-2

BES-1

17%

MG-3 DER-2

12%

MG-3 DER-4

11%

BES-2 17%

34%

(i2)MG-3

20%

6% 6%

9%

DER-4

(i2)

DER-4

14 18%

6% MG-2 DER-3

18%

MG-3 DER-3

MG-3

DER-2

DER-3 40%

MG-2 DER-2 MG-3 DER-1 MG-3 DER-2

MG-2 7% DER-3 MG-3 BES-1 8%

MG-3 DER-3

8% BES-1

MG-1 DER-2

7%

32%

22%

MG-2 DER-2 MG-2 DER-1

15 6%

MG-3 DER-1

7%

7% 17%

MG-2 BES-1

15%

6%

18%

MG-2 DER-3

MG-2 BES-1 9%

15%

MG-2 DER-2

19%

DER-1

DER-1

13

MG-3 DER-4

12 10%

DER-2

10%

10%

BES-2

25%

6%

10%

11

(i2)MG-2

49%

line-8 line-7 line-2

35%

line-3

18%

20%

20%

DER-2

MG-1

25%

line-4

20%

7%

25%

(i)MG-3 line-9 10%

line-2

10

line-3

line-6

line-3

20%

26%

9

(i)

line-1

BES-1

20%

line-1

10%

10%

(i)MG-2

(i2)MG-2 13% DER-1

41%

line-4

line-4

line-5

(i2)MG-2

DER-2

25%

line-7

DER-1

43%

(i2)MG-1

20%

11% (i2)MG-1 DER-2 (i2)MG-1

10%

8 10%

10%

8%

MG-3 DER-1

MG-3 11% DER-2

MG-3 BES-2 MG-3 BES-2 MG-2 BES-1 MG-3 BES-1

MG-3 DER-4 8% MG-3 DER-2

11%

8%

MG-3 BES-1

Fig. 13. Eigenvalue trajectory of CMG-C and participation factors of its eigenvalues.

of 1% was observed. Fig. 15 illustrates the time domain and statespace trajectory of active/reactive power of one of the DERs (i.e., DER-1 of MG-1) in CMG-C in each of the above 3 stability regions.

These results are captured from PSCAD/EMTDC for the same network and closely correspond to the expected stability performance of CMG-C from the developed MATLAB-based SSS analysis. The

56

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

Table 3 Characteristics of the critical eigenvalues of MG-2 for m and n variation.

k1,1 k2,20 k3,30 k4,40 0

r [1/s]

t [rad/s]

f [%]

xn [rad/s]

Tk [ms]

NoC50%

Critical for

14.699 8.624 5.966 314.159

±309.887 ±309.895 ±309.967 ±314.159

4.74 2.78 1.92 70.71

310.236 310.015 310.024 444.288

20 20 20 20

2.32 3.96 5.73 0.11

m, n m, n m, n m

Table 4 Characteristics of the critical eigenvalues of MG-3 for m and n variation.

k1,10 k2,20 k3,30 k4,40 k5,50 k6,60 k7,70

r [1/s]

t [rad/s]

f [%]

xn [rad/s]

Tk [ms]

NoC50%

Critical for

22.772 12.617 7.746 3.059 6.941 314.159 15.551

±309.933 ±309.948 ±309.913 ±310.093 ±309.883 ±314.159 ±44.959

7.33 4.07 2.50 0.99 2.24 70.71 32.69

310.769 310.205 310.009 310.108 309.961 444.288 47.572

20 20 20 20 20 20 139

1.50 2.70 4.41 11.18 4.92 0.11 0.31

m, m, m, m, m, m, m

n n n n n n

Table 5 Characteristics of the critical eigenvalues of CMG-A for m and n variation.

k1,10 k2,20 k3,30 k4,40 k5,50

r [1/s]

t [rad/s]

f [%]

xn [rad/s]

Tk [s]

NoC50%

Critical for

14.108 8.230 8.470 5.9486 314.159

±309.887 ±309.862 ±309.866 ±309.969 ±314.159

4.55 2.66 2.73 1.92 70.71

310.208 309.972 309.981 310.026 444.288

20 20 20 20 20

2.42 4.15 4.03 5.74 0.11

m, m, m, m, m

n n n n

Table 6 Characteristics of the critical eigenvalues of CMG-B for m and n variation.

k1,10 k2,20 k3,30 k4,40 k5,50 k6,60 k7,70 k8,80

r [1/s]

t [rad/s]

f [%]

xn [rad/s]

Tk [ms]

NoC50%

Critical for

20.926 3.051 11.735 8.359 6.903 7.729 314.159 314.159

±309.931 ±310.090 ±309.923 ±309.852 ±309.877 ±309.912 ±314.159 ±314.159

6.74 0.98 3.78 2.70 2.23 2.49 70.71 70.71

310.636 310.105 310.144 309.965 309.954 310.009 444.288 444.288

20 20 20 20 20 20 20 20

1.63 11.21 2.91 4.08 4.95 4.42 0.11 0.11

m, m, m, m, m, m, m m

n n n n n n

Table 7 Characteristics of the critical eigenvalues of CMG-C for m and n variation.

k1,1 k2,20 k3,30 k4,40 k5,50 k6,60 k7,70 k8,80 k9,90 k10, 100 k11,110 k12,120 k13,130 k14,140 k15,150 0

r [1/s]

t [rad/s]

f [%]

xn [rad/s]

Tk [ms]

NoC50%

Critical for

21.203 13.978 3.050 11.728 8.398 8.159 314.159 314.159 314.159 5.940 6.900 7.725 14.370 14.992 15.586

±309.930 ±309.895 ±310.090 ±309.923 ±309.858 ±309.878 ±314.159 ±314.159 ±314.159 ±309.976 ±309.877 ±309.911 ±38.141 ±42.616 ±45.194

6.83 4.51 0.98 3.78 2.71 2.63 70.71 70.71 70.71 1.92 2.23 2.49 35.26 33.19 32.60

310.655 310.210 310.105 310.145 309.972 309.985 444.288 444.288 444.288 310.033 309.954 310.007 40.758 45.176 47.807

20 20 20 20 20 20 20 20 20 20 20 20 164 147 139

1.61 2.44 11.21 2.91 4.07 4.18 0.11 0.11 0.11 5.75 4.95 4.42 0.29 0.313 0.319

m, m, m, m, m, m, m, m, m, m, m, m, m m m

time-domain active power response signals captured from PSCAD/ EMTDC have been also evaluated using Prony analysis to determine the eigenvalues of the corresponding system [57]. Among the defined eigenvalues, the dominant eigenvalues which have a stronger impact on the system stability have been illustrated in Fig. 15 and validate the defined stability status of each analyzed system. Similar results are also observed for the other DERs in each stability scenario in the other CMGs. Likewise, another study is car-

n n n n n n n n n n n n

ried out on n variations which also verifies the accuracy of the developed SSS analysis. Another study is also carried out to compare the outcome of the developed tool in MATLAB with the same system modeled in PSCAD/EMTDC during load variations. To this end, the above considered stable CMG-C is assumed to experience a demand increase of 17% at t1 = 0.33 s and a demand decrease of 31% at t2 = 0.66 s. Fig. 16an illustrates the time domain output active power of one

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

57

Fig. 14. Detection of (a) a stable CMG, (b) an unstable CMG, using the developed SSS analysis when both MGs of the CMG are either stable or marginally stable.

Fig. 15. Comparison of the PSCAD/EMTDC-based time-domain results of the active/reactive power for one of the DERs in the considered CMG (left figures) for different values of droop coefficient m with the corresponding state space trajectory results from MATLAB-based developed SSS tool (middle figures) and their corresponding Prony technique-based eigenvalues (right figures).

58

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

Fig. 16. (a) Active/reactive power for a DER in a considered stable CMG from time-domain analysis, (b) The corresponding state space trajectory results from the developed tool in MATLAB.

of the DERs in this CMG (i.e., DER-1 of MG-1) captured from the developed tool in MATLAB (denoted by a solid line) and the corresponding PSCAD/EMTDC-based result (denoted by dashed lines) around the steady-state points only. As it is seen from this figure, the two results match significantly at steady-state conditions while a small difference is observed in the system transients which is due to the fact that the conducted modeling in MATLAB does not consider the high-frequency dynamics of the system caused by the VSCs. A similar difference can also be observed in the transients of the output reactive power of the considered DER, which is not illustrated. Fig. 16b shows the state-space trajectory of active/reactive power of the considered DER in the CMG within this period. The time-domain results captured from PSCAD/EMTDC closely correspond to the state space trajectory given by the developed MATLAB-based analysis. Similar results are also observed for the other DERs.

SoC, RoD, energy capacity and internal parameters of BESes on the stability of a CMG will be studied. The developed tool will also be used in another future research to develop a mechanism to stabilize a new system which is determined unstable or determined to have small stability margins. This study focused on sustainable MGs of future remote areas, in which all electricity demand is generated by the help of converter-interfaced resources. The developed technique can be expanded to include the presence of inertial sources such as diesel/gas-based generators if they share a portion of the load. Furthermore, this research did not consider a cost-prioritized droop control, which can be incorporated in the developed technique in future research. Appendix A. Expansion of equations Eq. (14) can be expanded as

6. Conclusions Neighboring MGs, operating isolated and independently, may interconnect temporarily to support each other. This can occur when an MG is unexpectedly overloaded or when the generation level of one or more of its non-dispatchable DERs is reduced. Furthermore, the interconnection may be financially beneficial when a neighboring MG sells electricity with a lower price. Before the formation of a CMG, the stability of the new system needs to be examined to avoid an interconnection that can cause instability for the system. An eigenanalysis-based SSS evaluation tool was developed in this research which can be coded on a digital signal processor and used within the network controller. Through some numerical analyses, the performance of the developed tool was illustrated. The accuracy of the results was verified by comparing the timedomain performance of the considered examples with the findings of the developed tool. By the help of this developed tool, in the subsequent publications, the impact of different number, internal parameters and ratings of the DERs and loads will be evaluated for the stability of a CMG. The research will be further expanded to investigate the influence of the number of the MGs as well as the number, length and X/R ratio of lines. Moreover, the effect of different number,

2

3

2 Rf =Lf 6 q 7 6 x 6 i1 7 6 6 7 6 6 d 7 6 0 6 i2 7 D6 q 7 ¼ 6 6 0 6i 7 6 6 2 7 6 6 d7 4 1=C f 4 vc 5 0 v qc 2 d 3 2 1 i1 6 q 7 60 6 i1 7 6 6 7 6 d 7 aV 6 0 6 i2 7 dc 6 6 D6 q 7 þ 6i 7 Lf 6 60 6 2 7 6 6 d7 40 4v 5 d

i1

c

v qc

0

x

0

0

1=Lf

Rf =Lf

0

0

0

0

0

x

1=LT

0

x

0

0

0

1=C f

0

0

1=C f x 3 2 0 0 0 60 07 17 7 " # 7 " # 6 7 7 6 d d 61 07 07 7:D u  1 6 7:D v t 7 6 q L 07 0 1 u v qt T6 7 7 7 7 6 5 5 4 0 0 0 1=C f 3

0

0

0

0

3

1=Lf 7 7 7 0 7 7: 1=LT 7 7 7 x 5 0

ðA1Þ

0

Matrices ZDER, G, H, and K in (16) are expressed as T

ZDER ¼ ½i2 i2 ic ic v dc v qc   k2 0 k2  k1 G¼ 0 0 k2 d q d q

0 k2  k1

k3 0

0 k3



K ¼ ½ k1 k1 k2 k2 k3 k3    k1 0 k2 0 k3 0 H¼ 0 k1 0 k2 0 k3

ðA2Þ

59

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60 Table B1 Parameters of the network under consideration in Figs. 6 and 7. MG-1 and MG-2 parameters: 3-phase, 4-wire, multiple earthed neutral (MEN) type, Vbase = 230 V, Vmax = 235.75 V, Vmin = 224.25 V, fnominal = 50 Hz, f f min = 49.5 Hz, Zload = 50.57 + j 16.6504 X/phase, Zline = 0.1 + j 1 X/phase cap DER parameters: P cap DER = 3 kW, SDER = 3.9 kVA, LT = 13.6 mH

max

= 50.5 Hz,

cap cap BES parameters: Ecap BES ¼ 5 kW h, RoDBES = 2 kW, SBES = 2.6 kVA CMG parameters: Ztie-line MG-1,2 = 0.1 + j 1 X/phase, Ztie-line MG-2,3 = 0.1 + j 1 X/phase VSC Structure, filter and control parameters: Vdc = 350 V, a = 2, Rf = 0.1 X, Lf = 0.37 mH, Cf = 50 lF, k1 = 2.2968, k2 = 6.1802, k3 = 22.968, xc = 31.4159 rad (i.e. 5 Hz)

Eq. (18) and (19), (23), (24), and (35) are expanded respectively in (A3)–(A7) below.

"

D

i

d

#

 ¼

q

i

Rload =Lload x

" #  " d# x vd I2
2 i :D q þ :D tq L Rload =Lload vt load i

"

d

i q i

#

 " d# x Rline =Lline i :D q x Rline =Lline i 2 3 busð1Þ to busðkk1Þ busðkkÞ busðkkþ1Þ to busðll1Þ busðllÞ busðllþ1Þ to busðN bus Þ 1 6 zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ z}|{ zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ zfflffl}|fflffl{ zfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{ 7 þ 4 O2
2 . . . O2
2 I2
2 O2
2 . . . O2
2 I2
2 O2
2 . . . O2
2 5 Lline 

¼

2

busðkk1Þ

busð1Þ

busðll1Þ

The technical data of the studied network in Figs. 6 and 7 with the VSC structure and closed-loop control of Figs. 3 and 4 are given in Table B1.

busðllÞ

ðv qc Þo

_ ¼ DSoC

3 ½ 0 0 ðv dc Þo 2Ecap

:D½ id1

q

i1

d

i2

q

i2

ðA4Þ

ði2 Þo 

d

q

ði2 Þo

ðA5Þ

3 DQ_ ¼ xc DQ þ xc ½ 0 0 ðv qc Þo 2 T :D½ id1 iq1 id2 iq2 v dc v qc  ðv qc Þo

v dc v qc 

ðv dc Þo

q

ði2 Þ0

ði2 Þo  d

ðA6Þ d

ði2 Þo

ði2 Þo  q

T

ðA7Þ

where I2
2 represents the identity matrix of order 2. Matrix MiDER in (38) is expressed as

2

c11 6c 6 21

6 6 0 MiDER ¼ 6 6 0 6 6 4 0

3 c12 c13 7 c22 c23 7 7 C f nx cos do  C f V cmo x sin do 7 7and  C f nx sin do C f V cmo x cos do 7 7 n sin do

if DER  j is connected to bus  jj otherwise if BES  j is connected to bus  jj and is being discharged if BES  j is connected to bus  jj and is being charged otherwise ðA9Þ if load  j is connected to bus  jj otherwise if line  k enters bus  jj if BES  k exits bus  jj otherwise

Appendix B. Technical data

zfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflffl}|fflfflfflfflfflfflffl{ kkþ1 ll1 ll kkþ1 ll1 ll
ðv dt Þ ðv qt Þ . . . ðv dt Þ ðv qt Þ ðv dt Þ ðv qt Þ 3 T busðllþ1Þ busðN bus Þ zfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflffl{ 7 q llþ1 q N bus 7 d llþ1 d N bus
ðv t Þ ðv t Þ . . . ðv t Þ ðv t Þ 5

3 DP_ ¼ xc DP þ xc ½ 0 0 ðv dc Þo 2 T :D½ id1 iq1 id2 iq2 v dc v qc 

I2
2 O2
2 8 > < I2
2 lBESj ¼ I2
2 > : O2
2  I2
2 loadj l ¼ O2
2 8 > < I2
2 llinek ¼ I2
2 > : O2
2

busðkkÞ

zfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{ 6 d 1 q 1 kk kk1 kk d kk1 ðv qt Þ ðv dt Þ ðv qt Þ
D6 4ðv t Þ ðv t Þ . . . ðv t Þ busðkkþ1Þ



lDERj ¼

ðA3Þ D

Also, lDERj , lBESj , lloadj , and llinek in (44) are expressed as

V cmo cos do

7 5

0 n cos do V cmo sin do 8 c ¼ 2 sin d =3V > o cmo > > 11 > > > c ¼ 2½V cmo cos do þ nðQ o cos do  P o sin do Þ=3V 2cmo > 12 > > < c ¼ 2ðQ sin d þ P cos d Þ=3V o o o cmo o 13 ðA8Þ > c21 ¼ 2 cos do =3V cmo > > > > > c ¼ 2½V cmo sin do þ nðQ o sin do þ Po cos do Þ=3V 2 > 22 cmo > > : c23 ¼ 2ðQ o cos do  Po sin do Þ=3V cmo where subscript o denotes the operating point values for each state.

References [1] Lasseter RH. MicroGrids. IEEE Power Eng Soc Winter Meeting 2002;1:305–8. [2] Xu Z, Nthontho M, Chowdhury S. Rural electrification implementation strategies through microgrid approach in South African context. Int J Electr Power Energy Syst 2016;82:452–65. [3] Guerrero JM, Vasquez JC, Matas J, et al. Hierarchical control of droop-controlled ac and dc microgrids-a general approach toward standardization. IEEE Trans Ind Electr 2011;58(1):158–72. [4] Guerrero JM, Chandorkar M, Lee TL, Loh PC. Advanced control architectures for intelligent microgrids-part I: decentralized and hierarchical control. IEEE Trans Ind Electr 2013;60(4):1254–62. [5] Yazdanian M, Mehrizi-Sani A. Distributed control techniques in microgrids. IEEE Trans Smart Grid 2014;5(6):2901–9. [6] Hajipour E, Bozorg M, Fotuhi-Firuzabad M. Stochastic capacity expansion planning of remote microgrids with wind farms and energy storage. IEEE Trans Sust Energy 2015;6(2):491–8. [7] Lasseter RH. Smart distribution: coupled microgrids. Proc IEEE 2011;99 (6):1074–82. [8] Wang Z, Chen B, Wang J, Chen C. Networked microgrids for self-healing power systems. IEEE Trans Smart Grid 2016;7(1):310–9. [9] Pashajavid E, Shahnia F, Ghosh A. Development of a self-healing strategy to enhance the overloading resilience of islanded microgrids. IEEE Trans Smart Grid http://dx.doi.org/10.1109/TSG.2015.2477601. [10] Wang Z, Chen B, Wang J, et al. Coordinated energy management of networked microgrids in distribution systems. IEEE Trans Smart Grid 2015;6(1):45–53. [11] Wang Z, Wang J. Self-healing resilient distribution systems based on sectionalization into microgrids. IEEE Trans Power Syst 2015;30(6):3139–49. [12] Shahnia F, Bourbour S, Ghosh A. Coupling neighboring microgrids for overload management based on dynamic multicriteria decision-making. IEEE Trans Smart Grid http://dx.doi.org/10.1109/TSG.2015.2477845. [13] Zhang Y, Xie L, Ding Q. Interactive control of coupled microgrids for guaranteed system-wide small signal stability. IEEE Trans Smart Grid 2016;7 (2):1088–96. [14] Shahnia F, Chandrasena RPS, Rajakaruna S, Ghosh A. Primary control level of parallel distributed energy resources converters in system of multiple interconnected autonomous microgrids within self-healing networks. IET Gener Transm Distrib 2014;8:203–22. [15] Zhang Y, Xie L. Online dynamic security assessment of microgrid interconnections in smart distribution systems. IEEE Trans Power Syst 2015;30(6):3246–54. [16] Nikolakakos IP, Zeineldin HH, El-Moursi MS, Hatziargyriou ND. Stability evaluation of interconnected multi-inverter microgrids through critical clusters. IEEE Trans Power Syst 2016;31(4):3060–72. [17] Nikmehr N, Ravadanegh SN. Reliability evaluation of multi-microgrids considering optimal operation of small scale energy zones under loadgeneration uncertainties. Int J Electr Power Energy Syst 2016;78:80–7.

60

F. Shahnia, A. Arefi / Electrical Power and Energy Systems 91 (2017) 42–60

[18] Arefifar SA, Ordonez M, Mohamed Y. Voltage and current controllability in multi-microgrid smart distribution systems. IEEE Trans Smart Grid http:// dx.doi.org/10.1109/TSG.2016.2568999. [19] Majumder R, Bag G. Parallel operation of converter interfaced multiple microgrids. Int J Electr Power Energy Syst 2014;55:486–96. [20] Che L, Zhang X, Shahidehpour M, et al. Optimal interconnection planning of community microgrids with renewable energy sources. IEEE Trans Smart Grid http://dx.doi.org/10.1109/TSG.2015.2456834. [21] Coelho EAA, Cortizo PC, Garcia PFD. Small-signal stability for parallelconnected inverters in stand-alone ac supply systems. IEEE Trans Ind Appl 2002;38(2):533–42. [22] Middlebrook RD. Input filter consideration in design and application of switching regulators. IEEE Ind Appl Soc Annu Meeting 1976:366–82. [23] Belkhayat M. Stability criteria for ac power systems with regulated loads [Ph. D. thesis]. USA: Purdue University; 1997. [24] Haddadi A, Boulet B, Yazdani A, Joos G. A l-based approach to small-signal stability analysis of an interconnected distributed energy resource unit and load. IEEE Trans Power Delivery 2015;30(4):1715–26. [25] Sun J. Small-signal methods for ac distributed power systems–a review. IEEE Trans Power Electr 2009;24(11):2545–54. [26] Majumder R. Some aspects of stability in microgrids. IEEE Trans Power Syst 2013;28(3):3243–52. [27] Pogaku N, Prodanovic M, Green TC. Modeling, analysis and testing of autonomous operation of an inverter-based microgrid. IEEE Trans Power Electr 2007;22(2):613–25. [28] Marwali MN, Jung JW, Keyhani A. Stability analysis of load sharing control for distributed generation systems. IEEE Trans Energy Convers 2007;22 (3):737–45. [29] Barklund E, Pogaku N, Prodanovic M, et al. Energy management in autonomous microgrid using stability-constrained droop control of inverters. IEEE Trans Power Electr 2008;23(5):2346–52. [30] D’Arco S, Suul JA, Fosso OB. Small-signal modeling and parametric sensitivity of a virtual synchronous machine in islanded operation. Int J Electr Power Energy Syst 2015;72:3–15. [31] Majumder R, Chaudhuri B, Ghosh A, et al. Improvement of stability and load sharing in an autonomous microgrid using supplementary droop control loop. IEEE Trans Power Syst 2010;25(2):796–808. [32] Rasheduzzaman M, Mueller JA, Kimball JW. An accurate small-signal model of inverter-dominated islanded microgrids using dq reference frame. IEEE J Emerg Selected Topics Power Electr 2014;2(4):1070–80. [33] Shahnia F. Stability and eigenanalysis of a sustainable remote area microgrid with a transforming structure. Sust Energy Grids Networks 2016;8:37–50. [34] Katiraei F, Iravani MR, Lehn PW. Small-signal dynamic model of a micro-grid including conventional and electronically interfaced distributed resources. IET Gener Transm Distrib 2007;1(3):369–78. [35] Tang X, Deng W, Qi Z. Investigation of the dynamic stability of microgrid. IEEE Trans Power Syst 2014;29(2):698–706. [36] Guan Y, Vasquez JC, Guerrero JM, et al. Frequency stability of hierarchically controlled hybrid photovoltaic-battery-hydropower microgrids. IEEE Trans Ind Appl 2015;51(6):4729–42.

[37] Mishra S, Ramasubramanian D. Improving the small signal stability of a PVDE-dynamic load-based microgrid using an auxiliary signal in the PV control Loop. IEEE Trans Power Syst 2015;30(1):166–76. [38] Jamehbozorg A, Radman G. Small signal analysis of power systems with wind and energy storage units. IEEE Trans Power Syst 2015;30(1):298–305. [39] Jayawardena AV, Meegahapola LG, Robinson D, Perera S. Representation of a grid-tied microgrid as a reduced order entity for distribution system dynamics and stability studies. Int J Electr Power Energy Syst 2015;73:591–600. [40] Mariani V, Vasca F, Vasquez JC, Guerrero JM. Model order reductions for stability analysis of islanded microgrids with droop control. IEEE Trans Ind Electr 2015;62(7):4344–54. [41] Ahumada C, Cardenas R, Saez D, Guerrero JM. Secondary control strategies for frequency restoration in islanded microgrids with consideration of communication delays. IEEE Trans Smart Grid 2016;7(3):1430–41. [42] Bottrell N, Prodanovic M, Green TC. Dynamic stability of a microgrid with an active load. IEEE Trans Power Electr 2013;28(11):5107–19. [43] Radwan AAA, Mohamed YARI. Stabilization of medium-frequency modes in isolated microgrids supplying direct online induction motor loads. IEEE Trans Smart Grid 2014;5(1):358–70. [44] Uriarte FM, Smith C, VanBroekhoven S, Hebner RE. Microgrid ramp rates and the inertial stability margin. IEEE Trans Power Syst 2015;30(6):3209–16. [45] Eid BM, Rahim NA, Selvaraj J, Khateb AHE. Control methods and objectives for electronically coupled distributed energy resources in microgrids: a review. IEEE Syst J 2016;10(2):446–58. [46] Bevrani H, Shokoohi S. An intelligent droop control for simultaneous voltage and frequency regulation in islanded microgrids. IEEE Trans Smart Grid 2013;4 (3):1505–13. [47] Nutkani IU, Loh PC, Blaabjerg F. Cost-prioritized droop schemes for autonomous ac microgrids. IEEE Trans Power Electr 2015;30(2):1109–19. [48] Hosseinimehr T, Ghosh A, Shahnia F. Cooperative control of battery energy storage systems in microgrids. Int J Electr Power Energy Syst 2017;87:109–20. [49] Kim YS, Kim ES, Moon SI. Frequency and voltage control strategy of standalone microgrids with high penetration of intermittent renewable systems. IEEE Trans Power Syst 2016;31(1):718–28. [50] Yazdani A, Iravani R. Voltage-sourced converters in power systems. Wiley; 2010. [51] Ghosh A, Ledwich G. Power quality enhancement using custom power devices. Kluwer Academic; 2002. [52] Gopal M. Digital control and state variable methods. Tata McGraw-Hill; 2009. [53] Tewari A. Modern control design. Wiley; 2002. [54] Kundur P. Power system stability and control. McGraw-Hill; 1994. [55] Shahnia F. Eigenvalues and participation factors of a system of coupled adjacent microgrids, In: 2nd IEEE annual Southern Power Electronics Conference (SPEC), Auckland; 2016. [56] Kelly SG. Mechanical vibrations. Cengage Learning; 2013. [57] Hauer JF, Demeure CJ, Scharf LL. Initial results in Prony analysis of power system response signals. IEEE Trans Power Syst 1990;5(1):80–9.