Analysis of Harmonic Resonance Characteristics in Grid-Connected LCL Virtual Synchronous Generator

: The virtual synchronous generator (VSG), which emulates the essential behavior of the conventional synchronous generator, has attracted great attention. This paper proposes to analyze the harmonic resonance characteristics in VSG using the state-space model. The analysis is based on a full-order state-space small-signal model that fully considers the dynamic of the inner loops and the VSG-based outer power control loop. Participation analysis is used to point out the contributions of different states to the eigenvalues. Moreover, eigenvalue locus and singular value decomposition (SVD) are applied together to evaluate the impact of the inner loop parameters on the harmonic resonance characteristics around the LCL ﬁlter resonance frequency. The analysis indicates that the harmonic resonance instability is mainly caused by decreasing the proportional gains of the current loop and the voltage loop. Finally, extensive numerical simulation and experimental results are given to verify the validity of the theoretical analysis. Both the simulation and experimental results indicate that the voltage of the common coupling point is unstable after decreasing the proportional gains of the current and voltage controllers. As K pc decreases from 5 to 0.4 or K pv decreases from 0.6 to 0.2, the harmonic distortion factor (HDF) around the LCL ﬁlter resonance frequency increases. Furthermore, the consistency of simulation results, experimental results, and the theoretical analysis results is validated.


Introduction
The increasing integration of renewable energy sources (RESs) has received significant attention [1]. Moreover, most RESs are interfaced to the grid via power electronics converters, which provide power without contributing any inertia support. This further leads to serious frequency stability issues. To tackle this issue, the concept of the virtual synchronous generator (VSG) has been proposed. The idea of operating an inverter to mimic a synchronous generator (SG) is first motivated and developed in [2]. The robustness [3], dynamic response [4], and inertial control [5] have been deeply researched to improve the dynamic characteristics of VSG. However, only the outer power loops are taken into consideration in [2][3][4][5], and the impacts of the internal loops parameters are overlooked. Because of the grid-friendly characteristic of VSG, the research of VSG has been carried out in islanded microgrids (IMGs) [6], photovoltaic power generation [7], variable speed pumped storage hydropower with full-size converter [8], high-voltage direct current (HVDC) transmission [9], wind turbines and static VAR (volt-ampere reactive) compensator (SVC) [10], etc.
One of the critical concerns in the grid-connected virtual synchronous generator system is the harmonic resonance issue, which would arise from the interaction between the network and the voltage source converter, thus causing the harmonic resonance phenomena from hundreds of Hz to multiple kilohertz (kHz). LCL filters can weaken the switchingfrequency harmonics effectively, but the additional LCL filter resonance may interact with the control loops of VSG, leading to LCL filter harmonic resonance instability phenomena. With the increasing popularity of power electronic systems, this problem is becoming more and more serious. Many incidents about the grid integration of renewables [11] and high-speed trains have been reported [12], where the undesired distorted harmonic voltage seriously affects the stability of the control loops and even destroys the equipment in the system.
Consequently, effective modeling and analysis methods are urged to tackle the harmonic resonance instability issue. The transfer function model is used in studies [13,14] to analyze the LCL filter resonance characteristics. In one study [13], an enhanced current resonant controller for grid-tied converters with LCL filter is presented, and the developed method is based on the transfer function. The transfer function model is used to reveal the principle and relations of the resonant controllers in another study [14]. The impedancebased method is the most commonly applied approach to evaluate the stability of power electronics-based power systems. An impedance-based model is developed for the rectifier AC system and the inverter AC system for a VSC-HVDC system to analyze the resonance instability, and then the Nyquist stability criterion and impedance frequency responses are applied to detect resonances [15]. Besides, a great deal of impedance analysis concentrates on the low-frequency oscillations due to the phase-locked loop (PLL) dynamic [16] and the grid strength under weak-grid conditions [17]. The state-space model is another effective tool for evaluating the stability of the power electronics-based power system. To get the state-space model, the system is linearized around the operation points. Furthermore, by analyzing the state matrix, the stability of the system can be predicted. A detailed statespace model of VSG is presented in [18], and the influence of system controller parameters on the poles of the state-space model is also investigated by the parametric sensitivities of the system eigenvalues. In another study [19], a novel harmonic extraction algorithm based on the state-space model is presented to analyze the harmonic characteristics of VSG under a non-ideal grid. Resonance mitigation of VSG control under the weak-grid connected microgrid is analyzed in [20]. However, the current analysis of the harmonic resonance characteristics in VSG is concentrated on the impact of the non-ideal grid, and the impact of the internal controller parameters of VSG on LCL filter harmonic resonance is often overlooked.
Even though various kinds of modeling methods have been developed for several decades, they have their own advantages and disadvantages by themselves according to the purpose of usage. Table 1 outlines the performance comparison of modeling methods. The main advantages of the model built in this paper lie in the identification of dynamic modes and the calculation of the participation factor. In addition to the tools described above, the singular value decomposition (SVD) analysis provides an effective way to point out the contribution of each component on the resonance for the multiple input and multiple output (MIMO) systems [21]. This method is used together with the eigenvalue-based stability analysis to identify critical power systems nodes [22,23]. A robustness study of the power system is carried out by Rodríguez-Cabero  [24]. In [25], SVD is applied to evaluate the dynamic of VSG when the grid model and the virtual impedance varies. The influence of VSG numbers and outer loop parameters on the VSG-based multi-terminal direct current (VSG-MTDC) systems is intensively studied by SVD and modal analysis in [26]. However, this method has not been presented for VSG to uncover the constraints in inner parameters.
In order to avoid the undesired distorted voltage of the common coupling point and further improve the practical engineering applicability of VSG, this paper proposes to analyze the harmonic resonance characteristics around the LCL filter resonance frequency of VSG. The main contributions can be concluded as follows. A full-order state-space model that fully considers the dynamic of the inner loops is built step by step. The impacts of the parameters such as the proportional gains, integral coefficients, and the virtual inductance are analyzed in detail. Participation analysis is applied to assess and identify the contributions of different component on the eigenvalues. Specifically, SVD analysis is adopted to assess the impact of the inner loops parameters on LCL filter harmonic resonance oscillation phenomena. The variations of the inner loop control parameters are considered and the obtained results are compared critically. Finally, the consistency of simulation results, experimental results, and the theoretical analysis results is verified.
The rest of this paper is organized as follows: Section 2 provides the state-space model that fully considers the dynamic of the inner loops and the VSG-based outer power control loop. In Section 3, firstly, participation analysis is used to evaluate the contributions of different states to the eigenvalues, and then, the impacts of the controller parameters on the harmonic resonance are analyzed by eigenvalue-based analysis and SVD analysis. Through time-domain simulation and experimental results, Section 4 respectively validates the analysis. Finally, concluding remarks are provided in Section 5.

Small-Signal State-Space Modeling of VSG
VSG is a power electronic converter that has a SG mathematical model embedded into its power loop controller to simulate the SG rotor dynamic and provide inertia to the power grid. The VSG comprises active-and reactive-power loops and voltage and current inner loops. However, the inner loops are often omitted for simplification as shown in [27,28]. The VSG control in References [27,28] contains inner loops, but only the outer power loops are taken into consideration when building the small-signal model and conducting the small-signal analysis. The inner loops are omitted for a faster dynamic than the outer power loops. However, the inner loops have an essential influence on the harmonic stability analysis, and thus, this section builds a detailed state-space model that fully considers the dynamic of the inner loops. First, the overview of the traditional VSG is provided. Figure 1 shows a VSG-connected system, which consists of the topology of the gridconnected converter and the control diagram of VSG [29]. Table 2 provides the parameters of the VSG-connected system. The control diagram of VSG is composed of the VSG-based power controller, virtual impedance, voltage controller, and current controller. The VSG DC-bus voltage U dc is assumed to be constant. L c is the converter side filter inductance. L g is the grid inductance. C is the filter capacitance. L c , L g , and C form the LCL filter. Let u tdq , u gdq , i Ldq , and i gdq denote the capacitance voltage, the grid voltage, the converter side current, and the grid side current in the dq reference frame, respectively. The u tdqref is the reference value for u tdq , while the i Ldqref is the reference value for i Ldq . Let U g , U t , and E denote the root mean square (RMS) value of the grid voltage u g , the capacitance voltage u t , and the inner electric potential e, respectively. Let v pcc denote the line-to-line voltage of the common coupling point. Let U n denote the rated voltage value and ω n denote the rated angular speed value. Below, each control component of Figure 1 is explained in detail and the state-space model is developed step by step according to the procedure in Figure 2.   Figure 2. Procedure of the state-space modeling.

Overview of the Conventional VSG
Power caculation and filtering     Figure 1. Topology of the grid-connected converter and control diagram of the virtual synchronous generator (VSG) [29].  Figure 2. Procedure of the state-space modeling.
Power caculation and filtering   As depicted in Figure 1, the instantaneous active power p and reactive power q can be computed by [30]: The first-order low-pass filters (LPFs) are implemented to smooth the p and q. The smoothed signals are P f and Q f respectively, as follows [30]: where τ f is the time constant of the LPFs. All the filter time constants in this paper are assumed to be equal. Furthermore, P f and Q f are used as the feedback signals, as they reflect the actual active and reactive power (after being filtered) injected into the grid.

VSG-based Power Controller
The active power loop (APL) imitates the dynamic of the SG rotor and does not rely on a PLL during its normal operation. The standard swing equation can be adopted as [31]: where J is the inertia constant, P set is the active power reference value for P f , ω is the rotating speed of the VSG, and D p is often determined based on the grid code. According to EN 50438, it is required that the change of 100% active power corresponds to the change of 2% angular speed [32]. Then, the frequency droop coefficient, D p , is determined by: It can be seen from Figure 1 that the reactive power loop (RPL) of VSG regulates Q f and U t by adjusting E. Then, the dynamic of E is described by: where K is the integral coefficient that adjusts the dynamic of the RPL and U t can be expressed as: Let Q set denote the reference value of Q f . Being similar to D p , D q is also determined by the grid code. According to EN 50438, it is required that the change of 100% reactive power corresponds to the change of 10% grid nominal voltage [32]. Thus, the voltage droop coefficient, D q , is represented as: Commonly, the frequency droop coefficient D p and the voltage droop coefficient D q are chosen according to the application requirements, so that VSG can meet the required steady-state frequency and voltage droop characteristics and once set, generally no longer needs to be modified [31].
Denote the phase angle difference between the inner electric potential e and the grid voltage u g by δ, which obeys: where ω g is the angular speed of the grid voltage.

Virtual Impedance
The voltage amplitude reference generated by the reactive power loop in Figure 1 is passed through a virtual inductor before being used as a reference to control the capacitance voltage u t . To make the impedance inductive, a virtual inductor is introduced to be in series with the grid impedance through the control strategy [33]. The virtual impedance can be defined as follows: where L v is the virtual inductance. Then, the u tdqref can be calculated as:

The Inner Current Controller and Voltage Controller
The current controller and voltage controller consist of the classical PI controllers on the dq reference frame in Figure 1 [31]. Let K pc and K ic denote the proportional coefficient and integral coefficient of the PI in the current controller, respectively, and let K pv and K iv denote the proportional coefficient and integral coefficient of the PI in the voltage controller, respectively. The state equations and output equations of the current controller and the voltage controller are [34]: where ϕ d and ϕ q are the state variables of PI in the voltage controller, and γ d and γ q are the state variables of PI in the current controller.

LCL Filter and Coupling Inductance
According to Figure 1, the dynamic model of the LCL filter and coupling inductance can be expressed by the state equations on the dq reference frame as follows: where u sdq and u gdq can be calculated as:

State-Space Modeling of the Overall System
In this section, the full-order state-space small-signal model of the LCL grid-connected VSG system is developed. The dynamic behavior of the VSG-connected system is fully described by (11), and (12). The steady-state operating point can be calculated through the dynamic equations. The small-signal state-space model is obtained by linearizing these dynamic equations around the equilibrium point x o .
where the state vector ∆x and input vector ∆u are: When the state variable appears in the matrices, the values are denoted by subscript "o". For convenience of notation, the state matrix A is represented through a 15 × 8 sub-matrix A 11 and a 15 × 7 sub-matrix A 12 according to: The state matrix A and the input matrix B are depicted in Appendix A. The dynamic of the linearized model developed in Equation (14) is compared with that of the nonlinear model in MATLAB/SIMULINK in the next section to validate the effectiveness of the state-space model.

State-space Model Validation
The small-signal state-space model is validated in this section using a similar method to that in [4]. P f and Q f are the output variables in the output vector y = [P f , Q f ] T . Let ∆y denote the sufficiently small perturbations in y: where ∆x and ∆u are as described in Equation (15). The output matrix C and the feedforward matrix D are depicted in Appendix A.
Firstly, the linearized model consists of Equations (14) and (17) and is obtained by linearizing the nonlinear system around the equilibrium point corresponding to P set = 3000 W, Q set = 0 Var. The nonlinear grid-connected VSG system is modeled in MAT-LAB/SIMULINK. Both models use the parameters reported in Table 2. Time-domain simulations are conducted and the dynamic responses of the two models are compared. Figure 3a shows the simulation results of P f . The dynamic responses of the linearized and nonlinear models are compared. The active reference P set changes from 0 to 3000 W. It can be seen from Figure 3a that the dynamic responses of the linearized and nonlinear models match well. Secondly, the linearized model consists of Equations (14) and (17) and is obtained by linearizing the nonlinear system around the equilibrium point corresponding to P set = 0 W, Q set = 500 Var. Figure 3b shows the simulation results of the Q f . The active reference Q set changes from 0 to 500 Var. It can be seen from Figure 3b that the dynamic responses of the linearized and nonlinear models match well. The dynamic behavior of the nonlinear model is sufficiently emulated by the linearized model formed by Equations (14) and (17), and thus validate the effectiveness of the linearized model. Thirdly, the emulation of virtual inertia of the linearized model is verified. The power regulation process with the variation of J when grid frequency steps up by 0.1Hz is demonstrated in Figure 3c. It can be seen that the overshoot is small in active power regulation process when J = 0.01, and the adjusting time and overshoot become larger with the increase of J. Thus, the impact of the inertia on the grid-connected active power is validated, which is consistent with the conclusion in [35].
W, Qset = 0 Var. The nonlinear grid-connected VSG system is modeled in MATLAB/SIM-ULINK. Both models use the parameters reported in Table 2. Time-domain simulations are conducted and the dynamic responses of the two models are compared. Figure 3a shows the simulation results of Pf. The dynamic responses of the linearized and nonlinear models are compared. The active reference Pset changes from 0 to 3000 W. It can be seen from Figure 3a that the dynamic responses of the linearized and nonlinear models match well. Secondly, the linearized model consists of Equations (14) and (17) and is obtained by linearizing the nonlinear system around the equilibrium point corresponding to Pset = 0 W, Qset = 500 Var. Figure 3b shows the simulation results of the Qf. The active reference Qset changes from 0 to 500 Var. It can be seen from Figure 3b that the dynamic responses of the linearized and nonlinear models match well. The dynamic behavior of the nonlinear model is sufficiently emulated by the linearized model formed by Equations (14) and (17), and thus validate the effectiveness of the linearized model. Thirdly, the emulation of virtual inertia of the linearized model is verified. The power regulation process with the variation of J when grid frequency steps up by 0.1Hz is demonstrated in Figure 3c. It can be seen that the overshoot is small in active power regulation process when J = 0.01, and the adjusting time and overshoot become larger with the increase of J. Thus, the impact of the inertia on the grid-connected active power is validated, which is consistent with the conclusion in [35].

Analysis of Harmonic Resonance Characteristics
In Section 2, the full-order small-signal state-space model of the VSG-connected system is developed. In this section, on the basis of the small-signal model, the participation analysis, the eigenvalue locus, and the singular value decomposition analysis are conducted to study the impacts of the parameters on the harmonic resonance characteristics. The parameter values in this section are shown in Table 2 unless otherwise noted. The active power reference is Pset = 3000 W and the reactive power reference is Qset = 0 Var.

Participation Analysis
Participation analysis is used to uncover the contribution of each component on the small-signal stability [36]. The participation factor hik of the ith state variable and the kth eigenvalue [36] is defined as: where pik is the left eigenvalue element of the left eigenvector matrix and qik is the right eigenvalue element of the right eigenvector matrix. The matrix A in this paper has 15 eigenvalues. Table 3 shows the participation analysis results for harmonic eigenvalues. It can be seen from Table 3 that the state utd, utq, ild, ilq, igd, and igq are closely related to the high-frequency eigenvalues (λ1-6), while the lowfrequency eigenvalues (λ7-15) are not insensitive to those states. It is worth noting that the low-frequency oscillation analyzed by the method in this paper is consistent with the conclusions of other papers [4,37] and not described. This paper focuses on the harmonic res-

Analysis of Harmonic Resonance Characteristics
In Section 2, the full-order small-signal state-space model of the VSG-connected system is developed. In this section, on the basis of the small-signal model, the participation analysis, the eigenvalue locus, and the singular value decomposition analysis are conducted to study the impacts of the parameters on the harmonic resonance characteristics. The parameter values in this section are shown in Table 2 unless otherwise noted. The active power reference is P set = 3000 W and the reactive power reference is Q set = 0 Var.

Participation Analysis
Participation analysis is used to uncover the contribution of each component on the small-signal stability [36]. The participation factor h ik of the ith state variable and the kth eigenvalue [36] is defined as: where p ik is the left eigenvalue element of the left eigenvector matrix and q ik is the right eigenvalue element of the right eigenvector matrix. The matrix A in this paper has 15 eigenvalues. Table 3 shows the participation analysis results for harmonic eigenvalues. It can be seen from Table 3 that the state u td , u tq , i ld , i lq , i gd , and i gq are closely related to the high-frequency eigenvalues (λ 1-6 ), while the lowfrequency eigenvalues (λ 7-15 ) are not insensitive to those states. It is worth noting that the low-frequency oscillation analyzed by the method in this paper is consistent with the conclusions of other papers [4,37] and not described. This paper focuses on the harmonic resonance phenomena related to the inner loops controller parameters, where oscillations arise from 100 Hz to the LCL filter resonance frequency. Thus, the impacts of the inner loop controllers on λ 1-6 are analyzed in the further analysis below.

Singular Value Decomposition and Eigenvalue-Based Stability Analysis
Singular value composition is an effective way to identify the contribution of each component on the harmonic resonance characteristics. Let us consider a constant m × n complex matrix M, and then M can be decomposed into its SVD as shown in [26]. Besides, M can be expressed by the column vector elements of U and V [26], i.e., u i is the output singular vector and v i is the input singular vector. If the direction of an input vector is v i , the direction of the output vector is u i , and the corresponding gain of the vector length is σ i . The maximum gain in any input direction is equal to the maximum singular value σ 1 . Thus, Bodeplots of the max singular value can investigate the dynamic of the system [26]. For multiple input multiple output systems (MIMO), SVD can locate the frequencies at which the system is prone to be dynamically unstable [38][39][40]. The eigenvalue-based analysis is performed first to investigate the impact of the inner controller parameters on the harmonic resonances characteristics. The values of parameters are listed in Table 2, except the variables in each figure. Note that the minimum values of K pc and K pv are obtained through the eigenvalue-based analysis. The maximum values of K pc and K pv and the ranges of K ic and K iv can be calculated on the basis of the method in [41]. The range of L v can be obtained according to [25]. Next, the impacts of K pc and K pv on the LCL filter harmonic resonance were analyzed for a very low, medium, and a large value of K ic and K iv , respectively, while the impacts of K ic and K iv on the LCL filter harmonic resonance are analyzed for a very low, medium, and a large value of K pc and K pv , respectively. Figure 4 shows the harmonic-frequency eigenvalue trajectories and frequency response of the maximum singular value (SV) obtained when K pc decreases with different K ic . As shown in Figure 4(a1,b1,c1), the system is unstable when k pc is lower than 0.15, 0.2, and 1.4, respectively. It can be seen from Figure 4(a1,b1,c1) that four conjugate pairs are dominant in the harmonic-frequency oscillation characteristics. When K pc decreases, the harmonic-frequency eigenvalues (λ 1-4 ) move towards the unstable region. Based on the eigenvalue trajectories analysis, the SVD analysis is further performed in Figure 4(a2,b2,c2), respectively. The harmonic resonance around the LCL filter resonance frequency keeps amplifying with a decreases frequency as K pc decreases, which indicates that K pc has an essential effect on harmonic resonance instability.
Sustainability 2021, 13, x FOR PEER REVIEW 11 of 29 Figure 4 shows the harmonic-frequency eigenvalue trajectories and frequency response of the maximum singular value (SV) obtained when Kpc decreases with different Kic. As shown in Figure 4a1,b1,c1, the system is unstable when kpc is lower than 0.15, 0.2, and 1.4, respectively. It can be seen from Figure 4a1,b1,c1 that four conjugate pairs are dominant in the harmonic-frequency oscillation characteristics. When Kpc decreases, the harmonic-frequency eigenvalues (λ1-4) move towards the unstable region. Based on the eigenvalue trajectories analysis, the SVD analysis is further performed in Figure 4a2,b2,c2, respectively. The harmonic resonance around the LCL filter resonance frequency keeps amplifying with a decreases frequency as Kpc decreases, which indicates that Kpc has an essential effect on harmonic resonance instability.   Figure 5 shows the harmonic-frequency eigenvalue trajectories and frequency response of the maximum SV obtained when K ic decreases with different K pc . It can be seen that both the four harmonic-frequency eigenvalues and the LCL filter harmonic resonance around the LCL filter resonance frequency almost remain the same, which means that K ic has little influence on the harmonic resonance.
Sustainability 2021, 13, x FOR PEER REVIEW 12 of 29 Figure 5 shows the harmonic-frequency eigenvalue trajectories and frequency response of the maximum SV obtained when Kic decreases with different Kpc. It can be seen that both the four harmonic-frequency eigenvalues and the LCL filter harmonic resonance around the LCL filter resonance frequency almost remain the same, which means that Kic has little influence on the harmonic resonance.  Figure 6 depicts the harmonic-frequency eigenvalue trajectories and frequency response of the maximum SV obtained when the proportional gain of voltage controllers is in a range from 0.005 to 0.8 with different Kiv. As shown in Figure 6a1,b1,c1, the system is unstable when kpv is lower than 0.015, 0.17, and 0.35, respectively, and λ1-4 moves towards the virtual axis when Kpv decreases. It is worth noting that the peak of the harmonic resonance around the LCL filter resonance frequency keeps amplifying with the frequency decreases, as shown in Figure 6a2,b2,c2.  Figure 6 depicts the harmonic-frequency eigenvalue trajectories and frequency response of the maximum SV obtained when the proportional gain of voltage controllers is in a range from 0.005 to 0.8 with different K iv . As shown in Figure 6(a1,b1,c1), the system is unstable when k pv is lower than 0.015, 0.17, and 0.35, respectively, and λ 1-4 moves towards the virtual axis when K pv decreases. It is worth noting that the peak of the harmonic resonance around the LCL filter resonance frequency keeps amplifying with the frequency decreases, as shown in Figure 6(a2,b2,c2).  Figure 7 shows the harmonic-frequency eigenvalue trajectories and frequency response of the maximum SV obtained when Kiv decreases with different Kpv. Figure 8 shows the harmonic-frequency eigenvalue trajectories and the frequency response of the maximum SV obtained when LV varies. It can be seen that both the LCL filter harmonic resonance around the LCL filter resonance frequency and the four harmonic-frequency eigenvalues are insensitive to Kiv and LV.  Figure 7 shows the harmonic-frequency eigenvalue trajectories and frequency response of the maximum SV obtained when K iv decreases with different K pv . Figure 8 shows the harmonic-frequency eigenvalue trajectories and the frequency response of the maximum SV obtained when L V varies. It can be seen that both the LCL filter harmonic resonance around the LCL filter resonance frequency and the four harmonic-frequency eigenvalues are insensitive to K iv and L V.

Simulation and Experimental Verification
In this section, the impacts of the proportional gains on the harmonic resonance characteristics are validated via time-domain simulations and experiments.

Simulation Verification
The grid-connected VSG system is modeled in MATLAB/SIMULINK. The parameter values are listed in Table 2 unless otherwise noted. The active power reference is P set = 1500 W and the reactive power reference is Q set = 0 Var. Moreover, the frequency of the LCL filter resonance is 831 Hz, according to Table 2. Simulations with different parameters are performed and the key simulation results are given below.

Simulation Results in Unstable Case
Figures 9 and 10 depict the voltage waveform when the proportional gains are in an unstable case. Figure 9 depicts the voltage waveform of v pcc with K pc = 0.05 and K ic = 3. As shown in Figure 9, the voltage is unstable when K pc is lower than 0.2, which agrees with the theoretical analysis in Figure 4(b1) and Figure 4(b2). Figure 10 shows the voltage waveform of v pcc with K pv = 0.01 and K iv = 1. It can be seen from Figure 10 that the voltage is unstable when K pv is lower than 0.17, which agrees with the theoretical analysis in Figure 6(b1,b2).

Simulation and Experimental Verification
In this section, the impacts of the proportional gains on the harmonic resonance characteristics are validated via time-domain simulations and experiments.

Simulation Verification
The grid-connected VSG system is modeled in MATLAB/SIMULINK. The parameter values are listed in Table 2 unless otherwise noted. The active power reference is Pset = 1500 W and the reactive power reference is Qset = 0 Var. Moreover, the frequency of the LCL filter resonance is 831 Hz, according to Table 2. Simulations with different parameters are performed and the key simulation results are given below.

Simulation Results in Unstable Case
Figures 9 and 10 depict the voltage waveform when the proportional gains are in an unstable case. Figure 9 depicts the voltage waveform of vpcc with Kpc = 0.05 and Kic = 3. As shown in Figure 9, the voltage is unstable when Kpc is lower than 0.2, which agrees with the theoretical analysis in Figure 4b1 and Figure 4b2. Figure 10 shows the voltage waveform of vpcc with Kpv = 0.01 and Kiv = 1. It can be seen from Figure 10 that the voltage is unstable when Kpv is lower than 0.17, which agrees with the theoretical analysis in Figures 6b1 and 6b2.    Table 4. Figure 12 depicts the voltage waveform of vpcc when Kpc = 5 and Kpv = 0.2, while the HDF of Figure 12 is shown in Table 5.

Simulation and Experimental Verification
In this section, the impacts of the proportional gains on the harmonic resonance characteristics are validated via time-domain simulations and experiments.

Simulation Verification
The grid-connected VSG system is modeled in MATLAB/SIMULINK. The parameter values are listed in Table 2 unless otherwise noted. The active power reference is Pset = 1500 W and the reactive power reference is Qset = 0 Var. Moreover, the frequency of the LCL filter resonance is 831 Hz, according to Table 2. Simulations with different parameters are performed and the key simulation results are given below.

Simulation Results in Unstable Case
Figures 9 and 10 depict the voltage waveform when the proportional gains are in an unstable case. Figure 9 depicts the voltage waveform of vpcc with Kpc = 0.05 and Kic = 3. As shown in Figure 9, the voltage is unstable when Kpc is lower than 0.2, which agrees with the theoretical analysis in Figure 4b1 and Figure 4b2. Figure 10 shows the voltage waveform of vpcc with Kpv = 0.01 and Kiv = 1. It can be seen from Figure 10 that the voltage is unstable when Kpv is lower than 0.17, which agrees with the theoretical analysis in Figures 6b1 and 6b2.    Table 4. Figure 12 depicts the voltage waveform of vpcc when Kpc = 5 and Kpv = 0.2, while the HDF of Figure 12 is shown in Table 5.  Figure 11 shows the voltage waveform of v pcc when K pc = 0.4 and K pv = 0.6, of which the harmonic distortion factor (HDF) is shown in Table 4. Figure 12 depicts the voltage waveform of v pcc when K pc = 5 and K pv = 0.2, while the HDF of Figure 12 is shown in Table 5. Increasing K pc to 5 and K pv to 0.6, the voltage waveform of v pcc is shown in Figure 13, and the distribution of harmonic content is shown in Table 6. As shown in Figure 13, v pcc is stabilized after increasing the proportional gains of the current and voltage controllers. Increasing Kpc to 5 and Kpv to 0.6, the voltage waveform of vpcc is shown in Figure 13, and the distribution of harmonic content is shown in Table 6. As shown in Figure 13, vpcc is stabilized after increasing the proportional gains of the current and voltage controllers.      Increasing Kpc to 5 and Kpv to 0.6, the voltage waveform of vpcc is shown in Figure 13, and the distribution of harmonic content is shown in Table 6. As shown in Figure 13, vpcc is stabilized after increasing the proportional gains of the current and voltage controllers.          Figure 14 shows the comparison result of the HDF in Tables 4 and 6 with different proportional gains of the current controller. It can be seen from Figure 14 that the maximum harmonic content is around the LCL filter resonance frequency when K pc = 0.4. The HDF around the LCL filter resonance frequency increases when K pc decreases, which agrees with the analytical results in Figure 4.   Figure 14 shows the comparison result of the HDF in Tables 4 and 6 with different proportional gains of the current controller. It can be seen from Figure 14 that the maximum harmonic content is around the LCL filter resonance frequency when Kpc = 0.4. The HDF around the LCL filter resonance frequency increases when Kpc decreases, which agrees with the analytical results in Figure 4.   Tables 4 and 6 when Kpv = 0.6. Figure 15 shows the comparison result of the HDF in Tables 5 and 6 with different proportional gains of the voltage controller. As shown in Figure 15, the maximum harmonic content is around the LCL filter resonance frequency when Kpv = 0.2. The HDF around the LCL filter resonance increases when Kpv decreases, which agrees with the analytical results in Figure 6.  Tables 4 and 6 when K pv = 0.6. Figure 15 shows the comparison result of the HDF in Tables 5 and 6 with different proportional gains of the voltage controller. As shown in Figure 15, the maximum harmonic content is around the LCL filter resonance frequency when K pv = 0.2. The HDF around the LCL filter resonance increases when K pv decreases, which agrees with the analytical results in Figure 6.

Experimental Verification
In this section, the analysis is verified experimentally via the experimental setup. The diagram of the experimental verification platform is shown in Figure 16. The same configuration and parameters as the grid-connected system given in Figure 1 and Table 2 are used, unless otherwise noted. Because of safety considerations, the active power reference Pset is set to 1500 W, while the reactive power reference Qset is set to 0 Var. A three-phase two-level voltage source inverter with LCL filter instantiates VSG. Furthermore, the LCL filter resonance frequency was 831 Hz. The control strategy is implemented in the dSPACE DS1103 processor board. The Chroma 61,830 grid simulator emulates the grid voltage. The voltage waveform is measured using the Tektronic MDO 3034 oscilloscope. The distribution of harmonic content is analyzed through Yokogawa WT1804E power analyzer. Being similar to the simulations conducted in Section 4.1, the experiments described in this section are conducted with different parameters. Note that vab is the line-to-line voltage of the common coupling point. Figure 17 shows the voltage waveform of vab when Kpc = 0.4 and Kpv = 0.6, of which the HDF is shown in Table 7. Figure 18 depicts the voltage waveform of vab when Kpc = 5 and Kpv = 0.2, while the HDF is shown in Table 8. Increasing Kpc to 5 and Kpv to 0.6, the voltage waveform is shown in Figure 19 and the distribution of harmonic content is shown in Table 9. As shown in Figure 19, vab is stabilized after increasing the proportional gains of the current and voltage controllers, which is consistent with the results in the simulation.

Experimental Verification
In this section, the analysis is verified experimentally via the experimental setup. The diagram of the experimental verification platform is shown in Figure 16. The same configuration and parameters as the grid-connected system given in Figure 1 and Table  2 are used, unless otherwise noted. Because of safety considerations, the active power reference P set is set to 1500 W, while the reactive power reference Q set is set to 0 Var. A three-phase two-level voltage source inverter with LCL filter instantiates VSG. Furthermore, the LCL filter resonance frequency was 831 Hz. The control strategy is implemented in the dSPACE DS1103 processor board. The Chroma 61,830 grid simulator emulates the grid voltage. The voltage waveform is measured using the Tektronic MDO 3034 oscilloscope. The distribution of harmonic content is analyzed through Yokogawa WT1804E power analyzer.

Experimental Verification
In this section, the analysis is verified experimentally via the experimental setup. The diagram of the experimental verification platform is shown in Figure 16. The same configuration and parameters as the grid-connected system given in Figure 1 and Table 2 are used, unless otherwise noted. Because of safety considerations, the active power reference Pset is set to 1500 W, while the reactive power reference Qset is set to 0 Var. A three-phase two-level voltage source inverter with LCL filter instantiates VSG. Furthermore, the LCL filter resonance frequency was 831 Hz. The control strategy is implemented in the dSPACE DS1103 processor board. The Chroma 61,830 grid simulator emulates the grid voltage. The voltage waveform is measured using the Tektronic MDO 3034 oscilloscope. The distribution of harmonic content is analyzed through Yokogawa WT1804E power analyzer. Being similar to the simulations conducted in Section 4.1, the experiments described in this section are conducted with different parameters. Note that vab is the line-to-line voltage of the common coupling point. Figure 17 shows the voltage waveform of vab when Kpc = 0.4 and Kpv = 0.6, of which the HDF is shown in Table 7. Figure 18 depicts the voltage waveform of vab when Kpc = 5 and Kpv = 0.2, while the HDF is shown in Table 8. Increasing Kpc to 5 and Kpv to 0.6, the voltage waveform is shown in Figure 19 and the distribution of harmonic content is shown in Table 9. As shown in Figure 19, vab is stabilized after increasing the proportional gains of the current and voltage controllers, which is consistent with the results in the simulation. Being similar to the simulations conducted in Section 4.1, the experiments described in this section are conducted with different parameters. Note that v ab is the line-to-line voltage of the common coupling point. Figure 17 shows the voltage waveform of v ab when K pc = 0.4 and K pv = 0.6, of which the HDF is shown in Table 7. Figure 18 depicts the voltage waveform of v ab when K pc = 5 and K pv = 0.2, while the HDF is shown in Table 8. Increasing K pc to 5 and K pv to 0.6, the voltage waveform is shown in Figure 19 and the distribution of harmonic content is shown in Table 9. As shown in Figure 19, v ab is stabilized after increasing the proportional gains of the current and voltage controllers, which is consistent with the results in the simulation.             Figure 20 shows the comparison result of the HDF in Tables 7 and 9 with different proportional gains of the current controller. As can be seen from Figure 20, the maximum harmonic content is around the LCL filter resonance frequency when K pc = 0.4. The HDF around the LCL filter resonance frequency increases when K pc decreases, which agrees with the results in Section 4.1.   Figure 20 shows the comparison result of the HDF in Tables 7 and 9 with different proportional gains of the current controller. As can be seen from Figure 20, the maximum harmonic content is around the LCL filter resonance frequency when Kpc = 0.4. The HDF around the LCL filter resonance frequency increases when Kpc decreases, which agrees with the results in Section 4.1.   Tables 7 and 9 when Kpv = 0.6. Figure 21 shows the comparison result of the HDF in Tables 8 and 9 with different proportional gains of the voltage controller. As shown in Figure 21, the maximum harmonic content is around the LCL filter resonance frequency when Kpv = 0.2. The HDF around the LCL filter resonance frequency increases when Kpv decreases, which agrees with the simulation results in Section 4.2.  Tables 7 and 9 when K pv = 0.6. Figure 21 shows the comparison result of the HDF in Tables 8 and 9 with different proportional gains of the voltage controller. As shown in Figure 21, the maximum harmonic content is around the LCL filter resonance frequency when K pv = 0.2. The HDF around the LCL filter resonance frequency increases when K pv decreases, which agrees with the simulation results in Section 4.  Tables 8 and 9 when Kpc = 5.

Results Comparison
Further, the simulation results and the experimental results are compared in this section to verify the consistency of them. Figure 22 depicts the comparison result of the simulation and experiment when Kpc = 0.4 and Kpv = 0.6. Figure 23 depicts the comparison result of the simulation and experiment when Kpc = 5 and Kpv = 0.2. It can be seen from Figure  22 that when Kpc is 0.4, the peak difference between the simulation and experiment is 0.31%, and the frequency difference corresponding to the peak value in the simulation and experiment is 50 Hz. As shown in Figure 23, the peak values of the simulation and experiment are consistent, and the frequency difference corresponding to the peak value is 50 Hz.   Tables 8 and 9 when K pc = 5.

Results Comparison
Further, the simulation results and the experimental results are compared in this section to verify the consistency of them. Figure 22 depicts the comparison result of the simulation and experiment when K pc = 0.4 and K pv = 0.6. Figure 23 depicts the comparison result of the simulation and experiment when K pc = 5 and K pv = 0.2. It can be seen from Figure 22 that when K pc is 0.4, the peak difference between the simulation and experiment is 0.31%, and the frequency difference corresponding to the peak value in the simulation and experiment is 50 Hz. As shown in Figure 23, the peak values of the simulation and experiment are consistent, and the frequency difference corresponding to the peak value is 50 Hz.  Tables 8 and 9 when Kpc = 5.

Results Comparison
Further, the simulation results and the experimental results are compared in this section to verify the consistency of them. Figure 22 depicts the comparison result of the simulation and experiment when Kpc = 0.4 and Kpv = 0.6. Figure 23 depicts the comparison result of the simulation and experiment when Kpc = 5 and Kpv = 0.2. It can be seen from Figure  22 that when Kpc is 0.4, the peak difference between the simulation and experiment is 0.31%, and the frequency difference corresponding to the peak value in the simulation and experiment is 50 Hz. As shown in Figure 23, the peak values of the simulation and experiment are consistent, and the frequency difference corresponding to the peak value is 50 Hz.   Tables 8 and 9 when Kpc = 5.

Results Comparison
Further, the simulation results and the experimental results are compared in this section to verify the consistency of them. Figure 22 depicts the comparison result of the simulation and experiment when Kpc = 0.4 and Kpv = 0.6. Figure 23 depicts the comparison result of the simulation and experiment when Kpc = 5 and Kpv = 0.2. It can be seen from Figure  22 that when Kpc is 0.4, the peak difference between the simulation and experiment is 0.31%, and the frequency difference corresponding to the peak value in the simulation and experiment is 50 Hz. As shown in Figure 23, the peak values of the simulation and experiment are consistent, and the frequency difference corresponding to the peak value is 50 Hz.  The impact of line stray inductance as well as the internal impedance of the Chroma grid simulator on the difference of the simulation and experimental results is further analyzed. Note that the change in inductance of the line and Chroma grid simulator is equivalent to the variation in L g . Figures 24 and 25 show the resultant maximum SV frequency responses following the increase of L g . It can be seen that the harmonic resonance around the LCL filter resonance frequency varies with the decreases frequency as L g increases. Let R g denote the total resistance of the line and the Chroma grid simulator. In Figures 26 and 27, the frequency responses of the maximum SV with increasing R g are shown. As shown in Figure 26, increasing R g reduces the magnitude of the harmonic resonance around the LCL filter resonance frequency.
Based on the analysis above, the difference of the simulation and experiment is probably caused by the line stray inductance and the impedance of Chroma grid simulator. At the same time, the frequency difference between the simulation and experimental results is 50 Hz, and the accuracy of the power analyzer is 50 Hz, which means that the actual error is within 50 Hz. By further considering these factors, the difference between the simulation and experimental results is within a reasonable and acceptable range.
In brief, the experimental conclusions are in accord with the simulation conclusions in Section 4.1, and they validate the impact of the proportional gains on the harmonic resonance characteristics. Both the time-domain simulation and experimental results agrees with the analysis results from SVD and eigenvalue analysis. The experimental results, together with analysis results shown in Figures 4-7 and simulation results shown in Figures 14 and 15, demonstrate that decreasing the proportional gains of the inner loops causes the harmonic resonance oscillations around the LCL filter resonance frequency. The impact of line stray inductance as well as the internal impedance of the Chroma grid simulator on the difference of the simulation and experimental results is further analyzed. Note that the change in inductance of the line and Chroma grid simulator is equivalent to the variation in Lg. Figures 24 and 25 show the resultant maximum SV frequency responses following the increase of Lg. It can be seen that the harmonic resonance around the LCL filter resonance frequency varies with the decreases frequency as Lg increases. Let Rg denote the total resistance of the line and the Chroma grid simulator. In Figures 26 and  27, the frequency responses of the maximum SV with increasing Rg are shown. As shown in Figure 26, increasing Rg reduces the magnitude of the harmonic resonance around the LCL filter resonance frequency.
L g increases L g increases  The impact of line stray inductance as well as the internal impedance of the Chroma grid simulator on the difference of the simulation and experimental results is further analyzed. Note that the change in inductance of the line and Chroma grid simulator is equivalent to the variation in Lg. Figures 24 and 25 show the resultant maximum SV frequency responses following the increase of Lg. It can be seen that the harmonic resonance around the LCL filter resonance frequency varies with the decreases frequency as Lg increases. Let Rg denote the total resistance of the line and the Chroma grid simulator. In Figures 26 and  27, the frequency responses of the maximum SV with increasing Rg are shown. As shown in Figure 26, increasing Rg reduces the magnitude of the harmonic resonance around the LCL filter resonance frequency.
L g increases   Based on the analysis above, the difference of the simulation and experiment is probably caused by the line stray inductance and the impedance of Chroma grid simulator. At the same time, the frequency difference between the simulation and experimental results is 50 Hz, and the accuracy of the power analyzer is 50 Hz, which means that the actual error is within 50 Hz. By further considering these factors, the difference between the simulation and experimental results is within a reasonable and acceptable range.
In brief, the experimental conclusions are in accord with the simulation conclusions in Section 4.1, and they validate the impact of the proportional gains on the harmonic resonance characteristics. Both the time-domain simulation and experimental results agrees with the analysis results from SVD and eigenvalue analysis. The experimental results, together with analysis results shown in Figures 4-7 and simulation results shown in Figures 14 and 15, demonstrate that decreasing the proportional gains of the inner loops causes the harmonic resonance oscillations around the LCL filter resonance frequency.  Based on the analysis above, the difference of the simulation and experiment is probably caused by the line stray inductance and the impedance of Chroma grid simulator. At the same time, the frequency difference between the simulation and experimental results is 50 Hz, and the accuracy of the power analyzer is 50 Hz, which means that the actual error is within 50 Hz. By further considering these factors, the difference between the simulation and experimental results is within a reasonable and acceptable range.
In brief, the experimental conclusions are in accord with the simulation conclusions in Section 4.1, and they validate the impact of the proportional gains on the harmonic resonance characteristics. Both the time-domain simulation and experimental results agrees with the analysis results from SVD and eigenvalue analysis. The experimental results, together with analysis results shown in Figures 4-7 and simulation results shown in Figures 14 and 15, demonstrate that decreasing the proportional gains of the inner loops causes the harmonic resonance oscillations around the LCL filter resonance frequency.

Conclusions
In this paper, the harmonic resonance characteristics are analyzed in VSG using the full-order state-space model. First, the procedure to build the small-signal model step by step is summarized, and then, the state-space small-signal model that fully considers the dynamic of the current and voltage loops is built.
Secondly, the contribution of each state on harmonic instability is evaluated through participation analysis. The harmonic resonance instability is evaluated by performing the eigenvalue locus. Specifically, the SVD analysis is leveraged to study the harmonic resonance under varying parameters. The analysis results show that the harmonic resonance characteristics are mainly affected by proportional gains of the inner loops.
Thirdly, extensive numerical simulation and experimental results are given to verify the validity of the theoretical analysis. The simulation and experiments, both with different proportional gains of the inner loops, are respectively conducted. Both results indicate that the voltage of the common coupling point is unstable after decreasing the proportional gains of the current and voltage controllers. As K pc decreases from 5 to 0.4 or K pv decreases from 0.6 to 0.2, the HDF around the LCL filter frequency increases. Furthermore, the consistency of simulation and experimental results is validated.
In brief, the harmonic resonance characteristics around the LCL filter resonance frequency are analyzed to avoid the undesired distorted voltage of the common coupling point in this paper. The analytical methods and procedures are provided to guide the design and debugging of the inner loops parameters of VSG, which will further improve the practical engineering applicability of VSG.
Further work needs to be done to reveal the mechanism of LCL filter harmonic resonance instability when there are multiple VSGs in the system. Author Contributions: J.J. wrote the paper and carried out the simulation and experiments. W.W., X.W., F.T., and X.L. revised the whole paper. Z.Y. discussed and contributed to the writing of the paper. All authors have read and agreed to the published version of the manuscript.