Modeling and stability analysis of multi-time scale DC microgrid

G Model

ARTICLE IN PRESS

EPSR-4671; No. of Pages 11

Electric Power Systems Research xxx (2016) xxx–xxx

Contents lists available at ScienceDirect

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

Modeling and stability analysis of multi-time scale DC microgrid Nanfang Yang a,∗ , Babak Nahid-Mobarakeh b , Fei Gao a , Damien Paire a, Abdellatif Miraoui a , Weiguo Liu c a b c

Université de Technologie de Belfort-Montbéliard, 90100 Belfort, France Université de Lorraine, 54500 Vandœuvre-lès-Nancy, France Northwestern Polytechnical University, 710069 Xi’an, China

a r t i c l e

i n f o

Article history: Received 21 January 2016 Received in revised form 20 March 2016 Accepted 14 April 2016 Available online xxx Keywords: DC microgrid Stability analysis Multi-time scale system

a b s t r a c t The diverse dynamic characteristics of distributed generators (DGs) make the DC microgrid to be a multitime scale system. This paper investigates the dynamic modeling of multi-time scale DC microgrid, and a reduced-order multi-scale model (RMM) is proposed to reduce the system model complexity as well as to conserve major time scale information. The DGs with similar time constants are grouped together to form an equivalent DG, the equivalent DGs are then combined to construct RMM. The effectiveness of the proposed RMM is verified by stability analysis, numerical simulations and experimental tests. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Microgrids are proposed to supply reliable, clean and efficient power to local users, in which multiple distributed generators (DGs) are usually connected in parallel to improve reliability instead of using the central power plant. It avoids the long distance power transportation, thus can obtain higher efficiency. From the view of generation, many DGs, e.g., photovoltaic (PV), fuel cell (FC), and wind power are inherently DC type, or need to be converted in to DC at first. From the view of consumption, common loads, e.g., electronic loads, variable drives in air conditioners, refrigerators and washing machines all require DC power. Therefore, the adoption of DC microgrid can enhance overall system efficiency due to the reduction power conversion stags. Besides, it can avoid the complicated synchronization in AC grid [1]. The control of DC microgrid with multiple DGs falls into two popular methods: master–slave control and droop based control. Master–salve control is a single master method based on fast communication, where the master controls the bus voltage and sets current references to slaves. The reliability highly depends on the master unit and the communication line. However, droop control is a multi-master method (decentralized method), where each unit participates into bus voltage control. Thus system reliability is largely enhanced and fast communication is not vital.

∗ Corresponding author. Tel.: +33 384582021. E-mail address: [email protected] (N. Yang).

In the analysis of DC microgrid under droop control, usually DGs’ dynamics are not considered or assumed to be equal. The DG under droop control is modeled in Thévenin form to be an imperfect voltage source (a perfect voltage source and a resistor in series) [2], to investigate the DC-bus voltage control and load sharing issues, considering the connecting cables, communication, and voltage reference offsets [3–6]. Then the DC microgrid becomes the parallel of imperfect voltage sources [7]. This model is adopted in most research without considering the dynamics of DGs. However, in real situations, DGs may have largely diverse dynamic responses, referred as time scales or frequency scales. For example, the grid converter connecting the microgrid and utility grid requires a smooth power exchange, which indicates a small frequency scale [8]. Fuel cells are limited to slower dynamics to benefit long life-span; while super-capacitors can be used to absorb/supply high frequency power during limited time range, which have fast dynamics [9]. Although the power distribution in frequency spectrum has been investigated by some authors [10,11,8], the detailed dynamic modeling and stability analysis of DC microgrid with multi-time scales are still uncovered. These are the focus of this paper, and the main contributions are list as follows: 1 A comprehensive model (CM) is constructed by introducing a virtual inductor in series with droop resistor to represent DG’s dynamics, and then a reduced 2nd-order model (R2M) is deduced by using arithmetic mean values to replace distributed parameters.

http://dx.doi.org/10.1016/j.epsr.2016.04.014 0378-7796/© 2016 Elsevier B.V. All rights reserved.

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014

G Model

ARTICLE IN PRESS

EPSR-4671; No. of Pages 11

N. Yang et al. / Electric Power Systems Research xxx (2016) xxx–xxx

2

Fig. 3. The equivalent circuit of the DG under time scale droop control. Fig. 1. A general DC microgrid with three DGs, Z1, Z2 and Z3 are the line impedances.

2 A novel reduced-order multi-scale model (RMM) is proposed, in which DGs with similar time constants are grouped together, such that it reduces model complexity and keeps major dynamics.

With a unit voltage feedback, the input–output voltage closed-loop transfer function Gvc (s) is: Gvc (s) =

Dj ωj

(2)

Cj s2 + Cj ωj s + Dj ωj

Then, the voltage error to source current transfer function is: 2. DC microgrid dynamic modeling

Yj (s) =

A general DC microgrid with three DGs is shown in Fig. 1, in which the power converter controlled DGs are connected separately through cables to the common load point, where loads are connected. Every DG participates into voltage control as well as load sharing under droop control. 2.1. Equivalent circuit of DGs under droop control The traditional droop control only considers power scale, i.e., power-voltage droop in DC grids or real power-frequency droop and reactive power-voltage droop in AC grids. DG’s dynamics can be introduced as another freedom, referred as the frequency scale or time scale [10], to form time-scale droop control. It can be realized by a forward path low-pass filter (LPF), feedback LPF added to the traditional droop control [9], or directly by using a PI controller in the voltage control with a virtual droop resistor. The analysis of these three implementations are similar. Let’s take the first method with a forward path LPF as an example, as shown in Fig. 2. Cascaded control structure is adopted, where the inner current loop is much faster than other loops, and it can be assumed to be unit Gcur (s) = 1 [12]. Then the input–output voltage open-loop transfer function Gvo (s) is given by: Gvo (s) =

Dj ωj sC j (s + ωj )

(1)

=

Dj ωj s + ωj

=

1/Rdj j s + 1

(3)

where Voj is the output voltage of the DG as well as the voltage of the output capacitor; Rdj = 1/Dj is the virtual droop resistance; Ij is the current supplied by the DG;  j = 1/ωj is the dominant time constant of the DG, which is used to represent the DG’s time scale. According to these analyses, the equivalent circuit of droop controlled DG with time scale is deduced and shown in Fig. 3, in which the time scale is represented by the droop resistor and virtual inductor in series. A parallel capacitor is used to represent the output voltage dynamics. The time constant  j is the ratio of virtual inductance Ldj over the droop resistance Rdj . j =

Ldj Rdj

=

1 ωj

(4)

This equivalent circuit not only considers the effect of droop control in power scale but also time scale, thus more precise than the traditional model using only voltage source and resistor in series [13] or current source and capacitor in parallel [14]. An alternative equivalent circuit can be deduced by using Norton form: a perfect current source, resistor, inductor and capacitor connecting in parallel. The dynamic model of the DG in state-space form is:

⎧ dI Rdj 1 j ⎪ ⎨ dt = − L Ij + L (Vj − Voj ) dj

dj

⎪ ⎩ dV oj = 1 (Ij − Ioj ) dt

where the subscript j denotes the jth DG; Dj is the droop constant; ωj is the cut-off frequency of the LPF; and Cj is the output capacitance.

Ij Vj − Voj

(5)

Cj

where Ioj is the current injected by the DG to common DC-bus.

Fig. 2. Local voltage control using forward path LPF to realize time scale droop control.

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014

G Model EPSR-4671; No. of Pages 11

ARTICLE IN PRESS N. Yang et al. / Electric Power Systems Research xxx (2016) xxx–xxx

3

According to KCL the model of the input capacitor at CPL is written as: dV L 1 = (IL1 − IL11 ) CL dt Fig. 4. The equivalent circuit of the connecting cable between the jth DG to the common load point.

2.2. Equivalent circuit of connecting cables The lamped model of connecting cables can be represented by  or  type equivalent circuit. In the latter, the capacitors can be viewed distributed to the output of DGs and the input of loads, and then the equivalent circuit becomes the series connection of a resistor and inductor. Considering DGs in the microgrid are connected to the common load point through separate cables. The equivalent circuit of the connecting cable between the jth DG and the common load point is given by Fig. 4, and the dynamic cable model is: dI cj dt

=−

Rcj Lcj

Icj +

1 (V − VL ) Lcj oj

(10)

where IL1 denotes the current supplying to CPL; and CL is the equivalent input capacitance of loads. A non-idea CVL or a combination of CRL and CCL can also be represented by a resistor and a current sink connecting in parallel, as shown in Fig. 5b. The current flowing into this load is: IL2 =

VL + ICCL RCRL

(11)

where RCRL denotes the equivalent resistance of CRL; ICCL denotes the equivalent current sink of CCL. The total load current is the sum of all the load currents: IL = IL1 + IL2

(12)

where IL denotes the total current flowing to loads. (6)

2.4. Comprehensive model

where Icj denotes the current flowing from DG to loads through the connecting cable, and it equals to the output current of the DG; Rcj , Lcj are equivalent cable resistance and inductance, respectively; and VL is the voltage of loads at the common load point.

Considering a single bus DC microgrid with n DGs, they are connected to the common load point through separate cables. Then the model of DGs can be expressed in matrix notion:

2.3. Equivalent circuit of general loads

d I = −W d I + Y d (V − V o ) dt

Neglecting fast dynamics, a general DC microgrid load can be constant resistive load (CRL), constant current load (CCL), constant power load (CPL), or constant voltage load (CVL). Thus a general load can be written as: PL = PCPL + PCRL + PCCL + PCVL

(7)

where PCPL , PCRL , PCCL and PCVL are the powers of CPL, CRL, CCL and CVL at rated DC-bus voltage, respectively. The nonlinear relationship of load current and load voltage in an ideal CPL is written by: IL11

PCPL = VL

(8)

where IL11 is the current absorbed by CPL. The linearized model can be obtained by using Taylor expansion at the operating load voltage Ve . The equivalent circuit is given by Fig. 5a, and the load current is: IL11 ≈ ICPL +

VL PCPL PCPL =2 − 2 VL RCPL Ve Ve

(9)

The equivalent circuit is composed of an equivalent negative resistor RCPL = −Ve2 /PCPL and an equivalent current sink ICPL = 2PCPL /Ve .

(13)

where the state variables vector is I = [I1 , I2 , . . ., In ]T ; the voltage references vector is V = [V1 , V2 , . . ., Vn ]T and the output voltages vector is V o = [Vo1 , Vo2 , . . ., Von ]T . The superscript T indicates the transpose of vector. The state and input matrices are: W d = diag

Y d = diag

R

d1

Ld1

 1

,

Ld1

,

Rd2 R , . . ., dn Ld2 Ldn

1 1 , . . ., Ld2 Ldn



n×n

 n×n

where the diagonal components of state matrix W d are the inverses of DGs’ time constants. The model of the output capacitance is rewritten by: d V o = G o (I − I o ) dt

(14)

where the output current vector is I o = [Io1 , Io2 , . . ., Ion ]T . The state matrix is given: G o = diag

1

C1

,

1 1 , · · ·, C2 Cn



n×n

Since all DGs are viewed connecting through separate cables to the common load point, the output current equals to the current goes through the corresponding connecting cable. Icj = Ioj

or

Ic = Io

(15)

where I c = [Ic1 , Ic2 , . . ., Icn ]T denotes the connecting cable currents flowing into the common load point. The differential equations of n connecting cables are: d I o = −W c I o + Y c (V o − EVL ) dt

(16)

where E = [1, 1, . . ., 1]Tn×1 ; and the state and input matrices are: W c = diag

Fig. 5. The equivalent circuit of the (a) constant power load, (b) combined constant current load and constant resistance load or constant voltage load.

Y c = diag

R

c1

Lc1

 1 Lc1

,

,

Rc2 Rcn , · · ·, Lc2 Lcn

1 1 , · · ·, Lc2 Lcn



n×n

 n×n

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014

G Model

ARTICLE IN PRESS

EPSR-4671; No. of Pages 11

N. Yang et al. / Electric Power Systems Research xxx (2016) xxx–xxx

4

Fig. 6. The equivalent circuit of CM.

The order of CM is 3n + 1, with n the number of DGs. The system order increases linearly with the number of DGs, thus model reduction technologies need to be performed to reduce its complexity.

The total source current Is = IL supplied by DGs is: Is =

n 

Ioj = E T I o

(17)

i=1

3. Model reduction

Substituting (11), (12) and (17) into (10), the load model can be rewritten as: dV L 1 T 1 1 = E Ic − VL − ICC CL RL CL CL dt

3.1. Reduced 2nd-order model

(18)

A popular model reduction method is using lumped mean value to replace the distributed values, when they are similar [15]. In particular, this method is applicable to small-scale microgrids with similar DGs, thus:

where RL denotes the equivalent total load resistance: RL = RCRL RCPL

(19)

⎧ R ∼ Rd2 ∼ Rd ⎪ = = ··· ∼ ⎨ L d1 = Ld2 Ld d1 1 1 1 ⎪ ⎩ ∼ ∼ ∼

and ICC denotes the equivalent total current sink: ICC = ICPL + ICCL

(20)

Then CM can be obtained by combining the differential equations of DGs (13), output capacitors (14), connecting cables (16), and the load (18). The equivalent circuit of CM is shown in Fig. 6 and the linearized model in matrix notation is rewritten as: d X = AX + BU dt

−W d

⎢0 ⎢ n×n ⎢ A=⎢ ⎢ Go ⎣

01×n

0n×n

−Y d

0n×1

−W c

Yc

−Y c E

−G o

0n×n

0n×1

1 T E CL

01×n

−

1 RL CL

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

(3n+1)×(3n+1)

⎡

Yd

0n×1

Ld2

= ··· =

(22)

Ld

where Rd and Ld are the mean values of Rdj and Ldj , respectively. Assume all DGs have similar voltage references as:

(21)

where the state variables vector is X = [ I, I c , V o , VL ]T ; and the input variables vector is U = [ V , ICC ]T . The state and input matrices are:

⎡

Ld1

=

V1 ∼ = V2 ∼ = Vn ∼ = ··· ∼ = VN

(23)

⎤

⎢0 ⎥ ⎢ n×n 0n×1 ⎥ ⎢ ⎥ B=⎢ ⎥ ⎢ 0n×n 0n×1 ⎥ ⎣ ⎦ 01×n

−

1 CL

(3n+1)×(n+1)

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014

G Model

ARTICLE IN PRESS

EPSR-4671; No. of Pages 11

N. Yang et al. / Electric Power Systems Research xxx (2016) xxx–xxx

5

leads to R2M, as shown in Fig. 7, and given as:

d dt





I VL

⎡

−

Rs

−

⎢ Ld =⎣ 1

−

C dc

1

⎤

Ld 1

⎥ ⎦



I



VL

⎡ 1 ⎢L +⎣ d 0

C dc RL

⎤ 0 −

1 C dc

  ⎥ VN ⎦ I CC

(30)

where Rs = Rd + Rc denotes the equivalent source resistance; C dc = C + C L is the equivalent DC-bus capacitance. 3.2. Reduced-order multi-scale model

Fig. 7. Equivalent circuit of R2M.

where VN is the common voltage reference of the DGs. Then the dynamic model of DGs (13) can be reduced to be: dI R 1 =− dI+ (VN − V o ) dt Ld Ld

(24)

where I and V o are the mean values of Ij and Voj , respectively. The model of output capacitance (14) is reduced to be: dV o dt

=

1 C

(I − I c )

(25)

where C is the mean value of Cj . Supposing all the connecting cables between DGs and the common load are the same type, then the impedance per length is the same as well as the ratio of resistance to inductance [15]. Assuming close distances, it yields:

Although R2M can reduce system model complexity significantly, the time scale information is not well conserved. This may lead to large errors during the analysis of multi-scale systems. To reduce the model order and keep major time scale information, DGs with similar time constant can be grouped together to form an equivalent DG, then the equivalent DGs are combined to built RMM. The order of RMM is m + 1 with m the number of time constant groups. If the DGs are divided into two groups: one group has slow dynamics (large time constant) while the other has fast dynamics (small time constant), then the (m+1)-order model becomes a 3rdorder model as (31). The grouping of multiple DGs will be discussed in Section 4.5.

⎡ ⎡

⎤

L1L

L2L

Lc

L1L

L2L



C dc

⎡ 1 ⎢ Ld ⎢ 1 ⎢ + ⎢  ⎢ Ld ⎣

(26)

0

Lc

where Rc and Lc denote the mean values of Rcj and Lcj , respectively. Then the model of cables in (16) can be combined and represented by: dI c Rc 1 = − I c + (V o − VL ) dt Lc Lc

(27)

where the state variable I c denotes the average value of Icj . The corresponding load equation for the equivalent mean DG is: 1 dV L = dt CL



Ic −

VL RL



− I CC

(28)

⎧ R = n × RL ⎪ ⎪ L ⎪ ⎨ ICCL ICPL I CCL = CL =

n

+

n



−

0

C dc

− −

⎤

1 

Ld 1 

Ld 1

⎤ ⎥ ⎡ I ⎤ ⎥ ⎥ ⎢  ⎥ ⎥⎣ I ⎦ ⎥ ⎥ ⎦ VL

C dc RL

0

−

⎥  ⎥ VN ⎥ 0 ⎥ ⎥ I CC 1 ⎦

(29)

CL n

Most often, in small-scale DC microgrids, the droop control time constants  j are much higher than the connecting cables time constants. The higher bandwidth filter introduced by the connecting cable has very limit influence thus it can be neglected. The cable resistor is moved to fuse with virtual droop resistor, meanwhile output capacitors and input capacitors are combined together. This

(31)

C dc



where Rs and Rs denote source resistances of the equivalent slow  DG (ESDG) and the equivalent fast DG (EFDG), respectively; Ld and  Ld denote the equivalent virtual inductances of ESDG and EFDG,   respectively; I and I denote the source currents of ESDG and EFDG. The equivalent parameters are calculated by:

⎧ p 1  ⎪ ⎪ = (Rdj + Rcj ) R ⎪ s ⎪ p ⎪ ⎪ j=1 ⎪ ⎪ ⎪ n ⎪  ⎪ 1  ⎪ ⎪ Rs = (Rdj + Rcj ) ⎪ n−p ⎨ j=p+1

where the equivalent load parameters are given as:

⎪ ⎪ ⎪ ⎩

Rs

⎢ Ld ⎢  ⎢ d ⎢  ⎥ ⎢ 0 − Rs = ⎣ I ⎦ ⎢  dt ⎢ Ld ⎣ 1 VL 1 I

⎧ R2L ∼ R Rc ⎪ ⎨ 1L ∼ = = = ··· ∼ 1 ∼ 1 ⎪ ⎩ 1 ∼ = = ··· ∼ =



−

p ⎪ 1  ⎪ ⎪ L = Ldj ⎪ s ⎪ p ⎪ ⎪ j=1 ⎪ ⎪ n ⎪  ⎪  ⎪ 1 ⎪ L = Ldj ⎪ ⎩ s n−p

(32)

j=p+1

where the jth DG belongs to ESDG when j = 1, 2, . . ., p, and belongs to EFDG when j = p + 1, . . ., n. 3.3. Remarks of the reduced-order models The primary comparisons of the three models in frequency domain are conducted by examining the locations of eigenvalues

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014

G Model EPSR-4671; No. of Pages 11

ARTICLE IN PRESS N. Yang et al. / Electric Power Systems Research xxx (2016) xxx–xxx

6 Table 1 The parameters of the investigated DC microgrid.

DG1 DG2 DG3

Power (kW)

Resistance ()

Inductance (mH)

Time constant (s)

1.0 0.5 1.0

7.22* 14.44 7.22

72.2* 144 72.2

0.01* 0.01 0.01

0.25 1.00 1.50

Cable1 Cable2 Cable3 *

0.015 0.050 0.090

Length (m)

50 200 300

can be adjusted.

of state matrices to confirm the validity of the reduced-order models. The eigenvalues of the state matrices of CM, R2M and RMM are shown in Fig. 8, for a typical dc Microgrid with three DGs. The parameters of the examined DC microgrid is given by Table 1. CM has 10 eigenvalues, with the order of 10. The eigenvalues of CM with large real parts disappears in the reduced models. Those with small real parts (dominate poles) are shown in the second graph of Fig. 8. Both the reduced models have similar eigenvalues as that of CM, which proves the effectiveness of the reduced models to be used to represent the original system. Both the reduced-order models can reduce system model complexity significantly. R2M is deduced based on the assumption of similar time scales, thus it may become ineffective when the dynamics of the DGs are largely different. While the proposed RMM combines equivalent DGs of various dynamics together, it can reduce model complexity from the order of 3n + 1 to m + 1 as well as conserve major time scale information. 4. Stability analysis 4.1. Stability test methods The stability of the DC microgrid can be analyzed by either the location of the eigenvalues of the state matrix using linearized

model, or estimated domain of attraction (DoA), or numerical simulation using original nonlinear model. The stability of linear time-invariant (LTI) system depends on locations of the eigenvalues of state matrix. If all eigenvalues have negative real parts, i.e., located at the left half-plane (LHP), the system is asymptotically stable. However it requires the linearization of nonlinear components. The stable operation range of a nonlinear system needs to be examined using the large-signal stability analysis based on the second Lyapunov theorem [16]. Takagi–Sugeno (TS) multi-modeling method [17] can be used to determine the largest estimation of DoA for nonlinear electric systems [18]. A set of linear local models are deduced from the nonlinear system and interconnected by the nonlinear activation functions verifying the property of convex sum. It uses “if-then” rule to represent the input-output linear local relationships. Consider k distinct nonlinearities in the nonlinear model, and each nonlinearity can admit a maximum and a minimum in the studied domain. Then, replacing it by these two extremes, the nonlinear model can be represented by 2k local models, each one under the following LTI form [18]:

⎧ ⎨ d x(t) = A x(t) + B u(t) i i dt ⎩

(33)

y(t) = Hi x(t)

where Ai , Bi and Hi are constant matrices; the subscript i = 1, . . ., 2k . To achieve the nonlinear model, a normalized weight wi (x) is attributed to each local linear model. Then the following convex sum can represent the nonlinear model [18]:

⎧  ⎨ d x(t) = wi (x)(Ai x(t) + Bi u(t)) dt  ⎩ y(t) =

(34)

wi (x)Hi x(t)

The local models can be converted to autonomous models using coordination transformation to move the equilibrium point to the origin. Then the nonlinear system is stable if:



M = MT > 0 ATi M + MAi < 0, ∀i = 1, . . ., 2k

(35)

The existence of a common positive-definite matrix M satisfying the Lyapunov inequality for all the 2k local models is sufficient, not necessary, to guarantee stability of the nonlinear system. In this case, the Lyapunov candidate function is: F(x) = xT Mx

Fig. 8. Eigenvalues of state matrices in CM, R2M and RMM.

(36)

Then the domain in which all the 2k + 1 inequalities in (35) hold is an estimation of DoA. In the examined DC microgrid, only the CPL

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014

G Model

ARTICLE IN PRESS

EPSR-4671; No. of Pages 11

N. Yang et al. / Electric Power Systems Research xxx (2016) xxx–xxx

7

introduces the nonlinearity. Considering R2M with only CPL, the nonlinear model can be expressed by:

⎧ d Rs 1 1 ⎪ ⎪ ⎨ dt I = − I − VL + VN Ld

Ld

Ld

(37)

1 1 PCPL d ⎪ ⎪ I− ⎩ dt VL = C dc C dc VL

To analyze large-signal stability around its equilibrium point (Ie , Ve ), we introduce xI = I − Ie and xV = VL − Ve to move the equilibrium point to the origin. Ie and Ve are the equilibrium source current and DC-bus voltage for a given load, respectively. Then, the nonlinear model of (37) can be rewritten as:

 d dt

xI xV



⎡

−

Rs

⎢ Ld =⎣ 1 C dc

−

1 Ld

⎤

  ⎥ xI ⎦

f (xV )

(38)

xV

where the nonlinearity element, i.e., the last diagonal element in the state matrix, is expressed by: f (xV ) =

PCPL C dc Ve (xV + Ve )

(39)

Then, the boundary of the estimated DoA can been obtained by solving (35) as described in [18]. Numerical simulation can be performed to confirm the above results, and they are easily confirmed by experiments. However, numerical simulation can only test a few points, test of the overall domain is really time-consuming. Fig. 9. Eigenvalue traces when  1 varies from 0.01 s to 0.2 s.

4.2. Sensitivity of time scales To verify whether the stability properties are conserved in the reduced-order models, both small-signal and large-signal stability tests are evaluated by using CM, R2M and RMM. A three DG system is considered to perform the comparison, the time constant of DG1 ( 1 ) varies between 0.01 s and 0.2 s while the others keep constant 0.01 s. The capacitance ratio is chosen 146 ␮F/kW. Other parameters are listed in Table 1. The eigenvalue traces of the state matrices with various  1 in different models are presented in Fig. 9. The eigenvalues are approaching the right half-plane (RHP) with the increase of  1 . The eigenvalues of R2M enter into RHP when  1 is greater than 0.1 s, which indicates the system becomes unstable; While the eigenvalues of RMM are almost identical as those of CM, and the system can operate stably during the whole range when  1 increases from 0.01 s to 0.2 s. Therefore, the adoption of R2M in small-signal stability analysis gives conservative results, while RMM can represent CM pretty well and give more precise results. Estimated DoAs around the equilibrium point are given in Fig. 10 for different models using the multi-modeling method. It can be seen that when  1 goes from 0.01 s to 0.2 s, the DoA of R2M shrinks more quickly than that of CM and RMM with the increase of  1 . This confirms that R2M is more conservative, while RMM gives similar results as CM.

RMM can match those of CM pretty well, which confirms RMM can be adopted to replace CM and to conduct stability analysis. The estimated DoAs with variable capacitance ratio projected on the surface of the current of DG1 and load voltage are presented in Fig. 12. Similarly, the domain of R2M shrinks much faster than that of CM, indicating the results of R2M are conservative. The proposed RMM gives similar domains in different capacitance ratios, which matches CM properly. This confirms the effectiveness of RMM in large-signal stability analysis.

4.3. Sensitivity of capacitances Similar system configuration is adopted to investigate the influence of capacitances on system stability. Eigenvalues traces with various capacitance ratios in different models are shown in Fig. 11. The capacitance ratio decreases from the base value i.e., 146 ␮F with the step −2%, until 10% of the base value. Although R2M presents similar tendency as CM, the eigenvalues enter into RHP much earlier than those of CM. The eigenvalues of

Fig. 10. Estimated domains of attraction when  1 varies from 0.01 s to 0.2 s, in CM (solid line), R2M (line with mark) and RMM (dash line).

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014

G Model EPSR-4671; No. of Pages 11

ARTICLE IN PRESS N. Yang et al. / Electric Power Systems Research xxx (2016) xxx–xxx

8

Fig. 13. Numerical simulation results of CM with CPL step.

Fig. 11. Eigenvalue traces with various the capacitance ratio decreases from 146 ␮F/kW to 14.6 ␮F/kW.

4.4. Numerical simulation To further demonstrate the dynamic differences of the three models, numerical simulations are performed in Matlab. DG2 and DG3 have the time constant 0.01 s; the time constant of DG1 is 20 times of others’, i.e., 0.2 s. The CPL steps from 1 kW to 2 kW at t = 2.0 s. The DC-bus capacitance ratio is selected to be 86 ␮F/kW. The voltage and current responses of the three models are shown in Fig. 13, 14 and 15. The system operates stably in CM and RMM during the CPL step, while the system begins to oscillate in R2M when the load steps.

Fig. 12. Estimated domains of attraction with various capacitance ratios, in CM (solid line), R2M (line with mark) and RMM (dash line).

ESDG in RMM, indicating the slow DG (DG1), shows similar dynamics as that of DG1 in CM; EFDG, representing fast DGs, gives similar dynamics as that of DG2 and DG3 in CM. EFDG takes the high frequency power to maintain the voltage stability while ESDG only takes smooth power. It can be noticed that, ESDG and EFDG share proportional load according to their rated powers during the steady-state, i.e., 40% of the load is taken by ESDG while 60% is taken by EFDG. The average computation times of the three methods are compared in Table 2. The ordinary differential equations are solved by

Fig. 14. Numerical simulation results of R2M with CPL step.

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014

G Model

ARTICLE IN PRESS

EPSR-4671; No. of Pages 11

N. Yang et al. / Electric Power Systems Research xxx (2016) xxx–xxx

9

Fig. 16. Eigenvalue traces in CM and RMM when  2 increases from 0.01 s to 0.1 s, DG2 and DG3 are grouped together.

Fig. 15. Numerical simulation results of the proposed RMM with CPL step. Table 2 Average computation times of the three models. Model

Model order

Computation time ( s)

CM R2M RMM

10 2 3

59.175 57.660 58.072

using classical 4th order Runge-Kutta method, using the Intel CPU i7-2720M @ 2.20 GHz. The average computation time is obtained from the mean value of 150,000 execution cycles. Although only 2.56% and 1.86% of the computation time in CM is reduced in R2M and RMM compared to CM, respectively. The computation time of CM will increase with the number of DGs, while that of RMM depends on the number of groups. The proposed RMM will benefit more when adopted in a DC microgrid with higher number of DGs. Fig. 17. The diagram of the experimental test bench.

4.5. Discussions about the grouping effect The general procedure to construct RMM can be divided into three steps: (1) define the groups according to the time constants of DGs; (2) calculate the parameters of equivalent DGs; (3) construct RMM by connecting the equivalent DGs in parallel. The second and third steps have been presented in Section 3.1 and 3.2; the grouping of DGs is discussed as follows using a three DG microgrid. Considering the previous DC microgrid with three DGs divided into two groups, a slow group (DG1) has the time constant 0.1 s while a fast group (DG2 and DG3) has the time constant 0.01 s. The time constant of DG2  2 (which is grouped with DG3) varies from 0.01 s to 0.1 s. The three eigenvalues near the origin in CM and RMM

are shown in Fig. 16. The center eigenvalues give similar results with the various  2 ; the differences of complex eigenvalues at top and bottom for the two models increase with  2 . When  2 ≥ 4 3 the eigenvalues in RMM enter into RHP. Therefore, an approximated threshold to keep DGs in the same group is around three times. In real application the grouping also needs to consider the balance of accuracy and complexity. 5. Experimental verification A laboratory scale test bench is constructed to verify the proposed RMM and the previous analysis. As shown in Fig. 17, two DC power sources are used to represent DGs and they are

Table 3 The parameters of the microgrid test bench.

DG1 DG2 CPL

Power (W)

Voltage (V)

Current (A)

Time constant (s)

300 300 ∼500

60 60

5 5

0.05 various

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014

G Model EPSR-4671; No. of Pages 11 10

ARTICLE IN PRESS N. Yang et al. / Electric Power Systems Research xxx (2016) xxx–xxx

Fig. 18. Eigenvalue traces of R2M and RMM with  1 increases from 0.05 s to 0.5 s. The arrows show the change tendencies.

connected via DC/DC converters to the common DC-bus. An active load is connected to the DC-bus to represent the CPL. Connecting cables are inserted between the DGs and the load. Due to the limit of the hardware, the DC-bus voltage is set to be 100 V. The control algorithms are realized in dSPACE DS1104. Parameters of the DGs and the active load (referred as CPL) are listed in Table 3. Basic DC-bus capacitance is calculated by (A.1). Considering a time constant of 0.05 s, the DGs’ required output capacitance ratio is 10.55 mF/kW. Thus the required DC-bus capacitance is 12.66 mF. To verify the stability margin, a much aggressive DC-bus capacitance 2.8 mF is adopted. The time constant of DG2  2 keeps constant 0.05 s, while that of DG1  1 varies from 0.05 s to 0.5 s. The eigenvalues traces of R2M and RMM are shown in Fig. 18. The result of RMM indicates that the multi-scale system is stable even when the time constant of DG2 is as high as 10 times of DG1’s. While the R2M indicates that the system would become unstable when DG1’s time constant is 5 times of DG2’s. To compare RMM and R2M from the large-signal stability point of view, two scenarios are considered: a single-time scale microgrid and a multi-time scale microgrid. The estimated DoAs of RMM and R2M with single-time scale ( 1 =  2 = 0.05 s) and multi-time scale conditions ( 1 = 10 2 = 0.5 s) are shown in Fig. 19. Although, R2M gives same DoA as that of RMM for the single-time scale system, the results of R2M for multi-time scale system is too conservative

Fig. 19. Estimated DOAs of RMM and R2M in single and multi-time scales situations.

Fig. 20. Experimental results of DC microgrid with single time scale.

that no stable operation range is found. The RMM indicates the system can operate stably and this will be confirmed by experimental tests. To verify the previous analysis, some experimental tests are conducted. CPL steps from 300 W to 400 W at t = 3 s and steps back at

Fig. 21. Experimental results of DC microgrid with multi-time scale.

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014

G Model EPSR-4671; No. of Pages 11

ARTICLE IN PRESS N. Yang et al. / Electric Power Systems Research xxx (2016) xxx–xxx

t = 7 s. When single time scale system is considered, the DC-bus voltage, the source and load currents are shown in Fig. 20. The system operates stably; DG1 and DG2 show similar dynamic characteristics as the previous presented analysis. When multi-time scales are concerned, DG1’s time constant is 10 times of DG2’s. The experimental results are given in Fig. 21. With the step of CPL, the system still operates stable and DG2 responses quickly to absorb the high frequency term, while DG1 shows a slow dynamics. This confirms that this operation point is in the DoA and proves the effectiveness of the proposed RMM. 6. Conclusion This paper discusses the dynamic modeling of droop controlled multi-time scale DC microgrids. The complexity of the full order dynamic model CM increases with the number of DGs. The traditional reduced order model R2M can reduce the complexity but it cannot maintain similar performance in the multi-scale environment. The proposed RMM, which groups DGs with similar time constants together, can reduce the model complexity as well as keep major time scale information. The analysis using eigenvalues and DoAs have shown that it can represent CM very well even under various time scales, DC-bus capacitances and loads. Numerical simulations and experimental tests in the laboratory scale test bench also have confirmed the effectiveness of the proposed RMM. Appendix A. Basic bus capacitance calculation. √ Assuming that the damping ratio of voltage loop is Nj = 2/2 and the voltage tolerance ı = 5% with the nominal voltage 380 V and the frequency scale ωj = 100 rad/s, then the basic bus capacitance can be calculated by: Cj =

2000PNj 2 ωj ı(1 − ı)VNj

= 2916 × 10−6 PNj

(A.1)

If the DC-bus capacitor is equally distributed into the sources and loads, the output capacitance ratio of DGs is half of the DC-bus capacitance ratio 2916 ␮F/kW, i.e., 1458 ␮F/kW. It is referred as the basic capacitance ratio.

11

References [1] H. Kakigano, M. Nomura, T. Ise, Loss evaluation of DC distribution for residential houses compared with AC system, in: International Power Electronics Conference, Sapporo, 2010, pp. 480–486. [2] S. Luo, Z. Ye, R.-L. Lin, F. Lee, A classification and evaluation of paralleling methods for power supply modules, in: IEEE Power Electronics Specialists Conference, vol. 2, 1999, pp. 901–908. [3] J.M. Guerrero, J.C. Vasquez, J. Matas, L.G.D. Vicu na, M. Castilla, Hierarchical control of droop-controlled AC and DC microgrids a general approach toward standardization, IEEE Trans. Ind. Electron. 58 (1) (2011) 158–172. [4] S. Anand, B. Fernandes, Modified droop controller for paralleling of dc–dc converters in standalone dc system, IET Power Electron. 5 (6) (2012) 782–789. [5] X. Lu, J. Guerrero, K. Sun, J. Vasquez, An improved control method for DC microgrids based on low bandwidth communication with DC bus voltage restoration and enhanced current sharing accuracy, IEEE Trans. Power Electron. 29 (4) (2014) 1800–1812. [6] S. Augustine, M.K. Mishra, N. Lakshminarasamma, Adaptive droop control strategy for load sharing and circulating current minimization in standalone DC low-voltage microgrid, IEEE Trans. Energy Convers. 6 (1) (2015) 132–141. [7] Y. Huang, C. Tse, Circuit theoretic classification of parallel connected DC–DC converters, IEEE Trans. Circuits Syst. I: Reg. Pap. 54 (5) (2007) 1099–1108. [8] D. Chen, L. Xu, Autonomous DC voltage control of a DC microgrid with multiple slack terminals, IEEE Trans. Power Syst. 27 (4) (2012) 1897–1905. [9] N. Yang, D. Paire, F. Gao, A. Miraoui, Distributed control of DC microgrid considering dynamic responses of multiple generation units, in: IEEE IECON, 2014, pp. 1–6. [10] Y. Gu, W. Li, X. He, Frequency-coordinating virtual impedance for autonomous power management of DC microgrid, IEEE Trans. Power Electron. 30 (4) (2015) 2328–2337. [11] G. Xu, L. Xu, D. Morrow, D. Chen, Coordinated DC voltage control of wind turbine with embedded energy storage system, IEEE Trans. Energy Convers. 27 (4) (2012) 1036–1045. [12] P. Karlsson, J. Svensson, DC bus voltage control for a distributed power system, IEEE Trans. Power Electron. 18 (6) (2003) 1405–1412. [13] S. Anand, B.G. Fernandes, Reduced-order model and stability analysis of lowvoltage DC microgrid, IEEE Trans. Ind. Electron. 60 (11) (2013) 5040–5049. [14] E. Prieto-Araujo, F.D. Bianchi, A. Junyent-Ferré, O. Gomis-Bellmunt, Methodology for droop control dynamic analysis of multiterminal VSC-HVDC grids for offshore wind farms, IEEE Trans. Power Deliv. 26 (4) (2011) 2476–2485. [15] A.P. Nobrega Tahim, D.J. Pagano, E. Stramosk, V. Lenz, Modeling and stability analysis of islanded DC microgrids under droop control, IEEE Trans. Power Electron. 30 (8) (2015) 4597–4607. [16] S.F. Marinosson, Stability Analysis of Nonlinear Systems with Linear Programming – A Lyapunov Functions Based Approach, Gerhard-Mercator-Universität Duisburg, 2002 (Ph. D. thesis). [17] T. Takagi, M. Sugeno, Fuzzy Identification of Systems and Its Application to Modeling and Control, 1985, http://dx.doi.org/10.1109/TSMC.1985.6313399. [18] D. Marx, P. Magne, B. Nahid-Mobarakeh, S. Pierfederici, B. Davat, Large signal stability analysis tools in DC power systems with constant power loads and variable power loads – a review, IEEE Trans. Power Electron. 27 (4) (2012) 1773–1787.

Please cite this article in press as: N. Yang, et al., Modeling and stability analysis of multi-time scale DC microgrid, Electr. Power Syst. Res. (2016), http://dx.doi.org/10.1016/j.epsr.2016.04.014