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

Abstract

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 filter 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 filter resonance frequency increases. Furthermore, the consistency of simulation results, experimental results, and the theoretical analysis results is validated.

1. 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 switching-frequency 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 impedance-based 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 state-space 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. Performance comparison of different modeling methods.

	Transfer Function Model	The Impedance-Based Model	Model in Section 2
Identification of dynamic modes	−	−	+
Participation factor of the state variables	−	−	+
Scalability	−	+	−

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 et al. using SVD [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.

2. 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.

2.1. Overview of the Conventional VSG

Figure 1 shows a VSG-connected system, which consists of the topology of the grid-connected converter and the control diagram of VSG [29].

Table 2. Parameter values of gird-connected VSG system.

Parameters	Values	Parameters	Values
Rated voltage (phase-to-ground)	110 (V)	Rated angular speed value, ωn	314.159 (rad/s)
DC-link voltage, Udc	325 (V)	Time constant of the LPFs, τf	0.01s
Switching frequency	10 (kHz)	Filter inductance, Lc	8 (mH)
Rated capacity	3 (kW)	Filter capacitance, C	9.6 (μF)
Frequency droop coefficient, Dp	1.52	Voltage droop coefficient, Dq	96.4
Inertia, J	0.01	Integral coefficient, K	10
Proportional coefficient, Kpc	5	Integral coefficient, Kic	3
Proportional coefficient, Kpv	0.6	Integral coefficient, Kiv	1
Grid inductance, Lg	7.3 (mH)	Virtual inductance, Lv	6 (mH)

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.

2.2. The Dynamic Equations of Each Component

2.2.1. Power Calculation and Filtering

As depicted in Figure 1, the instantaneous active power p and reactive power q can be computed by [30]:

$$\begin{array}{l}p=\frac{3}{2}(u_{\mathrm{td}}{i}_{\mathrm{gd}}+u_{\mathrm{tq}}{i}_{\mathrm{gq}})\\ q=\frac{3}{2} (-u_{\mathrm{td}}{i}_{\mathrm{gq}}+u_{\mathrm{tq}}{i}_{\mathrm{gd}}) \end{array}$$

(1)

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]:

$$\begin{array}{l}{\tau }_{\mathrm{f}}\frac{d{P}_{\mathrm{f}}}{dt}=-P_{\mathrm{f}}+p\\ {\tau }_{\mathrm{f}}\frac{d{Q}_{\mathrm{f}}}{dt}=-Q_{\mathrm{f}}+q\end{array}$$

(2)

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.

2.2.2. 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]:

$$J{\omega }_{\mathrm{n}}\frac{\mathrm{d}\omega }{\mathrm{d}t}=P_{\mathrm{set}}-P_{\mathrm{f}}-D_{\mathrm{p}}{\omega }_{\mathrm{n}}(\omega -{\omega }_{\mathrm{n}})$$

(3)

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:

$$D_{\mathrm{p}}=\frac{\Delta T_{\max }}{\Delta {\omega }_{\max }}=\frac{\Delta P_{\max }}{{\omega }_{\mathrm{n}}\Delta {\omega }_{\max }}.$$

(4)

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:

$$\sqrt{2}K\frac{\mathrm{d}E}{\mathrm{d}t}=Q_{\mathrm{set}}-Q_{\mathrm{f}}-\sqrt{2}{D}_{\mathrm{q}}(U_{\mathrm{t}}-U_{\mathrm{n}})$$

(5)

where K is the integral coefficient that adjusts the dynamic of the RPL and U_t can be expressed as:

$$U_{\mathrm{t}}=\frac{\sqrt{{u_{\mathrm{td}}}^2+{u_{\mathrm{tq}}}^2}}{\sqrt{2}}.$$

(6)

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:

$$D_{\mathrm{q}}=\frac{\Delta Q_{\max }}{\sqrt{2}\Delta U_{\max }}.$$

(7)

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:

$$\frac{d\delta }{dt}=\omega -{\omega }_{\mathrm{g}}$$

(8)

where ω_g is the angular speed of the grid voltage.

2.2.3. 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:

$$Z_{\mathrm{v}}=L_{\mathrm{v}}s$$

(9)

where L_v is the virtual inductance. Then, the u_tdqref can be calculated as:

$$\begin{array}{l}{u}_{\mathrm{tdref}}=\sqrt{2}E+\omega L_{\mathrm{v}}{i}_{\mathrm{gq}}\\ u_{\mathrm{tqref}}=-\omega L_{\mathrm{v}}{i}_{\mathrm{gd}}\end{array}$$

(10)

2.2.4. 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]:

$$\begin{array}{l}\frac{d{\varphi }_{\mathrm{d}}}{dt}=u_{\mathrm{tdref}}-u_{\mathrm{td}}\\ \frac{d{\varphi }_{\mathrm{q}}}{dt}=u_{\mathrm{tqref}}-u_{\mathrm{tq}}\\ \frac{d{\gamma }_{\mathrm{d}}}{dt}=i_{\mathrm{Ldref}}-i_{\mathrm{Ld}}\\ \frac{d{\gamma }_{\mathrm{q}}}{dt}=i_{\mathrm{Lqref}}-i_{\mathrm{Lq}}\\ i_{\mathrm{Ldref}}=-\omega C{u}_{\mathrm{tq}}+K_{\mathrm{pv}}(u_{\mathrm{tdref}}-u_{\mathrm{td}})+K_{\mathrm{iv}}{\varphi }_{\mathrm{d}}\\ i_{\mathrm{Lqref}}=\omega C{u}_{\mathrm{td}}+K_{\mathrm{pv}}(u_{\mathrm{tqref}}-u_{\mathrm{tq}})+K_{\mathrm{iv}}{\varphi }_{\mathrm{q}}\end{array}$$

(11)

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.

2.2.5. 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:

$$\begin{array}{l}\frac{d{u}_{\mathrm{td}}}{dt}=\frac{1}{C}(\omega C{u}_{\mathrm{tq}}+i_{\mathrm{Ld}}-i_{\mathrm{gd}})\\ \frac{d{u}_{\mathrm{tq}}}{dt}=\frac{1}{C}(-\omega C{u}_{\mathrm{td}}+i_{\mathrm{Lq}}-i_{\mathrm{gq}})\\ \frac{d{i}_{\mathrm{Ld}}}{dt}=\frac{1}{L_{\mathrm{c}}}(\omega L_{\mathrm{c}}{i}_{\mathrm{Lq}}+u_{\mathrm{sd}}-u_{\mathrm{td}})\\ \frac{d{i}_{\mathrm{Lq}}}{dt}=\frac{1}{L_{\mathrm{c}}}(-\omega L_{\mathrm{c}}{i}_{\mathrm{Ld}}+u_{\mathrm{sq}}-u_{\mathrm{tq}})\\ \frac{d{i}_{\mathrm{gd}}}{dt}=\frac{1}{L_{\mathrm{g}}}(\omega L_{\mathrm{g}}{i}_{\mathrm{gq}}+u_{\mathrm{td}}-u_{\mathrm{gd}})\\ \frac{d{i}_{\mathrm{gq}}}{dt}=\frac{1}{L_{\mathrm{g}}}(-\omega L_{\mathrm{g}}{i}_{\mathrm{gd}}+u_{\mathrm{tq}}-u_{\mathrm{gq}})\end{array}$$

(12)

where u_sdq and u_gdq can be calculated as:

$$\begin{array}{l}{u}_{\mathrm{sd}}=u_{\mathrm{td}}-\omega L_{\mathrm{c}}{i}_{\mathrm{Lq}}+K_{\mathrm{pc}}(i_{\mathrm{Ldref}}-i_{\mathrm{Ld}})+K_{\mathrm{ic}}{\gamma }_{\mathrm{d}}\\ u_{\mathrm{sq}}=u_{\mathrm{tq}}+\omega L_{\mathrm{c}}{i}_{\mathrm{Ld}}+K_{\mathrm{pc}}(i_{\mathrm{Lqref}}-i_{\mathrm{Lq}})+K_{\mathrm{ic}}{\gamma }_{\mathrm{q}}\\ u_{\mathrm{gd}}=\sqrt{2}{U}_{\mathrm{g}}\cos \delta \\ u_{\mathrm{gq}}=-\sqrt{2}{U}_{\mathrm{g}}\sin \delta \end{array}$$

(13)

2.3. 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 Equations (2), (3), (5), (8), (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.

$$\frac{d\Delta x}{dt}=A\Delta x+{\rm B}\Delta u$$

(14)

where the state vector Δx and input vector Δu are:

$$\begin{array}{l}\Delta x=[\Delta P_{\mathrm{f}},\Delta Q_{\mathrm{f}},\Delta \omega ,\Delta E,\Delta \delta ,\Delta {\varphi }_{\mathrm{d}},\Delta {\varphi }_{\mathrm{q}},\\ \Delta {\gamma }_{\mathrm{d}},\Delta {\gamma }_{\mathrm{q}},\Delta u_{\mathrm{td}},\Delta u_{\mathrm{tq}},\Delta i_{\mathrm{Ld}},\Delta i_{\mathrm{Lq}},\Delta i_{\mathrm{gd}},\Delta i_{\mathrm{gq}}{]}^{\mathrm{T}}\\ \Delta u={[\Delta P_{\mathrm{set}},\Delta Q_{\mathrm{set}},\Delta {\omega }_{\mathrm{g}},\Delta U_{\mathrm{g}}]}^{\mathrm{T}}\end{array}$$

(15)

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₁₁ and a 15 × 7 sub-matrix A₁₂ according to:

$$\left[\begin{array}{c}\frac{d\Delta x_1}{dt}\\ \frac{d\Delta x_2}{dt}\end{array}\right]=\left[\begin{array}{cc}{A}_{11}& A_{12}\end{array}\right]\left[\begin{array}{c}\Delta x_1\\ \Delta x_2\end{array}\right]+B\Delta u$$

(16)

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.

2.4. 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:

$$\Delta y=C\Delta x+D\Delta u$$

(17)

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 MATLAB/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].

3. 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.

3.1. 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:

$$h_{\mathrm{ik}}=p_{\mathrm{ik}}{q}_{\mathrm{ik}}$$

(18)

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. Participation analysis for VSG.

	Value of λk, k = 1, …, 15	State	Participation
λ1	−231.7 + j7397.8	utd/utq	0.998/1
		ild/ilq	0.68/0.68
λ2	−231.7 − j7397.8	igd/igq	0.26/0.26
		Remaining states	<0.03
λ3	−220.9 + j7190.2	utd/utq	1/0.998
		ild/ilq	0.8/0.8
λ4	−220.9 − j7190.2	igd/igq	0.3/0.3
		Remaining states	<0.03
λ5	−169.2 + j418.7	ild/ilq	0.37/0.37
λ6	−169.2 − j418.7	igd/igq	1/1
		Remaining states	<0.03
λ7	−175.1	Pf/ω	0.3/1
		Remaining states	<0.03
λ8	−90.7	Qf/E	1/0.1
		Remaining states	<0.03
λ9	−62.3	Pf/ω	1/0.4
		Qf/E	0.22/0.16
		δ	0.56
		Remaining states	<0.03
λ10	−13 + j11	Pf/ω	0.19/0.12
		Qf/E	0.1/0.77
λ11	−13 − j11	δ	1
		Remaining states	<0.03
λ12	−1.6	φd/φq	1/0.99
λ13	−1.6	Remaining states	<0.03
λ14	−0.6	γd/γq	1/1
λ15	−0.6	Remaining states	<0.03

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 low-frequency 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.

3.2. 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 σ₁. 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.

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.

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 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.

4. 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.

4.1. 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.

4.1.1. Simulation Results in Unstable Case

Figure 9 and Figure 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).

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. Simulation result: distribution of harmonic content of v_ab.

Frequency	HDF (%)	Frequency	HDF (%)
450	0.09	500	0.25
550	0.09	600	0.09
650	0.82	700	2.36
750	0.5	800	0.2
850	0.14	900	0.17
950	0.02	1000	0.07
1050	0.02	1100	0.02

. 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. Simulation result: distribution of harmonic content of v_ab.

Frequency	HDF (%)	Frequency	HDF (%)
450	0.02	500	0.06
550	0.13	600	0.17
650	0.41	700	0.09
750	0.1	800	0.41
850	1.7	900	0.13
950	0.64	1000	0.04
1050	0.03	1100	0.02

. 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. Simulation result: distribution of harmonic content of v_ab.

Frequency	HDF (%)	Frequency	HDF (%)
450	0.04	500	0.07
550	0.13	600	0.11
650	0.05	700	0.02
750	0.03	800	0.11
850	0.22	900	0.21
950	0.07	1000	0.13
1050	0.19	1100	0.1

. As shown in Figure 13, v_pcc is stabilized after increasing the proportional gains of the current and voltage controllers.

4.1.2. LCL Filter Harmonic Resonance Associated with K_pc

Figure 14 shows the comparison result of the HDF in Table 4 and Table 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.

4.1.3. LCL Filter Harmonic Resonance Associated with K_pv

Figure 15 shows the comparison result of the HDF in Table 5 and Table 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.

4.2. 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.

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. Experimental result: distribution of harmonic content of v_ab.

Order	Frequency	HDF (%)	Order	Frequency	HDF (%)
9	450	0.72	10	500	0.69
11	550	0.31	12	600	0.78
13	650	0.35	14	700	1.26
15	750	2.67	16	800	1.93
17	850	0.85	18	900	0.33
19	950	0.51	20	1000	0.17
21	1050	0.10	22	1100	0.15

. 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. Experimental result: distribution of harmonic content of v_ab.

Order	Frequency	HDF (%)	Order	Frequency	HDF (%)
9	450	0.55	10	500	0.54
11	550	0.38	12	600	0.61
13	650	0.88	14	700	0.45
15	750	1.43	16	800	1.70
17	850	1.02	18	900	0.29
19	950	0.30	20	1000	0.15
21	1050	0.15	22	1100	0.20

. 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. Experimental result: distribution of harmonic content of v_ab.

Order	Frequency	HDF (%)	Order	Frequency	HDF (%)
9	450	0.21	10	500	0.23
11	550	0.26	12	600	0.39
13	650	0.27	14	700	0.15
15	750	0.12	16	800	0.25
17	850	0.17	18	900	0.09
19	950	0.05	20	1000	0.06
21	1050	0.06	22	1100	0.04

. 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.

4.2.1. LCL Filter Harmonic Resonance Associated with K_pc

Figure 20 shows the comparison result of the HDF in Table 7 and Table 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.

4.2.2. LCL Filter Harmonic Resonance Associated with K_pv

Figure 21 shows the comparison result of the HDF in Table 8 and Table 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.2.

4.3. 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.

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. Figure 24 and Figure 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 Figure 26 and Figure 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 Figure 4, Figure 5, Figure 6 and Figure 7 and simulation results shown in Figure 14 and Figure 15, demonstrate that decreasing the proportional gains of the inner loops causes the harmonic resonance oscillations around the LCL filter resonance frequency.

5. 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.