General Dynamic Equivalent Modeling of Microgrid Based on Physical Background

Abstract

Microgrid is a new power system concept consisting of small-scale distributed energy resources; storage devices and loads. It is necessary to employ a simplified model of microgrid in the simulation of a distribution network integrating large-scale microgrids. Based on the detailed model of the components, an equivalent model of microgrid is proposed in this paper. The equivalent model comprises two parts: namely, equivalent machine component and equivalent static component. Equivalent machine component describes the dynamics of synchronous generator, asynchronous wind turbine and induction motor, equivalent static component describes the dynamics of photovoltaic, storage and static load. The trajectory sensitivities of the equivalent model parameters with respect to the output variables are analyzed. The key parameters that play important roles in the dynamics of the output variables of the equivalent model are identified and included in further parameter estimation. Particle Swarm Optimization (PSO) is improved for the parameter estimation of the equivalent model. Simulations are performed in different microgrid operation conditions to evaluate the effectiveness of the equivalent model of microgrid.

1. Introduction

In order to increase the reliability of the electricity supply to the sensitive load, microgrid concept was proposed and developed in recent years [1,2]. Microgrid normally consists of distributed energy resources (DER), energy storage devices and loads. Most of the time, microgrid can be regarded as a self-controlled system that separates and isolates itself from the utility when a severe disturbance nearby occurs, and reconnects itself to the grid automatically when the disturbance is cleared. Obviously, the operational characteristics of the microgrid are quite different from those of the traditional electrical equivalent. Hence, the increasing penetration of the microgrid will have significant impact on the dynamic performances of the distribution network.

To investigate the interactive effect between microgrid and distribution network, a suitable microgrid model is needed. The detailed model of the microgrid comprises dozens of differential equations of all dynamic and static components [2,3,4]. In a simple distribution network with a small number of microgrids, the detailed model of the microgrid is suitable for the dynamic simulation of the distribution network when the microgrid under connected operation mode [5]. However, with increasing penetration of the microgrid into the distribution network, the simulation of a large-scale distribution network becomes very difficult. Under this condition, if an equivalent model of the microgrid is employed, the simulation of the distribution network can be simplified.

Microgrid should operate under connected mode most time to take full advantages of distributed generator. Compared with the distribution network, microgrid can be seen as a controlled load or a controlled electric source under this operation mode. In connected mode, the interactions between loads and distributed generations can be ignored, and the microgrid synthesized dynamic characteristics will be considered in the simulation of distribution network.

Distributed generation is the basis of microgrid; if the equivalent model of the distributed generation is utilized, dynamic simulation of the microgrid can be simplified. An equivalent model compared with the detailed model of the photovoltaic was discussed in [6], and the equivalent model could well describe the dynamic characteristics of the photovoltaic under different faults in power grid. The authors of [7,8] proposed a photovoltaic source dynamic model, the parameters of which were identified based on a least-squares regression-based data processing algorithm. The singular perturbations theory was applied to reduce the model order of the wind farm in [9], and the dynamics of the reduced-order model matched well with those of the detailed model under different operational conditions. Aggregate modeling and detailed modeling for the transient interaction between a large wind farm and a power system were discussed in [10], and the aggregate modeling decreased the simulation time without significantly compromising the accuracy in different conditions. In [11,12], an equivalent method was proposed for integrating wind power generation system in power flow and transient simulation, the unit plants equivalent method and the multiply equivalent method were used for the power flow calculation and transient dynamics simulation, respectively. A probabilistic clustering concept for aggregate modeling of wind farms was proposed in [13], the support vector clustering technique was used to cluster wind turbines based on wind farm layout and incoming wind. Due to the short distances of the electric circuits in the microgrid, there is a strong electromagnetic coupling between the electrical components. These characteristics increase the difficulty in the microgrid analysis. A generalized homology equivalence theory based on differential geometry was used for the microgrid equivalent modeling in [14], and the mathematical analysis of its reduced-order nature was discussed.

Parameter estimation method is a very difficult and challenging task in system modeling. Recently, global optimization techniques such as genetic algorithm [15], evolutionary algorithm [16] and differential evolution [17] have been proposed to solve the parameter estimation problems. Though the genetic algorithm (GA) was employed successfully to solve complex non-linear optimization problems, some deficiencies of GA have been identified in recent research [18]. This degradation in efficiency is apparent when the parameters being optimized are highly correlated and the premature convergence of the GA degrades its performance in terms of reducing the search capability.

Particle swarm optimization (PSO) is an evolutionary computation technique in nature motivated by the simulation of social behaviors. In searching the optimal solution of a problem, information of the best position of each individual particle and the best position among the whole swarm are used to direct the searching. Due to the simple concept, easy implementation and quick convergence, nowadays PSO has gained much attention and wide applications in different fields. Authors of papers [19,20,21,22,23,24] showed that PSO is a feasible approach to parameter estimation of nonlinear systems. In [19], PSO was applied in harmonic estimation. A modified PSO was utilized in the maximum power point tracking for the photovoltaic system in [20]. In the field of parameter estimation, PSO-based parameter estimation technique of proton exchange membrane fuel cell models was proposed in [21], and PSO with quantum was introduced successfully in synchronous generator offline and online parameters estimation problem. Parameter estimation of an induction machine using PSO was shown in [22], and the dynamic PSO and chaos PSO were better than the standard PSO. PSO was used for jointly estimating both the parameters and states of the lateral flow immunoassay model in [23]. Diffusion particle swarm optimization was proposed to optimize the maximum likelihood function in [24], and the PSO technique has been shown to provide a good solution to bearing estimation as it alleviates the effects of multi-modality.

Research has been carried out in the fields of detailed modeling of microgrid and equivalent modeling of the distributed generation. With the increasing penetration of microgrids, the interaction between the microgrid and the distribution network should not be ignored in the power system real-time simulation. However, this will increase the complexity of the simulation with the detailed model of the microgrid components. Hence, a simplified equivalent model of the microgrid is extremely urgent for the simulation analysis of the distribution network. Based on the component detailed models and synthetically dynamic characteristics of the microgrid, an equivalent model of microgrid is proposed in this paper. The proposed equivalent model contains two parts: equivalent static component and equivalent machine component. In order to increase the accuracy of the parameters estimation, trajectory sensitivity is used to identify the key parameters for the further steps of parameter estimation. Particle Swarm Optimization (PSO) is improved and employed to estimate the parameters of the equivalent model. The presented equivalent model and modeling method are shown to be effective by the simulation study on a microgrid connected into distribution network.

2. Microgrid Equivalent Model

A microgrid is made up of a large number of distribution generations, electrical loads and storage devices. Typically, there are two types of components: static components and rotating machines [1]. Static components contain photovoltaic (PV) and static loads. It is common that PV connects to the microgrid through power electronics equipment. Maximum Point Power Tracking (MPPT) and constant power control strategy are applied to the power electronics equipment when the microgrid is connected in the grid-connected operation mode [25]. The output power of these distribution generations is controllable and the dynamic characteristics of them are similar with the static load. In principle, the static load is represented by an exponential of the voltage and frequency [26]. Hence, the output power of the static components can be described by an exponential of the voltage and frequency.

Rotating machine components contain induction motor load, asynchronous induction wind generator and synchronous generator. The structure and the mathematical equations of the asynchronous induction wind generator are similar with those of the induction motor load [27]. The synchronous machine generator and the asynchronous wind generator have similar dynamic characteristics during faults, and the only difference between them is the modeling reference frames. The synchronous machine rotor angular velocity is constant and the velocity voltage is zero in steady-state conditions [28]. However, the rotor angular velocity of synchronous will deviate slightly from the synchronous velocity during a fault since the synchronous machine capability is small in most microgrids. The rotor angular velocity is not equal to the system synchronous velocity and its electrical structure is similar with that of the asynchronous induction wind generator, so the synchronous machine generator can be regarded as an asynchronous machine generator during a fault. Furthermore, the synchronous machine generator, the asynchronous induction wind generator and the induction motor load can be described with a unified mathematical model in the transient dynamic analysis.

As shown in Figure 1, the equivalent model of the microgrid is comprised of an equivalent static component and an equivalent machine component. The equivalent static component is parallel to the equivalent machine component, and they are connected to the distribution network through the Point of Coupling Common (PCC).

2.1. Equivalent Machine Component

The stator and rotor circuits of the equivalent machine component are shown in Figure 2.

Rotor angular velocity ${\text{ω}}_r$ is different from the stator angular velocity ${\text{ω}}_s$. Applying $dq0$ transformation [29], stator voltage equations can be written as

$$\begin{array}{l}{u}_{ds}=p{\text{ψ}}_{ds}-{\text{ω}}_s{\text{ψ}}_{qs}-r_s{i}_{ds}\\ u_{qs}=p{\text{ψ}}_{qs}+{\text{ω}}_s{\text{ψ}}_{ds}-r_s{i}_{qs}\end{array}$$

(1)

Figure 1. The equivalent model of microgrid.

Figure 1. The equivalent model of microgrid.

Figure 2. Equivalent circuit of the equivalent machine component.

Figure 2. Equivalent circuit of the equivalent machine component.

Rotor voltage equations can be written as

$$\begin{array}{l}{u}_{dr}=p{\text{ψ}}_{dr}-{\text{ω}}_{s}s{\text{ψ}}_{qr}+r_{dr}{i}_{dr}\\ u_{qr}=p{\text{ψ}}_{qr}+{\text{ω}}_{s}s{\text{ψ}}_{dr}+r_{qr}{i}_{qr}\end{array}$$

(2)

where $u_{ds}$ and $u_{qs}$ are the stator voltages; ${\text{ψ}}_{ds}$ and ${\text{ψ}}_{qs}$ are the stator flux linkages; $r_s$ is stator resistance; $i_{ds}$ and $i_{qs}$ are stator currents; $u_{dr}$ and $u_{qr}$ are rotor voltages; ${\text{ψ}}_{dr}$ and ${\text{ψ}}_{qr}$ are rotor flux linkages; $r_{dr}$ and $r_{qr}$ are rotor resistance; $i_{dr}$ and $i_{qr}$ are rotor currents; $s=\frac{{\text{ω}}_s-{\text{ω}}_r}{{\text{ω}}_s}$ is the rotor slip; and $p=\frac{d}{dt}$ is the per time derivative.

The stator flux linkage is

$$\begin{array}{l}{\text{ψ}}_{ds}=L_{ds}{i}_{ds}+L_{ad}{i}_{dr}\\ {\text{ψ}}_{qs}=L_{qs}{i}_{qs}+L_{aq}{i}_{qr}\end{array}$$

(3)

The rotor flux linkage is

$$\begin{array}{l}{\text{ψ}}_{dr}=L_{ad}{i}_{ds}+L_{dr}{i}_{dr}\\ {\text{ψ}}_{qr}=L_{aq}{i}_{qs}+L_{qr}{i}_{qr}\end{array}$$

(4)

where $L_{ds}$ and $L_{qs}$ are the stator inductances; $L_{dr}$ and $L_{qr}$ are the rotor inductances; and $L_{ad}$ and $L_{aq}$ are mutual inductances.

Using the definitions below:

$$\begin{array}{l}{L}_d^{\prime }=L_{ds}-\frac{L_{ad}^2}{L_{dr}},L_q^{\prime }=L_{qs}-\frac{L_{aq}^2}{L_{qr}}\\ T_{d0}^{\prime }=\frac{L_{dr}}{r_{dr}},T_{q0}^{\prime }=\frac{L_{qr}}{r_{qr}}\\ E_d^{\prime }=-\frac{L_{aq}}{L_{dr}}{\text{ψ}}_{dr},E_q^{\prime }=\frac{L_{ad}}{L_{qr}}{\text{ψ}}_{qr}\\ E_{dr}=\frac{L_{ad}}{L_{dr}}{u}_{dr},E_{qr}=-\frac{L_{aq}}{L_{qr}}{u}_{qr}\end{array}$$

(5)

rotor voltage equations may be rewritten as follows:

$$\begin{array}{l}\frac{d{E}_d^{\prime }}{dt}=E_{dr}+{\text{ω}}_{s}s\frac{L_{dr}}{L_{ad}}\frac{L_{aq}}{L_{qr}}{E}_q^{\prime }-\frac{1}{T_{q0}^{\prime }}(E_d^{\prime }+(L_{qs}-L_q^{\prime })i_{qs})\\ \frac{d{E}_q^{\prime }}{dt}=E_{qr}-{\text{ω}}_{s}s\frac{L_{qr}}{L_{aq}}\frac{L_{ad}}{L_{dr}}{E}_d^{\prime }-\frac{1}{{T^{\prime }}_{d0}}(E_q^{\prime }-(L_{ds}-L_d^{\prime })i_{ds})\end{array}$$

(6)

When representing power system stability studies, $p{\text{ψ}}_{ds}$ and $p{\text{ψ}}_{qs}$ are neglected in the stator voltage relations. Their neglect corresponds to ignoring the dc component in the stator transient currents, permitting representation of only fundamental frequency components [28]. With the stator transients neglected, stator voltage equations may be rewritten as:

$$\begin{array}{l}{u}_{ds}=-{\text{ω}}_s[(L_{qs}-\frac{L_{aq}^2}{L_{qr}})i_{qs}+\frac{L_{aq}}{L_{qr}}{\text{ψ}}_{qr}]+r_s{i}_{ds}\\ u_{qs}={\text{ω}}_s[(L_{ds}-\frac{L_{ad}^2}{L_{dr}})i_{ds}+\frac{L_{ad}}{L_{dr}}{\text{ψ}}_{dr}]+r_s{i}_{qs}\end{array}$$

(7)

From Equations (5) and (7), we have

$$\begin{array}{l}{u}_{ds}=r_s{i}_{ds}-{\text{ω}}_s{L}_q^{\prime }{i}_{qs}+{\text{ω}}_s{E}_d^{\prime }\\ u_{qs}=r_s{i}_{qs}+{\text{ω}}_s{L}_d^{\prime }{i}_{ds}+{\text{ω}}_s{E}_q^{\prime }\end{array}$$

(8)

The rotor acceleration equation, with time expressed in seconds, is

$$\frac{d{\text{ω}}_r}{dt}=\frac{1}{T_j}(T_M-T_e)$$

(9)

where $T_e$ is the electromagnetic torque, $T_M$ is the mechanical torque, and $T_j$ is the inertia constant of the rotor. Eliminating the rotor currents by expressing them in terms of the stator currents and rotor flux linkages, we find that the per unit electromagnetic torque is

$$T_e={\text{ψ}}_{qr}{i}_{dr}-{\text{ψ}}_{dr}{i}_{qr}=-E_q^{\prime }{i}_{qs}-E_d^{\prime }{i}_{ds}-({\text{ω}}_s{L}_d^{\prime }-{\text{ω}}_s{L}_q^{\prime })i_{ds}{i}_{qs}$$

(10)

The system frequency of the microgrid is constant when the microgrid operation in connected mode. Thus, rotor acceleration equation with $f={\text{ω}}_s=1$ pu can be written as:

$$\frac{d{\text{ω}}_r}{dt}=\frac{1}{T_j}[T_M-(-E_q^{\prime }{i}_{qs}-E_d^{\prime }{i}_{ds}-(L_d^{\prime }-L_q^{\prime })i_{ds}{i}_{qs})]$$

(11)

The output active power and reactive power of the equivalent machine component may be written as:

$$\begin{array}{l}{P}_m=-u_{ds}{i}_{ds}-u_{qs}{i}_{qs}\\ Q_m=u_{qs}{i}_{ds}-u_{ds}{i}_{qs}\end{array}$$

(12)

2.2. Equivalent Static Component

Static load model represents the load characteristics as an algebraic function of the bus voltage magnitude and frequency [26]. Equivalent static component, including static load and distribution generation such as PV in microgrid, is described using algebraic equations. The active and reactive power of the equivalent static component model are related to the system voltage and frequency in the following form:

$$\begin{array}{c}{P}_s=P_{s0}{(U/U_0)}^{p_u}{(f/f_0)}^{p_f}\\ Q_s=Q_{s0}{(U/U_0)}^{q_u}{(f/f_0)}^{q_f}\end{array}$$

(13)

where $P_s$ and $Q_s$ are the active and reactive power of the equivalent static component when the voltage magnitude is $U$ and frequency is $f$, respectively. The subscript $0$ identifies the values of the respective variables at the initial operating condition of PCC. The parameters of this model are the exponents $p_u$, $q_u$, $p_f$ and $q_f$, where $p_u$ is the coefficient of the active power and voltage, $p_f$ is the coefficient of active power and frequency, $q_u$ is the coefficient of reactive power and voltage, and $q_f$ is the coefficient of reactive power and frequency.

The system frequency of the microgrid is constant when the microgrid operates in grid-connected mode. Thus, with $f=1$ pu, the model of the equivalent static component can be written as:

$$\begin{array}{c}{P}_s=P_{s0}{(U/U_0)}^{p_u}\\ Q_s=Q_{s0}{(U/U_0)}^{q_u}\end{array}$$

(14)

2.3. Parameters of the Equivalent Model

From model Equations (6), (8), (11), (12) and (14), it can be seen that the equivalent model parameters include $r_s$, $L_{ds}$, $L_{qs}$, $L_d^{\prime }$, $L_q^{\prime }$, $L_{dr}$, $L_{qr}$, $L_{ad}$, $L_{aq}$, $T_{d0}^{\prime }$, $T_{q0}^{\prime }$, $T_j$, $p_u$, and $q_u$. In order to describe the physical characteristics of the equivalent model, the equivalent model parameters are initialized with corresponding fundamental parameters and will be further estimated in the microgrid modeling. Based on the definitions of the parameters in Equation (5), the fundamental parameters of the equivalent machine component are $r_s$, $L_{sl}$, $L_{ad}$, $L_{aq}$, $r_{dr}$, $L_{drl}$, $r_{qr}$, $L_{qrl}$, $T_j$, $p_u$, and $q_u$, where $L_{sl}$ is the stator leakage inductance; $L_{ad}$ and $L_{aq}$ are mutual inductances of $d$ and $q$ axis; $r_{dr}$ and $r_{qr}$ are rotor resistances of $d$ and $q$ axis; and $L_{drl}$ and $L_{qrl}$ are rotor leakage inductances of $d$ and $q$ axis.

There are other two important parameters, namely $s_0$ and $K_{mp}$. Where $s_0$ is the initial slip of the equivalent machine component, and $s_0$ presents the type of equivalent machine. If $s_0>0$, the equivalent machine component absorbs power from the distribution network, and has the characteristics of induction motor load. Oppositely, $s_0<0$ means that the equivalent machine component injects power into the distribution network, and has the characteristics of asynchronous generator. $K_{mp}$ is the fraction of the equivalent machine component active power with respect to the total initial active power $P$. The active power flow between the microgrid and the distribution network is bidirectional. Thus, $P>0$ indicates that the microgrid absorbs power from the distribution network, and $P<0$ indicates that the microgrid injects power into the distribution network. The similar definition is also applied to reactive power Q.

As a result, the 13 parameters, namely $K_{mp}$, $s_0$, $r_s$, $L_{sl}$, $L_{ad}$, $L_{aq}$, $r_{dr}$, $L_{drl}$, $r_{qr}$, $L_{qrl}$, $T_j$, $p_u$, and $q_u$, in the equivalent microgrid model need to be estimated in the microgrid equivalent dynamic modeling.

3. Equivalent Model Parameter Sensitivity Analysis

The number of the parameters to be estimated has significant impact on the accuracy of the parameters estimation. Parameter sensitivity analysis is an efficient method to determine the key parameters of the equivalent model. Parameter trajectory sensitivity is defined as:

$$\frac{\partial y(\text{θ},k)}{\partial {\text{θ}}_j}=\underset{\text{Δ}{\text{θ}}_j\to 0}{\text{lim}}\frac{y({\text{θ}}_1,\cdots ,{\text{θ}}_j+\text{Δ}{\text{θ}}_j,\cdots ,{\text{θ}}_m,k)-y({\text{θ}}_1,\cdots ,{\text{θ}}_j,\cdots ,{\text{θ}}_m,k)}{\text{Δ}{\text{θ}}_j}$$

(15)

where $y(\text{θ},k)$ is the time domain trajectory of the output variable, $\text{θ}$ is the vector of the parameters of the equivalent model, ${\text{θ}}_j$ is the jth parameter of the equivalent model, $m$ is the number of the parameters, and $k$ is the sampling sequence.

In order to improve the accuracy of parameter trajectory sensitivity, median method is used to calculate the trajectory sensitivity when $\text{Δ}{\text{θ}}_j$ is small enough, which is shown as:

$$\frac{\partial \left[y(\text{θ},k)/y_0\right]}{\partial \left[{\text{θ}}_j/{\text{θ}}_{j0}\right]}=\frac{\left[y({\text{θ}}_1,\cdots ,{\text{θ}}_j+\text{Δ}{\text{θ}}_j,\cdots ,{\text{θ}}_m,k)-y({\text{θ}}_1,\cdots ,{\text{θ}}_j-\text{Δ}{\text{θ}}_j,\cdots ,{\text{θ}}_m,k)\right]/y_0}{2\text{Δ}{\text{θ}}_j/{\text{θ}}_{j0}}$$

(16)

where ${\text{θ}}_{j0}$ is the initial value of the parameter ${\text{θ}}_j$, $\text{Δ}{\text{θ}}_j$ is the variation of ${\text{θ}}_j$, and $y_0$ is the initial value of the output variable in steady-state.

3.1. Trajectory Sensitivity

Time domain parameter trajectory sensitivity curve can describe the behaviors of the output variable. For convenience of comparison, the average sensitivity can be calculated as

$$A_j=\frac{1}{N}{\displaystyle \sum _{k=1}^N\left|\frac{\partial \left[y(\text{θ},k)/y_0\right]}{\partial \left[{\text{θ}}_j/{\text{θ}}_{j0}\right]}\right|}$$

(17)

where $A_j$ is the average of the j-th parameter with respect to the trajectory of the output variable $N$ is the number of sample points, $y_0$ is the initial value of the output variable.

Trajectory sensitivity demonstrates the impact of the variation of the parameter on that of the output variable’s trajectory. If the trajectory sensitivity of a parameter is larger than that of the other parameters, the parameter plays a more important role on the dynamics of the output variable; in other words, the parameter can be estimated easily by using the dynamics of the output variable. In contrast, if the trajectory sensitivity of a parameter is very small, e.g., even close to zero, it is difficult to estimate the parameter using the dynamics of the output variable.

3.2. Trajectory Sensitivity Phase

If a couple of parameters have an unknown relationship between each other, they are dependent on each other and unidentifiable as well. However, these unidentifiable parameters can be identified using trajectory sensitivity analysis [29].

Assuming that the parameters ${\text{θ}}_i$ and ${\text{θ}}_{i+1}$ are coupling with each other, the output of the power system can be written as

$$y=f[{\text{θ}}_1,{\text{θ}}_2,\cdots ,\text{φ}({\text{θ}}_i,{\text{θ}}_{i+1}),\cdots ,{\text{θ}}_n]$$

(18)

The sensitivities of the parameters ${\text{θ}}_i$ and ${\text{θ}}_{i+1}$ can be analyzed as [27]

$$\{\begin{array}{l}\frac{\partial y}{\partial {\text{θ}}_i}=\frac{dy}{\partial \text{φ}}\frac{\partial \text{φ}}{\partial {\text{θ}}_i}\\ \frac{\partial y}{\partial {\text{θ}}_{i+1}}=\frac{dy}{\partial \text{φ}}\frac{\partial \text{φ}}{\partial {\text{θ}}_{i+1}}\end{array}$$

(19)

Then

$$\frac{\partial y}{\partial {\text{θ}}_{i+1}}=\frac{dy}{\partial {\text{θ}}_i}[\frac{\partial \text{φ}/\partial {\text{θ}}_{i+1}}{\partial \text{φ}/\partial {\text{θ}}_i}]$$

(20)

It should be pointed out that $\partial y/\partial {\text{θ}}_{i+1}$ and $\partial y/\partial {\text{θ}}_i$ vary with time, while $\partial \text{φ}/\partial {\text{θ}}_{i+1}$ and $\partial \text{φ}/\partial {\text{θ}}_i$ are constant [29,30]. Hence, $\partial y/\partial {\text{θ}}_{i+1}$ and $\partial y/\partial {\text{θ}}_i$ reach zero at the same time. In other words, the trajectory sensitivities of these two couple unidentifiable parameters are either in phase or in reverse with each other. The trajectory sensitivity of these two parameters will pass zero at the same time. Oppositely, if the output variable curves of the two parameters do not pass zero at the same time approximately, they are independent and can be identified using the dynamics of output variables of the system [31].

4. Microgrid Parameter Estimation

4.1. Rotor Voltage Equivalence

Rotor voltage is an important operational variable in synchronous machine, and the parameters of the excitation system are estimated individually [32]. In the microgrid equivalent model, the excitation system is a virtual polymerization system, which maintains the terminal voltage. In [33], a synchronous machine with an excitation system is described by an asynchronous machine and a constant current load in load modeling. The dynamic of the synchronous machine is similar with the dynamic of the asynchronous machine and the constant current load. Hence, the rotor voltage of the equivalent microgrid model is replaced by a constant current load, which can be seen as a part of the equivalent static component in the modeling of microgrid.

4.2. Parameter Estimation Based on PSO

4.2.1. Review of Particle Swarm Optimization

Particle Swarm Optimization (PSO) is a heuristic optimization algorithm. It was first introduced by Kennedy and Eberhart, based on the observations of social behaviors of animals, such as bird flocking, fish schooling, and swarm theory [34]. This algorithm implements a global method that performs a search of parameters over a specified problem space. Like other evolutionary algorithms, PSO performs using a population of individuals that are updated iteratively. Swarm members communicate good position with each other and dynamically adjust their own position and velocity. Velocity adjustment is based upon the historical behaviors of the particles themselves as well as their neighbors. In each iteration, the velocity and position of each particle are updated according to the following equations.

$$v_i(k)=w{v}_i(k-1)+c_{1}ran{d}_1[pbes{t}_i-x_i(k-1)]+c_{2}ran{d}_2[gbest-x_i(k-1)]$$

(21)

$$x_i(k)=v_i(k)+x_i(k-1)$$

(22)

where $x_i$ is the position of the $i$th particle, $v_i$ is the velocity of the $i$-th particle, $w$ is the inertia weight which decreases linearly determined by Equation (23), $pbes{t}_i$ is the best position of the $i$-th particle, $gbest$ is the global optimal position of the current swarm with the best objective value, $c_1$ and $c_2$ are acceleration constants, and $ran{d}_1$ and $ran{d}_2$ are two independent random numbers uniformly distributed over [0 1].

$$w=w_{\text{min}}+\frac{w_{\text{max}}-w_{\text{min}}}{k_{\text{max}}}\times k$$

(23)

The algorithm achieves the optimal solution by two types of search memory, “cognitive” component and “social” component, which are shown as the second and third parts of Equation (21). The cognitive component makes the particle move toward its own best positions. The social component makes the particle toward the best position found by its neighbors, which means the collaborative behavior of particles.

4.2.2. Particle Swarm Optimization (PSO) with Chaos Neighborhood Searching

As optimization problem becomes more and more complex, the random characteristics of PSO will reduce the convergence velocity. In this paper, chaotic mutation theory is used to improve the global convergence velocity of the standard PSO. Chaos is a common phenomenon in non-linear systems that include infinite unstable period motions. Chaos-based neighborhood searching method with multiple different neighborhoods is designed and incorporated to enrich the searching behaviors, so as to avoid premature convergence [35]. In addition, an effective nonlinear adaptive inertia weight is employed to further enhance the exploitation ability after the chaos-based neighborhood searching.

Chaotic mapping is a discrete-time dynamical system, in which the chaotic sequences are considered as sources of random sequences, which can avoid getting stuck in a local optimum during the search process and overcomes the premature convergence phenomenon. Logistic mapping is a common chaos mapping, which is employed in the global optimal position neighborhood searching with the following steps.

(1)

Generating the initial variable $u_{0j}$ randomly based on the global optimal position of the particles.

(2)

Getting the chaotic sequences $u_{1j}$ using the logistic mapping $u_{1j}=4u_{0j}(1-u_{0j})$.

(3)

Generating a local neighborhood mutation variable $\text{Δ}{x}_j$: $\text{Δ}{x}_j=-\text{β}+2\text{β}{u}_{1j}$, where $\text{β}$ is the radius of the local neighborhood and is updated by Equation (24).

$$\text{β}=(x_{j\text{max}}-x_{j\text{min}})\text{cos}(\frac{\text{π}(t-1)}{2(t_{\text{max}}-1)})$$

(24)

where, $t$ and $t_{\text{max}}$ are the round of the iteration and the maximal iteration times of the neighborhood searching, respectively. $x_{j\max }$ and $x_{j\min }$ are the upper and lower bounds of the variable $x_j$, respectively.

(4)

Local neighborhood searching: A temporary global optimal position is defined as $gbes{t}^{\prime }=gbest+\text{Δ}X$, where $\text{Δ}X$ is the chaos mutation variables vector $[\text{Δ}{x}_1\text{\hspace{0.17em}Δ}{x}_j\cdots \text{Δ}{x}_N]$ and $gbest=[x_{g1}\cdots x_{gN}]$ is the current global optimal position. Comparing $gbes{t}^{\prime }$ with $gbest$, the current global optimal position is updated by the larger one.

4.2.3. Improved Particle Swarm Optimization(IPSO) Based Parameter Estimation

IPSO algorithm is employed to search the optimal parameters of the microgrid equivalent model to achieve an optimal matching between the detailed model and the equivalent model of the microgrid. The steps of the IPSO-based microgrid equivalent model parameter estimation are shown in Figure 3.

Figure 3. Chart of Improved Particle Swarm Optimization (IPSO)-based parameter estimation.

Figure 3. Chart of Improved Particle Swarm Optimization (IPSO)-based parameter estimation.

(1) Initialization: From Section 2, there are 13 parameters, $K_{mp}$, $s_0$, $r_s$, $L_{sl}$, $L_{ad}$, $L_{aq}$, $r_{dr}$, $L_{drl}$, $r_{qr}$, $L_{qrl}$, $T_j$, $p_u$ and $q_u$ for the optimal parameter of the microgrid model. The initial particles may be generated randomly with a specified upper and lower bounds, $X_{\text{max}}$ and $X_{\text{min}}$, respectively. The dimension of each particle is equal to the number of the parameters. Measured data include voltages ($u_x+j{u}_y$) and powers ($P+jQ$) of PCC.

The initial velocity of particles are generated randomly between the upper and lower bounds, $v_{i,\text{max}}$ and $v_{i,\text{min}}$. Where the upper and lower bounds is defined by

$$\begin{array}{l}{v}_{i,\text{max}}=\frac{x_{i,\text{max}}-x_{i,\text{min}}}{N}\\ v_{i,\text{min}}=-v_{i,\text{max}}\end{array}$$

(25)

where $x_{i,\text{max}}$ and $x_{i,\text{min}}$ are the upper and lower bounds of particles; and $v_{i,\text{max}}$ and $v_{i,\text{min}}$ are the maximum and minimum velocities of particles. $N$ is the interval of dimension, which is normally chosen to be between 5 and 10.

(2) Evaluation: The Objective function is employed to evaluate each particle, and it is usually defined as the error between the measured power and the output power of the microgrid equivalent model as follows:

$$E(\text{θ})=\text{min}{\displaystyle \sum _{k=1}^N\{{[P(k)-P_M(k,\text{θ})]}^2+{[Q(k)-Q_M(k,\text{θ})]}^2\}}$$

(26)

where $\text{θ}$ is the parameter vector of the equivalent model, $P(k)$ and $Q(k)$ are the measured power of sampling time $k$, and $P_M(k,\text{θ})$ and $Q_M(k,\text{θ})$ are the calculated power of the microgrid equivalent model. The aim of the objective function is that the output power of equivalent model can match the measured power.

(3) Updating movement velocities of the particles: Movement velocity updating is an important step in the process of evolution. The movement velocity of each particle is updated by Equation (21), and the velocity is limited in the upper and lower bounds as follows:

$$\begin{array}{l}if{v}_{ij}(t+1)>v_{j,\max }then{v}_{ij}(t+1)=v_{j,\max }\\ if{v}_{ij}(t+1)<v_{j,\min }then{v}_{ij}(t+1)=v_{j,\min }\end{array}$$

(27)

(4) Updating positions of particles: The position of each particle is updated by Equation (22). The position of each particle is limited as follows

$$\begin{array}{l}if{x}_{ij}(t+1)>x_{j,\max }then{x}_{ij}(t+1)=x_{j,\max }\\ if{x}_{ij}(t+1)<x_{j,\min }then{x}_{ij}(t+1)=x_{j,\min }\end{array}$$

(28)

(5) Chaos neighborhood searching: Following the steps described in the previous Section 4.2.2, the global optimal particle position is updated.

Firstly, an initial variable $u_{0j}$ is generated randomly, and the chaotic sequences $u_{1j}$ is generated by using the formula $u_{1j}=4u_{0j}(1-u_{0j})$. Secondly, a local neighborhood mutation variable $\text{Δ}{x}_j$ is generated by using the formula $\text{Δ}{x}_j=-\text{β}+2\text{β}{u}_{1j}$. Lastly, a new temporary global optimal position is generated which is the neighbor of the global optimal position. The fitness of the temporary global optimal position is compared with that of the current global optimal position, and the current global optimal position is replaced by the better one.

(6) Termination: The parameter estimation algorithm will be terminated if the iteration times exceeds the specified maximal iteration times or the fitness value of the global best is smaller than a given value.

5. Microgrid Equivalent Modeling and Discussion

5.1. Microgrid System

In this work, a microgrid system is built in DIgSILENT PowerFactory to test the effectiveness of the microgrid equivalent model and the feasibility of the modeling method. The structure of the microgrid is shown in Figure 4, and the detailed parameters and the operation mode can be found in [36]. The microgrid consists of wind generation, PV generation, micro-gas turbine, asynchronous induction motor and static load. The detailed parameters of microgrid are shown in Appendix A. The microgrid is radial and is connected to the distribution network through a 10 kv/20 kv transformer. The distribution network is a standard Benchmark medium voltage distribution network model in [37].

Figure 4. Microgrid system.

Figure 4. Microgrid system.

A single-phase short circuit is applied in the distribution network, and is cleared in 0.06 s. Using the dynamics of the exchange power between the microgrid and the distribution network, all the parameters of the equivalent model are estimated by the proposed PSO. The values of parameters used in PSO are depicted in Table 1. The algorithm has been implemented in C#.NET. The program runs on a 1.8 GHz, Intel Core2 Duo, with 2 GB RAM PC. The detail of the key parameters selected are: $w_{\text{max}}=0.9$, $w_{\text{min}}=0.4$, $c_1=c_2=2.0$, $t_{\text{max}}=10$, $k_{\text{max}}=500$. The bound of estimated parameters are shown in Appendix B. The estimated parameters are listed in Table 1, and the dynamics of the output power of the detailed model and the equivalent model are illustrated in Figure 5. It can be seen that the dynamics of the output power of the equivalent model is approximate to that of the detailed model.

Figure 5. (a) The active power dynamics of the detailed model and the equivalent model of the Microgrid; (b) The reactive power dynamics of the detailed model and the equivalent model of the Microgrid.

Figure 5. (a) The active power dynamics of the detailed model and the equivalent model of the Microgrid; (b) The reactive power dynamics of the detailed model and the equivalent model of the Microgrid.

Table 1. A group of feasible parameters for the equivalent model.

Parameter	Kmp	s0	rs	Lsl	Lad	Lap	rdr
Value	−8.70	−0.206	0.031	0.340	0.848	3.14	0.00023
Parameter	$L_{drl}$	$r_{qr}$	$L_{qrl}$	$T_j$	$p_u$	$q_u$	Error
Value	2.18	0.171	0.756	9.66	7.06	0.731	2.401

Table 1. A group of feasible parameters for the equivalent model.

Parameter	Kmp	s0	rs	Lsl	Lad	Lap	rdr
Value	−8.70	−0.206	0.031	0.340	0.848	3.14	0.00023
Parameter	$L_{drl}$	$r_{qr}$	$L_{qrl}$	$T_j$	$p_u$	$q_u$	Error
Value	2.18	0.171	0.756	9.66	7.06	0.731	2.401

The evolution of function objective is presented along the iterations in order to observe the behavior of improved PSO, and to verify if the PSO works as expected. Figure 6 shows the evolution of function objective along with the iterations. It can be seen from the graph that the objective function has converged to a constant value of 2.401 after 57 iterations. It is quite apparent that the improved PSO has the faster convergence than the standard PSO algorithm. Moreover, the proposed IPSO and standard PSO achieved in the same machine converged their corresponding optimal solutions in 25.3 s and 37.5 s, respectively.

Figure 6. Convergence of the improved particle swarm optimization (IPSO) and PSO Algorithms.

Figure 6. Convergence of the improved particle swarm optimization (IPSO) and PSO Algorithms.

5.2. Parameter Sensitivity Analysis and Error Analysis

The number of parameters to be estimated has a significant impact on the accuracy of the parameter estimation. If only the key parameters are estimated, the convergence of the PSO can be improved and the probability finding the true value of the parameters will be higher.

The sensitivities of the equivalent microgrid model parameters, $\text{θ}=[K_{mp},s_0,r_s,L_{sl},L_{ad},L_{aq},r_{dr},L_{drl},r_{qr},L_{qrl},T_j,p_u,q_u]$, with respect to the trajectory of the apparent power $S$, are calculated. $S$ is defined as:

$$S=\sqrt{P_M^2+Q_M^2}$$

(29)

where $P_M$ and $Q_M$ are the calculated active and reactive power based on the dynamic equivalent model of equal Equation (26).

The average trajectory sensitivities of the parameters are listed in Table 2, and the dynamics of the trajectory sensitivities are shown in Figure 7. It can be seen in Table 2 that the parameters $K_{mp}$, $s_0$, $L_{ad}$, $L_{aq}$, $r_{dr}$, $L_{drl}$, $L_{qrl}$ and $p_u$ have more significant impacts on the trajectory than the other parameters. Figure 7 shows that the phases of the sensitivity of $K_{mp}$, is almost reverse to that of $r_{dr}$, which means that $K_{mp}$ and $r_{dr}$ are related with each other and are unidentifiable by using the dynamics of $S$. As a result, only the seven parameters, i.e., $K_{mp}$, $s_0$, $L_{ad}$, $L_{aq}$, $L_{drl}$, $L_{qrl}$ and $p_u$, are included in the further steps of the parameter estimation, while the other parameters are set as default values.

Table 2. Parameter sensitivity absolute value.

Parameter	Kmp	s0	rs	Lsl	Lad	Lap	rdr
Sensitivity	0.0160	0.0352	0.0016	0.0040	0.1367	0.0377	0.0149
Parameter	$L_{drl}$	$r_{qr}$	$L_{qrl}$	$T_j$	$p_u$	$q_u$
Sensitivity	0.1543	0.00043	0.0583	0.0015	0.0120	0.0011

Table 2. Parameter sensitivity absolute value.

Parameter	Kmp	s0	rs	Lsl	Lad	Lap	rdr
Sensitivity	0.0160	0.0352	0.0016	0.0040	0.1367	0.0377	0.0149
Parameter	$L_{drl}$	$r_{qr}$	$L_{qrl}$	$T_j$	$p_u$	$q_u$
Sensitivity	0.1543	0.00043	0.0583	0.0015	0.0120	0.0011

Figure 7. The trajectory sensitivity of parameter $K_{mp}$ and $r_{dr}$.

Figure 7. The trajectory sensitivity of parameter $K_{mp}$ and $r_{dr}$.

5.3. Key Parameters Estimation and Error Analysis

By using the dynamics of the output power of the detailed model of microgrid, the key parameters are estimated and listed in Table 3. The dynamics of the output active power and reactive power are also illustrated in Figure 8a,b, respectively. As can be seen in Figure 8, the amplified windows show the differences of the two estimated method, and it is clear that the error between the dynamics of the output power of the equivalent model and the detailed model is smaller than that when all the parameters are estimated simultaneously, as shown in Table 1.

The simulation results reveal that the proposed dynamic equivalent model is accurate enough to describe the dynamic characteristics of the microgrid detailed model quickly during the simulation period. The equivalent model can track the fluctuation of active power and reactive power quickly and flexibly, especially in the moment of fault cutting. The estimated results of parameters are suitable according to the electrical parameter standards.

Figure 8. (a) The active power dynamics of microgrid; (b) The reactive power dynamics of microgrid.

Figure 8. (a) The active power dynamics of microgrid; (b) The reactive power dynamics of microgrid.

Table 3. Parameter estimated value.

Parameter	Kmp	s0	Lad	Lap	Ldrl	Lqrl	pu	Error
Value	−8.637	−0.2161	0.8042	3.417	2.120	0.7008	7.039	2.389

Table 3. Parameter estimated value.

Parameter	Kmp	s0	Lad	Lap	Ldrl	Lqrl	pu	Error
Value	−8.637	−0.2161	0.8042	3.417	2.120	0.7008	7.039	2.389

Because the proposed model is an equivalent model in electromechanical transient, the stator transient is neglected during the process of equivalent modeling, which will bring some error both in parameters and dynamic of the microgrid. In Figure 8a, it can be seen that at the time of the single-phase short fault occurs and cleared, the error between equivalent model and the detailed model is larger than other time relatively. At this moment, the active power of equivalent model cannot track the dynamic of microgrid completely. The components of microgrid are various and the dynamic of the output power of microgrid contains different components. In order to improve the representational capacity of the equivalent model, some elements are ignored, such as the stator transient of wind turbine and micro-gas turbine, the transient of PV inverter.

5.4. Comparing to Other Models

The proposed microgrid equivalent model is also compared with the equivalent machine component model and black-box model. The equivalent microgrid model contains two parts: equivalent machine component and equivalent static component. The equivalent machine component is used to describe the rotor characteristics of electrical equipment, and the equivalent static component is used to describe the characteristics of steady state of microgrid.

The equivalent machine model is compared with the comprehensive model in this paper. This means that the values of the parameters are $K_{mp}=1$, $p_u=0$ and $q_u=0$. Figure 9 shows the dynamics between different models of the microgrid. It can be seen that the error between the proposed equivalent microgrid model and detailed model is smaller than that of the equivalent machine model and detailed model. As can be seen in Figure 9b, the output reactive power of the proposed equivalent microgrid model can track the detailed model more accurately than that of the equivalent machine component model during the fault. It means that the equivalent static component is important in the proposed equivalent microgrid model. The equivalent static component can describe the dynamic characteristics of static component in the microgrid, such as: photovoltaic system, static load, and other static equipment.

Figure 9. (a) Dynamics of Active Power of different models; (b) Dynamics of Reactive Power of different models.

Figure 9. (a) Dynamics of Active Power of different models; (b) Dynamics of Reactive Power of different models.

The proposed microgrid equivalent model is based on the physics characteristics of the electrical components in microgrid, and can describe the electrical physical characteristics of microgrid components. Non-mechanism model is also used in the dynamic equivalent of power system. In the equivalent modeling of microgrid, a dynamic equivalent black-box model based on prony analysis is presented in [38], the proposed equivalent model is compared with the black-box equivalent model. Figure 10 shows the comparison the proposed model and the black-box model; it can be seen that the two models are in good agreement with the detailed model, and the proposed model is more accurate.

Figure 10. (a) Dynamics of Active Power of different models; (b) Dynamics of Reactive Power of different models.

Figure 10. (a) Dynamics of Active Power of different models; (b) Dynamics of Reactive Power of different models.

5.5. Parameters Estimation under Different Operational Conditions

The operational condition of the microgrid changes with the fluctuation of the distributed generation output power and the variation of loads, which leads to the power exchange between the microgrid and the distribution network changes accordingly. Therefore, the parameters of the microgrid equivalent model should be adjusted to adapt the changes of the microgrid operational condition. Three different typical operational conditions are studied in the paper, the detailed operational parameters are shown in Table 4 and the parameters estimated results are shown in Table 5.

Table 4. Describe of operational conditions.

Operational Condition	Description
Condition A	The output power of gas turbine reduces 40 percent with the output active power is 3 MW.
Condition B	The output power of gas turbine reduces 40 percent with the output active power is 3 MW, 40 percent of Static Load 2, which is 2 MW active power and 0.4 Mvar reactive power, is removed.
Condition C	The output power of gas turbine reduces 40 percent with the active power is 3 MW. The output active power of Wind Generation reduces a half, which become 0.75 MW. PV is removed. 40 percent of Static Load 1, which is 0.4 MW and 0.2 Mvar, is removed.

Table 4. Describe of operational conditions.

Operational Condition	Description
Condition A	The output power of gas turbine reduces 40 percent with the output active power is 3 MW.
Condition B	The output power of gas turbine reduces 40 percent with the output active power is 3 MW, 40 percent of Static Load 2, which is 2 MW active power and 0.4 Mvar reactive power, is removed.
Condition C	The output power of gas turbine reduces 40 percent with the active power is 3 MW. The output active power of Wind Generation reduces a half, which become 0.75 MW. PV is removed. 40 percent of Static Load 1, which is 0.4 MW and 0.2 Mvar, is removed.

As can be seen in Table 5, the equivalent machine component outputs active power as a generator with $K_{mp}<0$ and $s_0<0$, while the microgrid absorbs power from the distribution network in condition A. This means that the consume power of equivalent static component comes from distribution network and equivalent machine component. In conditions B and C, the equivalent machine component absorbs power as an asynchronous induction motor with $K_{mp}>0$ and $s_0>0$, while the microgrid absorbs power from the distribution network.

Table 5. Parameter value of different operational conditions.

Parameter	Kmp	s0	rs	Lsl	Lad	Laq	rdr	Ldrl	rqr	Lqrl	Tj	pu	qu
Condition A	−4.69	−0.17	0.0590	0.36	1.94	2.50	0.0003	3.04	0.2560	2.24	5.62	4.87	0.94
Condition B	1.64	0.07	0.0610	0.10	4.05	7.54	0.1186	1.49	0.0024	4.00	5.18	8.95	0.61
Condition C	0.99	0.12	0.0492	0.31	3.70	4.76	0.1450	0.65	0.0057	4.00	6.55	3.58	0.70

Table 5. Parameter value of different operational conditions.

Parameter	Kmp	s0	rs	Lsl	Lad	Laq	rdr	Ldrl	rqr	Lqrl	Tj	pu	qu
Condition A	−4.69	−0.17	0.0590	0.36	1.94	2.50	0.0003	3.04	0.2560	2.24	5.62	4.87	0.94
Condition B	1.64	0.07	0.0610	0.10	4.05	7.54	0.1186	1.49	0.0024	4.00	5.18	8.95	0.61
Condition C	0.99	0.12	0.0492	0.31	3.70	4.76	0.1450	0.65	0.0057	4.00	6.55	3.58	0.70

In condition A, with the power reduction of gas turbine, microgrid may absorb more power from distribution network. Figure 11 shows microgrid response to the single-phase short fault in distribution network. The exchange power between microgrid and distribution network fluctuates violently, but the results show the great agreement between proposed equivalent model and the detailed model, and the reactive power of equivalent model can track the detailed model with a small time delay.

Figure 12 and Figure 13 show the response of dynamic equivalent model during the fault of condition B and condition C. The microgrid construction of condition B and C is similar, and it can be seen that the dynamics of the equivalent model coincide with that of the detailed model well. Especially during the fault, the errors between equivalent model and detailed model are very small.

Furthermore, Table 5 shows the stability of estimated parameters of equivalent model in different microgrid operational conditions, and the equivalent model can descript the dynamic of microgrid under different operational condition, which means that the equivalent model is robust.

Figure 11. (a) The active power dynamics of the detailed model and the equivalent model under condition A; (b) The reactive power dynamics of the detailed model and the equivalent model under condition A.

Figure 11. (a) The active power dynamics of the detailed model and the equivalent model under condition A; (b) The reactive power dynamics of the detailed model and the equivalent model under condition A.

Figure 12. The active power dynamics of the detailed model and the equivalent model under condition B.

Figure 12. The active power dynamics of the detailed model and the equivalent model under condition B.

Figure 13. The active power dynamics of the detailed model and the equivalent model under condition C.

Figure 13. The active power dynamics of the detailed model and the equivalent model under condition C.

6. Conclusions

The equivalent model of the microgrid, based on electrical equipment physical characteristics, is proposed in this paper. The equivalent model consists of an equivalent machine component and an equivalent static component. Equivalent machine component is used to describe the dynamic of rotor machine characteristics and equivalent static component is used to describe the comprehensive characteristics of static component. In order to clear up the importance of each parameter toward the dynamics of the model, time-domain sensitivity of the parameters of the equivalent model with respect to the output power of the microgrid has been analyzed. The key parameters with important impact on the dynamics of the equivalent model have been detected. An improved particle with variable neighborhood searching is proposed where the radius of neighborhood is changed based on chaos iteration. The improved particle swarm optimization (IPSO) algorithm is used for parameters estimation, which improves the accuracy of the parameter estimation. The simulation results of microgrid under different operation conditions show that the simplified equivalent model is in good agreement with the detailed model under fault. The proposed model can be combined with other electrical component model in the power system simulation platform. With the proposed equivalent model, we can simplify the detailed model of microgrid obviously, and improved the simulation speed of power system.