Power Oscillation Analysis of PMSG Wind Power Generation System Considering Power Control Nonlinearity

High proportion of renewable energy generation, especially wind power generation, has changed the dynamic behavior of power systems and leads to emerging stability issues such as oscillation accidents. The commonly used stability analysis methods based on linear system theory cannot tackle the influence of nonlinear elements in the control loop. To fill this gap, the describing function method is used to analyze the stability and power oscillations of a grid-connected permanent magnet synchronous generator (PMSG) wind power generation system in this paper. A complete PMSG system model is taken into consideration, including a wind turbine, a generator, a machine-side converter, a grid-side converter, and a weak grid. Particularly, the nonlinear element introduced by maximum power point tracking control in the power control loop is modeled and analyzed, and the influences of critical system parameters are studied. With the high accuracy of the describing function method, it is revealed that the machine-side converter mainly influences the oscillation amplitude, while the grid-side converter mainly influences the oscillation frequency. Furthermore, the following results can be obtained: 1) High bandwidth of maximum power point tracking control and rotation speed control may enlarge the oscillation amplitude, while a fast current loop of the machine-side converter is beneficial to the stability; 2) Grid-side voltage loop can significantly affect the frequency of the oscillation, and improper high bandwidth of a phase-locked loop (PLL) is harmful to system stability; 3) A stronger grid is not always beneficial to the system stability, and the influence of grid strength should be carefully analyzed and examined when nonlinear elements are considered. All analyses are verified by simulations based on MATLAB/Simulink.


INTRODUCTION
Converter interfaced power generation technologies for renewable energy integration have significantly changed the dynamic behavior of power systems, leading to new types of stability problems such as power oscillations (Hatziargyriou et al., 2021). High proportion of renewable energy generation has brought a huge risk of power oscillations, and seriously threatens the stability of the power systems (Chi et al., 2019). In recent years, many power oscillation accidents have been reported in, for example, high-voltage direct current (HVDC) systems (Amin et al., 2017), offshore wind farms (Lyu et al., 2018), microgrids (Guo et al., 2017;Wang et al., 2017), and photovoltaic (PV) generators (Xia et al., 2019a;Zhao et al., 2019a). Some of the oscillations that are prone to occur in weak grid-connected wind power systems with the frequency from 2.5 to 50/60 Hz are classified as sub-synchronous oscillations (SSOs), and their harm is tremendous (Bi et al., 2017;Shi et al., 2020). It is reported that some SSO events that happened in Guyuan wind farms even involved more than a thousand wind turbine generators (Shair et al., 2019), and the oscillation frequency characteristics were strongly affected by the converter control parameters (Wang et al., 2015;Xie et al., 2017). Another severe SSO event that happened in Hami wind farms caused widely-propagated oscillations, leading to torsional oscillations of thermal power units 300 km away from the wind farms . Permanent magnet synchronous generator (PMSG)-based wind turbines were highly involved in these SSO events, and the role of PMSGs in these events is gradually gaining attention. Furthermore, the application of the PMSG-based wind power system is with rapid development due to its tremendous advantages (Mahmoud et al., 2021;Mahmoud et al., 2022), hence power oscillations of the PMSG-based wind power system deserve in-depth study.
Power oscillations generally start from small-signal instability, then rapidly grow due to the lack of damping in the systems, and finally evolve into sustained oscillations that can widely propagate until some of the devices participating in the oscillations are disconnected Liu and Xie, 2018;. The characteristics of the oscillations such as amplitude and frequency are strongly related to the control parameters of power converters, and are also influenced by the features of the interface grids, which can be described by grid strength (Wang and Blaabjerg, 2019).
To elaborate the mechanism of the oscillation phenomenon, a variety of linearized models and linear system methods are used for stability analyses. Eigenvalue analysis is a classic stability analysis method, based on which interactions among PMSGs are studied (Huang et al., 2019b;Zhao et al., 2019b). According to the root loci and participation factors of an SSO mode, it is found that a decrease of grid strength may drive the SSO mode towards the unstable region and this mode highly involves the DC capacitor, the AC network and the current control of grid-side converter (GSC) (Huang et al., 2019a). The eigenvalue analysis can be also combined with discrete-time state-space models so that it can be used for complicated high-order systems (Han et al., 2022). Another widely used method is impedance analysis, which has special advantages in "black-box" modeling (Sun, 2009). This method was originally designed for dc-dc converters and then was developed and applied to ac power systems (Liu et al., 2020). For single dc-ac converters, impedance models can be established in dq-domain or sequence-domain (Cespedes and Sun, 2014;Wen et al., 2015), and these two kinds of models are basically equivalent for stability analyses on the condition that there is no significant frequency coupling (Rygg et al., 2017). Based on dqdomain impedance analysis, it is pointed out that the phase-locked loop (PLL) with higher bandwidth may enlarge the frequency range of the negative incremental resistance characteristic introduced by converter controllers, and could make the system unstable (Wen et al., 2016). For wind power systems, impedance models of PMSG considering both machineside converter (MSC) and GSC are established, and it is proven that the parameters of the MSC controller have a noticeable effect on the impedance characteristic and significantly affect the stability margin (Xue et al., 2020;Liu et al., 2021). Furthermore, impedance analysis is also applied for large-scale power systems with wind farms as a result of being conveniently aggregated (Liu and Xie, 2018).
However, all aforementioned linearized methods are aimed at the small-signal instability stage of power oscillations, but are unsuitable for the stage of sustained oscillations which lasts for much more time with durative harm, since various nonlinear elements are participating in the sustained oscillations. The linearized methods also cannot predict whether the system can reach a new equilibrium point. Although some large-signal modeling and analysis methods such as the equal area method and the phase-plane method are applied to analyze whether new equilibrium points exist after a nonlinear transient process Wang et al., 2020), only extremely simplified systems can be analyzed and the oscillation characteristic cannot be reflected in these methods. To study the influence of nonlinear elements on the system stability, describing function (DF) method is introduced. The DF method, also called the harmonic linearization or the harmonic balance method, is a powerful frequency approximate method for nonlinear system analysis and design (Gasparyan, 2008). The main idea of this method is that a time-domain discontinuous nonlinear element can be linearized and replaced by a complex gain in the frequency domain. Hence, the nonlinearity which is commonly ignored by small-signal linearization methods can be analyzed. Furthermore, as a prediction of a limit cycle, or a steady-state periodical process with a certain frequency, the DF method gives a feasible way to calculate the amplitude and frequency of sustained oscillations. This feature can be utilized in the analysis of power oscillations in renewable energy systems, and provide valuable information for the identification and elimination of oscillation events (Xia et al., 2019b).
Power oscillations caused by the grid-connected PV generators are analyzed by the DF method . It is found that the power control nonlinearity has a significant influence on the system stability, and the improperly fast power control loop will cause sustained power oscillations. In (Xu et al., 2020) and (Wu et al., 2020), power oscillations for wind turbines are also studied based on the DF method. It is reported that the sustained SSO can be estimated more accurately compared with conventional eigenvalue analysis, and the nonlinearity may affect the frequency response of the system at sub-synchronous frequencies, which means the results of conventional smallsignal methods should be supplemented and carefully examined. However, the power control nonlinearity has not been considered in these articles, and more importantly, the MSCs and their controllers which have an important effect on system characteristics are omitted. Therefore, power oscillations involving a complete PMSG model with power control nonlinearity still need detailed analysis.
This paper adopts the DF method to analyze the power oscillations that occur in the weak grid-connected PMSG systems with the consideration of a complete PMSG model. First, the model of PMSG including the wind turbine, the generator, the MSC, the GSC, and the weak grid is established. In the model of MSC, the maximum power point tracking controller, i.e. the power control loop of the MSC, applies the perturbation and observation (P&O) algorithm, which is widely used in wind energy systems (Abdullah et al., 2012;Kumar and Chatterjee, 2016). The discontinuous nonlinear element of the P&O algorithm is modeled by a describing function, and other elements that can be linearized are modeled by a transfer function. Then, the stability of the system is analyzed by the DF method. The influence of MSC controller parameters, GSC controller parameters, grid strength, and PLL are studied respectively and the relationship between these parameters and the oscillation characteristics is revealed. The DF method is compared with the conventional linearized stability analysis method to show its advantages, and time-domain simulations are conducted to verify the accuracy of the method. The results show that the nonlinear element has a crucial effect on the system stability, and some of the results obtained by conventional linearized methods are incorrect under certain scenarios. A brief comparison with the existing research is shown in Table 1, and the main contributions of this paper are summarized as follows.
1) This paper adopts the describing function method to tackle the influence of power control nonlinearity, which cannot be tackled by the commonly used stability analysis method based on linear system theory.
2) The describing function-based stability analysis can calculate the amplitude and frequency of power oscillation accurately. 3) Some critical findings can be obtained by analyzing the influences of critical control parameters using describing function method. The different effects of GSC and MSC controllers on oscillation characteristics are distinguished.
The rest of this paper is organized as follows. In Section 2, the basic principles of the DF method are introduced. In Section 3, the complete model of the grid-connected PMSG system is established. In Section 4, stability and oscillations are analyzed by the DF method, and the simulation results are presented. In Section 5, the main conclusions of the analyses are drawn.

BRIEF INTRODUCTION TO DF METHOD
Describing function method is a kind of extended frequency response method that can be used to analyze the nonlinear behavior or predict potential limit cycle in nonlinear systems (Nassirharand, 2012). Figure 1 shows a nonlinear Single-Input-Single-Output (SISO) system with a single nonlinear element F(x), and the stability of the system is determined by the distribution of the roots of the closed-loop characteristic equation on the condition that F(x) can be properly linearized. The establishment of the describing function of a nonlinear element is a process of harmonic linearization (Atherton, 1975). Assume that a steady-state periodical process exists, i.e. the input φ(t) = 0 and x(t) = Asin(ωt), then the output y(t) of the nonlinear element is also a periodic function that can be expanded by Fourier series as To represent the nonlinear element by a describing function, there are several basic assumptions: Frontiers in Energy Research | www.frontiersin.org June 2022 | Volume 10 | Article 910716 1) The nonlinear element is odd, symmetrical, and memoryless, so that A 0 = 0. 2) Linear part G(s) is low-pass so that higher harmonics (beginning from the second) can be attenuated and neglected.
Based on these assumptions, the describing function is defined as the ratio of the fundamental term of the Fourier series of y(t) to the input signal x(t), i.e.

N(A)
A The describing function N(A) is a complex gain, but in the case of single-value nonlinearities, b 1 = 0, which means the nonlinear element does not cause any phase shift.
The nonlinear element F(x) is replaced by a linear element N(A) which depends on the unknown amplitude A, and the stability of the linearized system can be determined by the position of the roots of the closed-loop characteristic equation If the real part of all the roots is negative, the system will be stable. If the real part of any of the roots is positive, the system will be unstable. If the characteristic equation has a pure imaginary root (or a couple of complex-conjugate pure imaginary roots) and all other roots have negative real parts, the system will be critically stable and a steady oscillation with constant amplitude and frequency may occur. In the last case, the unknown amplitude A and the frequency ω can be calculated through two algebraic equations derived from (3), Once the closed-loop characteristic equation is established, the Nyquist criterion can be used, and the corresponding graphical analysis method will avoid repetitive solving of (3). Rewrite (3) as The stability of the whole system can be judged by the relative position of −1/N(A) and G(s) curves in the s plane. As shown in Figure 2, on the condition that G(s) does not have right-plane poles (the open-loop system is stable), the extended Nyquist criterion states that if the −1/N(A) curve is not surrounded by the G(s) curve, the system will be stable, and if −1/N(A) is surrounded by G(s), the system will be unstable. If the two curves intersect, the system will be critically stable, and at least one limit cycle may exist. Furthermore, the stability of the limit cycle depends on whether the system withstands small disturbances at the intersections. For example, provided that the system initially operates at the point P 1 in Figure 2, when it moves to P 1 ' due to a small disturbance, the amplitude A will continuously increase along the −1/N(A) curve because the system is unstable at P 1 ', until the operating point returns to P 1 . In the other direction, at the point P 1 ″, the amplitude A will decrease because the system is stable, and finally the operating point also returns to P 1 . Hence, the intersection P 1 corresponds to a stable limit cycle in the system, and the oscillation information can be derived by (4). Similarly, it can be inferred that P 2 corresponds to an unstable limit cycle, and any slight disturbance will cause the system to deviate from this operating point.

COMPLETE MODEL OF GRID-CONNECTED PMSG SYSTEM
In this section, a complete model of the grid-connected PMSG is derived considering the wind turbine, the generator, the MSC, the GSC, and the weak grid. Figure 3 shows the structure and the control strategy of the system. The wind turbine with a horizontal axis is directly coupled to the rotor of the PMSG, and the generator is symmetrical and cylindrical. The control methods are selected as voltage-oriented control under unity power coefficient for the gridside converter, zero d-axis current control for the machine-side converter, and P&O control for the maximum power point tracking. The output current of the generator is i a,b,cr , the mechanical angular velocity of the generator is ω g , the electrical angular velocity of the generator is ω e , and the electrical angle of the generator is θ e . A shunt capacitor C dc is connected between the converters, and the capacitor voltage is u dc . On the DC side, the current flowing out of the MSC is i dc2 , and the current flowing into the GSC is i dc1 . The output current of the GSC is i a,b,cg filtered by the inductor L f . The voltage at the point of common coupling (PCC) is u a,b,cg and the weak grid is represented as a grid inductor L g in series with an ideal voltage source u a,b,cs .

Model of the Wind Turbine and Drive Train
The wind turbine converts wind energy into mechanical motion and runs the generator. In general, the fraction of power extracted from the wind power is denoted by the power coefficient C p of the wind turbine, and the actual mechanical power P m provided by the wind turbine is where ρ denotes the air mass density, R denotes the blade radius, and v w denotes the wind velocity. C p is a function of tip velocity ratio λ (λ = ω g R/v w ) and pitch angle (in degrees) denoted by β, where Note that (7) and (8) are the mathematical approximation of a real curve of C p , and they have no physical meaning. When the PMSG works in maximum power point tracking mode, the wind velocity is smaller than the rated value, and hence only the condition that β = 0°is considered.
As the shaft of the wind turbine is directly connected to the rotor of the generator, their rotational speed is the same ω g . The drive train is modeled by a one-mass model, i.e.
where T m =P m /ω g is the mechanical torque provided by the wind turbine, T e is the electrical torque of the generator, and J is the equivalent inertia that contains both the turbine and the generator. Then the small-signal model (after Laplace transform) of the drive train is derived as where ΔT m B t Δω g , B t dTm dωg | ωg ω p g , ω p g denotes the steady-state value of ω g at equilibrium points. Besides, the mark Δ denotes small-signal disturbance of the corresponding variables, and capital letters with superscript * denote a steady-state value of the corresponding variables in this paper.

Model of the PMSG and MSC
Provided that the d-axis is aligned with the magnetic flux linkage of the rotor and the PMSG is a round rotor machine, the smallsignal model of PMSG can be expressed as ΔT e 3 2 n p ψ f Δi qr (12) where . L s and R s denote the armature inductance and stator resistance respectively, i dr and i qr denote the dq-axis stator currents respectively, and d dr and d qr denote the dq-axis duty ratios of the MSC controller respectively. ψ f is the permanent magnet flux and n p is the number of pole pairs. Based on (10) and (12), it is drawn that stator current disturbance can cause mechanical angular velocity disturbance, and their relationship is derived as ], and then Δω g G iqω Δi dqr , Δω e n p G iqω Δi dqr .
In the following, the model of the MSC controller is established. The MSC controller can be found in Figure 3, and it consists of the In Figure 3 MSC controller, P ref denotes the power reference value, P denotes the actual power output value, I drref and i qrref denote the dqaxis current reference value respectively, and ω ref g denotes the mechanical angular velocity reference value. Moreover, superscript "c" denotes variables in the controller dq frame, which should be distinguished from the system dq frame because the disturbance of electrical angular velocity can cause electrical angle disturbance, and hence affect the angle used for dq transformation. A typical P&O algorithm is shown in Figure 4, where T p is the control cycle and ε is the perturbation step. When needed this algorithm can make the PMSG output constant power in the maximum power point tracking mode.
In fact, by analyzing Figure 4, one can find that the accumulation of ε can be approximated as an integrator on the condition that T p is small enough. In each control cycle T p , ω ref g increases or decreases by one perturbation step ε, which means the controller is equivalent to an integrator whose integral coefficient is ±ε/T p , and the direction is determined by the sign of δP·δω g and P ref −P (note that P in Figure 3 has the same meaning as P n in Figure 4). Therefore, the algorithm shown in Figure 4 can be expressed as where sgn(x) is a sign function, when x ≥ 0, sgn(x) = 1, and when x < 0, sgn(x) = −1. Furthermore, the basic principle of the P&O algorithm states that the sign of δP·δω g reflects whether the system is currently running on the left or the right of the maximum power point, i.e. δP·δω g = 1 when ω g ≤ ω mpp (left), and δP·δω g = −1 when ω g > ω mpp (right). ω mpp denotes the mechanical angular velocity at the maximum power point of a certain wind velocity. If these two cases are discussed separately, (14) can be written as Form (15), it is shown that the power control loop contains the linear part and the nonlinear element. The linear part is composed of an integrator with coefficient ε/T p (or −ε/T p ), and the nonlinear element is the sign function. Based on the definition of the describing function introduced in Section 2, the describing function of the sign function can be derived as Then the approximate small-signal model of the power control loop in the frequency domain is The small-signal model of the rotational speed loop and the current control loop of the MSC controller is expressed as K ωi /s, K ωp and K ωi are the PI coefficients of the rotational speed loop respectively. H cr = K cpr + K cir /s, K cpr and K cir are the PI coefficients of the current control loop of the MSC controller respectively. In the MSC controller, the dq transformation is affected by the electrical angular disturbance, and hence the controller dq frame and the system dq frame should be distinguished. The superscript "c" denotes the variables in the controller dq frame. The same situation also occurs in the GSC model, but the difference is that angular disturbance in the GSC controller is introduced by PLL dynamics, which is widely used in converter modeling and stability analysis field. The small-signal disturbance of the electrical angular is derived as Δθ e n p Δω g s 3n 2 p ψ f 2s(B t − sJ) Δi qr H θe Δi qr where H θe 3n 2 p ψ f 2s(B t −sJ) , and the relationship of variables in the controller dq frame and the system dq frame is directly given as

Model of the GSC and Weak Grid
The small-signal model of the GSC is expressed as where and ω is the angular frequency of the power system. i dg and i qg denote the dq-axis output current of the GSC respectively, d dg and d qg denote the dq-axis duty ratios of the GSC controller respectively, and u dg and u qg denote the dq-axis PCC voltage respectively.
The small-signal model of the shunt capacitor is sC dc Δu dc 1.5 D p T dqr · Δi dqr + I p dqr · Δd dqr − 1.5 D p T dqg · Δi dqg The controller of GSC contains the commonly-used voltage loop and current loop, and the small-signal model is derived as K cig /s, K cpg and K cig are the PI coefficients of the current control loop of the GSC controller respectively. H v = K vp + K vi /s, K vp and K vi are the PI coefficients of the voltage control loop of the GSC controller respectively. U dcref is the reference value of dc voltage. As mentioned above, the superscript "c" denotes the controller dq frame. PLL dynamics lead to the angular difference between the controller and the system dq frame and are also significant to the stability of the system. Since the modeling process of PLL is well studied and the existence of angular difference is revealed in many papers (such as (Wen et al., 2016)), the small-signal model of PLL is directly given as (24) and the relationship of variables in the two different dq frames is given as (25).
where H pll = K ppll + K ipll /s, K ppll and K ipll are the proportional and the integral coefficient of the PI controller of PLL.
The small-signal model of the weak grid is expressed as where Z g sL g −ωL g ωL g sL g .
Moreover, the power output of the PMSG can be linearized as

A Complete Model of the Grid-Connected PMSG System
Through the above procedures, the model of each part of the PMSG has been established. It can be found that they are all linear except the model of the power control loop containing a nonlinear element. Combining these models and eliminating intermediate variables, the model of PMSG is reorganized as a Frontiers in Energy Research | www.frontiersin.org June 2022 | Volume 10 | Article 910716 typical SISO nonlinear system as shown in Figure 1. Note that ±ε/ (sT p ) is the linear part of the power control loop, and G 1 (s) expressed in (28) is a linear transfer function, and they form the linear part of the complete model. Finally, the linear model G(s) is expressed as (29).
In (29), 1/(1 + T f s) is a sampling filter of the power signal, and 1/(1 + 1.5T p s) is the PWM and controller delay. They are also not negligible for stability analyses.

POWER OSCILLATION ANALYSES BASED ON THE DF METHOD
In this section, the stability of the grid-connected PMSG system is analyzed. First, by comparing whether to consider the power control loop and its nonlinear element, the necessity and advantage of the DF method are clarified. Then several critical system parameters are chosen to show their influences on the system stability, and the accuracy of the DF method is presented. Finally, the influence of grid strength and PLL is studied. The basic system parameters are listed in  (Geng et al., 2011), which is used in all cases unless specifically mentioned.

Influence of MSC Controller Parameters
Commonly, the power control loop is ignored in the smallsignal analyses as the discontinuous nonlinear element cannot be handled. Then the system is considered linear and G 1 is the closed-loop transfer function of the system whose stability can be judged by the position of the poles of G 1 . However, when the power control loop is considered, a latent feedback path formed by the P&O algorithm is introduced, which makes G 1 become a part of the open-loop transfer function of the system whose stability is no longer determined by the poles of G 1 . In other words, the power control loop brings a new feedback path with a discontinuous nonlinear element, and the DF method is necessary for stability analyses.
To better explain the necessity of considering the power control loop, stability analyses and simulations when changing the rotational-speed-loop parameters K ωp and K ωi are shown in Figure 5. Figure 5A shows the stability analyses base on the DF method. When k changes from 0.5 to 1.5 (K ωp changes from 6,000 to 18,000, and K ωi changes from 1750 to 5,250), the intersection of −1/N(A) and G(s) moves in the direction that A increases and all these intersections correspond to stable limit cycles according to the principle introduced in Section 2. Therefore, oscillation occurs in the system, and the amplitude increases when k changes from 0.5 to 1.5. Figure 5B is based on the poles of G 1 which is the closed-loop transfer function when the power control loop is not considered. It is shown that when k changes from 0.5 to 1.5, the dominant pole moves away from the imaginary axis, which means the system should be stable and the stability should be strengthened as k increases. The conclusions of the two methods seem contradictory, however, simulation results presented in Figure 5C prove the correctness of the DF method. P ac is the output power of the PMSG at the PCC, and the power oscillation occurs. The amplitude expands when K ωp and K ωi increase, as predicted by the DF method.
If the intersection of −1/N(A) and G(s) curves corresponds to a stable limit cycle, the frequency and amplitude of the oscillations can be obtained. This information is valuable to the defense and Frontiers in Energy Research | www.frontiersin.org June 2022 | Volume 10 | Article 910716 8 identification of oscillations. It can be calculated that the oscillations at the point A 1 , A 2 , and A 3 are with the amplitude of 63, 124, and 185 kW respectively, and with the frequency of 7.3, 7.4, and 7.4 Hz respectively. The accuracy of the DF method is verified through FFT analyses, which are shown in Figure 5D. It is shown that the oscillations are with the amplitude of 63, 135, and 203 kW  respectively, and with the frequency of 7, 7, and 7 Hz respectively. These results are basically consistent with the theoretical analyses.
Moreover, the traditional method cannot analyze the influence of power control loop parameters on the system stability. Figure 6 shows the stability analyses and simulations when changing the perturbation step ε. The perturbation step is an important parameter because it determines how fast the power control loop can respond, or rather, the equivalent bandwidth of the power control loop. In Figure 6A, the intersection of −1/N(A) and G(s) also moves in the direction that A increases when ε changes from 0.1 × 10 −4 to 1 × 10 −4 , which means the increase of ε expands the amplitude of the power oscillation. The simulation result shown in Figure 6B verifies the analyses, and it can be seen that the power oscillation is enlarged as ε increases. The comparison between the analyses and the simulations is listed in Table 3, where FFT analyses are performed for the simulations. Figure 7 shows the effect of MSC current-loop parameters on the system stability. In Figure 7A, the intersection moves in the direction that A decreases when the bandwidth of the MSC current loop increases. It can be observed in the simulation that the sustained oscillation is suppressed with an incremental k, which is consistent with the analyses. Quantitative comparison results are also listed in Table 3.
This part shows the necessity and accuracy of the DF method. On the one hand, the latent feedback path with a nonlinear element can be considered in the stability analyses, which means the accuracy of the results is enhanced. On the other hand, some parameters that the conventional linear system method cannot handle can be analyzed by the DF method, which shows the superiority of the method. The higher bandwidth of the power control loop and rotational speed loop leads to larger oscillation amplitude, which should be avoided in parameter tuning. While higher bandwidth of the MSC current loop is beneficial to suppress the oscillation.

Influence of GSC Controller Parameters
The DF method can predict oscillation accidents precisely, and the influence of the GSC controller parameters will be analyzed in this part. Figure 8 shows the stability analyses and the simulations when changing the voltage-loop parameters K vp and K vi , and it can be observed that the amplitude is slightly changed but the frequency is significantly changed when K vp and K vi increase. According to the intersection information, the amplitude can be calculated as 139, 124, and 103 kW respectively when k is set to 0.2, 1, and 5 (K vp = 5k and K vi = 5k), and the frequency can be calculated as 3.9, 7.4, and 10.8 Hz respectively. Simulation results are shown in Figure 8B. The amplitude is 131, 122, and 100 kW respectively, and the frequency is 4, 7, and 9.3 Hz respectively. The simulation results are basically consistent with the analysis results.  Figure 9 shows the effect of the GSC current loop. The analyses and the simulations show that the amplitude and the frequency basically remain constant when the bandwidth of the GSC current loop increases, which means the GSC current-loop parameters have almost no effect on the oscillation characteristics. The quantitative comparison between the analyses and the simulations for changing K vp , K vi , K cpg , and K cig is shown in Table 4.
Although the GSC parameters have a smaller influence on the oscillation characteristics, the DF method can distinguish   these subtle differences, and accurate information on oscillations can be obtained. Furthermore, the bandwidth of the voltage loop has an evident impact on the oscillation frequency, but a minor impact on oscillation amplitude compared with MSC parameters. The bandwidth of the GSC current loop has almost no effect on the frequency and amplitude of the oscillation.

Influence of Grid Strength and PLL
Effects of grid strength and PLL on the system stability are strongly coupled, and hence they are discussed together. Grid strength is represented by the grid equivalent inductor L g , and a larger L g means a weaker grid strength. Figure 10 shows the stability analyses when changing L g and the PLL bandwidth, and Figure 10A is enlarged for better observation. When L g = 0.1 mH, the incremental bandwidth of PLL has almost no influence on the oscillation. When L g = 0.6 mH, the incremental bandwidth of PLL has a relatively evident influence on the oscillation. If k is set to 1 (K ppll = 2k, K ipll = 35k), the system is stable, but if k is set to 2, the system will be unstable because the open-loop transfer function G(s) has a pair of right-half plane poles, as shown in Figure 10B. Therefore, the excessively high bandwidth of PLL is harmful to the stability of the system, especially when the grid is weak.
Furthermore, a stronger grid is not always beneficial to keep the system stable, and it would enlarge the amplitude of oscillations under some specific system structure when nonlinear elements are considered. When L g is increased from 0.1 to 0.6 mH, the amplitude of the oscillation is decreased, in other words, oscillation accidents occurring in the weaker grid have less damage to the power system because the oscillation amplitude is suppressed, and the relationship between stability and the grid strength should be reconsidered and carefully analyzed when nonlinear elements are considered. This conclusion is proved by simulation results shown in Figure 11. When k = 1, the oscillation amplitude is smaller if L g is larger. The effect of the PLL is also verified in the simulation. When k = 2 and L g = 0.6 mH, the stable limit cycle vanishes and the system becomes unstable, which is consistent with the analysis results. The quantitative comparison for changing PLL and grid strength parameters is also shown in Table 4.

CONCLUSION
This paper presents the modeling and oscillation analyses of a grid-connected PMSG wind power generation system considering the discontinuous nonlinear element in the power controller. The complete model of the PMSG wind power generation system is established, including the wind turbine, the generator, the MSC, the GSC, and the weak grid. The MSC controller contains a nonlinear maximum power point tracking algorithm. The DF method is adopted to analyze the influence of the MSC controller parameters, the GSC controller parameters, and the grid strength. The analyses and simulations show that the DF method is more functional and accurate than the conventional linearized method. Based on the DF method,  several conclusions can be made. 1) Improper high bandwidth of a power control loop and a rotational speed loop is harmful to system stability, and it can enlarge the amplitude of oscillations. While a fast MSC current loop can suppress the oscillation.
2) The machine-side converter mainly influences the oscillation amplitude, while the grid-side converter mainly influences the oscillation frequency. More specifically, a higher bandwidth of a GSC voltage loop may lead to a higher oscillation frequency.
3) The effect of grid strength on stability should receive more attention, as a stronger grid is not always beneficial to stability when nonlinear elements are considered. In addition, PLLs with improper high bandwidth are also detrimental to the stability of the system. The DF method is with high practicability, and it can be further extended in some perspectives: 1) investigate various types of nonlinearity, such as saturation, nonlinear inductor, dead zone, and so on. 2) analyze the nonlinearity considering highorder harmonics. 3) consider multiple nonlinear elements and study their influence on stability.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.