Dynamic Analysis and Model Order Reduction of Virtual Synchronous Machine Based Microgrid

The concept of virtual synchronous machine (VSM) was proposed to deal with the shortcomings of low inertia and damping of traditional control strategies for power electronic converters. But what if all distributed energy resources and controllable loads in a microgrid adopt the VSM control strategy, and will it present better performance than conventional droop control-based microgrid (DMG)? In this paper, the VSM-based microgrid (VSMG) is analyzed. The small-signal modeling of the VSMG is studied at first. Then static stability and dynamic characteristics of the VSMG are analyzed and compared with the DMG in both frequency-domain and time-domain. With the growing scale of microgrids, their modeling and simulation are becoming significant computational burdens. Inspired by the participation factor analysis of the VSMG and the concept of coherency in power systems, the VSMG small-signal model is equivalent to a modified third-order synchronous generator (SG) model in this paper. The equivalencing involves gray-box system identification and is realized by estimating equivalent electrical parameters alternately and iteratively. The equivalent SG (EqSG) model is compared with three representative model order reduction methods to verify its effectiveness. Simulation results confirm the accuracy of the EqSG model substituting detailed VSMG model in time-domain simulations.


I. INTRODUCTION
With high penetration of distributed energy resources (DERs) and mass access of controllable loads, power systems are experiencing a paradigm shift from centralized and rotational generator-dominated systems to distributed and inverter-dominated systems [1]. In the modern power system, DERs and controllable loads are often integrated into a microgrid, and the hierarchical structure of microgrid -microgrid groups -active distribution network (ADN) is constructed. Microgrids not only promote the integration of renewable resources, but also improve the control flexibility of power systems [2]. However, the extensive introduction of power electronic devices reduces the overall inertia level and damages the frequency stability of power systems [3], which is more prominent in microgrids. This is because conventional The associate editor coordinating the review of this manuscript and approving it for publication was Junjian Qi . control strategies for grid-connected inverters, such as PQ control and droop control strategies, lack inertia and damping as synchronous generators (SGs). To deal with this dilemma, the technology of virtual synchronous machine (VSM) [4] was proposed to provide inertia support to power systems. Although control strategies emulating virtual inertia in existing literatures, such as VSM [5], virtual synchronous generator [6], [7], synchronverter [8] and demand side VSM [9], were slightly different from each other, the principles were similar in the aspect that all of them mimic the inertial characteristic of the SG by emulating its fundamental swing equation. Thus, it is widely recognized that control strategies providing virtual inertia by introducing the swing equation can be classified as VSM control.
But whether the introduction of VSM technology into microgrids will improve their dynamic responses and static stabilities still requires research and verification. VSM control was compared with frequency droop control in [10], and VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ it concluded that both schemes could be regarded as equivalence under certain conditions. But the low-pass filter (LPF) of the measured output power was ignored in the VSM control scheme, which was not in accord with the practical application. Besides, only the frequency-active power (f-P) response was analyzed and the coupling interaction between voltage and active power was dismissed, i.e., a low R/X ratio was considered as an implicit assumption, which was not suitable for the microgrid practice. Comparisons of dynamic characteristics between VSM control and droop control was discussed in [11], and it demonstrated that VSM control had better frequency stability than droop control. It was proved in [12] that frequency resonance between SGs could be suppressed with the introduction of a VSM. However, only single-VSM scenario and SG-connected scenario were considered in [11] and [12], which was too simple to meet the scale of an actual microgrid. Besides, few researches considered the demand side VSM in a microgrid.
To analyze dynamic characteristics of a microgrid, the first step is its small-signal modeling [13]. A practical microgrid can be described by numerous nonlinear differential algebraic equations (DAEs). The small-signal model is obtained by linearizing the nonlinear DAEs around a set of steady-state operating points. However, for microgrids, their small-signal models exhibit multi-time scale characteristics, which requires small time steps for fast dynamics and a long simulation time to capture slow dynamics [14]. With the growing scale of microgrids, their time-domain simulations are becoming time-consuming and computational burdens. The necessity of model order reduction (MOR) for microgrids has never been more obvious.
According to [15], MOR methods are generally classified into three categories, i.e., polynomial approximations, state truncations and parameter optimizations. Polynomial approximation methods are based on matching moments or Markov parameters between the original and reduced-order models, which are applied to transfer function models, e.g., Padé approximation (PA). State truncation methods usually involve transformations of the original state-space model to reconfigure states according to observability, controllability or response time, etc., and elimination of ''less significant'' states, e.g., balanced transformation (BT) method and singular perturbation (SP) method. Parameter optimization approaches are based on optimizing parameters of reduced-order models to minimize errors of response data between reduced and original models, e.g., gray-box system identification.
As small-signal models of microgrids are often described in state-space forms, polynomial approximation methods were seldom used in microgrids MOR, whereas state truncation methods, especially the SP method, were commonly used methods. In [16], the SP method was used to obtain reduced-order models for islanded microgrids by exploiting different dominant time constants. Stability properties of the reduced model turned out to be consistent with the original model, so that the stability analysis of islanded microgrids could be simplified. In [17], the SP method was used to obtain the reduced-order small-signal model of a microgrid in both grid-connected and islanded conditions. Time-domain dynamic responses of the reduced model presented slight difference against the full-order model. In [18], the SP method and Kron reduction were combined to reduce large-signal dynamic models of islanded microgrids in temporal and spatial aspects, respectively. In [19], a method properly eliminating fast states of network dynamics was proposed for microgrids to predict stable regions of droop coefficients more accurately than the quasi-stationary reduced model neglecting all network dynamics. However, the common disadvantage of these methods is that the states in reduced-order models lose their original physical meanings.
The parameter optimization approach can also be interpreted as developing the dynamic equivalent model. In [20], the evolutionary particle swarm optimization based gray-box system identification was used to obtain an equivalent model for f-P dynamics of a microgrid. In [21], [22], black-box models for microgrids using the Prony method were proposed based on measurements of voltage, current and output power at the point of common coupling (PCC). In [23], a novel structure gray-box model including basic electrical components was proposed for microgrid system identification. The identification procedure was carried out by four genetic algorithm-based optimization steps using measurement data at the PCC. The common disadvantage of these methods is that the accuracy of reduced models is affected by the scale of measurement data.
In this paper, the VSM-based microgrid (VSMG) is studied, in which both DERs (including energy storage systems) and controllable loads (with rectifiers at their front end) all adopt VSM control. In the frequency-domain, modal analysis is applied to the VSMG small-signal model to clarify the interaction between states and modes in the VSMG. Besides, bode plots, eigenvalue trajectories and time-domain dynamic responses of various scenarios are comparatively analyzed between VSMG and droop control-based microgrid (DMG) to figure out the influences of VSM control on dynamic characteristics of microgrids. Inspired by the conclusion of participation factor analysis and the thought of coherency in power systems, a new MOR method for VSMG is proposed by equivalencing the detailed small-signal model to a modified third-order SG model. The proposed equivalent SG (EqSG) method derives equivalent electrical parameters for the VSMG and remains certain physical meanings in the reduced states without relying on measurement data.
The remainder of this paper is organized as follows. The detailed small-signal model of the VSMG is built in Section II. Static stability and dynamic characteristics of the VSMG are analyzed and compared with the DMG in Section III. Then in Section IV, three representative MOR methods, i.e., BT, SP and poles clustering-based PA (PCPA) methods, are reviewed and applied to the VSMG model. In Section V, the proposed EqSG method is developed and introduced in detail. In Section VI, four MOR methods are compared to summarize their advantages and disadvantages, and the effectiveness of the proposed method is verified through simulation studies. Finally, the conclusions are drawn in Section VII.

II. SMALL-SIGNAL MODELING OF THE VSMG
In this Section, modeling of the VSMG is analyzed and expressed in terms of DAEs at first. These equations are nonlinear and need to be linearized around a set of equilibrium points, i.e., small-signal modeling, to study the system static stability and dynamic responses. Structure of the studied VSMG is shown as Fig. 1, where DER1 (SG) at Bus0, DER2 (VSM) at Bus1 and a controllable load (VSM) at Bus2 are connected to a local load at Bus3 through three distribution lines. This microgrid is connected to the upstream network from Bus3 through a grid-connection line. In view of existing researches on DMG modeling [13], [18] and [24], this Section focuses on the distinctions appeared in VSMG modeling.

A. MODELING OF A VSM
The block diagram of a VSM is shown in Fig. 2. The power part consists of a three-leg converter and a LC filter. Since the dc side of a VSM often contains an energy storage module to provide inertial energy support and responds fast enough to maintain dc voltage stable, it can be equivalent to an ideal dc source. The control part can be divided into three parts, i.e., VSM control loop, virtual impedance control and voltage current dual-loop control. The VSM model is implemented in its local rotational reference frame (dq-axis). Subscript d or q appeared in the following voltages and currents represents corresponding d or q axis component.

1) VSM CONTROL LOOP
VSM control loop generates the virtual internal electromotive force. The function of this loop is to control the output power of the converter. Active power control is realized by mimicking the swing equation: where J and D p are virtual moment of inertia and virtual damping factor, respectively; ω n and ω are nominal angular frequency and VSM's virtual angular frequency, respectively; P m and P e are active power reference and output active power, respectively. The dq-axis output voltage (v cd , v cq ) and current (i od , i oq ) measurements are used to calculate the instantaneous output active power, and a first-order LPF with the corner frequency ω c is used to obtain P e : Note that the filtered output power is selected as the feedback power value, so as to attenuate potential ripples in the instantaneous power, which will reflect in frequency and magnitude of the voltage reference [7]. Its reactive power control mainly consists of a droop link and an integral link, so as to mimic voltage-reactive power (V-Q) droop characteristic and excitation inertia.
where E is magnitude of the virtual internal electromotive force, U = v cd 2 + v cq 2 is the output voltage magnitude, U n denotes the nominal voltage magnitude, K and D q are virtual excitation inertia factor and V-Q damping factor, respectively; Q m and Q e are reactive power reference and output reactive power, respectively. Similar to P e , Q e is the filtered instantaneous reactive power:

2) VIRTUAL IMPEDANCE CONTROL
Virtual impedance control generates the virtual stator terminal voltage reference (v cd * , v cq * ). Its function is to simulate the stator impedance. It increases the output impedance to help inhibit harmonic circulation between paralleled VSMs and help the decoupling of active and reactive power control.
where r v and L v are virtual resistance and virtual inductance, respectively. Dynamics of the virtual inductance are not considered to avoid instability or poorly damped oscillations caused by the introduction of current derivative [25]. VOLUME 8, 2020 3) VOLTAGE CURRENT DUAL-LOOP CONTROL AND THE POWER PART Voltage current dual-loop control generates voltage command for the space vector pulse width modulation (SVPWM), in which PI controllers are adopted. Its function is to control v cd and v cq tracking v cd * and v cq * promptly and errorlessly. As the modeling of PI controllers, voltage source converters and LC filters is identical between VSMs and droop control-based inverters, equations of this part are listed in Appendix A.

B. COMPLETE MICROGRID MODEL
Network dynamics are generally neglected in small-signal modeling of conventional power systems. However, in the case of microgrids, network dynamics can greatly influence the slow states associated with inverter power controllers [19] and further influence the system stability under certain circumstances [13]. Here, dynamics of the network and the local load are considered, and corresponding DAEs are presented in Appendix A.
Modeling of the SG considers dynamics of the prime mover (p m ), excitation (v f ), rotor rotation (ω r ), stator (ϕ d , ϕ q ), field (ϕ fd ) and damper (ϕ kd , ϕ kq ) windings, which was discussed in detail in [26]. Then, combining all the DAEs of each component mentioned above, the nonlinear model of the studied VSMG can be established in the format as (6).
where y is the output vector, which is determined according to research needs, x and u denote the state vector and the input vector, respectively. For the islanded VSMG, u = [ r load L load ] T , where r load and L load are resistance and inductance of the local load, respectively. And for the grid-connected VSMG, u = [ ω g U g ] T , where ω g and U g denote angular frequency and voltage magnitude of the upstream network at PCC, respectively.

C. SMALL-SIGNAL MODEL OF THE SYSTEM
The small-signal model is obtained by linearizing (6) at a set of equilibrium points. There are two ways of finding equilibrium points. One method is to set the left part of the differential equations in (6) to zero, as 0 = f (x, u), and then solve these equations to determine the equilibrium points. The other is to simulate the nonlinear model and take a snapshot at a certain time when the system reaches steady-state. Because both methods yield same results, the simulation-based method was used in this paper. System parameters of the studied VSMG are displayed in Appendix A. The small-signal model represented in the state-space format as (7) can be generated using MATLAB Symbolic Math Toolbox.
where denotes small disturbance.
As the following dynamic analysis focuses on the influences introduced by VSM control in a microgrid and the comparison of dynamic characteristics between VSMG and DMG, only local control is considered in the microgrid modeling. Besides, the studied DMG maintains the same network structure and system parameters as the VSMG except for the power control loop. In the studied DMG, the power control loop of DER2 adopts the droop control strategy: where droop coefficients K f and K v are determined by realizing the same steady-state droop effect as D p and D q , respectively, according to (10).
Besides, the controllable load in the VSMG is substituted by a constant load with the same power consuming in the DMG.

III. DYNAMIC ANALYSIS OF THE VSMG A. FREQUENCY-DOMAIN ANALYSIS
Because the islanded operation of microgrids is more prone to instability, the frequency-domain analysis is mainly focused on the islanded VSMG. All eigenvalues (λ 1−42 ) of the islanded VSMG are listed in Appendix B and corresponding eigenvalue spectrum is shown in Fig. 3. As λ 39−42 are far from the imaginary axis and have little influence on the stability and dynamics of the system, they are not shown in Fig. 3. The eigenvalue spectrum can be divided into three groups according to its distribution in the s plane. The eigenvalues in group 1 are close to the imaginary axis and are critical to the small-signal stability. It is necessary to analyze which states are sensitive to these eigenvalues.

1) PARTICIPATION FACTOR ANALYSIS
Participation factor reflects the degree of each state participating in every mode (eigenvalue). Assume w i T and v i are respectively the left and right eigenvectors corresponding to λ i of matrix A. The components of v i indicate the relative activity of each state in the ith mode. The components of w i T weight the effect of initial conditions in exciting the ith mode. The normalized participation factor is defined as: where N is the number of eigenvalues.
Using (11), major participants of each mode can be found, which are listed in Appendix B. Eigenvalues in group 3 are mainly influenced by the states derived from the local load, LC filters and the network impedance. Eigenvalues λ 29−32,35−38 have relatively high natural frequency, which is determined by the inherent resonance characteristic of LC filters. Eigenvalues λ 33,34 indicate that dynamics of the network should not be ignored in the microgrid modeling, because this pair of conjugate eigenvalues is mainly stimulated by the interaction of LC filters and the network impedance. In terms of eigenvalues in group 2, they are mainly affected by the current controllers of VSMs and the stator inductance of the SG. Eigenvalues λ 21−24 indicate that PI parameters in the current controllers are well-designed as these modes are well-damped and have relatively small time constants. As for eigenvalues in group 1, the states that have relatively large participation factors are mainly from the SG and VSM control loop. Thus, these states are crucial in influencing the small-signal stability. Furthermore, eigenvalues in group 1 have relatively large time constants, which indicates that they are slow dynamics in the transient process. In other words, the quality of dynamic responses is mainly influenced by the states derived from the SG and VSM control loop.  case of VSMG than that of DMG, which indicates the voltage frequency and magnitude of VSMG are less disturbed when the local load changes.

3) SMALL-SIGNAL STABILITY COMPARISON
In order to test small-signal stabilities of islanded VSMG and DMG, trajectories of dominant eigenvalues are derived by varying the local load, damping factors (droop coefficients) and virtual impedance parameters in the small-signal models, as shown in Fig. 5, where arrows specify increasing of corresponding parameters. The influence of changing virtual inertia parameters is not analyzed here as it was discussed in [6] and [11] and cannot form a comparative analysis with the DMG. It can be observed form Fig. 5(a) and (b) that the variation of r load or L load has little effect on the dominant real eigenvalues for both microgrids, and their locations are very close in the s plane. As r load increases, the dominant conjugate eigenvalues of both microgrids move toward the left, indicating increase of small-signal stabilities, whereas the growing L load leads to the opposite situation. Besides, with the same r load or L load , dominant conjugate eigenvalues in the VSMG are distinctly farther from the imaginary axis than those in the DMG, indicating better damping ratios in the dominant oscillation modes of VSMG. Fig. 5(c) and (d) exhibit eigenvalue trajectories with damping factors D p and D q increasing (corresponding droop coefficients K f and K v are determined according to (10)). It can be seen from Fig. 5(c) that small-signal stabilities of both microgrids gradually increase with D p growing. Under the premise of achieving the same steady-state droop effect, D p in the VSMG has a larger stable region than K f in the DMG. In Fig. 5(d), when D q = 6.5, the DMG is small-signal unstable as there exists a pair of conjugate eigenvalues (15 ± i673) on the right half of the s plane (which are not shown in the figure since they are far from the origin), whereas the VSMG is small-signal stable. Then, as D q increases, the DMG becomes small-signal stable, and the dominant real eigenvalues in both systems travel toward the unstable region, which indicates that small-signal stabilities of both microgrids are degraded. When D q = 1e4, the DMG remains small-signal stable, whereas the VSMG becomes small-signal critical stable.
The reason why VSMG and DMG exhibit different stability with too small or too large D q can be explained by analyzing (3) and (9). According to (3), the reactive power control of VSM is realized by detecting the terminal voltage deviation and adjusting the virtual internal electromotive force E. Besides, with the additional integral link, its output reactive power can be accurately controlled as Q e = Q m + D q (U − U n ). Thus, a relatively large D q results in a large adjustment in E and a large deviation in Q e , which is prone to cause instability, whereas a relatively small D q can remain stable for the VSMG. According to (9), the reactive power control of droop control strategy is realized by detecting the output reactive power Q e and adjusting the voltage reference E. A relatively small D q results in a large adjustment in E, which is prone to cause instability, whereas a relatively large D q can remain stable for the DMG. Even though the stable regions of D q and K v cannot be compared as they perform opposite stability at both ends of the parameters interval, the dominant conjugate eigenvalues in the VSMG have lower natural oscillation frequency and greater damping ratio compared to the DMG, which indicates the VSMG has better dynamic performance than the DMG under the same steady-state reactive power droop.
It can be observed from Fig. 5(e) and (f) that the effects of increasing r v and L v are similar. For the DMG, the increasing of virtual impedance (no matter r v or L v ) introduces additional voltage drop and is equivalent to increasing K v [27], i.e., decreasing D q , so that eigenvalue trajectories of DMG are consistent with Fig. 5(d). However, for the VSMG, the additional voltage drop effect can be alleviated with the additional integral link in (3) by increasing E, so that the stability of VSMG is less sensitive to virtual impedance parameters, as indicated by the real dominant eigenvalues in Fig. 5(e) and (f). This also confirms that VSMG has wider stable regions for r v and L v than DMG. Hence, VSMG is significantly more robust against parameter variations and load uncertainties than DMG.

B. TIME-DOMAIN ANALYSIS
Time-domain responses (obtained by simulation of nonlinear DAEs) of both kinds of microgrids are compared in both islanded and grid-connected operation. In the islanded operation, dynamic responses of voltage frequency and magnitude at Bus3 are compared under the local load changing from 25 + j10 kVA to 50 + j10 kVA at time t = 5s, as shown in Fig. 6(a) and (b). Before the load changes, the voltage frequency and magnitude in the VSMG are closer to nominal values, i.e., 50Hz and 311V, respectively. From Fig. 6(a), it can be seen that the maximum frequency drop of the VSMG is f V ,max = 0.14Hz, while that of the DMG is f D,max = 0.21Hz. The steady-state frequency drop of the VSMG is f V = 0.07Hz, while that of the DMG is f D = 0.08Hz. The maximum absolute value of the rate of change of frequency (ROCOF) in the case of VSMG is 5.06Hz/s, while that of the DMG is 6.00Hz/s. From Fig. 6(b), it can be observed that the maximum voltage drop of the VSMG is U V ,max = 23.2V, while that of the DMS is U D,max = 28.1V. The steady-state voltage drop of the VSMG is U V = 9.3V, while that of the DMG is U D = 11.7V. The deviations of voltage frequency and magnitude are smaller in the VSMG than in the DMG in both dynamic process and steady-state, which agrees with the frequency-domain analysis.
In the grid-connected operation, dynamic responses of output active and reactive power are compared under disturbances of ω g and U g . Figure 6(c) depicts dynamic responses of output active power with ω g dropping 0.2πrad/s from the nominal value at t = 10s, where negative values mean the power is absorbed from the upstream network and vice versa. The variation of output active power in the VSMG ( P oV = 21.0kW) is greater than that in the DMG ( P oD = 10.6kW).
Note that before ω g drops, the steady-state output active power of the VSMG is slightly smaller than that of the DMG. This is because the voltage magnitude at Bus3 is higher in the case of VSMG, which leads to a larger absorbed active power by the local load in the VSMG. Figure 6(d) shows dynamic responses of output reactive power with U g dropping 0.1p.u. from the nominal value at t = 10s. The change of output reactive power in the VSMG ( Q oV = 18.8kvar) is greater than that in the DMG ( Q oD = 15.5kvar). It is clear that the VSMG provides more active and reactive power support than the DMG when voltage frequency and magnitude disturbances occur in the upstream network.
Time-domain comparisons indicate that VSMG performs better responses than DMG in both islanded and grid-connected operation. This is because, firstly, VSM control can be applied to both the supply side and the demand side. In the studied VSMG, both DER2 and the controllable load applying VSM control can automatically adjust their output power to reduce fluctuations of local voltage frequency and magnitude in the islanded operation and support the upstream network in the grid-connected operation, whereas only power supplies in the DMG can participate in the system regulation. Thus, VSMG have more satisfactory steady-state values than DMG. Secondly, in the case of reasonable design of virtual inertia parameters, the additional dynamics introduced by VSM control are coordinated with the inherent inertial dynamics of the SG, so as to reduce oscillations and smooth fluctuations in the dynamic process.

IV. MODEL ORDER REDUCTION OF THE VSMG
With the growing scale of VSMG, the order of its model will inevitably increase, so as the computational burdens of its analysis and simulation. In this Section, three representative MOR methods, i.e., BT, SP and PCPA methods are applied to the VSMG model.

A. MODEL ORDER REDUCTION BASED ON BALANCED TRANSFORMATION
Generally speaking, BT is a state coordinate transformation that makes the controllability and observability Grammians of the system identical and diagonal. Consider the small-signal VSMG model represented as (7), where A is assumed asymptotically stable, the pair (A, B) is assumed controllable and the pair (C, A) is assumed observable. The controllability and observability Grammians, denoted as P and Q, respectively, can be obtained by solving the Lyapunov equation: It is known that there exists an invertible matrix T to obtain the balanced realization [28] of (7) as ( is omitted for simple expression): where with corresponding controllability and observability Grammians P b and Q b satisfying the following equation: where σ i is the so called Hankel singular value (HSV) and σ i = √ λ i (PQ). λ i (PQ) denotes the ith eigenvalue of PQ. The HSVs are important criteria in system MOR because they measure the contribution of each state to the input-output behavior. Detailed steps of obtaining T is elaborated in [29].
Assuming that the ratio: is less than an acceptable small value, the states in (13) with indices from r + 1 to N can be regarded having little contribution to dynamic responses and can be eliminated by considering corresponding dynamics equal to 0. Then the balanced VSMG model (13) can be rewritten in the partitioned quasi-steady-state format as: where A 11 is of dimension r ×r, and A 22 is of (N −r)×(N −r), with remaining matrices having corresponding dimensions consistent with the system dimensions defined in (13). Even though the dynamics of x b2 are discarded, its steady-state gain can be remained by replacing (16) into (15).
Thus, the reduced VSMG model can be derived as: where Theoretically, BT method can reduce the system model to arbitrary order.

B. MODEL ORDER REDUCTION BASED ON SINGULAR PERTURBATION
The SP theory believes that the singular perturbed system can be decoupled into independent slow and fast subsystems, and the system reduction can be realized by neglecting the fast subsystem. For the VSMG with greater inertia than conventional microgrids, where slow dynamics are dominant, neglecting fast dynamics will have little influence in system responses. The VSMG model as (7) should be transformed VOLUME 8, 2020 into the singularly perturbed format before applying the SP method. According to the participation factor analysis, rearrange state vector of (7) such that slow states x 1 are placed at upper rows and fast states x 2 are placed at lower rows, and the singularly perturbed presentation of (7) can be denoted as ( is omitted for simple expression): where ε is a small positive singular perturbation parameter that indicates the separation between slow and fast states. Applying the Chang transformation, (18) can be decoupled into independent slow and fast subsystems as: where z 1 and z 2 are decoupled slow and fast states, respectively;

matrices L and M are introduced by
Chang transformation and can be solved by iterations detailed in [17]. Since small ε is associated with the decoupled fast subsystem, (19) can be reduced as: The SP method decouples the original system into independent fast and slow subsystems by calculating ''boundary layer'' corrections in separate time scales [30]. Generally speaking, this method allows for proper inclusion of possible effect that fast states have on the slow subsystem. The model proposed in [19] was essentially the first iteration of L and M and deriving corresponding reduced microgrid model. Additionally, SP method can only reduce the VSMG model to certain orders where the states present distinct separation of time scales, otherwise matrices L and M cannot be solved.

C. MODEL ORDER REDUCTION BASED ON POLES CLUSTERING
The PCPA method is suitable to reduce systems represented by transfer functions. Most MOR methods for transfer function models simplify denominator and numerator polynomials separately. The traditional approach of denominator simplification mainly considers preserving stability between the original and reduced systems [31], such as Mihailov stability criterion [32], stability equation method [33], Routh approximation [34], etc. Instead, the poles clustering method engages all properties of the original system poles by agglomerative clustering depending on their distributions in the s plane. The numerator simplification considers the matching of time moments or Markov parameters, e.g., PA, or regards minimizing response errors between the original and reduced systems by heuristic optimization algorithms.
Since the PCPA method is applied to the transfer function model, the VSMG model as (7) should be firstly transformed into several single-input single-output (SISO) transfer functions. The procedure of poles clustering consists of the following steps: 1) Divide all the left half poles of the original transfer function into a certain number of clusters by the Euclidean distance based hierarchical poles clustering algorithm. The number of possible clusters determines the order of reduced system. Note that real and complex poles should be clustered separately.
2) Estimate a cluster-center for each cluster.
3) Improve the cluster-center according to the dominant pole of this cluster.
The detailed procedure is explained in [35]. Note that poles on the left half and the right half of the s plane are clustered separately, whereas poles on the imaginary axis or the origin should be retained in the reduced model.
After the poles clustering procedure, coefficients for the numerator are determined through the PA procedure. Expand the original transfer function G(s) in Taylor series around s = 0 and s = ∞ respectively as: where T mi and M ki are the ith time moment and Markov parameter of G(s), respectively. Consider the rth-order reduced model as: G r (s) = c 0 + c 1 s + c 2 s 2 · · · + c r−1 s r−1 d 0 + d 1 s + d 2 s 2 · · · + d r s r (22) where c i is the ith numerator coefficient of the reduced model and d j is the jth denominator coefficient obtained by the poles clustering procedure. The coefficients of the reduced numerator are evaluated as: where γ is the number of time moments and β is the number of Markov parameters. Therefore, all coefficients required for the reduced model as (22) are known.

V. EQUIVALENT SG MODEL OF THE VSMG
The participation factor analysis confirms that modes in group 1, which are low frequency oscillation modes with relatively large time constants and highly influencing the system stability, are mainly derived from states related to the SG and VSM control loop. This inspires us that the VSMG model can be simplified into the structure of the SG model. The equivalent SG parameters of the VSMG can be obtained through the method of dynamic responses matching. Similarly, in the field of power system model reduction, one of the most commonly used methods is coherency and aggregation. Coherency means that some SGs in a wide area exhibit similar rotor angle swings after a disturbance. Then the coherent groups of SGs can be aggregated into a single equivalent SG so as to simplify the system model. As for the microgrid, voltage angle differences of each node are not distinct in steady-state or under small disturbances. Besides, both DERs and controllable loads in the VSMG inherit output characteristics of synchronous machines. Inspired by the concept of coherency and aggregation, it is supposed that the VSMG model can be equivalent as an SG model. The EqSG model can be used in time-domain simulations and dynamic analysis of ADNs containing multi-VSMGs by substituting detailed VSMG models.
Dynamic equations of the EqSG model are based on the modified third-order SG model: The third-order EqSG model considers the excitation transient and rotor dynamics, where the superscript eq denotes corresponding equivalent variable or parameter of the VSMG, P mg and Q mg are active and reactive power instructions derived by the microgrid secondary control, respectively; P eq e and Q eq e are internal active and reactive power of the EqSG, respectively, which can be calculated by the following equations: P eq e = 1.5[(E eq ) 2 cos α − E eq U g cos(δ eq + α)]/Z Q eq e = 1.5[(E eq ) 2 sin α − E eq U g sin(δ eq + α)]/Z (25) where Z is the equivalent impedance from the EqSG to the upstream network and α is the impedance angle. Note that each VSM in the VSMG cannot sample the PCC voltage and its output reactive power is directly influenced by the capacitor voltage of its LC filter, so the internal sampling voltage, denoted as U eq in (24), is introduced.
The EqSG model is proposed for grid-connected VSMG. Thus, ω g and U g are input signals. The output active and reactive power, denoted as P o and Q o , are defined as output signals: The structure of the EqSG model is defined through (24)- (27), and the complete EqSG model can be obtained by determining the following 7 equivalent parameters, i.e., J eq , D eq p , K eq , D eq q , Z , α, and r x . In this paper, these equivalent parameters are determined by matching dynamic responses of the EqSG model with the complete VSMG small-signal model, which can also be recognized as gray-box system identification. Since the equivalent model is effective under small disturbances, the linearized EqSG model is used in the equivalencing procedure, which can also speed up the gray-box system identification. The linearized EqSG model can be presented in the state-space format: x eq = A eq x eq + B eq u y = C eq x eq + D eq u (28) where state-space matrices A eq , B eq , C eq and D eq are displayed in Appendix C. The equivalencing procedure is illustrated by the flowchart in Fig. 7.

1) STEP 1: CREATE ESTIMATION DATA
Since the EqSG model is studied in the small-signal region, fluctuations in the input signals should satisfy the meaning of small disturbances in power systems. According to ''GB/T12325-2008'' and ''GB/T15945-2008'', fluctuations of voltage magnitude and frequency in power systems should not exceed ±5% and ±0.2Hz, respectively, which can be regarded as the range of small disturbances. In order to improve the applicability of estimated equivalent parameters, fluctuations in the input signals of estimation data reach their extreme values of the small-signal region in the form of step changes. The output signals of estimation data VOLUME 8, 2020 are dynamic responses derived from the small-signal gridconnected VSMG model according to input signals of estimation data. Thus, the estimation data can be collected into 3 sets as summarized in Table 1.

2) STEP 2: SET INITIAL VALUES OF EQUIVALENT PARAMETERS
The System Identification Toolbox in MATLAB was used to create the linear EqSG gray-box model (28). The initial values of equivalent parameters are obtained by an once Levenberg-Marquardt algorithm based parameters estimation with estimation dataset 1. The goal is to determine proper initial values of unknown equivalent parameters by minimizing the sum of squared errors between [ P oV 1 (t), Q oV 1 (t)] and time-domain responses of the gray-box model (28).

3) STEP 3: ESTIMATE EQUIVALENT PARAMETERS ALTERNATELY AND ITERATIVELY
Due to the coupling of output active and reactive power in the low-voltage VSMG, the estimated parameters derived from estimation dataset 1 are only effective within a narrow range of input signals. For the linearized EqSG model, dynamic responses of multi-inputs can be dealt with the sum of each SISO response. Thus, the estimation procedure is subdivided into two parts, i.e., the voltage-variation response estimation and the frequency-variation response estimation. These two subprocesses are executed alternately and iteratively until the unknown equivalent parameters convergence, so that the estimated parameters will satisfy a wide range of input signals.
In the voltage-variation response estimation, the estimation dataset 2 is used, in which step changes occur in U g while ω g remains unchanged in the input signals. To obtain output responses faster, the gray-box model (28) is separated into two SISO transfer functions, denoted as T UP (s) and T UQ (s), respectively. From (24), it's clear that J eq and D eq p mainly influence the response of output active power when ω g changes, whereas K eq , D eq q and r x primarily affect the response of output reactive power when U g changes. Thus, in the voltage-variation response estimation, K eq , D eq q , Z , α, and r x are allowed to change while other parameters (J eq and D eq p ) are fixed. Similarly, in the frequency-variation response estimation, the estimation dataset 3 is adopted. The gray-box model (28) is separated as T fP (s) and T fQ (s) where J eq , D eq p , Z and α are able to change while other parameters (K eq , D eq q and r x ) are constant.
In each iteration, parameters allowed to change are updated by minimizing the objective function J U (for voltage-variation response estimation) or J f (for frequency-variation response estimation) through interior point algorithm, which can be fulfilled with MATLAB Optimization Toolbox. The objective functions are weighted square average errors between time-domain responses of the estimated model and output signals of the estimation data: where y UP (t) and y UQ (t) are time-domain responses of T UP (s) and T UQ (s) with input signals of estimation dataset 2, respectively; y fP (t) and y fQ (t) are time-domain responses of T fP (s) and T fQ (s) with input signals of estimation dataset 3, respectively; N s is the number of sampling points; w UP , w UQ , w fP and w fQ are weight coefficients and w UP = w fQ = 0.2, w UQ = w fP = 0.8. After the parameters convergence, the EqSG model is obtained by updating the estimated parameters in (28).

VI. CASE STUDIES
Three MOR methods mentioned in Section IV and the proposed EqSG method are compared with each other in the studied grid-connected VSMG, so as to summarize their advantages and disadvantages in the application of microgrid model reduction. The input signals for testing different MOR methods are described in Table 2, where disturbances occur in ω g and U g respectively. The amplitudes of disturbances are 50% of their extreme values of the small-signal region. Since the proposed EqSG method fixes the reduced VSMG model into third-order, reduced models derived from BT and PCPA are also restricted to third-order for a fair comparison. Note that the reduced model derived from SP is a fourth-order model as the slow and fast subsystems of the studied VSMG cannot be decoupled at third-order. Reduced models derived by three MOR methods and equivalent parameters derived from the EqSG method are displayed in Appendix C.  Fig. 8 compares frequency-domain characteristics of f-P responses among the reduced-order models. It can be observed form Fig. 8(a) that all of the four MOR methods provide excellent matching with the full-order VSMG small-signal model in low frequency band (0-10rad/s). The EqSG method provides the widest frequency band (0-450rad/s) in properly approximating to the full-order model among four MOR methods. Both BT and SP methods exhibit obvious truncation characteristics because they both eliminate high frequency modes to some extent, but the properly matched frequency range of BT method (0-305rad/s) is larger than that of SP method (0-25rad/s). The PCPA method does not show truncation characteristics. The slope of PCPA method in high frequency band (>1e4 rad/s) is the same with that of the full-order model (−20dB/dec), but its amplitude response presents apparent difference with the original model from the corner frequency (10rad/s). This is because PCPA method is able to preserve high frequency poles, but the clustered high frequency poles are modified according to the dominant poles to preserve slow dynamics preferably, which leads to a slope match but amplitude mismatch in the high frequency domain. Figure 8(b) depicts phase-frequency characteristics of f-P responses, where BT method performs a better approximation to the full-order model than other methods while the SP method exhibits the worst matching.   Fig. 8(c) that PCPA method performs better approximation in high frequency band (>1e3 rad/s) than other MOR methods, because the slope of PCPA method is the same with the full-order model in high frequency band while other MOR methods exhibit high frequency truncation characteristics. There are mainly two reasons. Firstly, the V-Q response is generally faster than the f-P response. The elimination of some high frequency dynamics caused by BT, SP and EqSG methods will inevitably reduce the proximity in the V-Q response. Secondly, PCPA method is based on the transfer function model, which can provide the reduced model for the V-Q response independently, and the clustering of poles is not affected by other input-output responses. But in low-mid frequency band (0-360rad/s), BT and EqSG methods perform better approximation than SP and PCPA methods. From Fig. 8(d), it can be observed that BT method performs better approximation than other MOR methods in the phase-frequency characteristic.
Note that, in the f-P response, the full-order model is a minimum phase system, reduced models derived by PCPA and EqSG are also minimum phase systems, while reduced models derived by BT and SP are non-minimum phase systems. In the V-Q response, the full-order system is a non-minimum phase system and all reduced models are also non-minimum phase systems. This indicates that PCPA and EqSG methods are better than BT and SP in maintaining the structural properties of the original system.
Further, a comparison of time-domain dynamic responses is analyzed and shown in Fig. 9. It can be seen from Fig. 9(a) and (c) that all methods match well with the full-order model in the steady-state. Whereas in the transient process, for the f-P response, PCPA and SP methods show relatively larger response errors than BT and EqSG methods, as shown in Fig. 9(b). For the V-Q response, shown as Fig. 9(d), all methods fail to precisely match with the high frequency oscillations of the full-order model in the transient process. But the EqSG method performs better matching than other three MOR methods. To compare the matching degrees of different MOR methods in a quantitative way, the normalized root-mean-square error (NRMSE) is chosen as the index for evaluation of matching degrees:    Table 4 compares different properties of four MOR methods in VSMG model reduction, where the EqSG method shows more superiorities than other MOR methods. The EqSG method is suitable for both state-space model and transfer function model because it is based on dynamic equivalence of the original model. The EqSG method performs the best matching in low frequency of amplitude-frequency characteristics and has the best fitting degree in time-domain responses. Besides, this method maintains identical type of phase system with the original model. As for the stability consistency, BT and SP methods satisfy stability consistency because they require the original model to be asymptotically stable and their procedures are essentially linear transformations of state matrices which will not change system stability. The PCPA method preserves stability consistency because the poles on both sides of the s plane are clustered separately. The EqSG method also retains stability consistency because this method correctly fits output signals of the estimation data which reflects stability of the original system. Furthermore, the most significant advantage of the EqSG method is to keep equivalent physical meanings in the reduced states and provides equivalent electrical parameters which are instructive in power system analysis.
To further verify the effectiveness of EgSG model substituting detailed VSMG model in power system simulations, historical PCC frequency and voltage data of an actual microgrid, shown as  In this case study, fluctuations of ω g and U g are applied simultaneously to be consistent with the actual situation. It can be observed from Fig. 11 that responses of the EqSG model properly match with the detailed VSMG model (output active and reactive power values at nominal ω g and U g are subtracted in detailed VSMG's waveforms). Matching degrees of the EqSG model calculated according to (31) are 96.12% for output active power response and 97.93% for output reactive power response in this case study. Since step changes seldom occur in ω g and U g in actual operation, matching degrees in this case are better than those listed in Table 3. As the variation of U g inevitably influences active power in low-voltage microgrids, output active power of the EqSG model exhibits slightly error compared with the detailed VSMG when the variation of U g is relatively large, as shown in Fig. 11(a), so that the matching degree of active power response is slightly lower than the reactive power response in this case. Note that in order to ensure the accuracy of the EqSG model, disturbances in the input signals should not exceed the small-signal region (±0.2Hz in frequency and ±5% in voltage), as the equivalent parameters are obtained based on the small-signal model of the VSMG.

VII. CONCLUSION
This paper compares dynamic characteristics between VSMG and DMG. With additional inertia dynamics and proper inertia parameters, the VSMG exhibits better dynamic characteristics than DMG in both islanded and grid-connected operations. In the islanded operation, the VSMG is more robust against parameter variations, and performs less disturbance in transient process and faster recovery to steady-state when the local load changes. In the grid-connected operation, the VSMG provides more power support to the upstream network than DMG when voltage or frequency fluctuates at PCC.
The participation factor analysis indicates that static stability and dynamic characteristics of the VSMG are mainly influenced by states associated with the SG and VSM control loop. According to this, the third-order EqSG model is proposed to reduce the VSMG model. The EqSG method applies gray-box system identification and is effective in the small-signal region. To speed up the gray-box system identification, linearized EqSG model is used in the equivalencing procedure. To improve the matching degree of output responses, the equivalencing is subdivided into different input-output responses estimation and estimates equivalent parameters alternately and iteratively. The EqSG model shows more superiorities than BT, SP and PCPA methods in both frequency-domain and time-domain. It can be used in dynamic analysis and time-domain simulations of power systems containing multi-VSMGs, as its dynamic responses properly match with the detailed VSMG model. In addition, the EqSG model preserves equivalent physical meaning in reduced states and provides equivalent electrical parameters for power system analysis.

APPENDIX A
Equations for modeling the voltage current dual-loop control are expressed as: (35) where i fd * and i fq * are inductor current references of the LC filter; i fd and i fq are inductor current measurements of the LC filter; v d * and v d * are output voltage references of the converter; K pv , K iv , K pi and K ii are proportional gains and integral gains of the voltage controller and the current controller, respectively; x d , x q , y d and y q are states introduced by integral controllers.
Due to the realization of high switching frequencies, SVPWM and switching process of the converter can be neglected, and the converter can be regarded as a proportional link with gain K inv . For convenience, a gain K g = 1/K inv is introduced to make the output voltage of the converter (v d , v q ) equal to its reference.
The nonlinear model of the LC filter can be represented as: where r f is the parasitic resistance of inductor L f , and C f denotes the shunt capacitor.
Modeling of the network and the local load is represented in the common rotational reference frame (DQ-axis) of the whole microgrid system. For generality, the dq-axis frame of DER1, whose angular velocity is denoted as ω r , is chosen as the DQ-axis frame. Subscript D or Q appeared in the following voltages and currents represents corresponding D or Q axis component. The transformation equations between dq-axis frame and DQ-axis frame are expressed as: where the subscript i (or j) denotes the ith (or jth) Bus, the subscript ij denotes the connection between Busi and Busj, δ i is the angle of Busi's dq-axis frame with respect to the DQ-axis frame, v bD,i and v bQ,i are voltages of Busi, i lineD,ij and i lineQ,ij are currents of Lineij. The RL branch model is adopted for the network because the equivalent capacitance to ground of the line can be ignored in microgrids. where r line,ij and L line,ij are resistance and inductance of Lineij respectively.
The local load is considered as the combination of resistor r load and inductor L load so as to reflect the load dynamics.
i loadD = (v bD,i − r load i loadD )/L load + ω r i loadQ i loadQ = (v bQ,i − r load i loadQ )/L load − ω r i loadD (43) where i loadD and i loadQ are currents of the local load.
Since there is no corresponding output variable to provide voltage value for the node that is not directly connected to a capacitor, a virtual resistor r N is assumed between this type of node and ground to determine the voltage value. The value of r N should be sufficiently large such that its introduction has minimum influence on the system stability.
For the grid-connected VSMG, the input variable U g should also be transformed in the DQ-axis frame to integrate into the microgrid model: v gD = U g cos δ g v gQ = U g sin δ g where δ g is the angle of the upstream grid's rotational reference frame with respect to the DQ-axis frame, v gD and v gQ are grid voltages represented in the DQ-axis frame. Also, dynamics of the grid-connection line with resistance r g and inductance L g should be included in the grid-connected VSMG model. i gD = (v bD,i − v gD − r g i gD )/L g + ω r i gQ i gQ = (v bQ,i − v gQ − r g i gQ )/L g − ω r i gD (47) where i gD and i gQ are grid-connected currents of the VSMG. System parameters of the VSMG are listed in Table 5. For the convenience of analysis, control parameters of each VSM are the same.

APPENDIX B
Eigenvalues and corresponding participation factors are listed in Table 6.

APPENDIX C
The state-space matrices A eq , B eq , C eq and D eq are:  where U g0 is the measured initial value of U g , E eq 0 and δ eq 0 can be calculated according to (27) with measured initial values of output active and reactive power of the VSMG.