The DF has been widely used in nonlinear analysis and control in electronic power engineering (Vidal et al., 2001; Wang, 2017; Li et al., 2018; Liu et al., 2019; Xu et al., 2020; Wu et al., 2021). The DF of the nonlinear devices refers to the ratio of the fundamental component of the output signal to the input sinusoidal signal. It is essentially a method that uses the fundamental-frequency, linear characteristics to approximately describe the nonlinear characteristics of the devices. The nonlinear devices can be saturation limiter, dead zone, hysteresis loop, etc.

Generally, the DF can be expressed as N ( X ) = B 1 + j A 1 X (1) Where X denotes the amplitude of the input sinusoidal signal [x (t) = Xsin (ωt)], and B ₁+jA ₁ represents the output amplitude under the fundamental frequency (ω).

For an input sinusoidal signal, x (t) = Xsin (ωt), the output after the nonlinear device is usually non-sinusoidal. But we can express it by using the Fourier expansion as y ( t ) = A 0 2 + ∑ n = 1 ∞ [ A n ⁡ cos ( n ω t ) + B n ⁡ sin ( n ω t ) ]                     = A 0 2 + ∑ n = 1 ∞ Y n ⁡ sin ( n ω t + φ n ) (2) where Y n = A n 2 + B n 2 and φ n = a r c ⁡ tan ( A n / B n ) , and keep the fundamental-frequency terms only.

For the saturation limiter studied in the paper, its input and output waveforms are illustrated in Figure 1, and its nonlinear characteristics is described by y ( t ) = { a x ( t ) − a ( x ( t ) ≥ a ) ( − a ≤ x ( t ) < a ) ( x ( t ) < − a ) (3) where a denotes the limiter boundary.

Since the output signal is an odd function, A _i = 0 for all i = 0, 1, 2, …. The fundamental-frequency output signal can be simply expressed as y 1 ( t ) = 2 π [ X s i n − 1 a X + a 1 − ( a X ) 2 ] sin ( ω t ) (4)

Therefore, we obtain the DF of the saturation limiter explicitly N ( X ) = y 1 ( t ) x ( t ) = 2 π [ s i n − 1 a X + a X 1 − ( a X ) 2 ] ,    X ≥ a (5)

Here we can find that the DF is a real function, and we will see that this is convenient for the oscillation analysis.

First assume that the saturation limiter N(X), as the only one nonlinear component, can be separated from all linear components in the system [denoted by G(s) after linearization]. They form a unit negative feedback system, as shown in Figure 2, where letters r, e, and c represent the input, error, and output signals, respectively.

As the saturation limiter can be represented by a DF [N(X)], we obtain the transfer function of the closed-loop system c ( s , X ) r ( s , X ) = N ( X ) G ( s ) 1 + N ( X ) G ( s ) (6) With the closed-loop characteristic equation 1 + N ( X ) G ( s ) = 0 (7)

Therefore, the system stability should be determined by the relative position of the curve of N(X)G(s) and the point (−1, j0). As we know that the DF of the saturation limiter N(X) does not contain frequency variable s and G(s) does not contain any information of the saturation limiter, G(s) and N(X) can be separately analyzed. Based on the DFNC, we have.

1) If the curve of the transfer function G(s) surrounds the curve of −1/N(X), the closed-loop system is stable.

These three cases are schematically shown in Figures 3A–C, respectively. The dashed lines in Figures 3A,B represent the Nyquist curve of G(s) when s starts from j0⁻ and bypasses the origin to reach j0⁺. The arrow on the curve G(s) represents the increase direction of ω (s = jω). The arrow on the horizontal line N(X) represents the direction with the increase of X. Note that as the DF of the saturation limiter is always a positive, real function, −1/N(X) is exactly located on the negative real axis and it moves left with the increase of X, as shown in Figure 3.

Consequently, based on the intersection of these two functions [G(s) and −1/N(X)], we can obtain its oscillation frequency from arg [ G ( s ) ] = − 180 ° (8) And the oscillation amplitude from | N ( X ) G ( s ) | = 1 (9)

The next step should be separating the saturation limiter from other linear parts for our specific grid-tied VSC system.

Now we attempt to study the grid-tied VSC with the ACC and the phase-locked loop (PLL) and focus on the two different types of sustained oscillation caused by the saturation limiter: the DCO and SCO. The structure of the whole system is shown in Figure 4. Without losing generality, two saturation limiters are assumed to install at the outputs of the d-axis and q-axis PI controllers of the ACC. Here different with the previous relevant works (Xu et al., 2020; Wu et al., 2021), the reference value of the q-axis current i ^* _qg is not fixed as zero, and the dynamics of both the d-axis and q-axis currents in the ACC will be studied. This treatment may better reflect the fact that saturation limiters are omnipresent in any realistic converters.

In the main circuit in Figure 4A, L _f and L _g represent the filter inductance and the equivalent grid inductance, respectively. In our study, the capacitor of the main circuit is regarded as a constant voltage source, i.e., U _dc is a constant. We use u _gabc to denote the three-phase voltage at an infinite grid, u _tabc the voltage at the point of common coupling, and e _abc the terminal voltage of the VSC. The parameters k _p,pll and k _i,pll are the proportional and integral coefficients of the PLL, respectively. k _p,acc and k _i,acc are the proportional and integral coefficients of the ACC, respectively. Two saturation limiters denoted by N(X) are installed at the output of the PI controller of the d-axis and the q-axis to limit the overcurrent. N _d,in and N _q,in are the input signal of the d-axis limiter and the q-axis limiter, respectively. As usual, the controls rely on the classical vector-control method, and all the controllers work in the phase-locked dq frame. The angle relation between the phase-locked dq and the constant synchronous xy frames is schematically shown in Figure 4B. Hence the two coordinates rotate at frequencies ω _pll and ω ₀, respectively, and their angle mismatch in the dynamical process is △θ _pll . Since the focus of this paper is mainly on the sub-synchronous and super-synchronous oscillation, we assume that the outer-loop voltage control at slower timescales and the signal modulation at faster timescales are ignored (Ma et al., 2020; Ma et al., 2021). Thus, the modulator can ideally output the internal voltage of the VSC (e _abc ).

The difference between the two phenomena, i.e., the DCO and the SCO, is clear. Their output waveforms are different. In the DCO, both limiters show a sustained oscillation, whereas in the SCO, only one limiter shows a sustained oscillation, and the output of the other limiter keeps constant.

First let us analyze the DCO. Based on the Kirchhoff Voltage Law, we obtain the following equations for the main circuit { e x − u x t = s L f i x g − ω 0 L f i y g e y − u y t = s L f i y g + ω 0 L f i x g (10) { u x t − u x g = s L g i x g − ω 0 L g i y g u y t − u y g = s L g i y g + ω 0 L g i x g (11) where the subscripts x and y represent the electrical quantity within the synchronous xy frame. ω ₀ is the frequency of the system at the steady state.

Since the VSC is tied to an infinite grid, u _xg and u _yg are constant. Their corresponding small-signal equations are { Δ e x − Δ u x t = s L f Δ i x g − ω 0 L f Δ i y g Δ e y − Δ u y t = s L f Δ i y g + ω 0 L f Δ i x g (12) { Δ u x t = s L g Δ i x g − ω 0 L g Δ i y g Δ u y t = s L g Δ i y g + ω 0 L g Δ i x g (13)

As shown in Figure 4A, to improve stability, we consider the voltage feedforward and cross decoupling links in the ACC. When the DCO occurs, the output voltage for either the d-axis or the q-axis is composed of three parts: the output of the limiter, the feedforward voltage, and the cross-decoupling term. Therefore, we obtain the small-signal model for the ACC { Δ e d = G a c c N ( X ) ( Δ i d g ∗ − Δ i d g ) − ω 0 L f Δ i q g + Δ u d t Δ e q = G a c c N ( X ) ( Δ i q g ∗ − Δ i q g ) + ω 0 L f Δ i d g + Δ u q t (14) where the subscripts d and q represent the electrical quantity in the PLL dq frame, G _acc (s) = k _p,acc + k _i,acc /s for the transfer function of the PI controller, and N(X) denotes the DF of the saturation limiter.

Based on the angle relation between the two different coordinates in Figure 4B, we have { Δ f d = Δ f x + Δ θ p l l f y 0 Δ f q = Δ f y − Δ θ p l l f x 0 (15) where f _x0 and f _y0 represent the electrical quantity (voltage or current) in the steady state. Therefore, we obtain the q-axis voltage at the common point Δ u q t = Δ u y t − Δ θ p l l U x t 0 (16)

As the input signal of the PLL is the q-axis voltage at the common point, the small-signal model for the PLL is given by Δ θ p l l = Δ u q t G p l l 1 s (17) where G _pll (s) = k _p,pll + k _i,pll /s for the PI controller of the PLL. Combining Eqs 15–17 yields Δ θ p l l = G p l l s + U x t 0 G p l l Δ u y t = H p l l ( s ) Δ u y t (18)

Based on the above small-signal models of the main circuit, ACC, PLL, and saturation limiters in Eqs 12–18, we obtain { Δ i x g = G a c c ( s ) N ( X ) G a c c ( s ) N ( X ) + s L f Δ i d g ∗ − G a c c ( s ) N ( X ) G a c c ( s ) N ( X ) + s L f I y g 0 H p l l Δ u y t Δ i y g = G a c c ( s ) N ( X ) G a c c ( s ) N ( X ) + s L f Δ i q g ∗ + G a c c ( s ) N ( X ) G a c c ( s ) N ( X ) + s L f I x g 0 H p l l Δ u y t (19)

Finally, substituting Eq. 13 into Eq. 19 and eliminating △u _yt yield { Δ i x g = G 11 ( s ) 1 + G 0 _ d ( s ) N ( X ) Δ i d g ∗ + G 12 ( s ) 1 + G 0 _ d ( s ) N ( X ) Δ i q g ∗ Δ i y g = G 21 ( s ) 1 + G 0 _ d ( s ) N ( X ) Δ i d g ∗ + G 22 ( s ) 1 + G 0 _ d ( s ) N ( X ) Δ i q g ∗ (20) Where G 11 ( s ) = G a c c H 1 N ( X ) { 1 + G a c c s L f N ( X ) [ 1 − H p l l s L g I x g 0 ] } G 12 ( s ) = − G a c c 2 s L f H p l l H 1 N 2 ( X ) s L g I y g 0 G 21 ( s ) = G a c c 2 s L f H p l l H 1 N 2 ( X ) ω 0 L g I x g 0 G 22 ( s ) = G a c c H 1 N ( X ) { 1 + G a c c s L f N ( X ) [ 1 + H p l l ω 0 L g I y g 0 ] } G 0 _ d ( s ) = G a c c s L f [ 1 + H p l l ( ω 0 L g I y g 0 − s L g I x g 0 ) ] H 1 = 1 G a c c ( s ) N ( X ) + s L f N ( X ) = 2 π [ s i n − 1 ( a X ) + a X 1 − ( a X ) 2 ]

These MIMO relations between [△i _xg , △i _yg ] and [△i ^* _dg , △i ^* _yg ] can be expressed in a compact form G ′ 1 _ d ( s ) [ Δ i d g ∗ Δ i q g ∗ ] = [ I + G ′ 0 _ d ( s ) N ′ ( X ) ] [ Δ i x g Δ i y g ] (21) Where G ′ 1 _ d ( s ) = [ G 11 ( s ) G 12 ( s ) G 21 ( s ) G 22 ( s ) ] G ′ 0 _ d ( s ) = [ G 0 _ d ( s ) 0 0 G 0 _ d ( s ) ] N ′ ( X ) = [ N ( X ) 0 0 N ( X ) ] I = [ 1 1 ]

Figure 5A shows the corresponding MIMO model for the DCO, with the saturation limiter and other linear links well characterized and separated. Observing that both G ^’ _{0_d} (s) and N′(X) are diagonal matrices, we can further simplify the model into two independent, identical SISO models with the same structure and parameters, and the system stability should be solely determined by. L ( s ) = G 0 _ d ( s ) N ( X ) (22)

In addition, expanding the first element G ₁₁(s) in G ^’ _{0_d} (s), we obtain that its denominator contains s ² + U _xt0 k _p,pll s + U _xt0 k _i,pll , L _f s ², and L _f s ² + [k _p,acc N(X)]s + k _i,acc N(X). Therefore, all these three denominator terms should have no any positive, real solution, and they do not intersect the positive real axis. Thus, G ₁₁(s) has no right-half-plane pole and it has no contribution to the system stability. Similarly, we can prove that the other three elements of the G ^’ _{0_d} (s) are similar. Therefore, we deduce that the system stability should be determined by G _{0_d} (s) and N(X) only.

First for the SCO we assume that the output of the d-axis limiter is constant, i.e., it is fixed at the limiter boundary, but the output of the q-axis limiter oscillates. Therefore, the only difference between the DCO and SCO is the output of the d-axis in the ACC. Under this situation, the output of the q-axis is the same, but the linearized output of the d-axis should be composed of the feedforward voltage and cross-decoupling terms only. Hence correspondingly the small-signal dynamics of the ACC can be expressed as { Δ e d = − ω 0 L f Δ i q g + Δ u d t Δ e q = G a c c N ( X ) ( Δ i q g ∗ − Δ i q g ) + ω 0 L f Δ i d g + Δ u q t (23)

The parts of the main circuit and the PLL are the same as that for the DCO. Combining Eqs 12, 13, Eqs 15–18, and Eq. 23, we obtain the SISO model of the system { Δ i x g = 0 Δ i y g = G 1 _ s ( s ) 1 + N ( X ) G 0 _ s ( s ) Δ i q g ∗ (24) where G 1 _ s ( s ) = G a c c N ( X ) s L f , G 0 _ s ( s ) = G a c c ( 1 − H p l l s L g I x g 0 ) s L f .

Again, similar to the DCO, for G _{1_s} (s) we obtain that its denominator term is L _f s ². Therefore, G _{1_s} (s) has no right-half-plane pole and it has no effect on the system stability. The system stability should depend on G _{0_s} (s) and N(X) only. The corresponding SISO model is given in Figure 5B. From Figure 5, it can be seen that for both the DCO and the SCO, the nonlinear limiter has been well separated from the other parts. In particular, for the DCO, the original MIMO model can be simplified into a SISO model. Therefore, for the analyses of the DCO and SCO, they share the same control diagram, expect the different, explicit forms of G _{0_d} (s) for the DCO and G _{0_s} (s) for the SCO. Next it is ready to use the DFNC to calculate the oscillation amplitude and frequency.

We carry out the EMT simulation for the grid-tied VSC system based on MATLAB/Simulink. As one example, without losing generality, the grid inductance L _g = 1.2 is chosen to represent a weak grid, and the other parameters are set as L _f = 0.1, U _g = 1.0, I _xg0 = 0.8, I _yg0 = −0.21. All these parameters are listed in Table 1 (Yang et al., 2020a; Ma et al., 2021). In addition, in Supplementary AppendixS1, we give the details for how to obtain the state-space model and the associated Jacobian matrix. All these can be used for small-signal stability analysis. Accordingly, the parametric stability regions for the ACC and the PLL are calculated and shown in Figures 6A,B, respectively. Within the stable regions, each variable converges to a stable fixed point. While outside of the stable regions, complex dynamic behaviors such as limit cycles or chaos may occur and lead to the action of the saturation limiters. Next, we will focus on exploring self-sustained oscillations and studying the limiter effects within the unstable regions.

First, the typical controller parameters are chosen: k _p,acc = 0.6, k _i,acc = 160, k _p,pll = 310, and k _i,pll = 10,000. The input waveforms of the d-axis and q-axis limiters from the EMT simulation are shown in Figures 7A,B, respectively. The corresponding waveform of the VSC terminal voltage e _d is shown in Figure 7C. From Figures 7A–C it can be seen that the two limiters show sustained oscillations with a constant frequency and amplitude. Their frequency and amplitude are measured and presented in the panels. In addition, other variables also show a sustained oscillation with the same frequency. Therefore, the system exhibits the DCO for double-clipped oscillation for both d-axis and q-axis limiters.

Next, if we change the control parameters, the system dynamics could substantially change. For instance, when the parameters k _p,acc = 0.6, k _i,acc = 160, k _p,pll = 315, and k _i,pll = 20,000 are chosen, the typical pattern of SCO appears, as shown in Figures 7D–F. From Figure 7D, it can be seen that the input signal of the d-axis limiter does not show oscillation, and after it touches the limiter boundary (indicated by a dashed curve), it is fixed and remains constant. In a sharp contrast to this, in Figure 7E the q-axis limiter still exhibits a sustained oscillation with a constant amplitude. Again, the frequency is the same for all other variables, as shown in Figure 7F.

To further verify the accuracy of the theoretical model, the time-domain waveforms of i _d and i _q under these two types of oscillation are given in Figure 8. Figures 8A,B are the waveforms under the first set of parameters. Oppositely, Figures 8C,D show the waveforms under the second set of parameters. It can be seen from the EMT simulation that the corresponding currents in the PLL dq frame also oscillate at 55.73HZ and 116.55HZ, respectively. The oscillation frequencies are consistent with the two cases in Figures 7C,F. Therefore, we can see that under these two sets of parameters, all system variables show sustained oscillations, indicating that both limiters need to be considered.

Then these EMT simulation results are compared with the DF models for both the DCO and SCO. We bring these two sets of controller parameters into G _{0_d} (s) in Eq. 20 and G _{0_s} (s) in Eq. 24 and plot their corresponding DFNC analytical results in Figures 9A,B, respectively. From these results, we can see that indeed the Nyquist curve of G _{0_d} (s) and G _{0_s} (s) intersect the negative reciprocal curve of the DF, −1/N(X), for the DCO and SCO, respectively, indicative of an occurrence of sustained oscillation. The detailed quantitative results are summarized in Table 2, where f _model is the oscillation frequency from the DF model, f _sim is the oscillation frequency from the EMT simulation, correspondingly X _model is the oscillation amplitude from the DF model, X _sim is the oscillation amplitude from the EMT simulation, and E _f and E _X are their corresponding relative errors. Clearly for both the DCO and SCO, the mismatches are less than 10%, and comparatively the error for the SCO is much smaller. These comparisons well demonstrate the validity of the DF theoretical model.

It is well-known that control parameters have a significant influence on the oscillations (Xue et al., 2019; Ma et al., 2020). We take the SCO model as an example. We still calculate the oscillation amplitude and frequency based on the DFNC and compare them with the corresponding EMT simulation results.

In this subsection, we study the ACC parameters. The parameters of the PLL are k _p,pll = 315 and k _i,pll = 20,000, with other parameters keeping unchanged. As the first study, we select the proportional coefficient of the ACC as k _p,acc = 0.6 and gradually increase its integral coefficient k _i,acc as 150, 160, and 170. The corresponding Nyquist curves of G _{0_s} (s) for different k _i,acc and −1/N(X) are shown in Figure 10A. Similarly, we keep k _i,acc = 240 and gradually increase k _p,acc as 0.8, 0.9, and 1.0, with the corresponding results shown in Figure 10B.

From Figure 10 we can see that for different parameters, G _{0_s} (s) always intersects with −1/N(X), i.e., the system always has sustained oscillations. In addition, from Figure 10A one can find that with increasing k _i,acc , the oscillation amplitude gradually increases, as the intersection point moves left. Note that roughly N(X) is a monotonic-decrease function of X for a fixed a in Eq. 5. From Figure 10B, one can find that oppositely with increasing k _p,acc , the oscillation amplitude gradually decreases, as the intersection point moves right.

The corresponding time-domain waveforms of the q-axis limiter under the above six different parameters are illustrated in Figure 11, for the variations of k _i,acc and k _p,acc , respectively. From these plots, we can see that the system really shows the predicted tendency of oscillation amplitude and frequency with the parameters, consistent with the DF analysis. For quantitative comparisons, Table 3 lists the oscillation amplitude and frequency from the model analyses and the numerical simulations, accompanying with their mismatches by relative errors. Basically, the error is tiny. Therefore, they have well proved the correctness of the DF analysis.

Now we study the PLL parameters. We keep the parameters of the ACC and change the proportional and integral parameters of the PLL. The corresponding DF results are shown in Figure 12. In Figure 12A, k _i,pll = 15,000, 16,000, and 17,000. In Figure 12B, k _p,pll = 290, 300, and 310. The DFNC plots again show that the system has the sustained oscillations for all these six cases. In addition, we know that for smaller k _i,pll in Figure 12A and larger k _p,pll in Figure 12B, the oscillation amplitude should increase.

For these six parameters of the PLL, the corresponding EMT simulation results are shown in Figure 13. Clearly the tendencies of the oscillation amplitude and frequency are consistent with the theoretical model. Again, the quantitative comparisons in Table 4 demonstrates that the oscillation amplitude and frequency predicted by the theoretical model are perfect. Finally, if we compare the variation ranges for all these four different parameter sets k _i,acc , k _p,acc , k _i,pll , and k _p,pll in Tables 2, 3, we can find that comparatively k _p,pll has a greater influence than all other parameters. This might infer that k _p,pll is a key parameter for oscillation suppression and stability improvement.

For the system with the double saturation limiters, no matter whether the system has the DCO or SCO, we know that −1/N(X) is always located at the negative real axis, and G(s) is independent of the saturation limiter. In addition, according to the explicit DF for the saturation limiter in Eq. 5, the limiter boundary a is positively correlated with the oscillation amplitude X. Therefore, we may deduce that the size of the limiter boundary should affect the oscillation amplitude, but not the oscillation frequency, and in addition the oscillation amplitude should monotonically increase with the limiter boundary a.

These predictions have been well verified by the EMT simulation. For example, we keep all the main circuit parameters and the controller parameters unchanged, and set the limiter boundaries of the q-axis as 0.03, 0.028, and 0.026. Under these parameters, the Fourier spectrum analyses of the input signal of the q-axis limiter in Figure 14A shows that the oscillation frequencies are unchanged, but their amplitudes increase with the increase of a. Clearly this phenomenon is the same as what has been reported in the grid-tied VSC system with a single saturation limiter in the recent study (Wu et al., 2021).

In the end, let us study the extreme case of an infinite limiter boundary, i.e., the limiter is completely ignored. The EMT simulations are shown in Figure 14B, where clearly in the absence of the saturation limiter, the input of the q-axis of the ACC diverges after the initial exponential increase, whereas with the limiter, a sustained oscillation appears. Therefore, the limiter is a key factor for the appearance of sustained oscillation.